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.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(DimArray, "tha_pos")
tha_pos_gse = PySPEDAS.get_data(DimArray, "tha_pos_gse")
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{Dates.DateTime} [Dates.DateTime("2010-02-25T00:00:00"), …, Dates.DateTime("2010-02-25T23:59:00")] ForwardOrdered Irregular Points,
  → v_dim
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 3 entries:
  "data_att"     => PyDict{Any, Any}("coord_sys"=>"GSE", "units"=>"['km', 'km',…
  "CDF"          => PyDict{Any, Any}("VATT"=>PyDict{Any, Any}("CATDESC"=>"tha_p…
  "plot_options" => PyDict{Any, Any}("xaxis_opt"=>PyDict{Any, Any}("axis_label"…
└──────────────────────────────────────────────────────────────────────────────┘
 2010-02-25T00:00:00  -30914.5  3310.4   -12103.8
 2010-02-25T00:01:00  -31104.2  3185.34  -12111.8
 ⋮                                       
 2010-02-25T23:59:00  -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(DimArray, "tha_pos_new_geo")

@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(DimArray, "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(DimArray, "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")

(119.212 μs (3 allocs: 33.836 KiB), 2.374 ms (15 allocs: 90.383 KiB), 4.158 ms (20 allocs: 408 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")

(601.505 μs (3 allocs: 33.836 KiB), 4.658 ms (15 allocs: 90.383 KiB), 11.317 ms (20 allocs: 408 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.944 μs (3 allocs: 33.836 KiB), 4.678 ms (15 allocs: 90.383 KiB), 8.513 ms (20 allocs: 408 bytes))