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{Durations.Timestamp{Dates.Nanosecond}} [Durations.Timestamp{Dates.Nanosecond}("2010-02-25T00:00:00"), …, Durations.Timestamp{Dates.Nanosecond}("2010-02-25T23:59:00")] ForwardOrdered Irregular Points,
  → v_dim Sampled{Int32} [1, …, 3] ForwardOrdered Irregular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  SpaceDataModel.SchemaDict{PySPEDAS.PySPEDASSchema, String, Any, PythonCall.PyDict{String, Any}} with 3 entries:
  "CDF"          => PyDict{Any, Any}("VATT"=>PyDict{Any, Any}("CATDESC"=>"tha_p…
  "data_att"     => PyDict{Any, Any}("coord_sys"=>"GSE", "units"=>"km", "depend…
  "plot_options" => PyDict{Any, Any}("xaxis_opt"=>PyDict{Any, Any}("axis_label"…
└──────────────────────────────────────────────────────────────────────────────┘
 ↓ →                                                 1       2          3
  2010-02-25T00:00:00  -30914.5  3310.34  -12103.8
 ⋮                                                                 
  2010-02-25T23:59:00  -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")

(84.457 μs (29 allocs: 34.727 KiB), 2.393 ms (15 allocs: 90.383 KiB), 3.890 ms (16 allocs: 336 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")

(375.791 μs (29 allocs: 34.727 KiB), 4.697 ms (15 allocs: 90.383 KiB), 10.591 ms (16 allocs: 336 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")

(368.096 μs (29 allocs: 34.727 KiB), 4.695 ms (15 allocs: 90.383 KiB), 7.491 ms (16 allocs: 336 bytes))

Field line tracing and magnetic equator

find_magequator traces the field line through a point and locates the |B| minimum, matching IRBEM's find_magequator (IGRF only).

using OrdinaryDiffEqTsit5, Dates
using GeoCotrans: GEO
t = DateTime(2020, 1, 1)
pos = GEO(3.0, 0.5, 1.2)
x = collect(pos)
jl_eq = GeoCotrans.find_magequator(pos, t, Tsit5())
ir_eq = IRBEM.find_magequator(t, x, "GEO", (;); kext = 0)
@test jl_eq.pos ≈ ir_eq.XGEO rtol = 1e-3
@test jl_eq.Bmin ≈ ir_eq.Bmin rtol = 1e-3

Test Passed

Finding the equator costs about the same as IRBEM; tracing a whole field line is much faster in Julia.

@b GeoCotrans.find_magequator($pos, $t, Tsit5()), IRBEM.find_magequator($t, $x, "GEO", (;); kext = 0)

(58.229 μs (93 allocs: 6.531 KiB), 246.981 μs (35 allocs: 1.656 KiB))

@b trace($pos, $t, Tsit5()), IRBEM.trace_field_line($t, $x, "GEO", (;); kext = 0)

(38.812 μs (91 allocs: 8.625 KiB), 3.377 ms (46 allocs: 105.609 KiB))