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