improve diagnostic for filter and detection

Author: AntSimi

examples/02_eddy_identification/pet_eddy_detection.py

@@ -40,6 +40,8 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m)
 
 # %%
+# Get u/v
+# -------
 # U/V are deduced from ADT, this algortihm are not usable around equator (~+- 2°)
 g.add_uv("adt")
 ax = start_axes("U/V deduce from ADT (m)")
@@ -50,14 +52,18 @@ def update_axes(ax, mappable=None):
 update_axes(ax, m)
 
 # %%
+# Pre-processings
+# ---------------
 # Apply high filter to remove long scale to highlight mesoscale
-g.bessel_high_filter("adt", 500, order=2)
-ax = start_axes("ADT (m) filtered (500km, order 2)")
+g.bessel_high_filter("adt", 500)
+ax = start_axes("ADT (m) filtered (500km)")
 m = g.display(ax, "adt", vmin=-0.15, vmax=0.15)
 update_axes(ax, m)
 
 # %%
-#  Run identification with slice of 2 mm
+# Identification
+# --------------
+# run identification with slice of 2 mm
 date = datetime(2016, 5, 15)
 a, c = g.eddy_identification("adt", "u", "v", date, 0.002)
 
@@ -93,6 +99,8 @@ def update_axes(ax, mappable=None):
 update_axes(ax)
 
 # %%
+# Output
+# ------
 # Display detected eddies, dashed lines represent effective contour
 # and solid lines represent contour of maximum of speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1
 

examples/02_eddy_identification/pet_filter_and_detection.py

@@ -7,6 +7,7 @@
 from matplotlib import pyplot as plt
 from py_eddy_tracker.dataset.grid import RegularGridDataset
 from py_eddy_tracker import data
+from py_eddy_tracker.generic import reverse_index
 from numpy import arange
 
 
@@ -35,7 +36,7 @@ def update_axes(ax, mappable=None):
 g.add_uv("adt")
 g.copy("adt", "adt_high")
 wavelength = 400
-g.bessel_high_filter("adt_high", wavelength, order=3)
+g.bessel_high_filter("adt_high", wavelength)
 date = datetime(2016, 5, 15)
 
 # %%
@@ -60,7 +61,6 @@ def update_axes(ax, mappable=None):
 fig = plt.figure(figsize=(12, 5))
 ax_a = plt.subplot(121, xlabel="amplitdue(cm)")
 ax_r = plt.subplot(122, xlabel="speed radius (km)")
-ax_a.grid(), ax_r.grid()
 ax_a.hist(
     merge_f["amplitude"] * 100,
     bins=arange(0.0005, 100, 1),
@@ -86,15 +86,39 @@ def update_axes(ax, mappable=None):
 i, j, c = merge_f.match(merge_r)
 m = c > 0.1
 i_, j_ = i[m], j[m]
+
+
+# %%
+# where is lonely eddies
+kwargs_f = dict(lw=1.5, label="Lonely eddy from filtered grid", ref=-10, color="k")
+kwargs_r = dict(lw=1.5, label="Lonely eddy from raw grid", ref=-10, color="r")
+ax = start_axes("Eddies with no match, over filtered ADT")
+mappable = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
+merge_f.index(reverse_index(i_, len(merge_f))).display(ax, **kwargs_f)
+merge_r.index(reverse_index(j_, len(merge_r))).display(ax, **kwargs_r)
+ax.legend()
+update_axes(ax, mappable)
+
+ax = start_axes("Eddies with no match, over filtered ADT (zoom)")
+ax.set_xlim(25, 36), ax.set_ylim(31, 35.25)
+mappable = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
+u, v = g.grid("u").T, g.grid("v").T
+ax.quiver(g.x_c, g.y_c, u, v, scale=10, pivot="mid", color="gray")
+merge_f.index(reverse_index(i_, len(merge_f))).display(ax, **kwargs_f)
+merge_r.index(reverse_index(j_, len(merge_r))).display(ax, **kwargs_r)
+ax.legend()
+update_axes(ax, mappable)
+
+# %%
 fig = plt.figure(figsize=(12, 12))
-fig.suptitle(f"Scatter plot of speed_radius(km) ({m.sum()} matches)")
+fig.suptitle(f"Scatter plot ({m.sum()} matches)")
 
 for i, (label, field, factor, stop) in enumerate(
-    zip(
-        ("speed radius (km)", "outter radius (km)", "amplitude (cm)"),
-        ("radius_s", "radius_e", "amplitude"),
-        (0.001, 0.001, 100),
-        (80, 120, 25),
+    (
+        ("speed radius (km)", "radius_s", 0.001, 80),
+        ("outter radius (km)", "radius_e", 0.001, 120),
+        ("amplitude (cm)", "amplitude", 100, 25),
+        ("speed max (cm/s)", "speed_average", 100, 25),
     )
 ):
     ax = fig.add_subplot(

notebooks/python_module/02_eddy_identification/pet_eddy_detection.ipynb

@@ -62,7 +62,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "U/V are deduced from ADT, this algortihm are not usable around equator (~+- 2\u00b0)\n\n"
+        "Get u/v\n-------\nU/V are deduced from ADT, this algortihm are not usable around equator (~+- 2\u00b0)\n\n"
       ]
     },
     {
@@ -80,7 +80,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Apply high filter to remove long scale to highlight mesoscale\n\n"
+        "Pre-processings\n---------------\nApply high filter to remove long scale to highlight mesoscale\n\n"
       ]
     },
     {
@@ -91,14 +91,14 @@
       },
       "outputs": [],
       "source": [
-        "g.bessel_high_filter(\"adt\", 500, order=2)\nax = start_axes(\"ADT (m) filtered (500km, order 2)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
+        "g.bessel_high_filter(\"adt\", 500)\nax = start_axes(\"ADT (m) filtered (500km)\")\nm = g.display(ax, \"adt\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
       ]
     },
     {
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Run identification with slice of 2 mm\n\n"
+        "Identification\n--------------\nrun identification with slice of 2 mm\n\n"
       ]
     },
     {
@@ -188,7 +188,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Display detected eddies, dashed lines represent effective contour\nand solid lines represent contour of maximum of speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1\n\n"
+        "Output\n------\nDisplay detected eddies, dashed lines represent effective contour\nand solid lines represent contour of maximum of speed. See figure 1 of https://doi.org/10.1175/JTECH-D-14-00019.1\n\n"
       ]
     },
     {

notebooks/python_module/02_eddy_identification/pet_filter_and_detection.ipynb

@@ -26,7 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "from datetime import datetime\nfrom matplotlib import pyplot as plt\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom numpy import arange"
+        "from datetime import datetime\nfrom matplotlib import pyplot as plt\nfrom py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom py_eddy_tracker.generic import reverse_index\nfrom numpy import arange"
       ]
     },
     {
@@ -55,7 +55,7 @@
       },
       "outputs": [],
       "source": [
-        "g = RegularGridDataset(\n    data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\",\n)\ng.add_uv(\"adt\")\ng.copy(\"adt\", \"adt_high\")\nwavelength = 400\ng.bessel_high_filter(\"adt_high\", wavelength, order=3)\ndate = datetime(2016, 5, 15)"
+        "g = RegularGridDataset(\n    data.get_path(\"dt_med_allsat_phy_l4_20160515_20190101.nc\"), \"longitude\", \"latitude\",\n)\ng.add_uv(\"adt\")\ng.copy(\"adt\", \"adt_high\")\nwavelength = 400\ng.bessel_high_filter(\"adt_high\", wavelength)\ndate = datetime(2016, 5, 15)"
       ]
     },
     {
@@ -109,7 +109,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(12, 5))\nax_a = plt.subplot(121, xlabel=\"amplitdue(cm)\")\nax_r = plt.subplot(122, xlabel=\"speed radius (km)\")\nax_a.grid(), ax_r.grid()\nax_a.hist(\n    merge_f[\"amplitude\"] * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddy from filtered grid\",\n    histtype=\"step\",\n)\nax_a.hist(\n    merge_r[\"amplitude\"] * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddy from raw grid\",\n    histtype=\"step\",\n)\nax_a.set_xlim(0, 10)\nax_r.hist(merge_f[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.hist(merge_r[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.set_xlim(0, 100)\nax_a.legend()"
+        "fig = plt.figure(figsize=(12, 5))\nax_a = plt.subplot(121, xlabel=\"amplitdue(cm)\")\nax_r = plt.subplot(122, xlabel=\"speed radius (km)\")\nax_a.hist(\n    merge_f[\"amplitude\"] * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddy from filtered grid\",\n    histtype=\"step\",\n)\nax_a.hist(\n    merge_r[\"amplitude\"] * 100,\n    bins=arange(0.0005, 100, 1),\n    label=\"Eddy from raw grid\",\n    histtype=\"step\",\n)\nax_a.set_xlim(0, 10)\nax_r.hist(merge_f[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.hist(merge_r[\"radius_s\"] / 1000.0, bins=arange(0, 300, 5), histtype=\"step\")\nax_r.set_xlim(0, 100)\nax_a.legend()"
       ]
     },
     {
@@ -127,7 +127,36 @@
       },
       "outputs": [],
       "source": [
-        "i, j, c = merge_f.match(merge_r)\nm = c > 0.1\ni_, j_ = i[m], j[m]\nfig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot of speed_radius(km) ({m.sum()} matches)\")\n\nfor i, (label, field, factor, stop) in enumerate(\n    zip(\n        (\"speed radius (km)\", \"outter radius (km)\", \"amplitude (cm)\"),\n        (\"radius_s\", \"radius_e\", \"amplitude\"),\n        (0.001, 0.001, 100),\n        (80, 120, 25),\n    )\n):\n    ax = fig.add_subplot(\n        2, 2, i + 1, xlabel=\"filtered grid\", ylabel=\"raw grid\", title=label\n    )\n    ax.plot(merge_f[field][i_] * factor, merge_r[field][j_] * factor, \".\")\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.plot((0, 1000), (0, 1000), \"r\")\n    ax.set_xlim(0, stop), ax.set_ylim(0, stop)"
+        "i, j, c = merge_f.match(merge_r)\nm = c > 0.1\ni_, j_ = i[m], j[m]"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "where is lonely eddies\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "kwargs_f = dict(lw=1.5, label=\"Lonely eddy from filtered grid\", ref=-10, color=\"k\")\nkwargs_r = dict(lw=1.5, label=\"Lonely eddy from raw grid\", ref=-10, color=\"r\")\nax = start_axes(\"Eddies with no match, over filtered ADT\")\nmappable = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nmerge_f.index(reverse_index(i_, len(merge_f))).display(ax, **kwargs_f)\nmerge_r.index(reverse_index(j_, len(merge_r))).display(ax, **kwargs_r)\nax.legend()\nupdate_axes(ax, mappable)\n\nax = start_axes(\"Eddies with no match, over filtered ADT (zoom)\")\nax.set_xlim(25, 36), ax.set_ylim(31, 35.25)\nmappable = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nu, v = g.grid(\"u\").T, g.grid(\"v\").T\nax.quiver(g.x_c, g.y_c, u, v, scale=10, pivot=\"mid\", color=\"gray\")\nmerge_f.index(reverse_index(i_, len(merge_f))).display(ax, **kwargs_f)\nmerge_r.index(reverse_index(j_, len(merge_r))).display(ax, **kwargs_r)\nax.legend()\nupdate_axes(ax, mappable)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "fig = plt.figure(figsize=(12, 12))\nfig.suptitle(f\"Scatter plot ({m.sum()} 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(\n        2, 2, i + 1, xlabel=\"filtered grid\", ylabel=\"raw grid\", title=label\n    )\n    ax.plot(merge_f[field][i_] * factor, merge_r[field][j_] * factor, \".\")\n    ax.set_aspect(\"equal\"), ax.grid()\n    ax.plot((0, 1000), (0, 1000), \"r\")\n    ax.set_xlim(0, stop), ax.set_ylim(0, stop)"
       ]
     }
   ],

src/py_eddy_tracker/generic.py

@@ -37,10 +37,18 @@
     where,
     isnan,
 )
-from numba import njit, prange
+from numba import njit, prange, types as numba_types
 from numpy.linalg import lstsq
 
 
+@njit(cache=True)
+def reverse_index(index, nb):
+    m = ones(nb, dtype=numba_types.bool_)
+    for i in index:
+        m[i] = False
+    return where(m)[0]
+
+
 @njit(cache=True, fastmath=True, parallel=False)
 def distance_grid(lon0, lat0, lon1, lat1):
     """