{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eddy detection : Med\n\nScript will detect eddies on adt field, and compute u,v with method add_uv(which could use, only if equator is avoid)\n\nFigures will show different step to detect eddies.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\n\nfrom matplotlib import pyplot as plt\nfrom numpy import arange\n\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def start_axes(title):\n    fig = plt.figure(figsize=(13, 5))\n    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n    ax.set_aspect(\"equal\")\n    ax.set_title(title, weight=\"bold\")\n    return ax\n\n\ndef update_axes(ax, mappable=None):\n    ax.grid()\n    if mappable:\n        plt.colorbar(mappable, cax=ax.figure.add_axes([0.94, 0.05, 0.01, 0.9]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load Input grid, ADT is used to detect eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g = RegularGridDataset(\n    data.get_demo_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"),\n    \"longitude\",\n    \"latitude\",\n)\n\nax = start_axes(\"ADT (m)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15, cmap=\"RdBu_r\")\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Get geostrophic speed u,v\nU/V are deduced from ADT, this algortihm is not ok near the equator (~+- 2\u00b0)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g.add_uv(\"adt\")\nax = start_axes(\"U/V deduce from ADT (m)\")\nax.set_xlim(2.5, 9), ax.set_ylim(37.5, 40)\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15, cmap=\"RdBu_r\")\nu, v = g.grid(\"u\").T, g.grid(\"v\").T\nax.quiver(g.x_c, g.y_c, u, v, scale=10)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Pre-processings\nApply a high-pass filter to remove the large scale and highlight the mesoscale\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "g.bessel_high_filter(\"adt\", 500)\nax = start_axes(\"ADT (m) filtered (500km)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15, cmap=\"RdBu_r\")\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Identification\nRun the identification step with slices of 2 mm\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "date = datetime(2016, 5, 15)\na, c = g.eddy_identification(\"adt\", \"u\", \"v\", date, 0.002, shape_error=55)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display of all closed contours found in the grid (only 1 contour every 4)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT closed contours (only 1 / 4 levels)\")\ng.contours.display(ax, step=4)\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Contours included in eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT contours used as eddies\")\ng.contours.display(ax, only_used=True)\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Post analysis\nContours can be rejected for several reasons (shape error to high, several extremum in contour, ...)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT rejected contours\")\ng.contours.display(ax, only_unused=True)\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Criteria for rejecting a contour:\n 0. - Accepted (green)\n 1. - Rejection for shape error (red)\n 2. - Masked value within contour (blue)\n 3. - Under or over the pixel limit bounds (black)\n 4. - Amplitude criterion (yellow)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Contours' rejection criteria\")\ng.contours.display(ax, only_unused=True, lw=0.5, display_criterion=True)\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display the shape error of each tested contour, the limit of shape error is set to 55 %\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Contour shape error\")\nm = g.contours.display(\n    ax, lw=0.5, field=\"shape_error\", bins=arange(20, 90.1, 5), cmap=\"PRGn_r\"\n)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Some closed contours contains several eddies (aka, more than one extremum)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT rejected contours containing eddies\")\ng.contours.label_contour_unused_which_contain_eddies(a)\ng.contours.label_contour_unused_which_contain_eddies(c)\ng.contours.display(\n    ax,\n    only_contain_eddies=True,\n    color=\"k\",\n    lw=1,\n    label=\"Could be a contour of interaction\",\n)\na.display(ax, color=\"r\", linewidth=0.75, label=\"Anticyclonic\", ref=-10)\nc.display(ax, color=\"b\", linewidth=0.75, label=\"Cyclonic\", ref=-10)\nax.legend()\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Output\nWhen displaying the detected eddies, dashed lines are for effective contour, solide lines for the contour of\nthe maximum mean speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Detected Eddies\")\na.display(\n    ax, color=\"r\", linewidth=0.75, label=\"Anticyclonic ({nb_obs} eddies)\", ref=-10\n)\nc.display(ax, color=\"b\", linewidth=0.75, label=\"Cyclonic ({nb_obs} eddies)\", ref=-10)\nax.legend()\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display the speed radius of the detected eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Speed Radius (km)\")\nkwargs = dict(vmin=10, vmax=50, s=80, ref=-10, cmap=\"magma_r\", factor=0.001)\na.scatter(ax, \"radius_s\", **kwargs)\nm = c.scatter(ax, \"radius_s\", **kwargs)\nupdate_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Filling the effective radius contours with the effective radius values\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Effective Radius (km)\")\nkwargs = dict(vmin=10, vmax=80, cmap=\"magma_r\", factor=0.001, lut=14, ref=-10)\na.filled(ax, \"effective_radius\", **kwargs)\nm = c.filled(ax, \"radius_e\", **kwargs)\nupdate_axes(ax, m)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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",
      "version": "3.10.6"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}