Compare Coordinate Transformations with IRBEM and PySPEDAS

See Coordinate Systems for more information about integration with SPEDAS.jl.

References: test_cotrans.py - PySPEDAS

Setup

using GeoCotrans
using PySPEDAS
using PySPEDAS.PythonCall
using DimensionalData
using Chairmarks
using Test
using IRBEM

function irbem_cotrans(pos, in, out)
    IRBEM.transform(pos.dims[1].val, pos.data', in, out)'
end

irbem_cotrans (generic function with 1 method)

Using data from 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.34  -12103.8
 ⋮                                                  
  2010-02-25T23:59:00.000000000  -31397.2  3474.34  -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")

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

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")

@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")

@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")

(110.175 μs (4 allocs: 34.039 KiB), 2.379 ms (18 allocs: 101.703 KiB), 3.730 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")

(605.480 μs (4 allocs: 34.039 KiB), 4.655 ms (18 allocs: 101.703 KiB), 10.189 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")

(671.232 μs (4 allocs: 34.039 KiB), 4.669 ms (18 allocs: 101.703 KiB), 7.279 ms (14 allocs: 304 bytes))