SPEDAS.jl

Space data analysis toolbox in Julia

package
Julia
Published

April 16, 2025

Think of SPEDAS in Julia.

Installation

Code
using Pkg
Pkg.add(["Speasy", "SPEDAS", "DimensionalData", "Dates", "Unitful", "LinearAlgebra", "CairoMakie"])
Pkg.add(url="https://github.com/Beforerr/Beforerr.jl")
   Resolving package versions...
  No Changes to `~/projects/beforerr/docs/projects/coding/2025_SpaceTools.jl/Project.toml`
  No Changes to `~/projects/beforerr/docs/projects/coding/2025_SpaceTools.jl/Manifest.toml`
    Updating git-repo `https://github.com/Beforerr/Beforerr.jl`
   Resolving package versions...
    Updating `~/projects/beforerr/docs/projects/coding/2025_SpaceTools.jl/Project.toml`
  [49aed396] + Beforerr v0.1.1 `https://github.com/Beforerr/Beforerr.jl#main`
    Updating `~/projects/beforerr/docs/projects/coding/2025_SpaceTools.jl/Manifest.toml`
  [cbdf2221] + AlgebraOfGraphics v0.10.3
  [49aed396] + Beforerr v0.1.1 `https://github.com/Beforerr/Beforerr.jl#main`
  [8be319e6] + Chain v0.6.0
  [a93c6f00] + DataFrames v1.7.0
  [1313f7d8] + DataFramesMeta v0.15.4
  [85a47980] + Dictionaries v0.4.5
  [38e38edf] + GLM v1.9.0
  [313cdc1a] + Indexing v1.1.1
  [4345ca2d] + Loess v0.6.4
  [7eb4fadd] + Match v2.4.0
  [c020b1a1] + NaturalSort v1.0.0
  [2dfb63ee] + PooledArrays v1.4.3
  [91c51154] + SentinelArrays v1.4.8
  [1277b4bf] + ShiftedArrays v2.0.0
  [3eaba693] + StatsModels v0.7.4
  [9ce81f87] + TableMetadataTools v0.1.0
Precompiling project...
   1537.4 ms  ✓ StatsModels
   1516.3 ms  ✓ GLM
  18307.0 ms  ✓ AlgebraOfGraphics
   5639.8 ms  ✓ AlgebraOfGraphics → AlgebraOfGraphicsUnitfulExt
   5612.6 ms  ✓ DimensionalData → DimensionalDataAlgebraOfGraphicsExt
   6132.7 ms  ✓ Beforerr
    387.7 ms  ✓ Beforerr → LatexExt
  7 dependencies successfully precompiled in 34 seconds. 438 already precompiled.

Getting data is easy

Code
using Speasy
using Dates
using DimensionalData
using CairoMakie
using SPEDAS
using Unitful
using LinearAlgebra

const spz = speasy
Python: <module 'speasy' from '/Users/zijin/projects/beforerr/docs/projects/coding/2025_SpaceTools.jl/.CondaPkg/.pixi/envs/default/lib/python3.12/site-packages/speasy/__init__.py'>
Code
t0 = DateTime("2008-09-05T10:00:00")
t1 = DateTime("2008-09-05T22:00:00")
da = get_data("cda/THA_L2_FGM/tha_fgs_gse", t0, t1) |> DimArray
[ Info: Cannot parse tha_fgs_gse unit nT*GSE*(All*Qs)
14386×3 DimArray{Float32, 2} tha_fgs_gse
├──────────────────────────────────────────┴───────────────────────────── dims ┐
  Ti Sampled{Dates.DateTime} [2008-09-05T10:00:02.232, …, 2008-09-05T21:59:59.539] ForwardOrdered Irregular Points,
  Y  Sampled{Int64} 1:3 ForwardOrdered Regular Points
├──────────────────────────────────────────────────────────────────── metadata ┤
  Dict{Any, Any} with 26 entries:
  "FILLVAL"          => Any[-1.0e30]
  "DEPEND_EPOCH0"    => "tha_fgs_epoch0"
  "FIELDNAM"         => "BXYZ GSE Coordinate, nT units"
  "VALIDMAX"         => Any[25000.0, 25000.0, 25000.0]
  "SI_CONVERSION"    => "1e-9>T"
  "TENSOR_ORDER"     => "1"
  "SC_ID"            => "d"
  "DEPEND_1"         => "tha_fgs_compno"
  "SCALE_TYP"        => "linear"
  "CALIB_SOFTWARE"   => "fgm_calibrate, thm_cal_fgm, thm_cal_fit"
  "CATDESC"          => "---- FGS (spin-resolution/~3 sec) magnetic  field B in…
  "DEPEND_TIME"      => "tha_fgs_time"
  "AVG_TYPE"         => "standard"
  "DISPLAY_TYPE"     => "time_series"
  "VALIDMIN"         => Any[-25000.0, -25000.0, -25000.0]
  "UNITS"            => "nT GSE (All Qs)"
  "VAR_NOTES"        => "Units are in nanotesla"
  "REPRESENTATION_1" => "Rep_xyz"
  "DEPEND_0"         => "tha_fgs_epoch"
  ⋮                  => ⋮
└──────────────────────────────────────────────────────────────────────────────┘
 ⋮      ⋱  

Array-like Time Series

No more need to get_data and store_data. All the functions and operations that work on abstract array would also work on the time series.

Code
@info da_max = maximum(da)
@info norm.(eachrow(da))[1:10]
[ Info: 65.44341
[ Info: Float32[75.64841, 75.72411, 75.5134, 75.73276, 75.492195, 75.506645, 75.522064, 75.78088, 75.69144, 75.647]
Code
using Test
da1 = get_data("cda/OMNI_HRO_1MIN/flow_speed", t0, t1) |> DimArray
da2 = get_data("cda/OMNI_HRO_1MIN/BX_GSE", t0, t1) |> DimArray

try
    da1 + da2
catch e
    println(e)
    @test e isa Unitful.DimensionError
end

# in different times so add would not work

tspan3 = DateTime("2008-09-05T10:00:00"), DateTime("2008-09-05T11:00:00")
tspan4 = DateTime("2008-09-05T11:00:00"), DateTime("2008-09-05T12:00:00")
da3 = get_data("cda/OMNI_HRO_1MIN/flow_speed", tspan3...) |> DimArray
da4 = get_data("cda/OMNI_HRO_1MIN/flow_speed", tspan4...) |> DimArray
try
    da3 + da4
catch e
    println(e)
end
Unitful.DimensionError(522.0f0 km s⁻¹, 4.01f0 nT)
DimensionMismatch("Lookup values for DimensionalData.Dimensions.Ti of [Dates.DateTime(\"2008-09-05T10:00:00\"), Dates.DateTime(\"2008-09-05T10:01:00\"), Dates.DateTime(\"2008-09-05T10:02:00\"), Dates.DateTime(\"2008-09-05T10:03:00\"), Dates.DateTime(\"2008-09-05T10:04:00\"), Dates.DateTime(\"2008-09-05T10:05:00\"), Dates.DateTime(\"2008-09-05T10:06:00\"), Dates.DateTime(\"2008-09-05T10:07:00\"), Dates.DateTime(\"2008-09-05T10:08:00\"), Dates.DateTime(\"2008-09-05T10:09:00\"), Dates.DateTime(\"2008-09-05T10:10:00\"), Dates.DateTime(\"2008-09-05T10:11:00\"), Dates.DateTime(\"2008-09-05T10:12:00\"), Dates.DateTime(\"2008-09-05T10:13:00\"), Dates.DateTime(\"2008-09-05T10:14:00\"), Dates.DateTime(\"2008-09-05T10:15:00\"), Dates.DateTime(\"2008-09-05T10:16:00\"), Dates.DateTime(\"2008-09-05T10:17:00\"), Dates.DateTime(\"2008-09-05T10:18:00\"), Dates.DateTime(\"2008-09-05T10:19:00\"), Dates.DateTime(\"2008-09-05T10:20:00\"), Dates.DateTime(\"2008-09-05T10:21:00\"), Dates.DateTime(\"2008-09-05T10:22:00\"), Dates.DateTime(\"2008-09-05T10:23:00\"), Dates.DateTime(\"2008-09-05T10:24:00\"), Dates.DateTime(\"2008-09-05T10:25:00\"), Dates.DateTime(\"2008-09-05T10:26:00\"), Dates.DateTime(\"2008-09-05T10:27:00\"), Dates.DateTime(\"2008-09-05T10:28:00\"), Dates.DateTime(\"2008-09-05T10:29:00\"), Dates.DateTime(\"2008-09-05T10:30:00\"), Dates.DateTime(\"2008-09-05T10:31:00\"), Dates.DateTime(\"2008-09-05T10:32:00\"), Dates.DateTime(\"2008-09-05T10:33:00\"), Dates.DateTime(\"2008-09-05T10:34:00\"), Dates.DateTime(\"2008-09-05T10:35:00\"), Dates.DateTime(\"2008-09-05T10:36:00\"), Dates.DateTime(\"2008-09-05T10:37:00\"), Dates.DateTime(\"2008-09-05T10:38:00\"), Dates.DateTime(\"2008-09-05T10:39:00\"), Dates.DateTime(\"2008-09-05T10:40:00\"), Dates.DateTime(\"2008-09-05T10:41:00\"), Dates.DateTime(\"2008-09-05T10:42:00\"), Dates.DateTime(\"2008-09-05T10:43:00\"), Dates.DateTime(\"2008-09-05T10:44:00\"), Dates.DateTime(\"2008-09-05T10:45:00\"), Dates.DateTime(\"2008-09-05T10:46:00\"), Dates.DateTime(\"2008-09-05T10:47:00\"), Dates.DateTime(\"2008-09-05T10:48:00\"), Dates.DateTime(\"2008-09-05T10:49:00\"), Dates.DateTime(\"2008-09-05T10:50:00\"), Dates.DateTime(\"2008-09-05T10:51:00\"), Dates.DateTime(\"2008-09-05T10:52:00\"), Dates.DateTime(\"2008-09-05T10:53:00\"), Dates.DateTime(\"2008-09-05T10:54:00\"), Dates.DateTime(\"2008-09-05T10:55:00\"), Dates.DateTime(\"2008-09-05T10:56:00\"), Dates.DateTime(\"2008-09-05T10:57:00\"), Dates.DateTime(\"2008-09-05T10:58:00\"), Dates.DateTime(\"2008-09-05T10:59:00\")] and [Dates.DateTime(\"2008-09-05T11:00:00\"), Dates.DateTime(\"2008-09-05T11:01:00\"), Dates.DateTime(\"2008-09-05T11:02:00\"), Dates.DateTime(\"2008-09-05T11:03:00\"), Dates.DateTime(\"2008-09-05T11:04:00\"), Dates.DateTime(\"2008-09-05T11:05:00\"), Dates.DateTime(\"2008-09-05T11:06:00\"), Dates.DateTime(\"2008-09-05T11:07:00\"), Dates.DateTime(\"2008-09-05T11:08:00\"), Dates.DateTime(\"2008-09-05T11:09:00\"), Dates.DateTime(\"2008-09-05T11:10:00\"), Dates.DateTime(\"2008-09-05T11:11:00\"), Dates.DateTime(\"2008-09-05T11:12:00\"), Dates.DateTime(\"2008-09-05T11:13:00\"), Dates.DateTime(\"2008-09-05T11:14:00\"), Dates.DateTime(\"2008-09-05T11:15:00\"), Dates.DateTime(\"2008-09-05T11:16:00\"), Dates.DateTime(\"2008-09-05T11:17:00\"), Dates.DateTime(\"2008-09-05T11:18:00\"), Dates.DateTime(\"2008-09-05T11:19:00\"), Dates.DateTime(\"2008-09-05T11:20:00\"), Dates.DateTime(\"2008-09-05T11:21:00\"), Dates.DateTime(\"2008-09-05T11:22:00\"), Dates.DateTime(\"2008-09-05T11:23:00\"), Dates.DateTime(\"2008-09-05T11:24:00\"), Dates.DateTime(\"2008-09-05T11:25:00\"), Dates.DateTime(\"2008-09-05T11:26:00\"), Dates.DateTime(\"2008-09-05T11:27:00\"), Dates.DateTime(\"2008-09-05T11:28:00\"), Dates.DateTime(\"2008-09-05T11:29:00\"), Dates.DateTime(\"2008-09-05T11:30:00\"), Dates.DateTime(\"2008-09-05T11:31:00\"), Dates.DateTime(\"2008-09-05T11:32:00\"), Dates.DateTime(\"2008-09-05T11:33:00\"), Dates.DateTime(\"2008-09-05T11:34:00\"), Dates.DateTime(\"2008-09-05T11:35:00\"), Dates.DateTime(\"2008-09-05T11:36:00\"), Dates.DateTime(\"2008-09-05T11:37:00\"), Dates.DateTime(\"2008-09-05T11:38:00\"), Dates.DateTime(\"2008-09-05T11:39:00\"), Dates.DateTime(\"2008-09-05T11:40:00\"), Dates.DateTime(\"2008-09-05T11:41:00\"), Dates.DateTime(\"2008-09-05T11:42:00\"), Dates.DateTime(\"2008-09-05T11:43:00\"), Dates.DateTime(\"2008-09-05T11:44:00\"), Dates.DateTime(\"2008-09-05T11:45:00\"), Dates.DateTime(\"2008-09-05T11:46:00\"), Dates.DateTime(\"2008-09-05T11:47:00\"), Dates.DateTime(\"2008-09-05T11:48:00\"), Dates.DateTime(\"2008-09-05T11:49:00\"), Dates.DateTime(\"2008-09-05T11:50:00\"), Dates.DateTime(\"2008-09-05T11:51:00\"), Dates.DateTime(\"2008-09-05T11:52:00\"), Dates.DateTime(\"2008-09-05T11:53:00\"), Dates.DateTime(\"2008-09-05T11:54:00\"), Dates.DateTime(\"2008-09-05T11:55:00\"), Dates.DateTime(\"2008-09-05T11:56:00\"), Dates.DateTime(\"2008-09-05T11:57:00\"), Dates.DateTime(\"2008-09-05T11:58:00\"), Dates.DateTime(\"2008-09-05T11:59:00\")] do not match.")

Unit-aware Time Series

Unit is not just a plot label, but also has dimensions meaning. (Unitful.jl integration)

Code
@info da_dim = dimension(da.data[1])
[ Info: NoDims

Tplot

Code
f = tplot(da)

Code
figure = (; size=(1200, 800))

tvars = [
    "cda/THA_L2_FGM/tha_fgs_gsmQ",
    "cda/OMNI_HRO_1MIN/flow_speed",
    "cda/OMNI_HRO_1MIN/E",
    "cda/OMNI_HRO_1MIN/Pressure",
]

tplot(tvars, t0, t1; figure)
Can't get THA_L2_FGM/tha_fgs_gsmQ without web service, switching to web service
[ Info: Cannot parse tha_fgs_gsmQ unit nT*GSM
Can't get THA_L2_FGM/tha_fgs_gsmQ without web service, switching to web service
[ Info: Cannot parse tha_fgs_gsmQ unit nT*GSM
[ Info: Resampling 14386×3 DimArray{Float32, 2} tha_fgs_gsmQ along dimension 1 from 14386 to 6070 points

Spectrogram

https://cdaweb.gsfc.nasa.gov/misc/NotesM.html#MMS1_FPI_FAST_L2_DES-DIST

Code
ts0 = "2019-01-02T15" |> DateTime
ts1 = "2019-01-02T22" |> DateTime

tvars = [
    "cda/MMS1_FGM_SRVY_L2/mms1_fgm_b_gse_srvy_l2_clean",
    "cda/MMS1_FPI_FAST_L2_DES-MOMS/mms1_des_energyspectr_omni_fast",
]
tplot(tvars, ts0, ts1; add_title=true)
Can't get MMS1_FGM_SRVY_L2/mms1_fgm_b_gse_srvy_l2_clean without web service, switching to web service
Can't get MMS1_FGM_SRVY_L2/mms1_fgm_b_gse_srvy_l2_clean without web service, switching to web service
[ Info: Resampling 403194×4 DimArray{Unitful.Quantity{Float32, 𝐌 𝐈⁻¹ 𝐓⁻², Unitful.FreeUnits{(nT,), 𝐌 𝐈⁻¹ 𝐓⁻², nothing}}, 2} mms1_fgm_b_gse_srvy_l2_clean along dimension 1 from 403194 to 6070 points

Overlay plots

Code
SDS = SPEDAS.DataSet
ts0 = DateTime("2021-01-17")
ts1 = DateTime("2021-01-18")

density_tvar = SDS(
    "Density",
    [
        "cda/PSP_SWP_SPI_SF00_L3_MOM/DENS",
        "cda/PSP_FLD_L3_RFS_LFR_QTN/N_elec",
        "cda/PSP_FLD_L3_SQTN_RFS_V1V2/electron_density"
    ]
)

tvars = [density_tvar]

tplot(tvars, ts0, ts1; add_title=true)
[ Info: Resampling 24704×1 DimArray{Unitful.Quantity{Float32, 𝐋⁻³, Unitful.FreeUnits{(cm⁻³,), 𝐋⁻³, nothing}}, 2} DENS along dimension 1 from 24704 to 6070 points
[ Info: Resampling 24272×1 DimArray{Unitful.Quantity{Float32, 𝐋⁻³, Unitful.FreeUnits{(cm⁻³,), 𝐋⁻³, nothing}}, 2} N_elec along dimension 1 from 24272 to 6070 points
[ Info: Resampling 24282×1 DimArray{Unitful.Quantity{Float32, 𝐋⁻³, Unitful.FreeUnits{(cm⁻³,), 𝐋⁻³, nothing}}, 2} electron_density along dimension 1 from 24282 to 6070 points

Complex requests and flexible layout

Code
data = let intervals = ["2019-01-02T15", "2019-01-02T16"]
    products = [
        spz.inventories.tree.cda.MMS.MMS1.FGM.MMS1_FGM_SRVY_L2.mms1_fgm_b_gse_srvy_l2_clean,
        spz.inventories.tree.cda.MMS.MMS1.SCM.MMS1_SCM_SRVY_L2_SCSRVY.mms1_scm_acb_gse_scsrvy_srvy_l2,
        spz.inventories.tree.cda.MMS.MMS1.DES.MMS1_FPI_FAST_L2_DES_MOMS.mms1_des_bulkv_gse_fast,
        spz.inventories.tree.cda.MMS.MMS1.DES.MMS1_FPI_FAST_L2_DES_MOMS.mms1_des_temppara_fast,
        spz.inventories.tree.cda.MMS.MMS1.DES.MMS1_FPI_FAST_L2_DES_MOMS.mms1_des_tempperp_fast,
        spz.inventories.tree.cda.MMS.MMS1.DES.MMS1_FPI_FAST_L2_DES_MOMS.mms1_des_energyspectr_omni_fast,
        spz.inventories.tree.cda.MMS.MMS1.DIS.MMS1_FPI_FAST_L2_DIS_MOMS.mms1_dis_energyspectr_omni_fast
    ]
    get_data(products, intervals)
end

let figure = (; size=(1200, 1200)), add_title = true
    f = Figure(; figure...)
    tplot(f[1, 1], data[1:3]; add_title)
    tplot(f[1, 2], [data[4:5], data[6:7]...]; add_title)
    f
end
Can't get MMS1_FGM_SRVY_L2/mms1_fgm_b_gse_srvy_l2_clean without web service, switching to web service

Disk-cached Time Series

WebServer (CDAWeb) -> Files -> Data structure (tplot variable / numpy array)

However, loading cdf files is moderately slow, especially for large data. We can improve the performance by disk-caching the data structure (numpy array).

WebServer (CDAWeb) -> Files -> Disk-cached data structure -> Data structure

Interactive Time Series

Use GLMakie to zoom in and Zoom out to explore data.

Code
tplot(tvars, t0, t1; figure)

Multiple dispatch tplot

We are power Matlab users. Some of us are Lisp hackers. Some are Pythonistas, others Rubyists, still others Perl hackers. There are those of us who used Mathematica before we could grow facial hair. There are those who still can’t grow facial hair. We’ve generated more R plots than any sane person should. C is our desert island programming language. We love all of these languages; they are wonderful and powerful. For the work we do — scientific computing, machine learning, data mining, large-scale linear algebra, distributed and parallel computing — each one is perfect for some aspects of the work and terrible for others. Each one is a trade-off. We are greedy: we want more. We want a language that’s open source, with a liberal license. We want the speed of C with the dynamism of Ruby. We want a language that’s homoiconic, with true macros like Lisp, but with obvious, familiar mathematical notation like Matlab. We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell. Something that is dirt simple to learn, yet keeps the most serious hackers happy. We want it interactive and we want it compiled.

We want tplot to understand a time series represented as a vector-like / matrix-like (series/heatmap/spectrogram) / Function / String; We want tplot to plot multiple time series on the same figure; We want tplot to be interactive as well as publication-ready; and we want tplot to be composable.

Code
figure = (; size=(1600, 1200))
f = Figure(; figure...)

tvars1 = ["cda/OMNI_HRO_1MIN/flow_speed", "cda/OMNI_HRO_1MIN/E", "cda/OMNI_HRO_1MIN/Pressure"]
tvars2 = ["cda/THA_L2_FGM/tha_fgs_gse"]
tvars3 = ["cda/OMNI_HRO_1MIN/BX_GSE", "cda/OMNI_HRO_1MIN/BY_GSE"]

t0 = DateTime("2008-09-05T10:00:00")
t1 = DateTime("2008-09-05T22:00:00")

tplot(f[1, 1], tvars1, t0, t1)
tplot(f[1, 2], tvars2, t0, t1)
tplot(f[2, 1:2], tvars3, t0, t1)
f
[ Info: Cannot parse tha_fgs_gse unit nT*GSE*(All*Qs)
[ Info: Cannot parse tha_fgs_gse unit nT*GSE*(All*Qs)
[ Info: Resampling 14386×3 DimArray{Float32, 2} tha_fgs_gse along dimension 1 from 14386 to 6070 points