DSECS: Dipolar Spherical Elementary Current Systems
Contents
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
│ Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06 6.815e+06
│ Latitude (Timestamp) float64 -82.61 -82.55 -82.49 ... 70.61 70.55 70.48
│ QDLat (Timestamp) float64 -68.58 -68.54 -68.5 ... 65.61 65.55 65.48
│ B_NEC_Model (Timestamp, NEC) float64 1.411e+04 -1.231e+03 ... 4.676e+04
│ B_NEC (Timestamp, NEC) float64 1.413e+04 -1.228e+03 ... 4.675e+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
Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06 6.815e+06
Latitude (Timestamp) float64 -82.35 -82.29 -82.23 ... 70.34 70.27 70.21
QDLat (Timestamp) float64 -68.47 -68.43 -68.39 ... 65.36 65.3 65.23
B_NEC_Model (Timestamp, NEC) float64 1.41e+04 -1.655e+03 ... 4.662e+04
B_NEC (Timestamp, NEC) float64 1.412e+04 -1.64e+03 ... 4.661e+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
│ Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06
│ Latitude (Timestamp) float64 -82.61 -82.55 -82.49 ... 70.55 70.48
│ QDLat (Timestamp) float64 -68.58 -68.54 -68.5 ... 65.61 65.55 65.48
│ B_NEC_Model (Timestamp, NEC) float64 1.411e+04 -1.231e+03 ... 4.676e+04
│ B_NEC (Timestamp, NEC) float64 1.413e+04 -1.228e+03 ... 4.675e+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
Radius (Timestamp) float64 6.832e+06 6.832e+06 ... 6.815e+06
Latitude (Timestamp) float64 -82.35 -82.29 -82.23 ... 70.27 70.21
QDLat (Timestamp) float64 -68.47 -68.43 -68.39 ... 65.36 65.3 65.23
B_NEC_Model (Timestamp, NEC) float64 1.41e+04 -1.655e+03 ... 4.662e+04
B_NEC (Timestamp, NEC) float64 1.412e+04 -1.64e+03 ... 4.661e+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/pre-cleanup/lib/python3.8/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 18s, sys: 8.42 s, total: 3min 26s
Wall time: 2min 30s
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.18 13.42 13.85 14.3 ... 35.37 35.52 35.88
│ JNorthDf (x, y) float64 5.362 5.907 6.327 6.517 ... -8.926 -6.852 -3.939
│ Jr (x, y) float64 4.801 3.575 2.802 ... -1.815 -0.2686 -1.102
│ JEastCf (x, y) float64 3.845 3.72 3.581 3.379 ... -1.198 -1.677 -2.133
│ JNorthCf (x, y) float64 15.15 15.3 15.43 15.53 ... 4.118 4.485 4.85
│ JEastTotal (x, y) float64 17.03 17.14 17.43 17.68 ... 34.17 33.84 33.75
│ JNorthTotal (x, y) float64 20.51 21.21 21.76 22.05 ... -4.808 -2.366 0.911
│ 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.5839 10.71 7.145 ... -10.08 -9.709
│ B_NEC_data (Timestamp, NEC) float64 0.8435 11.94 7.262 ... -11.15 -9.344
│ 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.291 11.12 6.75 ... -10.12 -9.95
B_NEC_data (Timestamp, NEC) float64 0.602 12.92 6.793 ... -10.62 -9.493
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.18 13.42 13.85 14.3 ... 35.37 35.52 35.88
JNorthDf (x, y) float64 5.362 5.907 6.327 6.517 ... -8.926 -6.852 -3.939
Jr (x, y) float64 4.801 3.575 2.802 ... -1.815 -0.2686 -1.102
JEastCf (x, y) float64 3.845 3.72 3.581 3.379 ... -1.198 -1.677 -2.133
JNorthCf (x, y) float64 15.15 15.3 15.43 15.53 ... 4.118 4.485 4.85
JEastTotal (x, y) float64 17.03 17.14 17.43 17.68 ... 34.17 33.84 33.75
JNorthTotal (x, y) float64 20.51 21.21 21.76 22.05 ... -4.808 -2.366 0.911
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()
