Polar Region Plots

Polar Region Plots#

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

Visualising the data#

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

ds.plot.scatter(x="Longitude", y="Latitude", hue="FAC", s=1, vmin=-0.6, vmax=0.6, cmap="RdBu", linewidths=0);

/opt/conda/lib/python3.11/site-packages/matplotlib/colors.py:1370: RuntimeWarning: invalid value encountered in subtract
  resdat -= vmin

../_images/84b80af6a77c2f5305ae663684a3b0a85358bb55e90b2e88633583ce6aa0889c.png

ds.plot.scatter(x="Longitude", y="Latitude", hue="Ne", s=1, vmin=100_000, vmax=500_000, cmap="viridis", linewidths=0);

ds.plot.scatter(x="Longitude", y="Latitude", hue="Te", s=1, vmin=1500, vmax=5000, cmap="viridis", linewidths=0);

ds.plot.scatter(x="MLT", y="QDLat", hue="Te", s=1, vmin=1500, vmax=5000, cmap="viridis", linewidths=0);

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=1,
                    ax=ax, transform=ccrs.PlateCarree(),
                    vmin=1500, vmax=5000, cmap="viridis", linewidths=0)

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

/tmp/ipykernel_2534/3347982881.py:79: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  ax.set_xticklabels(["%2.2i" % (x*12/np.pi) for x in ax.get_xticks()])
/tmp/ipykernel_2534/3347982881.py:79: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  ax.set_xticklabels(["%2.2i" % (x*12/np.pi) for x in ax.get_xticks()])

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)

/tmp/ipykernel_2534/3347982881.py:79: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  ax.set_xticklabels(["%2.2i" % (x*12/np.pi) for x in ax.get_xticks()])
/tmp/ipykernel_2534/3347982881.py:79: UserWarning: set_ticklabels() should only be used with a fixed number of ticks, i.e. after set_ticks() or using a FixedLocator.
  ax.set_xticklabels(["%2.2i" % (x*12/np.pi) for x in ax.get_xticks()])

next steps:

visualisations of auroral oval boundaries etc, from AEBS products:
pacesm/jupyter_notebooks
use cartopy.feature.nightshade to show terminator (for shorter time periods)
https://scitools.org.uk/cartopy/docs/latest/gallery/aurora_forecast.html

Polar Region Plots

Contents

Polar Region Plots#

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

Why are some of the temperatures negative?#

Reorganise the data#

Visualising the data#

Geographic (Lon/Lat) and MLT/QDLat plots#

Functions to help plot onto the axes specified above#

Plot Te on both GEO and MLT views#

Plot both FAC and Ne on one figure#