example to compute statistics on raw identification

README.md

@@ -17,15 +17,17 @@ Method was described in :
 ### Use case ###
 
 Method is used in :
- 
+
 [Mason, E., A. Pascual, P. Gaube, S.Ruiz, J. Pelegrí, A. Delepoulle, 2017: Subregional characterization of mesoscale eddies across the Brazil-Malvinas Confluence](https://doi.org/10.1002/2016JC012611)
 
 ### How do I get set up? ###
 
 #### Short story ####
+
 ```bash
 pip install pyeddytracker
 ```
+
 #### Long story ####
 
 To avoid problems with installation, use of the virtualenv Python virtual environment is recommended.
@@ -36,12 +38,20 @@ Then use pip to install all dependencies (numpy, scipy, matplotlib, netCDF4, ...
 pip install numpy scipy netCDF4 matplotlib opencv-python pyyaml pint polygon3
 ```
 
-Then run the following to install the eddy tracker:
+Clone :
+
+```bash
+git clone https://github.com/AntSimi/py-eddy-tracker
+```
+
+Then run the following to install the eddy tracker :
 
 ```bash
 python setup.py install
 ```
+
 ### Tools gallery ###
+
 Several examples based on py eddy tracker module are [here](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html).
 
 [![](https://py-eddy-tracker.readthedocs.io/en/latest/_static/logo.png)](https://py-eddy-tracker.readthedocs.io/en/latest/python_module/index.html)

examples/02_eddy_identification/pet_statistics_on_identification.py

@@ -0,0 +1,105 @@
+"""
+Stastics on identification files
+================================
+
+Some statistics on raw identification without any tracking
+"""
+import numpy as np
+from matplotlib import pyplot as plt
+from matplotlib.dates import date2num
+
+from py_eddy_tracker import start_logger
+from py_eddy_tracker.data import get_remote_demo_sample
+from py_eddy_tracker.observations.observation import EddiesObservations
+
+start_logger().setLevel("ERROR")
+
+
+# %%
+def start_axes(title):
+    fig = plt.figure(figsize=(13, 5))
+    ax = fig.add_axes([0.03, 0.03, 0.90, 0.94])
+    ax.set_xlim(-6, 36.5), ax.set_ylim(30, 46)
+    ax.set_aspect("equal")
+    ax.set_title(title)
+    return ax
+
+
+def update_axes(ax, mappable=None):
+    ax.grid()
+    if mappable:
+        plt.colorbar(mappable, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
+
+
+# %%
+# We load demo sample and take only first year.
+#
+# Replace by a list of filename to apply on your own dataset.
+file_objects = get_remote_demo_sample(
+    "eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012"
+)[:365]
+
+# %%
+# Merge all identification datasets in one object
+all_a = EddiesObservations.concatenate(
+    [EddiesObservations.load_file(i) for i in file_objects]
+)
+
+# %%
+# We define polygon bound
+x0, x1, y0, y1 = 15, 20, 33, 38
+xs = np.array([[x0, x1, x1, x0, x0]], dtype="f8")
+ys = np.array([[y0, y0, y1, y1, y0]], dtype="f8")
+# Polygon object created for the match function use.
+polygon = dict(contour_lon_e=xs, contour_lat_e=ys, contour_lon_s=xs, contour_lat_s=ys)
+
+# %%
+# Geographic frequency of eddies
+step = 0.125
+ax = start_axes("")
+# Count pixel encompassed in each contour
+g_a = all_a.grid_count(bins=((-10, 37, step), (30, 46, step)), intern=True)
+m = g_a.display(
+    ax, cmap="terrain_r", vmin=0, vmax=0.75, factor=1 / all_a.nb_days, name="count"
+)
+ax.plot(polygon["contour_lon_e"][0], polygon["contour_lat_e"][0], "r")
+update_axes(ax, m)
+
+# %%
+# We use the match function to count the number of eddies that intersect the polygon defined previously
+# `p1_area` option allow to get in c_e/c_s output, precentage of area occupy by eddies in the polygon.
+i_e, j_e, c_e = all_a.match(polygon, p1_area=True, intern=False)
+i_s, j_s, c_s = all_a.match(polygon, p1_area=True, intern=True)
+
+# %%
+dt = np.datetime64("1970-01-01") - np.datetime64("1950-01-01")
+kw_hist = dict(
+    bins=date2num(np.arange(21900, 22300).astype("datetime64") - dt), histtype="step"
+)
+# translate julian day in datetime64
+t = all_a.time.astype("datetime64") - dt
+# %%
+# Number of eddies within a polygon
+ax = plt.figure(figsize=(12, 6)).add_subplot(111)
+ax.set_title("Different ways to count eddies within a polygon")
+ax.set_ylabel("Count")
+m = all_a.mask_from_polygons(((xs, ys),))
+ax.hist(t[m], label="Eddy Center in polygon", **kw_hist)
+ax.hist(t[i_s[c_s > 0]], label="Intersection Speed contour and polygon", **kw_hist)
+ax.hist(t[i_e[c_e > 0]], label="Intersection Effective contour and polygon", **kw_hist)
+ax.legend()
+ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
+ax.grid()
+
+# %%
+# Percent of the area of interest occupied by eddies.
+ax = plt.figure(figsize=(12, 6)).add_subplot(111)
+ax.set_title("Percent of polygon occupied by an anticyclonic eddy")
+ax.set_ylabel("Percent of polygon")
+ax.hist(t[i_s[c_s > 0]], weights=c_s[c_s > 0] * 100.0, label="speed contour", **kw_hist)
+ax.hist(
+    t[i_e[c_e > 0]], weights=c_e[c_e > 0] * 100.0, label="effective contour", **kw_hist
+)
+ax.legend(), ax.set_ylim(0, 25)
+ax.set_xlim(np.datetime64("2010"), np.datetime64("2011"))
+ax.grid()

notebooks/python_module/02_eddy_identification/pet_statistics_on_identification.ipynb

@@ -0,0 +1,202 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Stastics on identification files\n\nSome statistics on raw identification without any tracking\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import numpy as np\nfrom matplotlib import pyplot as plt\nfrom matplotlib.dates import date2num\n\nfrom py_eddy_tracker import start_logger\nfrom py_eddy_tracker.data import get_remote_demo_sample\nfrom py_eddy_tracker.observations.observation import EddiesObservations\n\nstart_logger().setLevel(\"ERROR\")"
+      ]
+    },
+    {
+      "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": [
+        "We load demo sample and take only first year.\n\nReplace by a list of filename to apply on your own dataset.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "file_objects = get_remote_demo_sample(\n    \"eddies_med_adt_allsat_dt2018/Anticyclonic_2010_2011_2012\"\n)[:365]"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Merge all identification dataset in one object\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "all_a = EddiesObservations.concatenate(\n    [EddiesObservations.load_file(i) for i in file_objects]\n)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We define polygon bound\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "x0, x1, y0, y1 = 15, 20, 33, 38\nxs = np.array([[x0, x1, x1, x0, x0]], dtype=\"f8\")\nys = np.array([[y0, y0, y1, y1, y0]], dtype=\"f8\")\n# Polygon object is create to be usable by match function.\npolygon = dict(contour_lon_e=xs, contour_lat_e=ys, contour_lon_s=xs, contour_lat_s=ys)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Geographic frequency of eddies\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "step = 0.125\nax = start_axes(\"\")\n# Count pixel used for each contour\ng_a = all_a.grid_count(bins=((-10, 37, step), (30, 46, step)), intern=True)\nm = g_a.display(\n    ax, cmap=\"terrain_r\", vmin=0, vmax=0.75, factor=1 / all_a.nb_days, name=\"count\"\n)\nax.plot(polygon[\"contour_lon_e\"][0], polygon[\"contour_lat_e\"][0], \"r\")\nupdate_axes(ax, m)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "We use match function to count number of eddies which intersect the polygon defined previously.\n`p1_area` option allow to get in c_e/c_s output, precentage of area occupy by eddies in the polygon.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "i_e, j_e, c_e = all_a.match(polygon, p1_area=True, intern=False)\ni_s, j_s, c_s = all_a.match(polygon, p1_area=True, intern=True)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "dt = np.datetime64(\"1970-01-01\") - np.datetime64(\"1950-01-01\")\nkw_hist = dict(\n    bins=date2num(np.arange(21900, 22300).astype(\"datetime64\") - dt), histtype=\"step\"\n)\n# translate julian day in datetime64\nt = all_a.time.astype(\"datetime64\") - dt"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Count how many are in polygon\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "ax = plt.figure(figsize=(12, 6)).add_subplot(111)\nax.set_title(\"Different way to count eddies presence in a polygon\")\nax.set_ylabel(\"Count\")\nm = all_a.mask_from_polygons(((xs, ys),))\nax.hist(t[m], label=\"center in polygon\", **kw_hist)\nax.hist(t[i_s[c_s > 0]], label=\"intersect speed contour with polygon\", **kw_hist)\nax.hist(t[i_e[c_e > 0]], label=\"intersect extern contour with polygon\", **kw_hist)\nax.legend()\nax.set_xlim(np.datetime64(\"2010\"), np.datetime64(\"2011\"))\nax.grid()"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Percent of are of interest occupy by eddies\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "ax = plt.figure(figsize=(12, 6)).add_subplot(111)\nax.set_title(\"Percent of polygon occupy by an anticyclonic eddy\")\nax.set_ylabel(\"Percent of polygon\")\nax.hist(t[i_s[c_s > 0]], weights=c_s[c_s > 0] * 100.0, label=\"speed contour\", **kw_hist)\nax.hist(t[i_e[c_e > 0]], weights=c_e[c_e > 0] * 100.0, label=\"effective contour\", **kw_hist)\nax.legend(), ax.set_ylim(0, 25)\nax.set_xlim(np.datetime64(\"2010\"), np.datetime64(\"2011\"))\nax.grid()"
+      ]
+    }
+  ],
+  "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.7"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

src/py_eddy_tracker/dataset/grid.py

@@ -807,7 +807,7 @@ def eddy_identification(
                         else:
                             centi = reset_centroid[0]
                         centj = reset_centroid[1]
-                        # To move in regular and unregular grid
+                        # FIXME : To move in regular and unregular grid
                         if len(x.shape) == 1:
                             centlon_e = x[centi]
                             centlat_e = y[centj]

src/py_eddy_tracker/observations/observation.py

@@ -8,7 +8,7 @@
 from tarfile import ExFileObject
 from tokenize import TokenError
 
-import packaging
+import packaging.version
 import zarr
 from matplotlib.cm import get_cmap
 from matplotlib.collections import PolyCollection