Contents

Polar Region Plots

Authors: Ashley Smith

Abstract: Demonstrates access to (FAC, Te, Ne) measurements, and visualisation of them on Lon/Lat and MLT/QDLat plots.

See also:

  • https://github.com/pacesm/jupyter_notebooks/blob/master/AEBS/AEBS_AOB_FAC.ipynb

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

2021-01-24T16:16:42+00:00

CPython 3.7.6
IPython 7.11.1

viresclient 0.7.1
pandas 0.25.3
xarray 0.15.0
matplotlib 3.1.2
cartopy 0.17.0

from viresclient import SwarmRequest
import datetime as dt
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

print("Collection names:")
print(SwarmRequest().available_collections("FAC", details=False))
print(SwarmRequest().available_collections("IPD", details=False), "\n")
print("FAC measurements:\n", SwarmRequest().available_measurements("FAC"))
print("IPD measurements:\n", SwarmRequest().available_measurements("IPD"))

Collection names:
{'FAC': ['SW_OPER_FACATMS_2F', 'SW_OPER_FACBTMS_2F', 'SW_OPER_FACCTMS_2F', 'SW_OPER_FAC_TMS_2F']}
{'IPD': ['SW_OPER_IPDAIRR_2F', 'SW_OPER_IPDBIRR_2F', 'SW_OPER_IPDCIRR_2F']} 

FAC measurements:
 ['IRC', 'IRC_Error', 'FAC', 'FAC_Error', 'Flags', 'Flags_F', 'Flags_B', 'Flags_q']
IPD measurements:
 ['Ne', 'Te', 'Background_Ne', 'Foreground_Ne', 'PCP_flag', 'Grad_Ne_at_100km', 'Grad_Ne_at_50km', 'Grad_Ne_at_20km', 'Grad_Ne_at_PCP_edge', 'ROD', 'RODI10s', 'RODI20s', 'delta_Ne10s', 'delta_Ne20s', 'delta_Ne40s', 'Num_GPS_satellites', 'mVTEC', 'mROT', 'mROTI10s', 'mROTI20s', 'IBI_flag', 'Ionosphere_region_flag', 'IPIR_index', 'Ne_quality_flag', 'TEC_STD']

Fetch data separately from FAC and IPD collections over the same time window

# Set by IS0-8601 time:
start = "2018-05-01T00:00:00Z"
end = "2018-05-01T12:00:00Z"
# Setting by datetime objects:
# start = dt.datetime(2018, 5, 1, 0)
# end = dt.datetime(2018, 5, 1, 12)

request = SwarmRequest()
request.set_collection("SW_OPER_FACATMS_2F")
request.set_products(
    measurements=["FAC"],
    auxiliaries=["MLT", "QDLat", "QDLon"])
data = request.get_between(start, end)
ds_fac = data.as_xarray()
ds_fac

[1/1] Processing:  100%|██████████|  [ Elapsed: 00:01, Remaining: 00:00 ]
      Downloading: 100%|██████████|  [ Elapsed: 00:00, Remaining: 00:00 ] (2.813MB)

<xarray.Dataset>
Dimensions:     (Timestamp: 43200)
Coordinates:
  * Timestamp   (Timestamp) datetime64[ns] 2018-05-01T00:00:00.500000 ... 2018-05-01T11:59:59.500000
Data variables:
    Spacecraft  (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A' 'A'
    QDLat       (Timestamp) float64 -16.37 -16.44 -16.5 ... 76.37 76.42 76.47
    Latitude    (Timestamp) float64 -4.918 -4.983 -5.047 ... 79.76 79.82 79.88
    FAC         (Timestamp) float64 0.006507 -0.006002 ... -0.007964 0.7068
    QDLon       (Timestamp) float64 89.94 89.93 89.92 ... 122.3 122.5 122.6
    MLT         (Timestamp) float64 1.343 1.342 1.342 1.342 ... 15.18 15.19 15.2
    Longitude   (Timestamp) float64 15.96 15.96 15.96 15.96 ... 29.91 30.0 30.09
    Radius      (Timestamp) float64 6.819e+06 6.819e+06 ... 6.806e+06 6.806e+06
Attributes:
    Sources:         ['SW_OPER_FACATMS_2F_20180501T000000_20180501T235959_0301']
    MagneticModels:  []
    RangeFilters:    []
request = SwarmRequest()
request.set_collection("SW_OPER_IPDAIRR_2F")
request.set_products(
    measurements=["Ne", "Te"],
    auxiliaries=["MLT", "QDLat", "QDLon"])
data = request.get_between(start, end)
ds_ipd = data.as_xarray()
ds_ipd

[1/1] Processing:  100%|██████████|  [ Elapsed: 00:01, Remaining: 00:00 ]
      Downloading: 100%|██████████|  [ Elapsed: 00:00, Remaining: 00:00 ] (3.16MB)

<xarray.Dataset>
Dimensions:     (Timestamp: 43200)
Coordinates:
  * Timestamp   (Timestamp) datetime64[ns] 2018-05-01T00:00:00.197000027 ... 2018-05-01T11:59:59.197000027
Data variables:
    Spacecraft  (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A' 'A'
    Te          (Timestamp) float64 1.747e+03 1.752e+03 ... 2.639e+03 2.873e+03
    QDLat       (Timestamp) float64 -16.27 -16.34 -16.4 ... 76.29 76.34 76.39
    Ne          (Timestamp) float64 8.814e+04 8.804e+04 ... 4.562e+04 4.832e+04
    Latitude    (Timestamp) float64 -4.822 -4.886 -4.951 ... 79.66 79.72 79.79
    QDLon       (Timestamp) float64 89.95 89.94 89.93 ... 122.0 122.2 122.4
    MLT         (Timestamp) float64 1.343 1.343 1.343 ... 15.16 15.17 15.18
    Longitude   (Timestamp) float64 15.96 15.96 15.96 ... 29.77 29.86 29.95
    Radius      (Timestamp) float64 6.819e+06 6.819e+06 ... 6.806e+06 6.806e+06
Attributes:
    Sources:         ['SW_OPER_IPDAIRR_2F_20180501T000000_20180501T235959_0301']
    MagneticModels:  []
    RangeFilters:    []

Why are some of the temperatures negative?

ds_ipd["Te"].where(ds_ipd["Te"] < 0, drop=True).head()  # NB only looking at first five entries

<xarray.DataArray 'Te' (Timestamp: 5)>
array([  -1657.87,  -34548.27,   -1694.56, -241857.02,   -7074.39])
Coordinates:
  * Timestamp  (Timestamp) datetime64[ns] 2018-05-01T00:19:36.197000027 ... 2018-05-01T00:24:35.197000027
Attributes:
    units:        K
    description:  Electron temperature, directly copied from the Langmuir pro...

Reorganise the data

Note that the data are recorded at different sampling points:

ds_fac["Timestamp"]

<xarray.DataArray 'Timestamp' (Timestamp: 43200)>
array(['2018-05-01T00:00:00.500000000', '2018-05-01T00:00:01.500000000',
       '2018-05-01T00:00:02.500000000', ..., '2018-05-01T11:59:57.500000000',
       '2018-05-01T11:59:58.500000000', '2018-05-01T11:59:59.500000000'],
      dtype='datetime64[ns]')
Coordinates:
  * Timestamp  (Timestamp) datetime64[ns] 2018-05-01T00:00:00.500000 ... 2018-05-01T11:59:59.500000
ds_ipd["Timestamp"]

<xarray.DataArray 'Timestamp' (Timestamp: 43200)>
array(['2018-05-01T00:00:00.197000027', '2018-05-01T00:00:01.197000027',
       '2018-05-01T00:00:02.197000027', ..., '2018-05-01T11:59:57.197000027',
       '2018-05-01T11:59:58.197000027', '2018-05-01T11:59:59.197000027'],
      dtype='datetime64[ns]')
Coordinates:
  * Timestamp  (Timestamp) datetime64[ns] 2018-05-01T00:00:00.197000027 ... 2018-05-01T11:59:59.197000027

We could have alternatively fetched the data together directly with request.set_collection("SW_OPER_FACATMS_2F", "SW_OPER_IPDAIRR_2F") ... - in this case the server resolves the time discrepancy by using just the sampling times (and rate) of the first collection given (i.e. "SW_OPER_FACATMS_2F") and interpolates the following collections onto the first time series with a nearest-neighbour method. This means that the IPD measurements would falsely be reported at the sampling times of the FAC measurements - this might not be a problem depending on your application, but we will avoid that issue here by accessing them as separate datasets.

We can perform an outer join to merge the datasets into one object, where nan’s fill the empty sampling times in each input dataset. Let’s also set the appropriate variables as coordinates:

ds = ds_ipd.merge(ds_fac, join="outer")
coords = ["Latitude", "Longitude", "Radius", "QDLat", "QDLon", "MLT", "Spacecraft"]
ds = ds.set_coords(coords)
# NB: xarray merge does not handle the attributes
#  so we must merge these manually
ds = ds.assign_attrs({"Sources": ds_fac.Sources + ds_ipd.Sources})
ds

<xarray.Dataset>
Dimensions:     (Timestamp: 86400)
Coordinates:
  * Timestamp   (Timestamp) datetime64[ns] 2018-05-01T00:00:00.197000027 ... 2018-05-01T11:59:59.500000
    Spacecraft  (Timestamp) object 'A' 'A' 'A' 'A' 'A' ... 'A' 'A' 'A' 'A' 'A'
    QDLat       (Timestamp) float64 -16.27 -16.37 -16.34 ... 76.42 76.39 76.47
    Latitude    (Timestamp) float64 -4.822 -4.918 -4.886 ... 79.82 79.79 79.88
    QDLon       (Timestamp) float64 89.95 89.94 89.94 ... 122.5 122.4 122.6
    MLT         (Timestamp) float64 1.343 1.343 1.343 1.342 ... 15.19 15.18 15.2
    Longitude   (Timestamp) float64 15.96 15.96 15.96 15.96 ... 30.0 29.95 30.09
    Radius      (Timestamp) float64 6.819e+06 6.819e+06 ... 6.806e+06 6.806e+06
Data variables:
    Te          (Timestamp) float64 1.747e+03 nan 1.752e+03 ... 2.873e+03 nan
    Ne          (Timestamp) float64 8.814e+04 nan 8.804e+04 ... 4.832e+04 nan
    FAC         (Timestamp) float64 nan 0.006507 nan ... -0.007964 nan 0.7068
Attributes:
    Sources:         ['SW_OPER_FACATMS_2F_20180501T000000_20180501T235959_030...
    MagneticModels:  []
    RangeFilters:    []

We use this single object, ds, to access data in the rest of the notebook.

Visualising the data

Let’s visualise the measurements using the shortcut xarray plotting methods:

ds.plot.scatter(x="Longitude", y="Latitude", hue="FAC", s=0.2, robust=True);

/opt/conda/lib/python3.7/site-packages/matplotlib/colors.py:972: RuntimeWarning: invalid value encountered in subtract
  resdat -= vmin
../_images/05a1_Polar-Region-Plots_19_1.png 
ds.plot.scatter(x="Longitude", y="Latitude", hue="Ne", s=0.2, robust=True);
../_images/05a1_Polar-Region-Plots_20_0.png 
ds.plot.scatter(x="Longitude", y="Latitude", hue="Te", s=0.2, robust=True);
../_images/05a1_Polar-Region-Plots_21_0.png 
ds.plot.scatter(x="MLT", y="QDLat", hue="Te", s=0.2, robust=True);
../_images/05a1_Polar-Region-Plots_22_0.png

Geographic (Lon/Lat) and MLT/QDLat plots

To make more complex figures, it is usually necessary to drop down to the lower level matplotlib interface.

First let’s set up figures onto which data can be plotted:

def _apply_circular_boundary(ax):
    """Make cartopy axes have round borders.
    See https://scitools.org.uk/cartopy/docs/v0.15/examples/always_circular_stereo.html
    
    Notes:
        Inline contour labels are still appearing outside the boundary
    """
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpl.path.Path(verts * radius + center)
    ax.set_boundary(circle, transform=ax.transAxes)

    
def create_axes_geo(title="", figsize=(10, 5)):
    """Create an empty geographic figure with North/South views
    
    Args:
        title (str)
        figsize (tuple): (width, height)
        
    Returns:
        matplotlib.figure.Figure, [GeoAxesSubplot, GeoAxesSubplot]
    """
    fig = plt.figure(figsize=figsize)
    fig.suptitle(title, fontsize=15)
    # Geographic (lat/lon) views
    ax_N_geo = plt.subplot2grid(
        (1, 2), (0, 0),
        projection=ccrs.AzimuthalEquidistant(
            central_longitude=0.0, central_latitude=90.0,
            false_easting=0.0, false_northing=0.0, globe=None
        )
    )
    ax_S_geo = plt.subplot2grid(
        (1, 2), (0, 1),
        projection=ccrs.AzimuthalEquidistant(
            central_longitude=0.0, central_latitude=-90.0,
            false_easting=0.0, false_northing=0.0, globe=None
        )
    )
    ax_N_geo.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())
    ax_S_geo.set_extent([-180, 180, -90, -50], ccrs.PlateCarree())
    for ax in [ax_N_geo, ax_S_geo]:
        _apply_circular_boundary(ax)
        ax.add_feature(cfeature.LAND, facecolor=(1.0, 1.0, 0.9))
        ax.add_feature(cfeature.OCEAN, facecolor=(0.9, 1.0, 1.0))
        ax.add_feature(cfeature.COASTLINE, edgecolor='silver')
    ax_N_geo.gridlines(ylocs=[50, 60, 70, 80, 90])
    ax_S_geo.gridlines(ylocs=[-90, -80, -70, -60, -50])
    ax_N_geo.set_title("North")
    ax_S_geo.set_title("South")
    return fig, [ax_N_geo, ax_S_geo]

def create_axes_mlt(title="", figsize=(10, 5)):
    """Create an empty MLT/QDLat figure with North/South views
    
    Args:
        title (str)
        figsize (tuple): (width, height)
        
    Returns:
        matplotlib.figure.Figure, [PolarAxesSubplot, PolarAxesSubplot]
    """
    fig = plt.figure(figsize=figsize)
    fig.suptitle(title, fontsize=15)
    # QDLat/MLT views
    ax_N_mlt = plt.subplot2grid(
        (1, 2), (0, 0),
        projection="polar"
    )
    ax_S_mlt = plt.subplot2grid(
        (1, 2), (0, 1),
        projection="polar"
    )
    for ax in [ax_N_mlt, ax_S_mlt]:
        ax.set_ylim(0, 40)
        ax.set_yticks([0, 10, 20, 30, 40])
        ax.set_xticklabels(["%2.2i" % (x*12/np.pi) for x in ax.get_xticks()])
        ax.set_theta_zero_location("S")
        ax.set_yticklabels([])
    ax_N_mlt.set_title("North")
    ax_S_mlt.set_title("South")
    return fig, [ax_N_mlt, ax_S_mlt]

These can be used together with the xarray plotting methods which will make some decisions for us automatically (like setting up the colorbars):

fig, axes = create_axes_geo()
for ax in axes:
    ds.plot.scatter(x="Longitude", y="Latitude", hue="Te", s=0.1,
                    ax=ax, transform=ccrs.PlateCarree(), robust=True)
../_images/05a1_Polar-Region-Plots_26_0.png

(the MLT plot case is more complicated because of the transformations needed to plot data onto the polar projections)

To give more control over the plotting, we use matplotlib directly:

Functions to help plot onto the axes specified above

def _subselect(ds, var, subselect_factor):
    _ds = ds.where(~np.isnan(ds[var]), drop=True)
    _ds = _ds.isel({"Timestamp": slice(0, -1, subselect_factor)})
    return _ds

def plot_geo(ax, ds, hemisphere="north", **kwargs):
    """
    kwargs must contain "var" to plot from ds
    
    Args:
        ax (matplotlib.axes)
        ds (xarray.Dataset)
        hemisphere (str): "north" or "south"
    """
    # Identify data variable to plot
    var = kwargs.pop("var")
    # Sub-select data by a given factor
    subselect = kwargs.pop("subselect", None)
    if subselect:
        _ds = _subselect(ds, var, subselect)
    else:
        _ds = ds
    # Select only either NH or SH data
    if hemisphere == "north":
        _ds = _ds.where(ds["Latitude"] > 50, drop=True)
    elif hemisphere == "south":
        _ds = _ds.where(ds["Latitude"] < -50, drop=True)
    # Do the actual plotting
    p = ax.scatter(_ds["Longitude"], _ds["Latitude"], c=_ds[var],
                   transform=ccrs.PlateCarree(), **kwargs)
    return p

def plot_mlt(ax, ds, hemisphere="north", **kwargs):
    """
    kwargs must contain "var" to plot from ds
    
    Args:
        ax (matplotlib.axes)
        ds (xarray.Dataset)
        hemisphere (str): "north" or "south"
    """
    # Identify data variable to plot
    var = kwargs.pop("var")
    # Sub-select data by a given factor
    subselect = kwargs.pop("subselect", None)
    if subselect:
        _ds = _subselect(ds, var, subselect)
    else:
        _ds = ds
    # Select only either NH or SH data
    if hemisphere == "north":
        _ds = _ds.where(ds["QDLat"] > 50, drop=True)
        h = 1
    elif hemisphere == "south":
        _ds = _ds.where(ds["QDLat"] < -50, drop=True)
        h = -1
    # Transformation to the polar coordinates
    def _scatter(x, y, *args, **kwargs):
        return ax.scatter(x*(np.pi/12), 90 - y*h, *args, **kwargs)
    p = _scatter(_ds["MLT"], _ds["QDLat"], c=_ds[var],
                 **kwargs)
    return p

Plot Te on both GEO and MLT views

# Create the bare figures
fig_geo, (ax_N_geo, ax_S_geo) = create_axes_geo()
fig_mlt, (ax_N_mlt, ax_S_mlt) = create_axes_mlt()

# Set the parameters to control the plotting
options = dict(
    var="Te",
    cmap=mpl.cm.plasma,
    norm=mpl.colors.Normalize(vmin=3e3, vmax=6e3),
    s=0.1,
)

# Plot onto each axes
plot_geo(ax_N_geo, ds, hemisphere="north", **options)
plot_geo(ax_S_geo, ds, hemisphere="south", **options)
plot_mlt(ax_N_mlt, ds, hemisphere="north", **options)
plot_mlt(ax_S_mlt, ds, hemisphere="south", **options)

# Add colorbars
for fig in [fig_geo, fig_mlt]:
    cax = fig.add_axes([0.4, 0.15, 0.2, 0.02])
    var = options["var"]
    label = f"{var} [{ds[var].units}]"
    cbar = mpl.colorbar.ColorbarBase(cax,
                                     cmap=options["cmap"],
                                     norm=options["norm"],
                                     orientation='horizontal',
    #                                  ticks=[norm.vmin, (norm.vmax+norm.vmin)/2, norm.vmax],
                                     label=label)
../_images/05a1_Polar-Region-Plots_30_0.png ../_images/05a1_Polar-Region-Plots_30_1.png

Plot both FAC and Ne on one figure

# set up title based on approx start/end times
time1 = ds["Timestamp"][0].values
time2 = ds["Timestamp"][-1].values
def format_dt64(dt64, form="%Y-%m-%d %H:%M"):
    # Convert numpy.datetime64 [ns] -> datetime
    time = dt.datetime.utcfromtimestamp(dt64.astype(int) * 1e-9)
    return time.strftime(form)
spacecraft = ds["Spacecraft"][0].values
title = f"{format_dt64(time1)}\n{format_dt64(time2)}\n(Swarm {spacecraft})"

fig, (ax_N_mlt, ax_S_mlt) = create_axes_mlt(title=title)

options_Ne = dict(
    var="Ne",
    cmap=mpl.cm.viridis,
    norm=mpl.colors.Normalize(vmin=10e3, vmax=100e3),
    s=0.1, zorder=2,
)
# - scatterplot point size (s) and data rate (subselect)
#     need to be balanced to make them legible
# - zorder controls the layers of plots (which one goes on top)
# - Try a horizontal line as the marker
#     (only works when the satellite tracks run vertically)
#     ref: https://matplotlib.org/3.1.0/api/markers_api.html
options_FAC = dict(
    var="FAC", subselect=10,
    cmap=mpl.cm.seismic,
    norm=mpl.colors.SymLogNorm(linthresh=1, linscale=1, vmin=-20,vmax=20),
    s=80, zorder=1, marker="_",
)

for options in [options_FAC, options_Ne]:
    plot_mlt(ax_N_mlt, ds, hemisphere="north", **options)
    plot_mlt(ax_S_mlt, ds, hemisphere="south", **options)

# Set colorbar locations
cax1 = fig.add_axes([0.4, 0.15, 0.2, 0.02])
cax2 = fig.add_axes([0.4, 0.0, 0.2, 0.02])
# Draw colorbars
for cax, options in zip([cax1, cax2], [options_FAC, options_Ne]):
    var = options["var"]
    label = f"{var} [{ds[var].units}]"
    cbar = mpl.colorbar.ColorbarBase(cax,
                                     cmap=options["cmap"],
                                     norm=options["norm"],
                                     orientation='horizontal',
                                     label=label)
../_images/05a1_Polar-Region-Plots_32_0.png

next steps:

  • visualisations of auroral oval boundaries etc, from AEBS products:

  • https://github.com/pacesm/jupyter_notebooks/blob/master/AEBS/AEBS_AOB_FAC.ipynb

  • use cartopy.feature.nightshade to show terminator (for shorter time periods)

  • https://scitools.org.uk/cartopy/docs/latest/gallery/aurora_forecast.html

Ground Observatory Data - VirES