{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# AEJxLPS (Auroral electrojets SECS)\n", "\n", "> Abstract: Access to the AEBS products, SECS type. This notebook uses code from the previous notebook to build a routine that is flexible to plot either the LC or SECS products - this demonstrates a prototype quicklook routine." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "%load_ext watermark\n", "%watermark -i -v -p viresclient,pandas,xarray,matplotlib" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "from viresclient import SwarmRequest\n", "import datetime as dt\n", "import numpy as np\n", "import pandas as pd\n", "import xarray as xr\n", "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "request = SwarmRequest()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## AEBS product information\n", "\n", "See previous notebook, \"Demo AEBS products (LC)\", for an introduction to these products." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function to request data from VirES and reshape it" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "def fetch_data(start_time=None, end_time=None, spacecraft=None, AEBS_type=\"L\"):\n", " \"\"\"DUPLICATED FROM PREVIOUS NOTEBOOK. TO BE REFACTORED\"\"\"\n", "\n", " # Fetch data from VirES\n", " auxiliaries = ['OrbitNumber', 'QDLat', 'QDOrbitDirection', 'OrbitDirection', 'MLT']\n", " if AEBS_type == \"L\":\n", " measurement_vars = [\"J_NE\"]\n", " elif AEBS_type == \"S\":\n", " measurement_vars = [\"J_CF_NE\", \"J_DF_NE\"]\n", " # Fetch LPL/LPS\n", " request.set_collection(f'SW_OPER_AEJ{spacecraft}LP{AEBS_type}_2F')\n", " request.set_products(\n", " measurements=measurement_vars,\n", " auxiliaries=auxiliaries,\n", " )\n", " data = request.get_between(start_time, end_time, asynchronous=False, show_progress=False)\n", " ds_lp = data.as_xarray()\n", " # Fetch LPL/LPS Quality\n", " request.set_collection(f'SW_OPER_AEJ{spacecraft}LP{AEBS_type}_2F:Quality')\n", " request.set_products(\n", " measurements=['RMS_misfit', 'Confidence'],\n", " )\n", " data = request.get_between(start_time, end_time, asynchronous=False, show_progress=False)\n", " ds_lpq = data.as_xarray()\n", " # Fetch PBL\n", " request.set_collection(f'SW_OPER_AEJ{spacecraft}PB{AEBS_type}_2F')\n", " request.set_products(\n", " measurements=['PointType', 'Flags'],\n", " auxiliaries=auxiliaries\n", " )\n", " data = request.get_between(start_time, end_time, asynchronous=False, show_progress=False)\n", " ds_pb = data.as_xarray()\n", "\n", " # Meaning of PointType\n", " PointType_meanings = {\n", " \"WEJ_peak\": 0, # minimum\n", " \"EEJ_peak\": 1, # maximum\n", " \"WEJ_eq_bound_s\": 2, # equatorward (pair start)\n", " \"EEJ_eq_bound_s\": 3,\n", " \"WEJ_po_bound_s\": 6, # poleward\n", " \"EEJ_po_bound_s\": 7,\n", " \"WEJ_eq_bound_e\": 10, # equatorward (pair end)\n", " \"EEJ_eq_bound_e\": 11,\n", " \"WEJ_po_bound_e\": 14, # poleward\n", " \"EEJ_po_bound_e\": 15,\n", " }\n", " # Add new data variables (boolean Type) according to the dictionary above\n", " ds_pb = ds_pb.assign(\n", " {name: ds_pb[\"PointType\"] == PointType_meanings[name]\n", " for name in PointType_meanings.keys()}\n", " )\n", "\n", " # Merge datasets together\n", " def drop_duplicate_times(_ds):\n", " _, index = np.unique(_ds['Timestamp'], return_index=True)\n", " return _ds.isel(Timestamp=index)\n", " def merge_attrs(_ds1, _ds2):\n", " attrs = {\"Sources\":[], \"MagneticModels\":[], \"AppliedFilters\":[]}\n", " for item in [\"Sources\", \"MagneticModels\", \"AppliedFilters\"]:\n", " attrs[item] = list(set(_ds1.attrs[item] + _ds2.attrs[item]))\n", " return attrs\n", " # Create new dataset from just the newly created PointType arrays\n", " # This is created on a non-repeating Timestamp coordinate\n", " ds = xr.Dataset(\n", " {name: ds_pb[name].where(ds_pb[name], drop=True)\n", " for name in PointType_meanings.keys()}\n", " )\n", " # Merge in the positional and auxiliary data\n", " data_vars = list(set(ds_pb.data_vars).difference(set(PointType_meanings.keys())))\n", " data_vars.remove(\"PointType\")\n", " ds = ds.merge(\n", " (ds_pb[data_vars]\n", " .pipe(drop_duplicate_times))\n", " )\n", " # Merge together with the LPL data\n", " # Note that the Timestamp coordinates aren't equal\n", "\n", " # Separately merge data with matching and missing time sample points in ds_lpl\n", " idx_present = list(set(ds[\"Timestamp\"].values).intersection(set(ds_lp[\"Timestamp\"].values)))\n", " idx_missing = list(set(ds[\"Timestamp\"].values).difference(set(ds_lp[\"Timestamp\"].values)))\n", " # Override prioritises the first dataset (ds_lpl) where there are conflicts\n", " ds2 = ds_lp.merge(ds.sel(Timestamp=idx_present), join=\"outer\", compat=\"override\")\n", " ds2 = ds2.merge(ds.sel(Timestamp=idx_missing), join=\"outer\")\n", " # Update the metadata\n", " ds2.attrs = merge_attrs(ds_lp, ds_pb)\n", "\n", " # Switch the point type arrays to uint8 or bool for performance?\n", " # But the .where operations later cast them back to float64 since gaps are filled with nan\n", " for name in PointType_meanings.keys():\n", " ds2[name] = ds2[name].astype(\"uint8\").fillna(False)\n", " # ds2[name] = ds2[name].fillna(False).astype(bool)\n", "\n", " ds = ds2\n", "\n", " # Append the PBL Flags information into the LPL:Quality dataset to use as a lookup table\n", " ds_lpq = ds_lpq.assign(\n", " Flags_PBL=\n", " ds_pb[\"Flags\"]\n", " .pipe(drop_duplicate_times)\n", " .reindex_like(ds_lpq, method=\"nearest\"),\n", " )\n", "\n", " return ds, ds_lpq\n", "\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting function" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Bit numbers which indicate non-nominal state\n", "# Check SW-DS-DTU-GS-003_AEBS_PDD for details\n", "BITS_PBL_FLAGS_EEJ_MINOR = (2, 3, 6)\n", "BITS_PBL_FLAGS_WEJ_MINOR = (4, 5, 6)\n", "BITS_PBL_FLAGS_EEJ_BAD = (0, 7, 8, 11)\n", "BITS_PBL_FLAGS_WEJ_BAD = (1, 9, 10, 12)\n", "\n", "def check_PBL_Flags(flags=0b0, EJ_type=\"WEJ\"):\n", " \"\"\"Return \"good\", \"poor\", or \"bad\" depending on status\"\"\"\n", " def _check_bits(bitno_set):\n", " return any(flags & (1 << bitno) for bitno in bitno_set)\n", " if EJ_type == \"WEJ\":\n", " if _check_bits(BITS_PBL_FLAGS_WEJ_BAD):\n", " return \"bad\"\n", " elif _check_bits(BITS_PBL_FLAGS_WEJ_MINOR):\n", " return \"poor\"\n", " else:\n", " return \"good\"\n", " elif EJ_type == \"EEJ\":\n", " if _check_bits(BITS_PBL_FLAGS_EEJ_BAD):\n", " return \"bad\"\n", " elif _check_bits(BITS_PBL_FLAGS_EEJ_MINOR):\n", " return \"poor\"\n", " else:\n", " return \"good\"\n", "\n", "glyphs = {\n", " \"WEJ_peak\": {\"marker\": 'v', \"color\":'tab:red'}, # minimum\n", " \"EEJ_peak\": {\"marker\": '^', \"color\":'tab:purple'}, # maximum\n", " \"WEJ_eq_bound_s\": {\"marker\": '>', \"color\":'black'}, # equatorward (pair start)\n", " \"EEJ_eq_bound_s\": {\"marker\": '>', \"color\":'black'},\n", " \"WEJ_po_bound_s\": {\"marker\": '>', \"color\":'black'}, # poleward\n", " \"EEJ_po_bound_s\": {\"marker\": '>', \"color\":'black'},\n", " \"WEJ_eq_bound_e\": {\"marker\": '<', \"color\":'black'}, # equatorward (pair end)\n", " \"EEJ_eq_bound_e\": {\"marker\": '<', \"color\":'black'},\n", " \"WEJ_po_bound_e\": {\"marker\": '<', \"color\":'black'}, # poleward\n", " \"EEJ_po_bound_e\": {\"marker\": '<', \"color\":'black'},\n", "}\n", "\n", "\n", "def plot_stack(ds, ds_lpq, hemisphere=\"North\", x_axis=\"Latitude\", AEBS_type=\"L\"):\n", " # Identify which variable to plot from dataset\n", " # If accessing the SECS (LPS) data, sum the DF & CF parts\n", " if \"J_CF_NE\" in ds.data_vars:\n", " ds[\"J_NE\"] = ds[\"J_DF_NE\"] + ds[\"J_CF_NE\"]\n", " plotvar = \"J_NE\"\n", " orbdir = \"OrbitDirection\" if x_axis==\"Latitude\" else \"QDOrbitDirection\"\n", " markersize = 1 if AEBS_type==\"S\" else 5\n", " # Select hemisphere\n", " if hemisphere == \"North\":\n", " ds = ds.where(ds[\"Latitude\"]>0, drop=True)\n", " elif hemisphere == \"South\":\n", " ds = ds.where(ds[\"Latitude\"]<0, drop=True)\n", " # Generate plot with split by columns: ascending/descending to/from pole\n", " # by rows: successive orbits\n", " # Skip when no data available\n", " if len(ds[\"Timestamp\"]) == 0:\n", " return None, None\n", " fig, axes = plt.subplots(\n", " nrows=len(ds.groupby(\"OrbitNumber\")), ncols=2, sharex=\"col\", sharey=\"all\",\n", " figsize=(10, 20)\n", " )\n", " max_ylim = np.max(np.abs(ds[plotvar].sel({\"NE\": \"E\"})))\n", " # Loop through each orbit\n", " for i, (_, ds_orbit) in enumerate(ds.groupby(\"OrbitNumber\")):\n", " if hemisphere == \"North\":\n", " ds_orb_asc = ds_orbit.where(ds_orbit[orbdir] == 1, drop=True)\n", " ds_orb_desc = ds_orbit.where(ds_orbit[orbdir] == -1, drop=True)\n", " if hemisphere == \"South\":\n", " ds_orb_asc = ds_orbit.where(ds_orbit[orbdir] == -1, drop=True)\n", " ds_orb_desc = ds_orbit.where(ds_orbit[orbdir] == 1, drop=True)\n", " # Loop through ascending and descending sections\n", " for j, _ds in enumerate((ds_orb_asc, ds_orb_desc)):\n", " if len(_ds.Timestamp) == 0:\n", " continue\n", " # Line plot of current strength\n", " axes[i, j].plot(\n", " _ds[x_axis], _ds[plotvar].sel({\"NE\": \"E\"}),\n", " color=\"tab:blue\", marker=\".\", markersize=markersize, linestyle=\"\"\n", " )\n", " axes[i, j].plot(\n", " _ds[x_axis], _ds[plotvar].sel({\"NE\": \"N\"}),\n", " color=\"tab:grey\", marker=\".\", markersize=markersize, linestyle=\"\"\n", " )\n", " # Plot glyphs at the peaks and boundaries locations\n", " for name in glyphs.keys():\n", " __ds = _ds.where(_ds[name], drop=True)\n", " try:\n", " for lat in __ds[x_axis]:\n", " axes[i, j].plot(\n", " lat, 0,\n", " marker=glyphs[name][\"marker\"], color=glyphs[name][\"color\"]\n", " )\n", " except Exception:\n", " pass\n", " # Identify Quality and Flags info\n", " # Use either the start time of the section or the end, depending on asc or desc\n", " index = 0 if j == 0 else -1\n", " t = _ds[\"Timestamp\"].isel(Timestamp=index).values\n", " _ds_qualflags = ds_lpq.sel(Timestamp=t, method=\"nearest\")\n", " pbl_flags = int(_ds_qualflags[\"Flags_PBL\"].values)\n", " lpl_rms_misfit = float(_ds_qualflags[\"RMS_misfit\"].values)\n", " lpl_confidence = float(_ds_qualflags[\"Confidence\"].values)\n", " # Shade WEJ and EEJ regions, only if well-defined\n", " # def _shade_EJ_region(_ds=None, EJ=\"WEJ\", color=\"tab:red\", alpha=0.3):\n", " wej_status = check_PBL_Flags(pbl_flags, \"WEJ\")\n", " eej_status = check_PBL_Flags(pbl_flags, \"EEJ\")\n", " if wej_status in [\"good\", \"poor\"]:\n", " alpha = 0.3 if wej_status == \"good\" else 0.1\n", " try:\n", " WEJ_left = _ds.where(\n", " (_ds[\"WEJ_eq_bound_s\"] == 1) | (_ds[\"WEJ_po_bound_s\"] == 1), drop=True)\n", " WEJ_right = _ds.where(\n", " (_ds[\"WEJ_eq_bound_e\"] == 1) | (_ds[\"WEJ_po_bound_e\"] == 1), drop=True)\n", " x1 = WEJ_left[x_axis][0]\n", " x2 = WEJ_right[x_axis][0]\n", " axes[i, j].fill_betweenx(\n", " [-max_ylim, max_ylim], [x1, x1], [x2, x2], color=\"tab:red\", alpha=alpha)\n", " except Exception:\n", " pass\n", " if eej_status in [\"good\", \"poor\"]:\n", " alpha = 0.3 if eej_status == \"good\" else 0.15\n", " try:\n", " EEJ_left = _ds.where(\n", " (_ds[\"EEJ_eq_bound_s\"] == 1) | (_ds[\"EEJ_po_bound_s\"] == 1), drop=True)\n", " EEJ_right = _ds.where(\n", " (_ds[\"EEJ_eq_bound_e\"] == 1) | (_ds[\"EEJ_po_bound_e\"] == 1), drop=True)\n", " x1 = EEJ_left[x_axis][0]\n", " x2 = EEJ_right[x_axis][0]\n", " axes[i, j].fill_betweenx(\n", " [-max_ylim, max_ylim], [x1, x1], [x2, x2], color=\"tab:purple\", alpha=alpha)\n", " except Exception:\n", " pass\n", " # Write the LPL:Quality and PBL Flags info\n", " ha = \"right\" if j == 0 else \"left\"\n", " textx = 0.98 if j == 0 else 0.02\n", " axes[i, j].text(\n", " textx, 0.95,\n", " f\"RMS Misfit {np.round(lpl_rms_misfit, 2)}; Confidence {np.round(lpl_confidence, 2)}\",\n", " transform=axes[i, j].transAxes, verticalalignment=\"top\", horizontalalignment=ha\n", " )\n", " axes[i, j].text(\n", " textx, 0.05,\n", " f\"PBL Flags {pbl_flags:013b}\",\n", " transform=axes[i, j].transAxes, verticalalignment=\"bottom\", horizontalalignment=ha\n", " )\n", " # Write the start/end time and MLT of the section, and the orbit number\n", " def _format_utc(t):\n", " return f\"UTC {t.strftime('%H:%M')}\"\n", " def _format_mlt(mlt):\n", " hour, fraction = divmod(mlt, 1)\n", " t = dt.time(int(hour), minute=int(60*fraction))\n", " return f\"MLT {t.strftime('%H:%M')}\"\n", " try:\n", " # Left part (section starting UTC, MLT, OrbitNumber)\n", " time_s = pd.to_datetime(ds_orb_asc[\"Timestamp\"].isel(Timestamp=0).data)\n", " mlt_s = ds_orb_asc[\"MLT\"].dropna(dim=\"Timestamp\").isel(Timestamp=0).data\n", " orbit_number = int(ds_orb_asc[\"OrbitNumber\"].isel(Timestamp=0).data)\n", " axes[i, 0].text(\n", " 0.01, 0.95, f\"{_format_utc(time_s)}\\n{_format_mlt(mlt_s)}\",\n", " transform=axes[i, 0].transAxes, verticalalignment=\"top\"\n", " )\n", " axes[i, 0].text(\n", " 0.01, 0.05, f\"Orbit {orbit_number}\",\n", " transform=axes[i, 0].transAxes, verticalalignment=\"bottom\"\n", " )\n", " except Exception:\n", " pass\n", " try:\n", " # Right part (section ending UTC, MLT)\n", " time_e = pd.to_datetime(ds_orb_desc[\"Timestamp\"].isel(Timestamp=-1).data)\n", " mlt_e = ds_orb_desc[\"MLT\"].dropna(dim=\"Timestamp\").isel(Timestamp=-1).data\n", " axes[i, 1].text(\n", " 0.99, 0.95, f\"{_format_utc(time_e)}\\n{_format_mlt(mlt_e)}\",\n", " transform=axes[i, 1].transAxes, verticalalignment=\"top\", horizontalalignment=\"right\"\n", " )\n", " except Exception:\n", " pass\n", " # Extra config of axes and figure text\n", " axes[0, 0].set_ylim(-max_ylim, max_ylim)\n", " if hemisphere == \"North\":\n", " axes[0, 0].set_xlim(50, 90)\n", " axes[0, 1].set_xlim(90, 50)\n", " elif hemisphere == \"South\":\n", " axes[0, 0].set_xlim(-50, -90)\n", " axes[0, 1].set_xlim(-90, -50)\n", " for ax in axes.flatten():\n", " ax.grid()\n", " axes[-1, 0].set_xlabel(x_axis)\n", " axes[-1, 0].set_ylabel(\"Horizontal currents\\n[ A.km$^{-1}$ ]\")\n", " time = pd.to_datetime(ds[\"Timestamp\"].isel(Timestamp=0).data)\n", " spacecraft = ds[\"Spacecraft\"].dropna(dim=\"Timestamp\").isel(Timestamp=0).data\n", " AEBS_type_name = \"LC\" if AEBS_type == \"L\" else \"SECS\"\n", " fig.text(\n", " 0.5, 0.9, f\"{time.strftime('%Y-%m-%d')}\\nSwarm {spacecraft}\\n{hemisphere}\\nAEBS: {AEBS_type_name}\",\n", " transform=fig.transFigure, horizontalalignment=\"center\",\n", " )\n", " fig.subplots_adjust(wspace=0, hspace=0)\n", " return fig, axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Fetching and plotting function\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def quicklook(day=\"2015-01-01\", hemisphere=\"North\", spacecraft=\"A\", AEBS_type=\"L\", xaxis=\"Latitude\"):\n", " start_time = dt.datetime.fromisoformat(day)\n", " end_time = start_time + dt.timedelta(days=1)\n", " ds, ds_lpq = fetch_data(start_time, end_time, spacecraft, AEBS_type)\n", " fig, axes = plot_stack(ds, ds_lpq, hemisphere, xaxis, AEBS_type)\n", " return ds, fig, axes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Consecutive orbits are shown in consecutive rows, centered over the pole. The starting and ending times (UTC and MLT) of the orbital section are shown at the left and right. Westward (WEJ) and Eastward (EEJ) electrojet extents and peak intensities are indicated:\n", "- Blue dots: Estimated current density in Eastward direction, J_NE (E)\n", "- Grey dots: Estimated current density in Northward direction, J_NE (N)\n", "- Red/Purple shaded region: WEJ/EEJ extent (boundaries marked by black triangles)\n", "- Red/Purple triangles: Locations of peak WEJ/EEJ intensity\n", "\n", "Select AEBS_type as S to get SECS results, L to get LC results \n", "SECS = spherical elementary current systems method \n", "LC = Line current method\n", "\n", "Notes: \n", "The code is currently quite fragile, so it is broken on some days. Sometimes the electrojet regions are not shaded correctly. Only the horizontal currents are currently shown." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "quicklook(day=\"2016-01-01\", hemisphere=\"North\", spacecraft=\"A\", AEBS_type=\"S\", xaxis=\"Latitude\");" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "quicklook(day=\"2016-01-01\", hemisphere=\"North\", spacecraft=\"A\", AEBS_type=\"L\", xaxis=\"Latitude\");" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3" } }, "nbformat": 4, "nbformat_minor": 4 }