{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eddy detection : Gulf stream\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\nfrom py_eddy_tracker.eddy_feature import Contours"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def start_axes(title):\n    fig = plt.figure(figsize=(13, 8))\n    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])\n    ax.set_xlim(279, 304), ax.set_ylim(29, 44)\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": [
        "margin = 30\ng = RegularGridDataset(\n    data.get_path(\"nrt_global_allsat_phy_l4_20190223_20190226.nc\"),\n    \"longitude\",\n    \"latitude\",\n    # Manual area subset\n    indexs=dict(\n        longitude=slice(1116 - margin, 1216 + margin),\n        latitude=slice(476 - margin, 536 + margin),\n    ),\n)\n\nax = start_axes(\"ADT (m)\")\nm = g.display(ax, \"adt\", vmin=-1, vmax=1, cmap=\"RdBu_r\")\n# Draw line on the gulf stream front\ngreat_current = Contours(g.x_c, g.y_c, g.grid(\"adt\"), levels=(0.35,), keep_unclose=True)\ngreat_current.display(ax, color=\"k\")\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\")"
      ]
    },
    {
      "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\", 700)\nax = start_axes(\"ADT (m) filtered (700km)\")\nm = g.display(ax, \"adt\", vmin=-0.4, vmax=0.4, cmap=\"RdBu_r\")\ngreat_current.display(ax, color=\"k\")\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 5)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT closed contours (only 1 / 5 levels)\")\ng.contours.display(ax, step=5, lw=1)\ngreat_current.display(ax, color=\"k\")\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, lw=0.25)\ngreat_current.display(ax, color=\"k\")\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, lw=0.25)\ngreat_current.display(ax, color=\"k\")\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 the\nmaximum 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(\"Eddies detected\")\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()\ngreat_current.display(ax, color=\"k\")\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display the effective radius of the detected eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Effective radius (km)\")\na.filled(ax, \"radius_e\", vmin=10, vmax=150, cmap=\"magma_r\", factor=0.001, lut=14)\nm = c.filled(ax, \"radius_e\", vmin=10, vmax=150, cmap=\"magma_r\", factor=0.001, lut=14)\ngreat_current.display(ax, color=\"k\")\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.7.9"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}