{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eddy detection on SLA and ADT\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from datetime import datetime\n\nfrom matplotlib import pyplot as plt\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)\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.95, 0.05, 0.01, 0.9]))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load Input grid, ADT will be 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)\ng.add_uv(\"adt\", \"ugos\", \"vgos\")\ng.add_uv(\"sla\", \"ugosa\", \"vgosa\")\nwavelength = 400\ng.copy(\"adt\", \"adt_raw\")\ng.copy(\"sla\", \"sla_raw\")\ng.bessel_high_filter(\"adt\", wavelength)\ng.bessel_high_filter(\"sla\", wavelength)\ndate = datetime(2016, 5, 15)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "kwargs_a_adt = dict(\n    lw=0.5, label=\"Anticyclonic ADT ({nb_obs} eddies)\", ref=-10, color=\"k\"\n)\nkwargs_c_adt = dict(lw=0.5, label=\"Cyclonic ADT ({nb_obs} eddies)\", ref=-10, color=\"r\")\nkwargs_a_sla = dict(\n    lw=0.5, label=\"Anticyclonic SLA ({nb_obs} eddies)\", ref=-10, color=\"g\"\n)\nkwargs_c_sla = dict(lw=0.5, label=\"Cyclonic SLA ({nb_obs} eddies)\", ref=-10, color=\"b\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run algorithm of detection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a_adt, c_adt = g.eddy_identification(\"adt\", \"ugos\", \"vgos\", date, 0.002)\na_sla, c_sla = g.eddy_identification(\"sla\", \"ugosa\", \"vgosa\", date, 0.002)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "over filtered\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(f\"ADT (m) filtered ({wavelength}km)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\na_adt.display(ax, **kwargs_a_adt), c_adt.display(ax, **kwargs_c_adt)\nax.legend(), update_axes(ax, m)\n\nax = start_axes(f\"SLA (m) filtered ({wavelength}km)\")\nm = g.display(ax, \"sla\", vmin=-0.15, vmax=0.15)\na_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)\nax.legend(), update_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "over raw\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"ADT (m)\")\nm = g.display(ax, \"adt_raw\", vmin=-0.15, vmax=0.15)\na_adt.display(ax, **kwargs_a_adt), c_adt.display(ax, **kwargs_c_adt)\nax.legend(), update_axes(ax, m)\n\nax = start_axes(\"SLA (m)\")\nm = g.display(ax, \"sla_raw\", vmin=-0.15, vmax=0.15)\na_sla.display(ax, **kwargs_a_sla), c_sla.display(ax, **kwargs_c_sla)\nax.legend(), update_axes(ax, m)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Display detection\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Eddies detected\")\na_adt.display(ax, **kwargs_a_adt)\na_sla.display(ax, **kwargs_a_sla)\nc_adt.display(ax, **kwargs_c_adt)\nc_sla.display(ax, **kwargs_c_sla)\nax.legend()\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Match\nWhere cyclone meet anticyclone\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_c_adt, i_a_sla, c = c_adt.match(a_sla, cmin=0.01)\ni_a_adt, i_c_sla, c = a_adt.match(c_sla, cmin=0.01)\n\nax = start_axes(\"Cyclone share area with anticyclone\")\na_adt.index(i_a_adt).display(ax, **kwargs_a_adt)\nc_adt.index(i_c_adt).display(ax, **kwargs_c_adt)\na_sla.index(i_a_sla).display(ax, **kwargs_a_sla)\nc_sla.index(i_c_sla).display(ax, **kwargs_c_sla)\nax.legend()\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Scatter plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "i_a_adt, i_a_sla, c = a_adt.match(a_sla, cmin=0.1)\ni_c_adt, i_c_sla, c = c_adt.match(c_sla, cmin=0.1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "where is lonely eddies\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = start_axes(\"Eddies with no match\")\na_adt.index(i_a_adt, reverse=True).display(ax, **kwargs_a_adt)\nc_adt.index(i_c_adt, reverse=True).display(ax, **kwargs_c_adt)\na_sla.index(i_a_sla, reverse=True).display(ax, **kwargs_a_sla)\nc_sla.index(i_c_sla, reverse=True).display(ax, **kwargs_c_sla)\nax.legend()\nupdate_axes(ax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot (A : {i_a_adt.shape[0]}, C : {i_c_adt.shape[0]} matches)\")\n\nfor i, (label, field, factor, stop) in enumerate(\n    (\n        (\"speed radius (km)\", \"radius_s\", 0.001, 80),\n        (\"outter radius (km)\", \"radius_e\", 0.001, 120),\n        (\"amplitude (cm)\", \"amplitude\", 100, 25),\n        (\"speed max (cm/s)\", \"speed_average\", 100, 25),\n    )\n):\n    ax = fig.add_subplot(2, 2, i + 1, title=label)\n    ax.set_xlabel(\"Absolute Dynamic Topography\")\n    ax.set_ylabel(\"Sea Level Anomaly\")\n\n    ax.plot(\n        a_adt[field][i_a_adt] * factor,\n        a_sla[field][i_a_sla] * factor,\n        \"r.\",\n        label=\"Anticyclonic\",\n    )\n    ax.plot(\n        c_adt[field][i_c_adt] * factor,\n        c_sla[field][i_c_sla] * factor,\n        \"b.\",\n        label=\"Cyclonic\",\n    )\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.plot((0, 1000), (0, 1000), \"g\")\n    ax.set_xlim(0, stop), ax.set_ylim(0, stop)\n    ax.legend()"
      ]
    }
  ],
  "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
}