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.585 μs (4 allocs: 34.039 KiB), 2.378 ms (18 allocs: 101.703 KiB), 3.745 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.579 μs (4 allocs: 34.039 KiB), 4.659 ms (18 allocs: 101.703 KiB), 10.303 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")

(707.688 μs (4 allocs: 34.039 KiB), 4.674 ms (18 allocs: 101.703 KiB), 7.286 ms (14 allocs: 304 bytes))