Data from the solar wind, from LEO, and from the ground

TODO:

preamble around options for accessing the different data, with links to other projects..

For this notebook we focus on:

Using hapiclient to access solar wind data from OMNI
Using viresclient to access FACs & AEJs from Swarm, and B-field from ground observatories
Using xarray and matplotlib for data wrangling and plotting
- xarray tutorial website

The notebook can be used to extend the PyHC example gallery and the Swarm notebook gallery

import numpy as np import pandas as pd import xarray as xr import matplotlib.pyplot as plt from viresclient import SwarmRequest from hapiclient import hapi, hapiplot, hapitime2datetime

Let's choose an interesting time period to study - the "St. Patrick's day storm" of 17th March 2015

You can look at the wider context of this event using the interactive Space Weather Data Portal from the University of Colorado

start_time = '2015-03-15T00:00:00Z' end_time = '2015-03-20T00:00:00Z'

Part 1: solar wind data

HAPI is an access protocol supported by a wide array of heliophysics datasets. We can use the Python package "hapiclient" to retrieve data from HAPI servers. In this case we will access the OMNI HRO2 dataset which provides consolidated solar wind data, and then we will show how we can load these data into pandas and xarray objects.

To generate code snippets to use, and to see what data are available:
http://hapi-server.org/servers/#server=CDAWeb&dataset=OMNI_HRO2_1MIN¶meters=flow_speed&start=2000-01-01T00:00:00Z&stop=2000-02-01T00:00:00Z&return=script&format=python

For more examples with hapiclient, take a look at the demonstration notebooks:
https://github.com/hapi-server/client-python-notebooks

print(hapi.__doc__)

def fetch_omni_data(start, stop): server = 'https://cdaweb.gsfc.nasa.gov/hapi'; dataset = 'OMNI_HRO2_1MIN'; # Use parameters='' to request all data. Multiple parameters # can be requested using a comma-separated list, e.g., parameters='IMF,PLS' parameters = 'BX_GSE,BY_GSM,BZ_GSM,flow_speed'; opts = {'usecache': True} data, meta = hapi(server, dataset, parameters, start, stop, **opts) return data, meta data, meta = fetch_omni_data(start_time, end_time)

Data retrieved from HAPI are labelled numpy arrays (data) accompanied by a metadata dictionary (meta):

data

meta

You could extract an array for a particular value like data["BZ_GSM"], and use the metadata to get full descriptions and units for the chosen parameter.

The metadata sometimes contains fill values used during data gaps (e.g. the 9999... values appearing above). Let's use those to replace the gaps with NaN values:

def fill2nan(hapidata,hapimeta): """Replace bad values (fill values given in metadata) with NaN""" #Hapi returns metadata for parameters as #a list of dictionaries for metavar in hapimeta['parameters']: varname = metavar['name'] fillvalstr = metavar['fill'] if fillvalstr is None: continue vardata = hapidata[varname] mask = vardata==float(fillvalstr) nbad = np.count_nonzero(mask) print('{}: {} fills NaNd'.format(varname,nbad)) vardata[mask]=np.nan return hapidata, hapimeta data, meta = fill2nan(data,meta) data

We can load the data into a pandas DataFrame to more readily use for analysis:

def to_pandas(hapidata): df = pd.DataFrame( columns=hapidata.dtype.names, data=hapidata, ).set_index("Time") # Convert from hapitime to standard datetime df.index = hapitime2datetime(df.index.values) # Rename to Timestamp to match viresclient df.index.name = "Timestamp" return df to_pandas(data)

How can we get the extra information like the units from the metadata? Let's construct dictionaries, units and description, that allow easier access to these:

def get_units_description(meta): units = {} description = {} for paramdict in meta["parameters"]: units[paramdict["name"]] = paramdict.get("units", None) description[paramdict["name"]] = paramdict.get("description", None) return units, description units, description = get_units_description(meta) units, description

The xarray.Dataset object has advantages for handling multi-dimensional data and for attaching of metadata like units. Let's convert the data to an xarray.Dataset:

def to_xarray(hapidata, hapimeta): # Here we will conveniently re-use the pandas function we just built, # and use the pandas API to build the xarray Dataset. # NB: if performance is important, it's better to build the Dataset directly ds = to_pandas(hapidata).to_xarray() units, description = get_units_description(hapimeta) for param in ds: ds[param].attrs = { "units": units[param], "description": description[param] } return ds ds_sw = to_xarray(data, meta) ds_sw

Now let's plot these data:

fig, axes = plt.subplots(nrows=2, ncols=1, sharex=True, figsize=(15, 5)) for IMF_var in ["BX_GSE", "BY_GSM", "BZ_GSM"]: ds_sw[IMF_var].plot.line( x="Timestamp", linewidth=1, alpha=0.8, ax=axes[0], label=IMF_var ) axes[0].legend() axes[0].set_ylabel("IMF [nT]") ds_sw["flow_speed"].plot.line( x="Timestamp", linewidth=1, alpha=0.8, ax=axes[1] )

Part 2: Swarm data

2a: Auroral electrojet peaks from Swarm orbits

Since spacecraft move, it is difficult to extract a simple time series that can be easily tracked. From the complex Swarm product portfolio, we will pick a particular derived parameter: the peak auroral electrojet intensities derived from each pass over the current system. This signal tracks reasonably well from one orbit to the next (when separated into four orbital segments - accounting for two passes over the auroral oval in different local time sectors, and over the northern and southern hemispheres).

To keep things a bit simpler, we will retrieve data only from Swarm Alpha over the nothern hemisphere

def fetch_Swarm_AEJ(start_time, end_time): request = SwarmRequest() # Meaning of AEJAPBL: (AEJ) Auroral electrojets # (A) Swarm Alpha # (PBL) Peaks and boundaries from LC method # J_QD is the current intensity along QD-latitude contours # QDOrbitDirection is a flag (1, -1) marking the direction of the # satellite (ascending, descending) relative to the QD pole request.set_collection("SW_OPER_AEJAPBL_2F") request.set_products( measurements=["J_QD", "PointType"], auxiliaries=["QDOrbitDirection", "MLT"] ) # PointType of 0 refers to WEJ (westward electrojet) peaks # PointType of 1 refers to EEJ (eastward electrojet) peaks # See https://nbviewer.jupyter.org/github/pacesm/jupyter_notebooks/blob/master/AEBS/AEBS_00_data_access.ipynb#AEJxPBL-product request.set_range_filter("Latitude", 0, 90) # Northern hemisphere request.set_range_filter("PointType", 0, 1) # Extract only peaks data = request.get_between(start_time, end_time) ds_AEJ_peaks = data.as_xarray() return ds_AEJ_peaks ds_AEJ_peaks = fetch_Swarm_AEJ(start_time, end_time) ds_AEJ_peaks

Now we need some complex logic to plot the eastward and westward electrojet intensities, separated for each local time sector:

# Masks to identify which sector the satellite is in # and which current type (WEJ/EEJ) is given mask_asc = ds_AEJ_peaks["QDOrbitDirection"] == 1 mask_desc = ds_AEJ_peaks["QDOrbitDirection"] == -1 mask_WEJ = ds_AEJ_peaks["PointType"] == 0 mask_EEJ = ds_AEJ_peaks["PointType"] == 1 fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True) # Select and plot from the ascending orbital segments # on axes 0 # Eastward electrojet: _ds = ds_AEJ_peaks.where(mask_EEJ & mask_asc, drop=True) _ds["J_QD"].plot.line(x="Timestamp", ax=axes[0], label="EEJ") # Westward electrojet: _ds = ds_AEJ_peaks.where(mask_WEJ & mask_asc, drop=True) _ds["J_QD"].plot.line(x="Timestamp", ax=axes[0], label="WEJ") # Identify approximate MLT of sector _ds = ds_AEJ_peaks.where(mask_asc, drop=True) mlt = round(float(_ds["MLT"].mean())) axes[0].set_ylabel(axes[0].get_ylabel() + f"\nMLT: ~{mlt}") # ... and for descending segments # on axes 1 # Eastward electrojet: _ds = ds_AEJ_peaks.where(mask_EEJ & mask_desc, drop=True) _ds["J_QD"].plot.line(x="Timestamp", ax=axes[1], label="EEJ") # Westward electrojet: _ds = ds_AEJ_peaks.where(mask_WEJ & mask_desc, drop=True) _ds["J_QD"].plot.line(x="Timestamp", ax=axes[1], label="WEJ") # Identify approximate MLT of sector _ds = ds_AEJ_peaks.where(mask_desc, drop=True) mlt = round(float(_ds["MLT"].mean())) axes[1].set_ylabel(axes[1].get_ylabel() + f"\nMLT: ~{mlt}") axes[1].legend();

2b: Estimates of peak ground magnetic disturbances below satellite tracks

There are also predictions made from Swarm data of the location and strength of the peak disturbance on the ground (along the satellite ground-track) caused by the auroral electrojets. Note that this is from the AEJ_PBS (using the SECS method) collection rather than the AEJ_PBL (using the LC method) used above.

request = SwarmRequest() # request.available_collections() request.available_measurements("SW_OPER_AEJAPBS_2F:GroundMagneticDisturbance")

request.set_collection("SW_OPER_AEJAPBS_2F:GroundMagneticDisturbance") request.set_products( measurements=["B_NE"], # auxiliaries=["MLT"] #Quasidipole Magnetic Local Time - not available because there is no Radius parameter associated with this collection auxiliaries=["OrbitNumber"] ) data = request.get_between(start_time, end_time) ds = data.as_xarray()

TODO: explain plots below

Expected ground magnetic disturbance below SWARM Alpha (as timeseries)
Expected ground magnetic disturbance below SWARM Alpha in northern hemisphere only
Average ground magnetic disturbance from each SWARM A orbit (rough estimate of 'global' disturbance)

fig, ax = plt.subplots(nrows=1) ax.plot()

ds["B_NE"].plot.line(x="Timestamp")

_ds_N = ds.where(ds["Latitude"] > 0) _ds_N["Latitude"].plot.line(x="Timestamp")

print(ds)

# ds.isel({"Timestamp": slice(0, -1, 4)})

orbav_ds = ds.groupby('OrbitNumber').mean() orbrsmp_ds = ds.resample({'Timestamp':'90Min'}).mean() #90 minutes approximates the SWARM orbital period print(orbrsmp_ds) orbrsmp_ds["B_NE"].plot.line(x='Timestamp')

Data from the solar wind, from LEO, and from the ground

Part 1: solar wind data

Part 2: Swarm data

2a: Auroral electrojet peaks from Swarm orbits

2b: Estimates of peak ground magnetic disturbances below satellite tracks

Part 3: Ground observatory data

Recommended on Deepnote

Stock market analysis

10 ways to create NumPy arrays

Wide Residual Networks