Add a method to get mean of grid inside a contour

Add docstrings

Author: AntSimi

doc/autodoc/featured_tracking.rst

@@ -4,4 +4,4 @@ Featured tracking
 .. automodule:: py_eddy_tracker.featured_tracking.old_tracker_reference
     :members:
     :undoc-members:
-    :show-inheritance:
+    :show-inheritance:

doc/autodoc/generic.rst

@@ -0,0 +1,7 @@
+Generic function
+================
+
+.. automodule:: py_eddy_tracker.generic
+    :members:
+    :undoc-members:
+    :show-inheritance:

doc/conf.py

@@ -15,7 +15,7 @@
 # import sys
 # import os
 import sphinx_rtd_theme
-
+import py_eddy_tracker
 
 # If extensions (or modules to document with autodoc) are in another directory,
 # add these directories to sys.path here. If the directory is relative to the
@@ -82,9 +82,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = u"3.0"
+version = py_eddy_tracker.__version__
 # The full version, including alpha/beta/rc tags.
-release = u"3.0"
+release = py_eddy_tracker.__version__
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

examples/02_eddy_identification/pet_interp_grid_on_dataset.py

@@ -0,0 +1,64 @@
+"""
+Get mean of grid in each eddies
+===============================
+
+"""
+
+from matplotlib import pyplot as plt
+from py_eddy_tracker.dataset.grid import RegularGridDataset
+from py_eddy_tracker.observations.observation import EddiesObservations
+from py_eddy_tracker import data
+
+
+# %%
+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()
+    ax.legend()
+    if mappable:
+        plt.colorbar(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))
+
+
+# %%
+# Load detection files and data to interp
+a = EddiesObservations.load_file(data.get_path("Anticyclonic_20160515.nc"))
+c = EddiesObservations.load_file(data.get_path("Cyclonic_20160515.nc"))
+
+aviso_map = RegularGridDataset(
+    data.get_path("dt_med_allsat_phy_l4_20160515_20190101.nc"), "longitude", "latitude"
+)
+aviso_map.add_uv("adt")
+
+# %%
+# Compute and store eke in cm²/s²
+aviso_map.add_grid(
+    "eke", (aviso_map.grid("u") ** 2 + aviso_map.grid("v") ** 2) * 0.5 * (100 ** 2)
+)
+
+eke_kwargs = dict(vmin=1, vmax=1000, cmap="magma_r")
+
+ax = start_axes("EKE (cm²/s²)")
+m = aviso_map.display(ax, "eke", **eke_kwargs)
+a.display(ax, color="r", linewidth=0.5, label="Anticyclonic", ref=-10)
+c.display(ax, color="b", linewidth=0.5, label="Cyclonic", ref=-10)
+update_axes(ax, m)
+
+# %%
+# Get mean of eke in each effective contour
+
+ax = start_axes("EKE (cm²/s²)")
+a.display(ax, color="r", linewidth=0.5, label="Anticyclonic", ref=-10)
+c.display(ax, color="b", linewidth=0.5, label="Cyclonic", ref=-10)
+eke = a.interp_grid(aviso_map, "eke", method="mean", intern=False)
+m = ax.scatter((a.longitude + 10) % 360 - 10, a.latitude, c=eke, **eke_kwargs)
+eke = c.interp_grid(aviso_map, "eke", method="mean", intern=False)
+m = ax.scatter((c.longitude + 10) % 360 - 10, c.latitude, c=eke, **eke_kwargs)
+update_axes(ax, m)

examples/10_tracking_diagnostics/pet_propagation.py

@@ -1,6 +1,6 @@
 """
 Propagation Histogram
-===================
+=====================
 
 """
 from matplotlib import pyplot as plt

notebooks/python_module/02_eddy_identification/pet_interp_grid_on_dataset.ipynb

@@ -0,0 +1,119 @@
+{
+  "cells": [
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "%matplotlib inline"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Get mean of grid in each eddies\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "from matplotlib import pyplot as plt\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker.observations.observation import EddiesObservations\nfrom py_eddy_tracker import data"
+      ]
+    },
+    {
+      "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    ax.legend()\n    if mappable:\n        plt.colorbar(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Load detection files and data to interp\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "a = EddiesObservations.load_file(data.get_path(\"Anticyclonic_20160515.nc\"))\nc = EddiesObservations.load_file(data.get_path(\"Cyclonic_20160515.nc\"))\n\naviso_map = RegularGridDataset(\n    data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\"\n)\n# aviso_map.add_uv(\"adt\")\naviso_map.add_uv(\"sla\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Compute and store eke in cm\u00b2/s\u00b2\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "aviso_map.add_grid(\n    \"eke\", (aviso_map.grid(\"u\") ** 2 + aviso_map.grid(\"v\") ** 2) * 0.5 * (100 ** 2)\n)\n\neke_kwargs = dict(vmin=1, vmax=1000, cmap=\"magma_r\")\n\nax = start_axes(\"EKE (cm\u00b2/s\u00b2)\")\nm = aviso_map.display(ax, \"eke\", **eke_kwargs)\na.display(ax, color=\"r\", linewidth=0.5, label=\"Anticyclonic\", ref=-10)\nc.display(ax, color=\"b\", linewidth=0.5, label=\"Cyclonic\", ref=-10)\nupdate_axes(ax, m)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Get mean of eke in each effective contour\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "ax = start_axes(\"EKE (cm\u00b2/s\u00b2)\")\na.display(ax, color=\"r\", linewidth=0.5, label=\"Anticyclonic\", ref=-10)\nc.display(ax, color=\"b\", linewidth=0.5, label=\"Cyclonic\", ref=-10)\neke = a.interp_grid(aviso_map, \"eke\", method=\"mean\", intern=False)\nm = ax.scatter((a.longitude + 10) % 360 - 10, a.latitude, c=eke, **eke_kwargs)\neke = c.interp_grid(aviso_map, \"eke\", method=\"mean\", intern=False)\nm = ax.scatter((c.longitude + 10) % 360 - 10, c.latitude, c=eke, **eke_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.7.0"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

src/py_eddy_tracker/data/Anticyclonic_20160515.nc

@@ -3,6 +3,10 @@
     nrt_global_allsat_phy_l4_20190223_20190226.nc \
     20190223 adt ugos vgos longitude latitude . \
     --cut 800 --fil 1
+EddyId \
+    dt_med_allsat_phy_l4_20160515_20190101.nc \
+    20160515 adt None None longitude latitude . \
+    --cut 800 --fil 1
 """
 from os import path
 import requests

src/py_eddy_tracker/data/Cyclonic_20160515.nc

@@ -54,10 +54,12 @@
     uniform_resample,
     coordinates_to_local,
     local_to_coordinates,
+    get_pixel_in_regular,
+    nearest_grd_indice,
+    bbox_indice_regular,
 )
 from ..poly import (
     poly_contain_poly,
-    winding_number_grid_in_poly,
     winding_number_poly,
     create_vertice,
     poly_area,
@@ -174,36 +176,6 @@ def _fit_circle_path(vertice):
     return centlon, centlat, eddy_radius, err
 
 
-@njit(cache=True, fastmath=True)
-def _get_pixel_in_regular(vertices, x_c, y_c, x_start, x_stop, y_start, y_stop):
-    if x_stop < x_start:
-        x_ref = vertices[0, 0]
-        x_array = (
-            (concatenate((x_c[x_start:], x_c[:x_stop])) - x_ref + 180) % 360
-            + x_ref
-            - 180
-        )
-        return winding_number_grid_in_poly(
-            x_array,
-            y_c[y_start:y_stop],
-            x_start,
-            x_stop,
-            x_c.shape[0],
-            y_start,
-            vertices,
-        )
-    else:
-        return winding_number_grid_in_poly(
-            x_c[x_start:x_stop],
-            y_c[y_start:y_stop],
-            x_start,
-            x_stop,
-            x_c.shape[0],
-            y_start,
-            vertices,
-        )
-
-
 @njit(cache=True, fastmath=True)
 def _get_pixel_in_unregular(vertices, x_c, y_c, x_start, x_stop, y_start, y_stop):
     nb_x, nb_y = x_stop - x_start, y_stop - y_start
@@ -454,6 +426,7 @@ def is_circular(self):
         return False
 
     def units(self, varname):
+        """Get unit from variable"""
         stored_units = self.variables_description[varname]["attrs"].get("units", None)
         if stored_units is not None:
             return stored_units
@@ -479,6 +452,15 @@ def copy(self, grid_in, grid_out):
         )
         self.vars[grid_out] = self.grid(grid_in).copy()
 
+    def add_grid(self, varname, grid):
+        """
+        Add a grid in handler
+
+        :param str varname: name of future grid
+        :param array grid: grid array
+        """
+        self.vars[varname] = grid
+
     def grid(self, varname, indexs=None):
         """give grid required
         """
@@ -1180,15 +1162,15 @@ def bbox_indice(self, vertices):
 
     def get_pixels_in(self, contour):
         (x_start, x_stop), (y_start, y_stop) = contour.bbox_slice
-        return _get_pixel_in_regular(
+        return get_pixel_in_regular(
             contour.vertices, self.x_c, self.y_c, x_start, x_stop, y_start, y_stop
         )
 
     def normalize_x_indice(self, indices):
         return indices % self.x_size
 
     def nearest_grd_indice(self, x, y):
-        return _nearest_grd_indice(
+        return nearest_grd_indice(
             x, y, self.x_bounds, self.y_bounds, self.xstep, self.ystep
         )
 
@@ -1846,29 +1828,6 @@ def compute_pixel_path(x0, y0, x1, y1, x_ori, y_ori, x_step, y_step, nb_x):
     return i_g, j_g, d_max
 
 
-@njit(cache=True)
-def bbox_indice_regular(vertices, x0, y0, xstep, ystep, N, circular, x_size):
-    lon, lat = vertices[:, 0], vertices[:, 1]
-    lon_min, lon_max = lon.min(), lon.max()
-    lat_min, lat_max = lat.min(), lat.max()
-    i_x0, i_y0 = _nearest_grd_indice(lon_min, lat_min, x0, y0, xstep, ystep)
-    i_x1, i_y1 = _nearest_grd_indice(lon_max, lat_max, x0, y0, xstep, ystep)
-    if circular:
-        slice_x = (i_x0 - N) % x_size, (i_x1 + N + 1) % x_size
-    else:
-        slice_x = max(i_x0 - N, 0), i_x1 + N + 1
-    slice_y = i_y0 - N, i_y1 + N + 1
-    return slice_x, slice_y
-
-
-@njit(cache=True, fastmath=True)
-def _nearest_grd_indice(x, y, x0, y0, xstep, ystep):
-    return (
-        numba_types.int32(round(((x - x0[0]) % 360.0) / xstep)),
-        numba_types.int32(round((y - y0[0]) / ystep)),
-    )
-
-
 @njit(cache=True)
 def has_masked_value(grid, i_x, i_y):
     for i, j in zip(i_x, i_y):