Compare Coordinate Transformations with IRBEM and PySPEDAS

See Coordinate Systems for more information.

Takeaway

Julia's implementation yields results very close to IRBEM's and PySPEDAS's one, and is an order of magnitude faster. (Julia's one uses finer interpolation than IRBEM's and PySPEDAS's one to determine IGRF coefficients and sun's direction, leading to more accurate transformations.)

References: cotrans, test_cotrans.py - PySPEDAS

Setup

using PySPEDAS
using SPEDAS
using SPEDAS: irbem_cotrans
using PythonCall
using DimensionalData
using Chairmarks
using Test

Setup using PySPEDAS test cases.

@py import pyspedas.cotrans_tools.tests.test_cotrans: CotransTestCases

pytest = CotransTestCases()
pytest.test_cotrans()

trange = ["2010-02-25/00:00:00", "2010-02-25/23:59:59"]
pyspedas.projects.themis.state(trange, probe="a", time_clip=true)

tha_pos = PySPEDAS.get_data("tha_pos") |> DimArray
tha_pos_gse = PySPEDAS.get_data("tha_pos_gse") |> DimArray
jl_tha_pos = set(tha_pos, Dim{:time}=>Ti)
jl_tha_pos_gse = set(tha_pos_gse, Dim{:time}=>Ti)

┌ 1440×3 DimArray{Float32, 2} tha_pos_gse ┐
├─────────────────────────────────────────┴────────────────────────────── dims ┐
  ↓ Ti Sampled{UnixTimes.UnixTime} [2010-02-25T00:00:00.000000000, …, 2010-02-25T23:59:00.000000000] ForwardOrdered Irregular Points,
  → v_dim Sampled{Int32} [1, …, 3] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 3 entries:
  "data_att"     => PyDict{String, Any}("coord_sys"=>"GSE", "units"=>"['km', 'k…
  "CDF"          => PyDict{String, Any}("VATT"=>PyDict{Any, Any}("CATDESC"=>"th…
  "plot_options" => PyDict{String, Any}("xaxis_opt"=>PyDict{Any, Any}("axis_lab…
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →                                  1       2         3
  2010-02-25T00:00:00.000000000  -30914.5  3310.4  -12103.8
 ⋮                                                 
  2010-02-25T23:59:00.000000000  -31397.2  3474.4  -12112.3

Validation

Transform coordinates using Julia native implementation, IRBEM, and PySPEDAS.

GEI <-> GEO

jl_tha_pos_geo = gei2geo(jl_tha_pos)
ir_tha_pos_geo = irbem_cotrans(jl_tha_pos', "GEI", "GEO")'
py_tha_pos_geo = PySPEDAS.get_data("tha_pos_new_geo") |> DimArray

@test jl_tha_pos_geo ≈ parent(py_tha_pos_geo)
@test jl_tha_pos_geo ≈ ir_tha_pos_geo

Test Passed

GEI <-> GSM

jl_tha_pos_gsm = gei2gsm(jl_tha_pos)
ir_tha_pos_gsm = irbem_cotrans(jl_tha_pos', "GEI", "GSM")'
pyspedas.cotrans("tha_pos", "tha_pos_new_gsm", coord_in="GEI", coord_out="GSM")
py_tha_pos_gsm = PySPEDAS.get_data("tha_pos_new_gsm") |> DimArray

@test isapprox(jl_tha_pos_gsm, parent(py_tha_pos_gsm), rtol=1e-5)
@test isapprox(jl_tha_pos_gsm, ir_tha_pos_gsm, rtol=1e-3)

Test Passed

GSE <-> GSM

jl_tha_pos_gsm = gse2gsm(jl_tha_pos_gse)
ir_tha_pos_gsm = irbem_cotrans(jl_tha_pos_gse', "GSE", "GSM")'
pyspedas.cotrans("tha_pos_gse", "tha_pos_new_gsm", coord_in="GSE", coord_out="GSM")
py_tha_pos_gsm = PySPEDAS.get_data("tha_pos_new_gsm") |> DimArray

@test isapprox(jl_tha_pos_gsm, parent(py_tha_pos_gsm), rtol=1e-5)
@test isapprox(jl_tha_pos_gsm, ir_tha_pos_gsm, rtol=1e-3)

Test Passed

Validate results: GEI/GEO transformations is quite accurate, while there are some differences in GSE/GSM transformations between Julia native implementation and IRBEM's one.

Benchmark

Depends on the transformation, Julia's implementation is about 10-40 times faster than IRBEM's (Fortran) implementation, and 20-50 times faster than PySPEDAS's (Python) implementation.

@b gei2geo($jl_tha_pos), irbem_cotrans($jl_tha_pos', "GEI", "GEO"), pyspedas.cotrans("tha_pos", "tha_pos_new_geo", coord_in="GEI", coord_out="GEO")

(115.175 μs (4 allocs: 34.039 KiB), 2.390 ms (18 allocs: 101.703 KiB), 3.777 ms (14 allocs: 304 bytes))

@b gei2gsm($jl_tha_pos), irbem_cotrans($jl_tha_pos', "GEI", "GSM"), pyspedas.cotrans("tha_pos", "tha_pos_new_gsm", coord_in="GEI", coord_out="GSM")

(595.421 μs (4 allocs: 34.039 KiB), 4.665 ms (18 allocs: 101.703 KiB), 10.508 ms (14 allocs: 304 bytes))

@b gse2gsm($jl_tha_pos_gse), irbem_cotrans($jl_tha_pos_gse', "GSE", "GSM"), pyspedas.cotrans("tha_pos_gse", "tha_pos_new_gsm", coord_in="GSE", coord_out="GSM")

(700.217 μs (4 allocs: 34.039 KiB), 4.685 ms (18 allocs: 101.703 KiB), 7.422 ms (14 allocs: 304 bytes))