split kernel function for filtering to simplify

Author: AntSimi

examples/06_grid_manipulation/pet_filter.py

@@ -10,6 +10,7 @@
 from py_eddy_tracker.dataset.grid import RegularGridDataset
 from py_eddy_tracker import data
 from matplotlib import pyplot as plt
+from numpy import arange
 
 
 def start_axes(title):
@@ -35,11 +36,11 @@ def update_axes(ax, mappable=None):
 # %%
 # Kernel
 # ------
-# Shape of kernel will increase in x when latitude increase
+# Shape of kernel will increase in x, when latitude increase
 fig = plt.figure(figsize=(12, 8))
 for i, latitude in enumerate((15, 35, 55, 75)):
     k = g.kernel_bessel(latitude, 500, order=3).T
-    ax0 = plt.subplot(
+    ax0 = fig.add_subplot(
         2,
         2,
         i + 1,
@@ -49,6 +50,26 @@ def update_axes(ax, mappable=None):
     m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap="viridis_r")
 plt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))
 
+# %%
+# Kernel along latitude
+
+fig = plt.figure(figsize=(12, 8))
+ax = fig.add_subplot(
+    111,
+    ylabel="Kernel weight",
+    xlabel="Latitude in °",
+    title="Kernel in latitude, centered at 0° of latitude ",
+)
+k = g.kernel_bessel(0, 500, order=3)
+k_lat = k[k.shape[0] // 2 + 1]
+nb = k_lat.shape[0] // 2
+ax.plot(
+    arange(-nb * g.xstep, (nb + 0.5) * g.xstep, g.xstep), k_lat, label="Bessel kernel"
+)
+
+ax.legend()
+ax.grid()
+
 # %%
 # Kernel applying
 # ---------------
@@ -62,28 +83,28 @@ def update_axes(ax, mappable=None):
 #
 # Low frequency
 ax = start_axes("ADT low frequency")
-g.copy("adt", "adt_low")
-g.bessel_low_filter("adt_low", 300, order=3)
-m = g.display(ax, "adt_low", vmin=-0.15, vmax=0.15)
+g.copy("adt", "adt_low_300")
+g.bessel_low_filter("adt_low_300", 300, order=3)
+m = g.display(ax, "adt_low_300", vmin=-0.15, vmax=0.15)
 update_axes(ax, m)
 
 # %%
 # High frequency
 ax = start_axes("ADT high frequency")
-g.copy("adt", "adt_high")
-g.bessel_high_filter("adt_high", 300, order=3)
-m = g.display(ax, "adt_high", vmin=-0.15, vmax=0.15)
+g.copy("adt", "adt_high_300")
+g.bessel_high_filter("adt_high_300", 300, order=3)
+m = g.display(ax, "adt_high_300", vmin=-0.15, vmax=0.15)
 update_axes(ax, m)
 
 # %%
 # Clues
 # -----
 # wavelength : 80km
 
-g.copy("adt", "adt_high_80")
-g.bessel_high_filter("adt_high_80", 80, order=3)
-g.copy("adt", "adt_low_80")
-g.bessel_low_filter("adt_low_80", 80, order=3)
+g.copy("adt", "adt_high_bessel")
+g.bessel_high_filter("adt_high_bessel", 80, order=3)
+g.copy("adt", "adt_low_bessel")
+g.bessel_low_filter("adt_low_bessel", 80, order=3)
 
 area = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)
 
@@ -94,17 +115,12 @@ def update_axes(ax, mappable=None):
 ax.set_title("Spectrum")
 ax.set_xlabel("km")
 
-lon_spec, lat_spec = g.spectrum_lonlat("adt", area=area)
-mappable = ax.loglog(*lat_spec, label="lat raw")[0]
-ax.loglog(*lon_spec, label="lon raw", color=mappable.get_color(), linestyle="--")
-
-lon_spec, lat_spec = g.spectrum_lonlat("adt_high_80", area=area)
-mappable = ax.loglog(*lat_spec, label="lat high")[0]
-ax.loglog(*lon_spec, label="lon high", color=mappable.get_color(), linestyle="--")
-
-lon_spec, lat_spec = g.spectrum_lonlat("adt_low_80", area=area)
-mappable = ax.loglog(*lat_spec, label="lat low")[0]
-ax.loglog(*lon_spec, label="lon low", color=mappable.get_color(), linestyle="--")
+for label in ("adt_high_bessel", "adt_low_bessel", "adt"):
+    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area)
+    mappable = ax.loglog(*lat_spec, label=f"lat {label}")[0]
+    ax.loglog(
+        *lon_spec, label=f"lon {label}", color=mappable.get_color(), linestyle="--"
+    )
 
 ax.set_xlim(10, 1000)
 ax.set_ylim(1e-6, 1)
@@ -119,17 +135,10 @@ def update_axes(ax, mappable=None):
 ax.set_title("Spectrum ratio")
 ax.set_xlabel("km")
 
-lon_spec, lat_spec = g.spectrum_lonlat(
-    "adt_high_80", area=area, ref=g, ref_grid_name="adt"
-)
-mappable = ax.plot(*lat_spec, label="lat high")[0]
-ax.plot(*lon_spec, label="lon high", color=mappable.get_color(), linestyle="--")
-
-lon_spec, lat_spec = g.spectrum_lonlat(
-    "adt_low_80", area=area, ref=g, ref_grid_name="adt"
-)
-mappable = ax.plot(*lat_spec, label="lat low")[0]
-ax.plot(*lon_spec, label="lon low", color=mappable.get_color(), linestyle="--")
+for label in ("adt_high_bessel", "adt_low_bessel"):
+    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area, ref=g, ref_grid_name="adt")
+    mappable = ax.plot(*lat_spec, label=f"lat {label}")[0]
+    ax.plot(*lon_spec, label=f"lon {label}", color=mappable.get_color(), linestyle="--")
 
 ax.set_xlim(10, 1000)
 ax.set_ylim(0, 1)
@@ -141,4 +150,3 @@ def update_axes(ax, mappable=None):
 # Old filter
 # ----------
 # To do ...
-

examples/08_tracking_manipulation/pet_run_a_tracking.py

@@ -10,22 +10,9 @@
 from py_eddy_tracker.data import get_remote_sample
 from py_eddy_tracker.tracking import Correspondances
 from py_eddy_tracker.featured_tracking.area_tracker import AreaTracker
-from numpy import where, empty
 from py_eddy_tracker.gui import GUI
 
 
-# %%
-# Function to have track with contiguous longitude
-def wrap_longitude(eddies):
-    lon = eddies.longitude
-    first = where(eddies.obs["n"] == 0)[0]
-    nb_obs = empty(first.shape, dtype="u4")
-    nb_obs[:-1] = first[1:] - first[:-1]
-    nb_obs[-1] = lon.shape[0] - first[-1]
-    lon0 = (lon[first] - 180).repeat(nb_obs)
-    lon[:] = (lon - lon0) % 360 + lon0
-
-
 # %%
 # Get remote data, we will keep only 180 first days
 file_objects = get_remote_sample(
@@ -34,20 +21,22 @@ def wrap_longitude(eddies):
 
 # %%
 # We run a traking with a tracker which use contour overlap
-c = Correspondances(datasets=file_objects, class_method=AreaTracker)
+c = Correspondances(datasets=file_objects, class_method=AreaTracker, virtual=3)
 c.track()
 c.prepare_merging()
 # We have now an eddy object
 eddies_area_tracker = c.merge(raw_data=False)
-wrap_longitude(eddies_area_tracker)
+eddies_area_tracker["virtual"][:] = eddies_area_tracker["time"] == 0
+eddies_area_tracker.filled_by_interpolation(eddies_area_tracker["virtual"] == 1)
 
 # %%
 # We run a traking with default tracker
-c = Correspondances(datasets=file_objects)
+c = Correspondances(datasets=file_objects, virtual=3)
 c.track()
 c.prepare_merging()
 eddies_default_tracker = c.merge(raw_data=False)
-wrap_longitude(eddies_default_tracker)
+eddies_default_tracker["virtual"][:] = eddies_default_tracker["time"] == 0
+eddies_default_tracker.filled_by_interpolation(eddies_default_tracker["virtual"] == 1)
 
 # %%
 # Start GUI to compare tracking

notebooks/python_module/06_grid_manipulation/pet_filter.ipynb

@@ -26,7 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "from py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom matplotlib import pyplot as plt\n\n\ndef 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(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
+        "from py_eddy_tracker.dataset.grid import RegularGridDataset\nfrom py_eddy_tracker import data\nfrom matplotlib import pyplot as plt\nfrom numpy import arange\n\n\ndef 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(m, cax=ax.figure.add_axes([0.95, 0.05, 0.01, 0.9]))"
       ]
     },
     {
@@ -51,7 +51,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "Kernel\n------\nShape of kernel will increase in x when latitude increase\n\n"
+        "Kernel\n------\nShape of kernel will increase in x, when latitude increase\n\n"
       ]
     },
     {
@@ -62,7 +62,25 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(12, 8))\nfor i, latitude in enumerate((15, 35, 55, 75)):\n    k = g.kernel_bessel(latitude, 500, order=3).T\n    ax0 = plt.subplot(\n        2,\n        2,\n        i + 1,\n        title=f\"Kernel at {latitude}\u00b0 of latitude\\nfor 1/8\u00b0 grid, shape : {k.shape}\",\n        aspect=\"equal\",\n    )\n    m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap=\"viridis_r\")\nplt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))"
+        "fig = plt.figure(figsize=(12, 8))\nfor i, latitude in enumerate((15, 35, 55, 75)):\n    k = g.kernel_bessel(latitude, 500, order=3).T\n    ax0 = fig.add_subplot(\n        2,\n        2,\n        i + 1,\n        title=f\"Kernel at {latitude}\u00b0 of latitude\\nfor 1/8\u00b0 grid, shape : {k.shape}\",\n        aspect=\"equal\",\n    )\n    m = ax0.pcolormesh(k, vmin=-0.5, vmax=2, cmap=\"viridis_r\")\nplt.colorbar(m, cax=fig.add_axes((0.92, 0.05, 0.01, 0.9)))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "Kernel along latitude\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "fig = plt.figure(figsize=(12, 8))\nax = fig.add_subplot(\n    111,\n    ylabel=\"Kernel weight\",\n    xlabel=\"Latitude in \u00b0\",\n    title=\"Kernel in latitude, centered at 0\u00b0 of latitude \",\n)\nk = g.kernel_bessel(0, 500, order=3)\nk_lat = k[k.shape[0] // 2 + 1]\nnb = k_lat.shape[0] // 2\nax.plot(\n    arange(-nb * g.xstep, (nb + 0.5) * g.xstep, g.xstep), k_lat, label=\"Bessel kernel\"\n)\n\nax.legend()\nax.grid()"
       ]
     },
     {
@@ -98,7 +116,7 @@
       },
       "outputs": [],
       "source": [
-        "ax = start_axes(\"ADT low frequency\")\ng.copy(\"adt\", \"adt_low\")\ng.bessel_low_filter(\"adt_low\", 300, order=3)\nm = g.display(ax, \"adt_low\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
+        "ax = start_axes(\"ADT low frequency\")\ng.copy(\"adt\", \"adt_low_300\")\ng.bessel_low_filter(\"adt_low_300\", 300, order=3)\nm = g.display(ax, \"adt_low_300\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
       ]
     },
     {
@@ -116,7 +134,7 @@
       },
       "outputs": [],
       "source": [
-        "ax = start_axes(\"ADT high frequency\")\ng.copy(\"adt\", \"adt_high\")\ng.bessel_high_filter(\"adt_high\", 300, order=3)\nm = g.display(ax, \"adt_high\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
+        "ax = start_axes(\"ADT high frequency\")\ng.copy(\"adt\", \"adt_high_300\")\ng.bessel_high_filter(\"adt_high_300\", 300, order=3)\nm = g.display(ax, \"adt_high_300\", vmin=-0.15, vmax=0.15)\nupdate_axes(ax, m)"
       ]
     },
     {
@@ -134,7 +152,7 @@
       },
       "outputs": [],
       "source": [
-        "g.copy(\"adt\", \"adt_high_80\")\ng.bessel_high_filter(\"adt_high_80\", 80, order=3)\ng.copy(\"adt\", \"adt_low_80\")\ng.bessel_low_filter(\"adt_low_80\", 80, order=3)\n\narea = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)"
+        "g.copy(\"adt\", \"adt_high_bessel\")\ng.bessel_high_filter(\"adt_high_bessel\", 80, order=3)\ng.copy(\"adt\", \"adt_low_bessel\")\ng.bessel_low_filter(\"adt_low_bessel\", 80, order=3)\n\narea = dict(llcrnrlon=11.75, urcrnrlon=21, llcrnrlat=33, urcrnrlat=36.75)"
       ]
     },
     {
@@ -152,7 +170,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum\")\nax.set_xlabel(\"km\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat raw\")[0]\nax.loglog(*lon_spec, label=\"lon raw\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt_high_80\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat high\")[0]\nax.loglog(*lon_spec, label=\"lon high\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\"adt_low_80\", area=area)\nmappable = ax.loglog(*lat_spec, label=\"lat low\")[0]\nax.loglog(*lon_spec, label=\"lon low\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(1e-6, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
+        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum\")\nax.set_xlabel(\"km\")\n\nfor label in (\"adt_high_bessel\", \"adt_low_bessel\", \"adt\"):\n    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area)\n    mappable = ax.loglog(*lat_spec, label=f\"lat {label}\")[0]\n    ax.loglog(\n        *lon_spec, label=f\"lon {label}\", color=mappable.get_color(), linestyle=\"--\"\n    )\n\nax.set_xlim(10, 1000)\nax.set_ylim(1e-6, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
       ]
     },
     {
@@ -170,7 +188,7 @@
       },
       "outputs": [],
       "source": [
-        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum ratio\")\nax.set_xlabel(\"km\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\n    \"adt_high_80\", area=area, ref=g, ref_grid_name=\"adt\"\n)\nmappable = ax.plot(*lat_spec, label=\"lat high\")[0]\nax.plot(*lon_spec, label=\"lon high\", color=mappable.get_color(), linestyle=\"--\")\n\nlon_spec, lat_spec = g.spectrum_lonlat(\n    \"adt_low_80\", area=area, ref=g, ref_grid_name=\"adt\"\n)\nmappable = ax.plot(*lat_spec, label=\"lat low\")[0]\nax.plot(*lon_spec, label=\"lon low\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(0, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
+        "fig = plt.figure(figsize=(10, 6))\nax = fig.add_subplot(111)\nax.set_title(\"Spectrum ratio\")\nax.set_xlabel(\"km\")\n\nfor label in (\"adt_high_bessel\", \"adt_low_bessel\"):\n    lon_spec, lat_spec = g.spectrum_lonlat(label, area=area, ref=g, ref_grid_name=\"adt\")\n    mappable = ax.plot(*lat_spec, label=f\"lat {label}\")[0]\n    ax.plot(*lon_spec, label=f\"lon {label}\", color=mappable.get_color(), linestyle=\"--\")\n\nax.set_xlim(10, 1000)\nax.set_ylim(0, 1)\nax.set_xscale(\"log\")\nax.legend()\nax.grid()"
       ]
     },
     {

notebooks/python_module/08_tracking_manipulation/pet_run_a_tracking.ipynb

@@ -26,25 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "from py_eddy_tracker.data import get_remote_sample\nfrom py_eddy_tracker.tracking import Correspondances\nfrom py_eddy_tracker.featured_tracking.area_tracker import AreaTracker\nfrom numpy import where, empty\nfrom py_eddy_tracker.gui import GUI"
-      ]
-    },
-    {
-      "cell_type": "markdown",
-      "metadata": {},
-      "source": [
-        "Function to have track with contiguous longitude\n\n"
-      ]
-    },
-    {
-      "cell_type": "code",
-      "execution_count": null,
-      "metadata": {
-        "collapsed": false
-      },
-      "outputs": [],
-      "source": [
-        "def wrap_longitude(eddies):\n    lon = eddies.longitude\n    first = where(eddies.obs[\"n\"] == 0)[0]\n    nb_obs = empty(first.shape, dtype=\"u4\")\n    nb_obs[:-1] = first[1:] - first[:-1]\n    nb_obs[-1] = lon.shape[0] - first[-1]\n    lon0 = (lon[first] - 180).repeat(nb_obs)\n    lon[:] = (lon - lon0) % 360 + lon0"
+        "from py_eddy_tracker.data import get_remote_sample\nfrom py_eddy_tracker.tracking import Correspondances\nfrom py_eddy_tracker.featured_tracking.area_tracker import AreaTracker\nfrom py_eddy_tracker.gui import GUI"
       ]
     },
     {
@@ -80,7 +62,7 @@
       },
       "outputs": [],
       "source": [
-        "c = Correspondances(datasets=file_objects, class_method=AreaTracker)\nc.track()\nc.prepare_merging()\n# We have now an eddy object\neddies_area_tracker = c.merge(raw_data=False)\nwrap_longitude(eddies_area_tracker)"
+        "c = Correspondances(datasets=file_objects, class_method=AreaTracker, virtual=3)\nc.track()\nc.prepare_merging()\n# We have now an eddy object\neddies_area_tracker = c.merge(raw_data=False)\neddies_area_tracker[\"virtual\"][:] = eddies_area_tracker[\"time\"] == 0\neddies_area_tracker.filled_by_interpolation(eddies_area_tracker[\"virtual\"] == 1)"
       ]
     },
     {
@@ -98,7 +80,7 @@
       },
       "outputs": [],
       "source": [
-        "c = Correspondances(datasets=file_objects)\nc.track()\nc.prepare_merging()\neddies_default_tracker = c.merge(raw_data=False)\nwrap_longitude(eddies_default_tracker)"
+        "c = Correspondances(datasets=file_objects, virtual=3)\nc.track()\nc.prepare_merging()\neddies_default_tracker = c.merge(raw_data=False)\neddies_default_tracker[\"virtual\"][:] = eddies_default_tracker[\"time\"] == 0\neddies_default_tracker.filled_by_interpolation(eddies_default_tracker[\"virtual\"] == 1)"
       ]
     },
     {

src/py_eddy_tracker/data/__init__.py

@@ -11,8 +11,9 @@ def get_path(name):
 
 def get_remote_sample(path):
     url = (
-    f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}.tar.xz"
+        f"https://github.com/AntSimi/py-eddy-tracker-sample-id/raw/master/{path}.tar.xz"
     )
+
     content = requests.get(url).content
 
     # Tar module could manage lzma tar, but it will apply un compress for each extractfile