DSECS: Dipolar Spherical Elementary Current Systems#
For more information about the project and the technique, see:
Vanhamäki, H., Maute, A., Alken, P. et al. Dipolar elementary current systems for ionospheric current reconstruction at low and middle latitudes. Earth Planets Space 72, 146 (2020). https://doi.org/10.1186/s40623-020-01284-1
import logging
import datetime as dt
from swarmpal.io import create_paldata, PalDataItem
from swarmpal.toolboxes import dsecs
import matplotlib.pyplot as plt
# To enable logging in the notebook, uncomment this line:
# logging.basicConfig(level=logging.INFO, force=True)
Fetching inputs to the toolbox#
def data_params(spacecraft="A"):
return dict(
server_url="https://vires.services/ows",
collection=f"SW_OPER_MAG{spacecraft}_LR_1B",
measurements=["B_NEC"],
models=["Model = CHAOS"], # currently must use name "Model"
auxiliaries=["QDLat"],
start_time="2016-03-18T11:00:00",
end_time="2016-03-18T14:00:00",
filters=["OrbitDirection == 1"], # Filters according to ascending passes
options=dict(asynchronous=False, show_progress=False),
)
data = create_paldata(
PalDataItem.from_vires(**data_params("A")),
PalDataItem.from_vires(**data_params("C")),
)
print(data)
<xarray.DataTree 'paldata'>
Group: /
├── Group: /SW_OPER_MAGA_LR_1B
│ Dimensions: (Timestamp: 5506, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 2016-...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 44kB 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A'
│ B_NEC (Timestamp, NEC) float64 132kB 1.413e+04 ... 4.687e+04
│ B_NEC_Model (Timestamp, NEC) float64 132kB 1.411e+04 ... 4.688e+04
│ Longitude (Timestamp) float64 44kB -26.01 -25.83 -25.66 ... 48.7 50.09
│ QDLat (Timestamp) float64 44kB -68.58 -68.54 -68.5 ... 82.04 82.01
│ Latitude (Timestamp) float64 44kB -82.61 -82.55 -82.49 ... 87.35 87.35
│ Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
│ Attributes:
│ Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGA_LR_1B_20160318T00...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: ['OrbitDirection == 1']
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
└── Group: /SW_OPER_MAGC_LR_1B
Dimensions: (Timestamp: 5500, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 2016-...
* NEC (NEC) <U1 12B 'N' 'E' 'C'
Data variables:
Spacecraft (Timestamp) object 44kB 'C' 'C' 'C' 'C' 'C' ... 'C' 'C' 'C' 'C'
B_NEC (Timestamp, NEC) float64 132kB 1.412e+04 ... 4.688e+04
B_NEC_Model (Timestamp, NEC) float64 132kB 1.41e+04 ... 4.688e+04
Longitude (Timestamp) float64 44kB -23.86 -23.69 -23.54 ... 49.16 50.54
QDLat (Timestamp) float64 44kB -68.47 -68.43 -68.39 ... 82.03 82.0
Latitude (Timestamp) float64 44kB -82.35 -82.29 -82.23 ... 87.35 87.35
Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
Attributes:
Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGC_LR_1B_20160318T00...
MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
AppliedFilters: ['OrbitDirection == 1']
PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
Applying the DSECS process#
Preprocess#
This initial process adds in Apex coordinates to the datasets.
p1 = dsecs.processes.Preprocess()
p1(data)
print(data)
<xarray.DataTree 'paldata'>
Group: /
│ Attributes:
│ PAL_meta: {"output_datasets": ["DSECS_output"]}
├── Group: /SW_OPER_MAGA_LR_1B
│ Dimensions: (Timestamp: 5506, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 201...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 44kB 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A'
│ B_NEC (Timestamp, NEC) float64 132kB 1.413e+04 ... 4.687e+04
│ B_NEC_Model (Timestamp, NEC) float64 132kB 1.411e+04 ... 4.688e+04
│ Longitude (Timestamp) float64 44kB -26.01 -25.83 -25.66 ... 48.7 50.09
│ QDLat (Timestamp) float64 44kB -68.58 -68.54 -68.5 ... 82.04 82.01
│ Latitude (Timestamp) float64 44kB -82.61 -82.55 -82.49 ... 87.35 87.35
│ Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
│ ApexLatitude (Timestamp) float64 44kB -69.37 -69.33 -69.29 ... 82.31 82.29
│ ApexLongitude (Timestamp) float64 44kB 25.92 26.02 26.12 ... 155.1 155.5
│ Attributes:
│ Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGA_LR_1B_20160318T00...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: ['OrbitDirection == 1']
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
├── Group: /SW_OPER_MAGC_LR_1B
│ Dimensions: (Timestamp: 5500, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 201...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 44kB 'C' 'C' 'C' 'C' ... 'C' 'C' 'C' 'C'
│ B_NEC (Timestamp, NEC) float64 132kB 1.412e+04 ... 4.688e+04
│ B_NEC_Model (Timestamp, NEC) float64 132kB 1.41e+04 ... 4.688e+04
│ Longitude (Timestamp) float64 44kB -23.86 -23.69 -23.54 ... 49.16 50.54
│ QDLat (Timestamp) float64 44kB -68.47 -68.43 -68.39 ... 82.03 82.0
│ Latitude (Timestamp) float64 44kB -82.35 -82.29 -82.23 ... 87.35 87.35
│ Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
│ ApexLatitude (Timestamp) float64 44kB -69.27 -69.23 -69.19 ... 82.3 82.28
│ ApexLongitude (Timestamp) float64 44kB 26.81 26.91 27.01 ... 155.2 155.6
│ Attributes:
│ Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGC_LR_1B_20160318T00...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: ['OrbitDirection == 1']
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
└── Group: /DSECS_output
Attributes:
PAL_meta: {"DSECS_Preprocess": {"output_dataset": "DSECS_output", "datas...
Analysis#
This process performs the DSECS analysis. It currently takes about 3 minutes to process one pass over the mid-latitudes.
%%time
p2 = dsecs.processes.Analysis()
p2(data)
CPU times: user 3min 41s, sys: 2.62 s, total: 3min 44s
Wall time: 2min 43s
<xarray.DataTree 'paldata'>
Group: /
│ Attributes:
│ PAL_meta: {"output_datasets": ["DSECS_output"]}
├── Group: /SW_OPER_MAGA_LR_1B
│ Dimensions: (Timestamp: 5506, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 201...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 44kB 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A'
│ B_NEC (Timestamp, NEC) float64 132kB 1.413e+04 ... 4.687e+04
│ B_NEC_Model (Timestamp, NEC) float64 132kB 1.411e+04 ... 4.688e+04
│ Longitude (Timestamp) float64 44kB -26.01 -25.83 -25.66 ... 48.7 50.09
│ QDLat (Timestamp) float64 44kB -68.58 -68.54 -68.5 ... 82.04 82.01
│ Latitude (Timestamp) float64 44kB -82.61 -82.55 -82.49 ... 87.35 87.35
│ Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
│ ApexLatitude (Timestamp) float64 44kB -69.37 -69.33 -69.29 ... 82.31 82.29
│ ApexLongitude (Timestamp) float64 44kB 25.92 26.02 26.12 ... 155.1 155.5
│ Attributes:
│ Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGA_LR_1B_20160318T00...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: ['OrbitDirection == 1']
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
├── Group: /SW_OPER_MAGC_LR_1B
│ Dimensions: (Timestamp: 5500, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 44kB 2016-03-18T11:00:00 ... 201...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 44kB 'C' 'C' 'C' 'C' ... 'C' 'C' 'C' 'C'
│ B_NEC (Timestamp, NEC) float64 132kB 1.412e+04 ... 4.688e+04
│ B_NEC_Model (Timestamp, NEC) float64 132kB 1.41e+04 ... 4.688e+04
│ Longitude (Timestamp) float64 44kB -23.86 -23.69 -23.54 ... 49.16 50.54
│ QDLat (Timestamp) float64 44kB -68.47 -68.43 -68.39 ... 82.03 82.0
│ Latitude (Timestamp) float64 44kB -82.35 -82.29 -82.23 ... 87.35 87.35
│ Radius (Timestamp) float64 44kB 6.832e+06 6.832e+06 ... 6.814e+06
│ ApexLatitude (Timestamp) float64 44kB -69.27 -69.23 -69.19 ... 82.3 82.28
│ ApexLongitude (Timestamp) float64 44kB 26.81 26.91 27.01 ... 155.2 155.6
│ Attributes:
│ Sources: ['CHAOS-8.1_static.shc', 'SW_OPER_MAGC_LR_1B_20160318T00...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: ['OrbitDirection == 1']
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
└── Group: /DSECS_output
│ Attributes:
│ PAL_meta: {"DSECS_Preprocess": {"output_dataset": "DSECS_output", "datas...
├── Group: /DSECS_output/0
│ ├── Group: /DSECS_output/0/currents
│ │ Dimensions: (x: 160, y: 7)
│ │ Coordinates:
│ │ Longitude (x, y) float64 9kB 348.4 349.1 349.8 ... 350.5 351.2 351.9
│ │ Latitude (x, y) float64 9kB -39.96 -39.96 -39.96 ... 39.54 39.54 39.54
│ │ Dimensions without coordinates: x, y
│ │ Data variables:
│ │ JEastDf (x, y) float64 9kB 12.72 12.93 13.3 13.73 ... 37.02 37.14 37.48
│ │ JNorthDf (x, y) float64 9kB 3.516 3.657 3.778 ... -2.793 -1.304 0.7451
│ │ Jr (x, y) float64 9kB 5.582 4.441 3.013 ... -1.077 0.4146 -0.4951
│ │ JEastCf (x, y) float64 9kB 3.991 3.817 3.626 ... -1.831 -2.357 -2.855
│ │ JNorthCf (x, y) float64 9kB 15.71 15.87 16.01 ... 3.147 3.519 3.884
│ │ JEastTotal (x, y) float64 9kB 16.71 16.75 16.93 17.2 ... 35.19 34.78 34.62
│ │ JNorthTotal (x, y) float64 9kB 19.22 19.52 19.78 ... 0.3541 2.215 4.629
│ │ Attributes:
│ │ description: DSECS result geographic coordinates
│ │ ionospheric height (km): 6501
│ │ Time interval: 2016-03-18T11:06:00 - 2016-03-18T11:37:13
│ │ Mean time: 2016-03-18T11:21:36
│ ├── Group: /DSECS_output/0/Fit_Alpha
│ │ Dimensions: (Timestamp: 1874, NEC: 3)
│ │ Coordinates:
│ │ * Timestamp (Timestamp) datetime64[s] 15kB 2016-03-18T11:06:00 ... 2016-0...
│ │ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ │ Data variables:
│ │ B_NEC_fit (Timestamp, NEC) float64 45kB 0.9512 10.68 ... -10.07 -9.467
│ │ B_NEC_data (Timestamp, NEC) float64 45kB 1.2 12.2 7.152 ... -11.31 -9.008
│ │ Longitude (Timestamp) float64 15kB -11.22 -11.21 -11.2 ... -9.878 -9.871
│ │ Latitude (Timestamp) float64 15kB -59.96 -59.9 -59.83 ... 59.92 59.99
│ │ Radius (Timestamp) float64 15kB 6.832e+06 6.832e+06 ... 6.816e+06
│ │ Attributes:
│ │ description: DSECS magnetic field fit.
│ │ satellite: Swarm Alpha
│ └── Group: /DSECS_output/0/Fit_Charlie
│ Dimensions: (Timestamp: 1874, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 15kB 2016-03-18T11:06:00 ... 2016-0...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ B_NEC_fit (Timestamp, NEC) float64 45kB -1.112 11.08 ... -10.11 -9.272
│ B_NEC_data (Timestamp, NEC) float64 45kB 0.683 13.04 ... -10.94 -8.747
│ Longitude (Timestamp) float64 15kB -9.758 -9.75 -9.743 ... -8.419 -8.411
│ Latitude (Timestamp) float64 15kB -59.69 -59.62 -59.56 ... 60.2 60.26
│ Radius (Timestamp) float64 15kB 6.831e+06 6.831e+06 ... 6.816e+06
│ Attributes:
│ description: DSECS magnetic field fit.
│ satellite: Swarm Charlie
└── Group: /DSECS_output/1
├── Group: /DSECS_output/1/currents
│ Dimensions: (x: 160, y: 7)
│ Coordinates:
│ Longitude (x, y) float64 9kB 324.9 325.7 326.4 ... 327.0 327.7 328.4
│ Latitude (x, y) float64 9kB -39.95 -39.95 -39.95 ... 39.55 39.55 39.55
│ Dimensions without coordinates: x, y
│ Data variables:
│ JEastDf (x, y) float64 9kB 17.71 17.8 18.02 17.99 ... 17.09 16.93 17.01
│ JNorthDf (x, y) float64 9kB 14.82 15.99 16.73 ... -17.06 -14.76 -11.38
│ Jr (x, y) float64 9kB -0.15 4.891 11.06 ... -19.09 -16.88 -14.24
│ JEastCf (x, y) float64 9kB 7.956 7.311 6.408 ... -4.29 -3.409 -2.703
│ JNorthCf (x, y) float64 9kB 32.09 31.83 31.77 ... -9.872 -9.959 -9.941
│ JEastTotal (x, y) float64 9kB 25.67 25.11 24.43 23.76 ... 12.8 13.52 14.31
│ JNorthTotal (x, y) float64 9kB 46.9 47.82 48.51 ... -26.93 -24.72 -21.32
│ Attributes:
│ description: DSECS result geographic coordinates
│ ionospheric height (km): 6501
│ Time interval: 2016-03-18T12:39:33 - 2016-03-18T13:10:46
│ Mean time: 2016-03-18T12:55:09
├── Group: /DSECS_output/1/Fit_Alpha
│ Dimensions: (Timestamp: 1874, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[s] 15kB 2016-03-18T12:39:33 ... 2016-0...
│ * NEC (NEC) <U1 12B 'N' 'E' 'C'
│ Data variables:
│ B_NEC_fit (Timestamp, NEC) float64 45kB -1.297 34.78 ... -32.73 -5.795
│ B_NEC_data (Timestamp, NEC) float64 45kB -3.149 40.9 4.93 ... -29.42 -2.885
│ Longitude (Timestamp) float64 15kB -34.69 -34.68 -34.67 ... -33.35 -33.34
│ Latitude (Timestamp) float64 15kB -59.95 -59.89 -59.82 ... 59.93 60.0
│ Radius (Timestamp) float64 15kB 6.832e+06 6.832e+06 ... 6.816e+06
│ Attributes:
│ description: DSECS magnetic field fit.
│ satellite: Swarm Alpha
└── Group: /DSECS_output/1/Fit_Charlie
Dimensions: (Timestamp: 1874, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[s] 15kB 2016-03-18T12:39:33 ... 2016-0...
* NEC (NEC) <U1 12B 'N' 'E' 'C'
Data variables:
B_NEC_fit (Timestamp, NEC) float64 45kB -3.886 35.33 ... -32.17 -6.292
B_NEC_data (Timestamp, NEC) float64 45kB -5.221 43.65 ... -23.76 -1.899
Longitude (Timestamp) float64 15kB -33.22 -33.21 -33.21 ... -31.88 -31.88
Latitude (Timestamp) float64 15kB -59.61 -59.55 -59.49 ... 60.27 60.33
Radius (Timestamp) float64 15kB 6.832e+06 6.832e+06 ... 6.816e+06
Attributes:
description: DSECS magnetic field fit.
satellite: Swarm CharlieOutputs are stored under a new “DSECS_output” branch of the datatree. This branch is further divided into N branches (“0”, “1”, “2”, …) depending on how many passes have been analysed, and under each branch into “currents” (the estimated current densities), “Fit_Alpha” and “Fit_Charlie” (the residuals for the fitted magnetic data).
data.groups
('/',
'/SW_OPER_MAGA_LR_1B',
'/SW_OPER_MAGC_LR_1B',
'/DSECS_output',
'/DSECS_output/0',
'/DSECS_output/1',
'/DSECS_output/0/currents',
'/DSECS_output/0/Fit_Alpha',
'/DSECS_output/0/Fit_Charlie',
'/DSECS_output/1/currents',
'/DSECS_output/1/Fit_Alpha',
'/DSECS_output/1/Fit_Charlie')
Plotting the outputs and fitted magnetic field#
Simple line plots of the currents from the central latitudinal slice#
The outputs in the datatree are enumerated, for each analyzed equatorial crossing included inside the timestamps. (“DSECS_output/0/currents”,”DSECS_output/1/currents” etc.)
print(data["DSECS_output/0/currents"].coords)
print(data["DSECS_output/0/currents"].data_vars)
Coordinates:
Longitude (x, y) float64 9kB 348.4 349.1 349.8 350.6 ... 350.5 351.2 351.9
Latitude (x, y) float64 9kB -39.96 -39.96 -39.96 ... 39.54 39.54 39.54
Data variables:
JEastDf (x, y) float64 9kB 12.72 12.93 13.3 13.73 ... 37.02 37.14 37.48
JNorthDf (x, y) float64 9kB 3.516 3.657 3.778 ... -2.793 -1.304 0.7451
Jr (x, y) float64 9kB 5.582 4.441 3.013 ... -1.077 0.4146 -0.4951
JEastCf (x, y) float64 9kB 3.991 3.817 3.626 ... -1.831 -2.357 -2.855
JNorthCf (x, y) float64 9kB 15.71 15.87 16.01 ... 3.147 3.519 3.884
JEastTotal (x, y) float64 9kB 16.71 16.75 16.93 17.2 ... 35.19 34.78 34.62
JNorthTotal (x, y) float64 9kB 19.22 19.52 19.78 ... 0.3541 2.215 4.629
latitudes = data["DSECS_output/0/currents"]["Latitude"][:, 3]
fig, ax = plt.subplots(4, 1, figsize=(10, 20))
lineE = ax[0].plot(
latitudes,
data["DSECS_output/0/currents"]["JEastTotal"][:, 3],
"r",
label="Eastward",
)
lineN = ax[0].plot(
latitudes,
data["DSECS_output/0/currents"]["JNorthTotal"][:, 3],
"k",
label="Northward",
)
ax[0].set_title("Total Current")
ax[0].set_ylabel("Current density (A/km)")
lineE = ax[1].plot(
latitudes, data["DSECS_output/0/currents"]["JEastCf"][:, 3], "r", label="Eastward"
)
lineN = ax[1].plot(
latitudes, data["DSECS_output/0/currents"]["JNorthCf"][:, 3], "k", label="Northward"
)
ax[1].set_title("Curl free current")
ax[1].set_ylabel("Current density (A/km)")
# ax[0].set_title('Total Current (DF + CF)')
lineE = ax[2].plot(
latitudes, data["DSECS_output/0/currents"]["JEastDf"][:, 3], "r", label="Eastward"
)
lineN = ax[2].plot(
latitudes, data["DSECS_output/0/currents"]["JNorthDf"][:, 3], "k", label="Northward"
)
ax[2].set_title("Divergence free current")
ax[2].set_ylabel("Current density (A/km)")
liner = ax[3].plot(
latitudes,
data["DSECS_output/0/currents"]["Jr"][:, 3],
"r",
label="Radial current density",
)
ax[3].set_title("Radial current density")
ax[3].set_ylabel("Current density (nA/m^2)")
for axv in ax:
axv.set_xlabel("Latitude", loc="left")
axv.legend()
(Preview) Quicklook plots#
from swarmpal.experimental import dsecs_plotting
Plot a specific pass
dsecs_plotting.plot_analysed_pass(data, "DSECS_output", pass_no=0)
/home/docs/checkouts/readthedocs.org/user_builds/swarmpal/envs/latest/lib/python3.11/site-packages/cartopy/io/__init__.py:242: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
warnings.warn(f'Downloading: {url}', DownloadWarning)
Automatically plot all passes
figs = dsecs_plotting.quicklook(data)
figs[0]
figs[1]
Create an animation of those quicklook plots (using ipwidgets)
(If you are viewing this on the web, it will not be interactive)
dsecs_plotting.quicklook_animated(data)
data.to_netcdf("dsecs_example.nc")