Abstract: Access to the field aligned currents evaluated by the dual satellite method (level 2 product). We also show an orbit-by-orbit plot using a periodic axis to display (centred over both poles) an overview of the FAC estimates over two weeks.
%load_ext watermark
%watermark -i -v -p viresclient,pandas,xarray,matplotlib
from viresclient import SwarmRequest
import datetime as dt
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
request = SwarmRequest()
This is derived from data from both Swarm Alpha and Charlie by the Ampère's integral method
Documentation:
NB: these are the same as in the FACxTMS_2F
single-satellite FAC product
request.available_collections("FAC", details=False)
request.available_measurements("FAC")
Also fetch the quasidipole (QD) coordinates at the same time.
request.set_collection("SW_OPER_FAC_TMS_2F")
request.set_products(
measurements=["FAC", "FAC_Error",
"Flags", "Flags_F", "Flags_B", "Flags_q"],
auxiliaries=["QDLat", "QDLon"],
)
data = request.get_between(
dt.datetime(2016,1,1),
dt.datetime(2016,1,2)
)
data.sources
Load as a pandas dataframe:
df = data.as_dataframe()
df.head()
Depending on your application, you should probably do some filtering according to each of the flags. This can be done on the dataframe here, or beforehand on the server using request.set_range_filter()
. See https://earth.esa.int/documents/10174/1514862/Swarm-L2-FAC-Dual-Product-Description for more about the data
Load as xarray dataset (we will use this instead of the pandas dataframe)
ds = data.as_xarray()
ds
ds["FAC"].plot(x="Timestamp");
fig, axes = plt.subplots(ncols=1, nrows=2, figsize=(15,5))
axes[0].plot(ds["Timestamp"], ds["FAC"])
axes[1].plot(ds["Timestamp"], ds["FAC_Error"], color="orange")
axes[0].set_ylabel("FAC\n[$\mu A / m^2$]");
axes[1].set_ylabel("Error\n[$\mu A / m^2$]");
axes[1].set_xlabel("Timestamp");
date_format = mdates.DateFormatter('%Y-%m-%d\n%H:%M')
axes[1].xaxis.set_major_formatter(date_format)
axes[0].set_ylim(-5, 5);
axes[1].set_ylim(0, 1);
axes[0].set_xticklabels([])
axes[0].grid(True)
axes[1].grid(True)
fig.subplots_adjust(hspace=0.1)
Note that the errors are lower than with the single satellite product
# Fetch Dst over this period
# NB this is just a convenient way to fetch Dst in this notebook
# Not really a good way to do it in general
start_time = dt.datetime(2016,1,1)
end_time = dt.datetime(2018,1,1)
request = SwarmRequest()
request.set_collection("SW_OPER_FAC_TMS_2F")
request.set_products(
measurements=[],
auxiliaries=["Dst"],
sampling_step="PT120M"
)
ds_dst = request.get_between(start_time, end_time).as_xarray()
# Count the occurences of Dst < -50 nT each month
ds_dst["Dst"].where(ds_dst["Dst"]<-50).resample({"Timestamp": "1m"}).count().plot();
# Plot Dst from that peak month to then choose a two-week period from it
ds_dst["Dst"].sel({"Timestamp": slice("2016-10-01","2016-10-30")}).plot();
start_time = dt.datetime(2016,10,9)
## If you want an exact number of days:
end_time = start_time + dt.timedelta(days=14)
request = SwarmRequest()
request.set_collection("SW_OPER_FAC_TMS_2F")
request.set_products(
measurements=["FAC", "Flags"],
auxiliaries=["QDLat", "QDLon", "MLT",
"OrbitNumber", "QDOrbitDirection"],
sampling_step="PT10S"
)
data = request.get_between(start_time, end_time)
ds = data.as_xarray()
# NB OrbitNumber is not available for the dual satellite product because it is ambiguous
ds
# Functions to produce periodic axes
# Source: https://github.com/pacesm/jupyter_notebooks/blob/master/Periodic%20Axis.ipynb
from matplotlib import pyplot as plt
from numpy import mod, arange, floor_divide, asarray, concatenate, empty, array
from itertools import chain
import matplotlib.cm as cm
from matplotlib.colors import Normalize, LogNorm, SymLogNorm
from matplotlib.colorbar import ColorbarBase
class PeriodicAxis(object):
def __init__(self, period=1.0, offset=0):
self.period = period
self.offset = offset
def _period_index(self, x):
return floor_divide(x - self.offset, self.period)
def periods(self, offset, size):
return self.period * arange(self._period_index(offset), self._period_index(offset + size) + 1)
def normalize(self, x):
return mod(x - self.offset, self.period) + self.offset
# def periodic_plot(pax, x, y, xmin, xmax, *args, **kwargs):
# xx = pax.normalize(x)
# for period in pax.periods(xmin, xmax - xmin):
# plt.plot(xx + period, y, *args, **kwargs)
# plt.xlim(xmin, xmax)
# def periodic_xticks(pax, xmin, xmax, ticks, labels=None):
# ticks = asarray(ticks)
# labels = labels or ticks
# ticks_locations = concatenate([
# ticks + period
# for period in pax.periods(xmin, xmax - xmin)
# ])
# ticks_labels = list(chain.from_iterable(
# labels for _ in pax.periods(xmin, xmax - xmin)
# ))
# plt.xticks(ticks_locations, ticks_labels)
# plt.xlim(xmin, xmax)
class PeriodicLatitudeAxis(PeriodicAxis):
def __init__(self, period=360, offset=-270):
super().__init__(period, offset)
def periods(self, offset, size):
return self.period * arange(self._period_index(offset), self._period_index(offset + size) + 1)
def normalize(self, x):
return mod(x - self.offset, self.period) + self.offset
def periodic_latscatter(pax, x, y, ymin, ymax, *args, **kwargs):
yy = pax.normalize(y)
xmin = kwargs.pop("xmin")
xmax = kwargs.pop("xmax")
for period in pax.periods(xmin, xmax - xmin):
plt.scatter(x, yy + period, *args, **kwargs)
plt.ylim(ymin, ymax)
def periodic_yticks(pax, ymin, ymax, ticks, labels=None):
ticks = asarray(ticks)
labels = labels or ticks
ticks_locations = concatenate([
ticks + period
for period in pax.periods(ymin, ymax - ymin)
])
ticks_labels = list(chain.from_iterable(
labels for _ in pax.periods(ymin, ymax - ymin)
))
plt.yticks(ticks_locations, ticks_labels)
plt.ylim(ymin, ymax)
def fac_overview_plot(ds):
fig, axes = plt.subplots(figsize=(16, 8), dpi=100)
tmp1 = array(ds['QDLat'][1:])
tmp0 = array(ds['QDLat'][:-1])
pass_flag = tmp1 - tmp0
latitudes = tmp0
times = array(ds['Timestamp'][:-1])
values = array(ds['FAC'][:-1])
# # Subsample data
# pass_flag = pass_flag[::10]
# latitudes = latitudes[::10]
# times = times[::10]
# values = values[::10]
plotted_latitudes = array(latitudes)
descending = pass_flag < 0
plotted_latitudes[descending] = -180 - latitudes[descending]
vmax = 20
plax = PeriodicLatitudeAxis()
periodic_latscatter(
plax, times, plotted_latitudes, -190, 190,
xmin=-360-90, xmax=360+90, c=values, s=1,
# https://matplotlib.org/tutorials/colors/colormapnorms.html#symmetric-logarithmic
cmap=cm.coolwarm, norm=SymLogNorm(linthresh=0.1, linscale=1, vmin=-vmax,vmax=vmax)
)
cax = plt.colorbar()
cax.ax.set_ylabel("FAC [$\mu A / m^2$]")
plt.xlim(times.min(), times.max())
plt.xlabel("time")
plt.ylabel("QD-latitude / deg")
periodic_yticks(plax, -190, +190, [-225, -180, -135, -90, -45, 0, 45, 90], labels=[
'+45\u2193', '0\u2193', '\u221245\u2193', '\u221290', '\u221245\u2191', '0\u2191', '+45\u2191', '+90'
])
return fig, axes
fac_overview_plot(ds);
# Maximum values for FAC indicate saturation of the colorbar
ds["FAC"].max(), ds["FAC"].min()
# Flags indicate something wrong around 2016-10-19
fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(20,5), dpi=200)
ds["FAC"].plot(ax=axes[0])
ds["Flags"].plot(ax=axes[1]);
ds.plot.scatter(y="FAC", x="QDLat", s=1, linewidths=0);
fig, ax = plt.subplots(figsize=(15,5), dpi=200)
ds.plot.scatter(
ax=ax,
x="MLT", y="QDLat", hue="FAC", cmap=cm.coolwarm, s=1, linewidths=0,
norm=SymLogNorm(linthresh=0.1, linscale=1, vmin=-20, vmax=20),
);
# Plotting the maximum FAC encountered on each orbit
# Do this with the single satellite product instead for now
# so that we can quickly use the OrbitNumber parameter
# which is not available for the dual-sat product
# You could achieve this by generating a flag based on
# zero-crossings of Latitude
start_time = dt.datetime(2016,10,9)
end_time = start_time + dt.timedelta(days=14)
request = SwarmRequest()
request.set_collection("SW_OPER_FACATMS_2F")
request.set_products(
measurements=["FAC", "Flags"],
auxiliaries=["QDLat", "QDLon", "MLT",
"OrbitNumber", "QDOrbitDirection"],
sampling_step="PT10S"
)
data = request.get_between(start_time, end_time)
ds = data.as_xarray()
fig, axes = plt.subplots(nrows=2, sharex=True)
ds.groupby("OrbitNumber").max()["FAC"].plot(ax=axes[0])
ds.groupby("OrbitNumber").min()["FAC"].plot(ax=axes[1]);