{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Count center\n\nDo Geo stat with center and compare with frequency method\nshow: `sphx_glr_python_module_10_tracking_diagnostics_pet_pixel_used.py`\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib import pyplot as plt\nfrom matplotlib.colors import LogNorm\nimport py_eddy_tracker_sample\n\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Load an experimental med atlas over a period of 26 years (1993-2019)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "a = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_demo_path(\n        \"eddies_med_adt_allsat_dt2018/Anticyclonic.zarr\"\n    )\n)\nc = TrackEddiesObservations.load_file(\n    py_eddy_tracker_sample.get_demo_path(\"eddies_med_adt_allsat_dt2018/Cyclonic.zarr\")\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Parameters\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "step = 0.125\nbins = ((-10, 37, step), (30, 46, step))\nkwargs_pcolormesh = dict(\n    cmap=\"terrain_r\", vmin=0, vmax=2, factor=1 / (a.nb_days * step**2), name=\"count\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 18.5))\nax_a = fig.add_axes([0.03, 0.75, 0.90, 0.25])\nax_a.set_title(\"Anticyclonic center frequency\")\nax_c = fig.add_axes([0.03, 0.5, 0.90, 0.25])\nax_c.set_title(\"Cyclonic center frequency\")\nax_all = fig.add_axes([0.03, 0.25, 0.90, 0.25])\nax_all.set_title(\"All eddies center frequency\")\nax_ratio = fig.add_axes([0.03, 0.0, 0.90, 0.25])\nax_ratio.set_title(\"Ratio cyclonic / Anticyclonic\")\n\n# Count pixel used for each center\ng_a = a.grid_count(bins, intern=True, center=True)\ng_a.display(ax_a, **kwargs_pcolormesh)\ng_c = c.grid_count(bins, intern=True, center=True)\ng_c.display(ax_c, **kwargs_pcolormesh)\n# Compute a ratio Cyclonic / Anticyclonic\nratio = g_c.vars[\"count\"] / g_a.vars[\"count\"]\n\n# Mask manipulation to be able to sum the 2 grids\nm_c = g_c.vars[\"count\"].mask\nm = m_c & g_a.vars[\"count\"].mask\ng_c.vars[\"count\"][m_c] = 0\ng_c.vars[\"count\"] += g_a.vars[\"count\"]\ng_c.vars[\"count\"].mask = m\n\nm = g_c.display(ax_all, **kwargs_pcolormesh)\ncb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.27, 0.01, 0.7]))\ncb.set_label(\"Eddies by 1\u00b0^2 by day\")\n\ng_c.vars[\"count\"] = ratio\nm = g_c.display(\n    ax_ratio, name=\"count\", norm=LogNorm(vmin=0.1, vmax=10), cmap=\"coolwarm_r\"\n)\nplt.colorbar(m, cax=fig.add_axes([0.94, 0.02, 0.01, 0.2]))\n\nfor ax in (ax_a, ax_c, ax_all, ax_ratio):\n    ax.set_aspect(\"equal\")\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\n    ax.grid()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Count Anticyclones as a function of lifetime\nCount at the center's position\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig = plt.figure(figsize=(12, 10))\nmask = a.lifetime >= 60\nax_long = fig.add_axes([0.03, 0.53, 0.90, 0.45])\ng_a = a.grid_count(bins, center=True, filter=mask)\ng_a.display(ax_long, **kwargs_pcolormesh)\nax_long.set_title(f\"Anticyclones with lifetime >= 60 days ({mask.sum()} Obs)\")\nax_short = fig.add_axes([0.03, 0.03, 0.90, 0.45])\ng_a = a.grid_count(bins, center=True, filter=~mask)\nm = g_a.display(ax_short, **kwargs_pcolormesh)\nax_short.set_title(f\"Anticyclones with lifetime < 60 days ({(~mask).sum()} Obs)\")\nfor ax in (ax_short, ax_long):\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)\ncb = plt.colorbar(m, cax=fig.add_axes([0.94, 0.05, 0.015, 0.9]))"
      ]
    }
  ],
  "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
}