From the solar wind to the ground

From the solar wind to the ground#

Abstract: We demonstrate a basic analysis of a geomagnetic storm using hapiclient & viresclient to access data from the solar wind (OMNI IMF), Low Earth Orbit (Swarm-derived auroral electrojet estimates), and the ground (INTERMAGNET observatory magnetic measurements).

Packages to use#

hapiclient to access solar wind data from OMNI (alternatively we could use pysat)
- For more examples with hapiclient, take a look at the demonstration notebooks
viresclient to access AEJs from Swarm, and B-field from ground observatories
xarray and matplotlib for data wrangling and plotting
- See the xarray tutorial website to learn more

%load_ext watermark
%watermark -i -v -p viresclient,hapiclient,pandas,xarray,matplotlib

Python implementation: CPython
Python version       : 3.11.6
IPython version      : 8.18.0

viresclient: 0.12.0
hapiclient : 0.2.6
pandas     : 2.1.3
xarray     : 2023.12.0
matplotlib : 3.8.2

from copy import deepcopy
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, hapitime2datetime

Time selection#

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

We will use the same time window to fetch data from the different sources:

START_TIME = '2015-03-15T00:00:00Z'
END_TIME = '2015-03-20T00:00:00Z'

Solar wind data (OMNI)#

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.

OMNI Combined, Definitive 1-minute IMF and Definitive Plasma Data Time-Shifted to the Nose of the Earth’s Bow Shock, plus Magnetic Indices - J.H. King, N. Papatashvilli (AdnetSystems, NASA GSFC)

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

Here we will access five-minute-resolution measurements of the Interplanetary Magnetic Field (IMF) vector and the bulk flow speed of the solar wind:

def fetch_omni_data(start, stop):
    server     = 'https://cdaweb.gsfc.nasa.gov/hapi'
    dataset    = 'OMNI_HRO2_5MIN'
    parameters = 'BX_GSE,BY_GSM,BZ_GSM,flow_speed';
    data, meta = hapi(server, dataset, parameters, start, stop)
    return data, meta

data, meta = fetch_omni_data(START_TIME, END_TIME)

Data are automatically loaded as a NumPy structured array and metadata as a dictionary:

data

array([(b'2015-03-15T00:00:00.000Z', 9.99999000e+03,  9.99999000e+03,  9.99999000e+03, 99999.9),
       (b'2015-03-15T00:05:00.000Z', 9.99999000e+03,  9.99999000e+03,  9.99999000e+03, 99999.9),
       (b'2015-03-15T00:10:00.000Z', 9.99999000e+03,  9.99999000e+03,  9.99999000e+03, 99999.9),
       ...,
       (b'2015-03-19T23:45:00.000Z', 3.03999996e+00, -4.92000008e+00, -5.99999987e-02, 99999.9),
       (b'2015-03-19T23:50:00.000Z', 2.91000009e+00, -4.63999987e+00, -1.02999997e+00, 99999.9),
       (b'2015-03-19T23:55:00.000Z', 3.20000005e+00, -5.13999987e+00, -1.04999995e+00, 99999.9)],
      dtype=[('Time', 'S24'), ('BX_GSE', '<f8'), ('BY_GSM', '<f8'), ('BZ_GSM', '<f8'), ('flow_speed', '<f8')])

meta

{'HAPI': '2.0',
 'resourceURL': 'https://cdaweb.gsfc.nasa.gov/misc/NotesO.html#OMNI_HRO2_5MIN',
 'contact': 'J.H. King, N. Papatashvilli @ AdnetSystems, NASA GSFC',
 'parameters': [{'name': 'Time',
   'length': 24,
   'units': 'UTC',
   'type': 'isotime',
   'fill': None},
  {'name': 'BX_GSE',
   'description': 'Bx (nT), GSE',
   'units': 'nT',
   'type': 'double',
   'fill': '9999.99'},
  {'name': 'BY_GSM',
   'description': 'By (nT), GSM, determined from post-shift GSE components',
   'units': 'nT',
   'type': 'double',
   'fill': '9999.99'},
  {'name': 'BZ_GSM',
   'description': 'Bz (nT), GSM, determined from post-shift GSE components',
   'units': 'nT',
   'type': 'double',
   'fill': '9999.99'},
  {'name': 'flow_speed',
   'description': 'Flow Speed (km/s), GSE',
   'units': 'km/s',
   'type': 'double',
   'fill': '99999.9'}],
 'startDate': '1995-01-01T00:00:00Z',
 'stopDate': '2024-08-31T23:55:00Z',
 'status': {'code': 1200, 'message': 'OK'},
 'x_server': 'https://cdaweb.gsfc.nasa.gov/hapi',
 'x_dataset': 'OMNI_HRO2_5MIN',
 'x_parameters': 'BX_GSE,BY_GSM,BZ_GSM,flow_speed',
 'x_time.min': '2015-03-15T00:00:00Z',
 'x_time.max': '2015-03-20T00:00:00Z',
 'x_requestDate': '2024-09-20T09:54:41',
 'x_cacheDir': '/tmp/hapi-data/cdaweb.gsfc.nasa.gov_hapi',
 'x_downloadTime': 15.295110940933228,
 'x_readTime': 0.0005857944488525391,
 'x_metaFileParsed': '/tmp/hapi-data/cdaweb.gsfc.nasa.gov_hapi/OMNI_HRO2_5MIN___.pkl',
 'x_dataFileParsed': '/tmp/hapi-data/cdaweb.gsfc.nasa.gov_hapi/OMNI_HRO2_5MIN_BX_GSE,BY_GSM,BZ_GSM,flow_speed_20150315T000000_20150320T000000.npy',
 'x_metaFile': '/tmp/hapi-data/cdaweb.gsfc.nasa.gov_hapi/OMNI_HRO2_5MIN___.json',
 'x_dataFile': '/tmp/hapi-data/cdaweb.gsfc.nasa.gov_hapi/OMNI_HRO2_5MIN_BX_GSE,BY_GSM,BZ_GSM,flow_speed_20150315T000000_20150320T000000.bin',
 'x_totalTime': 16.384485483169556}

We are now able to 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_in, hapimeta):
    """Replace bad values (fill values given in metadata) with NaN"""
    hapidata = deepcopy(hapidata_in)
    # HAPI returns metadata for parameters as a list of dictionaries
    # - Loop through them
    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

BX_GSE: 297 fills NaNd
BY_GSM: 297 fills NaNd
BZ_GSM: 297 fills NaNd
flow_speed: 401 fills NaNd

array([(b'2015-03-15T00:00:00.000Z',        nan,         nan,         nan, nan),
       (b'2015-03-15T00:05:00.000Z',        nan,         nan,         nan, nan),
       (b'2015-03-15T00:10:00.000Z',        nan,         nan,         nan, nan),
       ...,
       (b'2015-03-19T23:45:00.000Z', 3.03999996, -4.92000008, -0.06      , nan),
       (b'2015-03-19T23:50:00.000Z', 2.91000009, -4.63999987, -1.02999997, nan),
       (b'2015-03-19T23:55:00.000Z', 3.20000005, -5.13999987, -1.04999995, nan)],
      dtype=[('Time', 'S24'), ('BX_GSE', '<f8'), ('BY_GSM', '<f8'), ('BZ_GSM', '<f8'), ('flow_speed', '<f8')])

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 Python datetime
    df.index = hapitime2datetime(df.index.values.astype(str))
    # df.index = pd.DatetimeIndex(df.index.values.astype(str))
    # Remove timezone awareness
    df.index = df.index.tz_convert("UTC").tz_convert(None)
    # Rename to Timestamp to match viresclient
    df.index.name = "Timestamp"
    return df

df = to_pandas(data)
df

/opt/conda/lib/python3.11/site-packages/hapiclient/hapitime.py:284: UserWarning: The argument 'infer_datetime_format' is deprecated and will be removed in a future version. A strict version of it is now the default, see https://pandas.pydata.org/pdeps/0004-consistent-to-datetime-parsing.html. You can safely remove this argument.
  Time = pandas.to_datetime(Time, infer_datetime_format=True).tz_convert(tzinfo).to_pydatetime()

	BX_GSE	BY_GSM	BZ_GSM	flow_speed
Timestamp
2015-03-15 00:00:00	NaN	NaN	NaN	NaN
2015-03-15 00:05:00	NaN	NaN	NaN	NaN
2015-03-15 00:10:00	NaN	NaN	NaN	NaN
2015-03-15 00:15:00	NaN	NaN	NaN	NaN
2015-03-15 00:20:00	NaN	NaN	NaN	NaN
...	...	...	...	...
2015-03-19 23:35:00	4.30	-3.83	0.61	NaN
2015-03-19 23:40:00	2.26	-5.46	1.24	NaN
2015-03-19 23:45:00	3.04	-4.92	-0.06	NaN
2015-03-19 23:50:00	2.91	-4.64	-1.03	NaN
2015-03-19 23:55:00	3.20	-5.14	-1.05	NaN

1440 rows × 4 columns

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

({'Time': 'UTC',
  'BX_GSE': 'nT',
  'BY_GSM': 'nT',
  'BZ_GSM': 'nT',
  'flow_speed': 'km/s'},
 {'Time': None,
  'BX_GSE': 'Bx (nT), GSE',
  'BY_GSM': 'By (nT), GSM, determined from post-shift GSE components',
  'BZ_GSM': 'Bz (nT), GSM, determined from post-shift GSE components',
  'flow_speed': 'Flow Speed (km/s), GSE'})

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:

Now let’s plot these data:

def plot_solar_wind(ds_sw):
    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\n[nT]")
    axes[0].set_xlabel("")
    ds_sw["flow_speed"].plot.line(
        x="Timestamp", linewidth=1, alpha=0.8, ax=axes[1]
    )
    axes[1].set_ylabel("flow_speed\n[km/s]")
    axes[1].set_xlabel("")
    axes[0].grid()
    axes[1].grid()
    fig.suptitle("Interplanetary Magnetic Field and Solar Wind flow")
    return fig, axes

fig_sw, axes_sw = plot_solar_wind(ds_sw)

../_images/a2e56c37d9f6bb712afa9fe708a551ad3dde06f8f01ced2d97ba5f5cfe7b190c.png

From the solar wind to the ground

Contents

From the solar wind to the ground#

Packages to use#

Time selection#

Solar wind data (OMNI)#

Auroral electrojets as measured by Swarm#

Peak ground magnetic disturbances below satellite tracks#

Ground observatory data (INTERMAGNET)#