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__)

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

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

# 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()

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)

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')

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

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")