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

(120.474 μs (3 allocs: 33.836 KiB), 2.920 ms (15 allocs: 90.383 KiB), 4.261 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")

(485.563 μs (3 allocs: 33.836 KiB), 5.871 ms (15 allocs: 90.383 KiB), 11.481 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")

(578.245 μs (3 allocs: 33.836 KiB), 5.884 ms (15 allocs: 90.383 KiB), 8.729 ms (20 allocs: 408 bytes))