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
start_time = '2015-03-15T00:00:00Z'
end_time = '2015-03-20T00:00:00Z'
print(hapi.__doc__)
Request data from a HAPI server.
For additional documentation and demonstration, see
<https://github.com/hapi-server/client-python-notebooks/blob/master/hapi_demo.ipynb>
Version: 0.1.4
Parameters
----------
server : str
A string with the URL to a HAPI compliant server. (A HAPI URL
always ends with "/hapi").
dataset : str
A string specifying a dataset from a server
parameters: str
A Comma-separated list of parameters in dataset
start: str
The start time of the requested data
stop: str
The end time of the requested data; end times are exclusive - the
last data record returned by a HAPI server should have a timestamp
before end_time.
options : dict
The following options are available.
logging (False) - Log to console
cache (True) - Save responses and processed responses in cachedir
cachedir (./hapi-data)
usecache (True) - Use files in cachedir if found
serverlist (https://github.com/hapi-server/servers/raw/master/all.txt)
Returns
-------
result : various
Results depend on the input parameters.
Servers = hapi() returns a list of available HAPI server URLs from
<https://github.com/hapi-server/data-specification/blob/master/all.txt>
Dataset = hapi(Server) returns a dict of datasets available from a
URL given by the string Server. The dictionary structure follows the
HAPI JSON structure.
Parameters = hapi(Server, Dataset) returns a dictionary of parameters
in the string Dataset. The dictionary structure follows the HAPI JSON
structure.
Metadata = hapi(Server, Dataset, Parameters) returns metadata
associated each parameter in the comma-separated string Parameters. The
dictionary structure follows the HAPI JSON structure.
Data = hapi(Server, Dataset, Parameters, Start, Stop) returns a
dictionary with elements corresponding to Parameters, e.g., if
Parameters = 'scalar,vector' and the number of records in the time
range Start <= t < Stop returned is N, then
Data['scalar'] is a NumPy array of shape (N)
Data['vector'] is a NumPy array of shape (N,3)
Data['Time'] is a NumPy array of byte literals with shape (N).
Byte literal times can be converted to Python datetimes using
dtarray = hapitime2datetime(Data['Time'])
Data, Meta = hapi(Server, Dataset, Parameters, Start, Stop) returns
the metadata for parameters in Meta.
References
----------
* `HAPI Server Definition <https://github.com/hapi-server/data-specification>`
Examples
----------
See <https://github.com/hapi-server/client-python-notebooks>
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
meta
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
BX_GSE: 0 fills NaNd
BY_GSM: 0 fills NaNd
BZ_GSM: 0 fills NaNd
flow_speed: 0 fills NaNd
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)
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
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
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]
)
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
[1/1] Processing: 100%|██████████| [ Elapsed: 00:01, Remaining: 00:00 ]
Downloading: 100%|██████████| [ Elapsed: 00:00, Remaining: 00:00 ] (0.072MB)
# 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();
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()
[1/1] Processing: 100%|██████████| [ Elapsed: 00:01, Remaining: 00:00 ]
Downloading: 100%|██████████| [ Elapsed: 00:00, Remaining: 00:00 ] (0.056MB)
ds
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)
<xarray.Dataset>
Dimensions: (NE: 2, Timestamp: 613)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2015-03-15T00:00:15 ... 2015-03-1...
* NE (NE) <U1 'N' 'E'
Data variables:
Spacecraft (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A' 'A'
B_NE (Timestamp, NE) float64 23.81 -67.72 85.51 ... 28.03 -53.16
Longitude (Timestamp) float64 95.92 115.1 113.7 ... 125.2 -82.54 -144.0
OrbitNumber (Timestamp) int32 7314 7314 7314 7314 ... 7390 7390 7390 7390
Latitude (Timestamp) float64 83.5 -63.85 -55.77 ... -53.96 -81.24 -86.87
Attributes:
Sources: ['SW_OPER_AEJAPBS_2F_20150101T000000_20151231T235959_010...
MagneticModels: []
RangeFilters: []
# 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')
<xarray.Dataset>
Dimensions: (NE: 2, Timestamp: 80)
Coordinates:
* Timestamp (Timestamp) datetime64[ns] 2015-03-15 ... 2015-03-19T22:30:00
* NE (NE) <U1 'N' 'E'
Data variables:
B_NE (Timestamp, NE) float64 27.91 -19.99 22.15 ... 31.37 -40.68
Longitude (Timestamp) float64 -2.764 4.562 -22.25 ... 113.3 81.96 20.98
OrbitNumber (Timestamp) float64 7.314e+03 7.315e+03 ... 7.389e+03 7.39e+03
Latitude (Timestamp) float64 -9.733 1.443 -10.27 ... -12.09 1.797 2.364
def fetch_ground_obs(IAGA_code, start_time, end_time):
request = SwarmRequest()
request.set_collection(f"SW_OPER_AUX_OBSM2_:{IAGA_code}")
request.set_products(
measurements=["B_NEC"]
)
data = request.get_between(start_time, end_time)
ds = data.as_xarray()
ds.attrs["Latitude_GEO"] = ds["Latitude"].values[0]
ds.attrs["Longitude_GEO"] = ds["Longitude"].values[0]
ds.attrs["Radius_GEO"] = ds["Radius"].values[0]
ds = ds.drop_vars(["Spacecraft", "Latitude", "Longitude", "Radius"])
return ds
ds_ESK = fetch_ground_obs("ESK", start_time, end_time)
ds_ESK
[1/1] Processing: 0%| | [ Elapsed: 00:00, Remaining: ? ]
Accessing INTERMAGNET and/or WDC data
Check usage terms at ftp://ftp.nerc-murchison.ac.uk/geomag/Swarm/AUX_OBS/minute/README
[1/1] Processing: 100%|██████████| [ Elapsed: 00:01, Remaining: 00:00 ]
Downloading: 100%|██████████| [ Elapsed: 00:00, Remaining: 00:00 ] (0.417MB)
fig, axes = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=(15, 5))
for i, NEC_var in enumerate("NEC"):
ds_ESK["B_NEC"].sel(NEC=NEC_var).plot.line(ax=axes[i])
axes[i].set_xlabel("")
axes[i].set_ylabel(f"B_{NEC_var} [nT]")
axes[i].set_title("")
fig.suptitle("ESK")