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(
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-18T11:50:00",
server_url="https://vires.services/ows",
options=dict(asynchronous=False, show_progress=False),
)
data = create_paldata(
PalDataItem.from_vires(**data_params("A")),
PalDataItem.from_vires(**data_params("C")),
)
print(data)
DataTree('paldata', parent=None)
├── DataTree('SW_OPER_MAGA_LR_1B')
│ Dimensions: (Timestamp: 3000, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:00:00 ... 2016-03-1...
│ * NEC (NEC) <U1 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A' 'A'
│ Longitude (Timestamp) float64 -26.01 -25.83 -25.66 ... 154.8 154.8 154.8
│ Latitude (Timestamp) float64 -82.61 -82.55 -82.49 ... 70.61 70.55 70.48
│ B_NEC (Timestamp, NEC) float64 1.413e+04 -1.228e+03 ... 4.675e+04
│ QDLat (Timestamp) float64 -68.58 -68.54 -68.5 ... 65.61 65.55 65.48
│ Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06 6.815e+06
│ B_NEC_Model (Timestamp, NEC) float64 1.411e+04 -1.23e+03 ... 4.676e+04
│ Attributes:
│ Sources: ['CHAOS-7_static.shc', 'SW_OPER_MAGA_LR_1B_20160318T0000...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: []
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
└── DataTree('SW_OPER_MAGC_LR_1B')
Dimensions: (Timestamp: 3000, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:00:00 ... 2016-03-1...
* NEC (NEC) <U1 'N' 'E' 'C'
Data variables:
Spacecraft (Timestamp) object 'C' 'C' 'C' 'C' 'C' ... 'C' 'C' 'C' 'C' 'C'
Longitude (Timestamp) float64 -23.86 -23.69 -23.54 ... 156.3 156.3 156.4
Latitude (Timestamp) float64 -82.35 -82.29 -82.23 ... 70.34 70.27 70.21
B_NEC (Timestamp, NEC) float64 1.412e+04 -1.64e+03 ... 4.661e+04
QDLat (Timestamp) float64 -68.47 -68.43 -68.39 ... 65.36 65.3 65.23
Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06 6.815e+06
B_NEC_Model (Timestamp, NEC) float64 1.41e+04 -1.655e+03 ... 4.662e+04
Attributes:
Sources: ['CHAOS-7_static.shc', 'SW_OPER_MAGC_LR_1B_20160318T0000...
MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
AppliedFilters: []
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)
DataTree('paldata', parent=None)
│ Dimensions: ()
│ Data variables:
│ *empty*
│ Attributes:
│ PAL_meta: {"DSECS_Preprocess": {"dataset_alpha": "SW_OPER_MAGA_LR_1B", "...
├── DataTree('SW_OPER_MAGA_LR_1B')
│ Dimensions: (Timestamp: 3000, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:00:00 ... 2016-03...
│ * NEC (NEC) <U1 'N' 'E' 'C'
│ Data variables:
│ Spacecraft (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A'
│ Longitude (Timestamp) float64 -26.01 -25.83 -25.66 ... 154.8 154.8
│ Latitude (Timestamp) float64 -82.61 -82.55 -82.49 ... 70.55 70.48
│ B_NEC (Timestamp, NEC) float64 1.413e+04 -1.228e+03 ... 4.675e+04
│ QDLat (Timestamp) float64 -68.58 -68.54 -68.5 ... 65.61 65.55 65.48
│ Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06
│ B_NEC_Model (Timestamp, NEC) float64 1.411e+04 -1.23e+03 ... 4.676e+04
│ ApexLatitude (Timestamp) float64 -69.37 -69.33 -69.29 ... 66.43 66.36
│ ApexLongitude (Timestamp) float64 25.92 26.02 26.12 ... -140.5 -140.4
│ Attributes:
│ Sources: ['CHAOS-7_static.shc', 'SW_OPER_MAGA_LR_1B_20160318T0000...
│ MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
│ AppliedFilters: []
│ PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
└── DataTree('SW_OPER_MAGC_LR_1B')
Dimensions: (Timestamp: 3000, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:00:00 ... 2016-03...
* NEC (NEC) <U1 'N' 'E' 'C'
Data variables:
Spacecraft (Timestamp) object 'C' 'C' 'C' 'C' 'C' ... 'C' 'C' 'C' 'C'
Longitude (Timestamp) float64 -23.86 -23.69 -23.54 ... 156.3 156.4
Latitude (Timestamp) float64 -82.35 -82.29 -82.23 ... 70.27 70.21
B_NEC (Timestamp, NEC) float64 1.412e+04 -1.64e+03 ... 4.661e+04
QDLat (Timestamp) float64 -68.47 -68.43 -68.39 ... 65.36 65.3 65.23
Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06
B_NEC_Model (Timestamp, NEC) float64 1.41e+04 -1.655e+03 ... 4.662e+04
ApexLatitude (Timestamp) float64 -69.27 -69.23 -69.19 ... 66.19 66.12
ApexLongitude (Timestamp) float64 26.81 26.91 27.01 ... -139.2 -139.1
Attributes:
Sources: ['CHAOS-7_static.shc', 'SW_OPER_MAGC_LR_1B_20160318T0000...
MagneticModels: ["Model = 'CHAOS-Core'(max_degree=20,min_degree=1) + 'CH...
AppliedFilters: []
PAL_meta: {"analysis_window": ["2016-03-18T11:00:00", "2016-03-18T...
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()
data = p2(data)
/home/docs/checkouts/readthedocs.org/user_builds/swarmpal/envs/latest/lib/python3.11/site-packages/swarmpal/toolboxes/dsecs/dsecs_algorithms.py:1490: UserWarning: Discarding nonzero nanoseconds in conversion.
date = pandas.to_datetime(SwA.Timestamp).mean().to_pydatetime()
CPU times: user 3min 20s, sys: 8.81 s, total: 3min 29s
Wall time: 2min 33s
Outputs 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).
print(data["DSECS_output"])
DataTree('DSECS_output', parent="paldata")
└── DataTree('0')
├── DataTree('currents')
│ Dimensions: (x: 160, y: 7)
│ Coordinates:
│ Longitude (x, y) float64 348.4 349.1 349.8 350.6 ... 350.5 351.2 351.9
│ Latitude (x, y) float64 -39.96 -39.96 -39.96 ... 39.54 39.54 39.54
│ Dimensions without coordinates: x, y
│ Data variables:
│ JEastDf (x, y) float64 13.01 13.22 13.63 14.06 ... 34.95 35.13 35.52
│ JNorthDf (x, y) float64 5.236 5.76 6.163 6.347 ... -8.563 -6.535 -3.691
│ Jr (x, y) float64 4.651 3.4 2.643 2.164 ... -1.813 -0.2648 -1.101
│ JEastCf (x, y) float64 3.822 3.704 3.575 3.379 ... -1.196 -1.675 -2.13
│ JNorthCf (x, y) float64 15.12 15.27 15.41 15.51 ... 4.131 4.495 4.856
│ JEastTotal (x, y) float64 16.83 16.93 17.21 17.44 ... 33.75 33.45 33.39
│ JNorthTotal (x, y) float64 20.36 21.03 21.57 21.85 ... -4.432 -2.041 1.165
│ Attributes:
│ description: DSECS result geographic coordinates
│ ionospheric height (km): 6501
│ Time interval: 2016-03-18T11:06:00.000000000 - 2016-03-18T11:3...
│ Mean time: 2016-03-18T11:21:36.500000000
├── DataTree('Fit_Alpha')
│ Dimensions: (Timestamp: 1874, NEC: 3)
│ Coordinates:
│ * Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:06:00 ... 2016-03-18...
│ * NEC (NEC) <U1 'N' 'E' 'C'
│ Data variables:
│ B_NEC_fit (Timestamp, NEC) float64 0.7235 10.67 6.973 ... -10.25 -9.54
│ B_NEC_data (Timestamp, NEC) float64 1.002 11.88 7.123 ... -11.35 -9.178
│ Longitude (Timestamp) float64 -11.22 -11.21 -11.2 ... -9.886 -9.878 -9.871
│ Latitude (Timestamp) float64 -59.96 -59.9 -59.83 ... 59.86 59.92 59.99
│ Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.816e+06 6.816e+06
│ Attributes:
│ description: DSECS magnetic field fit.
│ satellite: Swarm Alpha
└── DataTree('Fit_Charlie')
Dimensions: (Timestamp: 1874, NEC: 3)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2016-03-18T11:06:00 ... 2016-03-18...
* NEC (NEC) <U1 'N' 'E' 'C'
Data variables:
B_NEC_fit (Timestamp, NEC) float64 -1.145 11.08 6.579 ... -10.29 -9.798
B_NEC_data (Timestamp, NEC) float64 0.764 12.86 6.657 ... -10.8 -9.344
Longitude (Timestamp) float64 -9.758 -9.75 -9.743 ... -8.426 -8.419 -8.411
Latitude (Timestamp) float64 -59.69 -59.62 -59.56 ... 60.13 60.2 60.26
Radius (Timestamp) float64 6.831e+06 6.831e+06 ... 6.816e+06 6.816e+06
Attributes:
description: DSECS magnetic field fit.
satellite: Swarm 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 348.4 349.1 349.8 350.6 ... 350.5 351.2 351.9
Latitude (x, y) float64 -39.96 -39.96 -39.96 -39.96 ... 39.54 39.54 39.54
Data variables:
JEastDf (x, y) float64 13.01 13.22 13.63 14.06 ... 34.95 35.13 35.52
JNorthDf (x, y) float64 5.236 5.76 6.163 6.347 ... -8.563 -6.535 -3.691
Jr (x, y) float64 4.651 3.4 2.643 2.164 ... -1.813 -0.2648 -1.101
JEastCf (x, y) float64 3.822 3.704 3.575 3.379 ... -1.196 -1.675 -2.13
JNorthCf (x, y) float64 15.12 15.27 15.41 15.51 ... 4.131 4.495 4.856
JEastTotal (x, y) float64 16.83 16.93 17.21 17.44 ... 33.75 33.45 33.39
JNorthTotal (x, y) float64 20.36 21.03 21.57 21.85 ... -4.432 -2.041 1.165
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()
