coarser resolution in FSLE example

Author: AntSimi

examples/07_cube_manipulation/pet_fsle_med.py

@@ -2,7 +2,7 @@
 FSLE experiment in med
 ======================
 
-Example to build FSLE, parameter values must be adapted for your case.
+Example to build Finite Size Lyapunov Exponents, parameter values must be adapted for your case.
 
 Example use a method similar to `AVISO flse`_
 
@@ -25,7 +25,7 @@
 
 # %%
 # ADT in med
-# ---------------------
+# ----------
 c = GridCollection.from_netcdf_cube(
     get_path("dt_med_allsat_phy_l4_2005T2.nc"),
     "longitude",
@@ -36,7 +36,7 @@
 
 
 # %%
-# Methods to compute fsle
+# Methods to compute FSLE
 # -----------------------
 @njit(cache=True, fastmath=True)
 def check_p(x, y, g, m, dt, dist_init=0.02, dist_max=0.6):
@@ -86,13 +86,33 @@ def build_triplet(x, y, step=0.02):
     return x_, y_
 
 
+# %%
+# Settings
+# --------
+
+# Step in degrees for ouput
+step_grid_out = 1 / 25.0
+# Initial separation in degrees
+dist_init = 1 / 50.0
+# Final separation in degrees
+dist_max = 1 / 5.0
+# Time of start
+t0 = 20268
+# Number of time step by days
+time_step_by_days = 5
+# Maximal time of advection
+# Here we limit because our data cube cover only 3 month
+nb_days = 85
+# Backward or forward
+backward = True
+
 # %%
 # Particles
 # ---------
-step = 0.02
-t0 = 20268
 x0_, y0_ = -5, 30
-lon_p, lat_p = arange(x0_, x0_ + 43, step), arange(y0_, y0_ + 16, step)
+lon_p, lat_p = arange(x0_, x0_ + 43, step_grid_out), arange(
+    y0_, y0_ + 16, step_grid_out
+)
 x0, y0 = meshgrid(lon_p, lat_p)
 grid_shape = x0.shape
 x0, y0 = x0.reshape(-1), y0.reshape(-1)
@@ -103,26 +123,25 @@ def build_triplet(x, y, step=0.02):
 # %%
 # FSLE
 # ----
-time_step_by_days = 5
+
 # Array to compute fsle
 fsle = zeros(x0.shape[0], dtype="f4")
-x, y = build_triplet(x0, y0)
+x, y = build_triplet(x0, y0, dist_init)
 used = zeros(x.shape[0], dtype="bool")
 
 # advection generator
-kw = dict(t_init=t0, nb_step=1, backward=True, mask_particule=used)
+kw = dict(t_init=t0, nb_step=1, backward=backward, mask_particule=used)
 p = c.advect(x, y, "u", "v", time_step=86400 / time_step_by_days, **kw)
 
-nb_days = 85
 # We check at each step of advection if particle distance is over `dist_max`
 for i in range(time_step_by_days * nb_days):
     t, xt, yt = p.__next__()
     dt = t / 86400.0 - t0
-    check_p(xt, yt, fsle, used, dt, dist_max=0.2, dist_init=step)
+    check_p(xt, yt, fsle, used, dt, dist_max=dist_max, dist_init=dist_init)
 
 # Get index with original_position
-i = ((x0 - x0_) / step).astype("i4")
-j = ((y0 - y0_) / step).astype("i4")
+i = ((x0 - x0_) / step_grid_out).astype("i4")
+j = ((y0 - y0_) / step_grid_out).astype("i4")
 fsle_ = empty(grid_shape, dtype="f4")
 used_ = zeros(grid_shape, dtype="bool")
 fsle_[j, i] = fsle

notebooks/python_module/07_cube_manipulation/pet_fsle_med.ipynb

@@ -15,7 +15,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "\n# FSLE experiment in med\n\nExample to build FSLE, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n    https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n    fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
+        "\n# FSLE experiment in med\n\nExample to build Finite Size Lyapunov Exponents, parameter values must be adapted for your case.\n\nExample use a method similar to `AVISO flse`_\n\n    https://www.aviso.altimetry.fr/en/data/products/value-added-products/\n    fsle-finite-size-lyapunov-exponents/fsle-description.html\n"
       ]
     },
     {
@@ -51,7 +51,7 @@
       "cell_type": "markdown",
       "metadata": {},
       "source": [
-        "## Methods to compute fsle\n\n"
+        "## Methods to compute FSLE\n\n"
       ]
     },
     {
@@ -65,6 +65,24 @@
         "@njit(cache=True, fastmath=True)\ndef check_p(x, y, g, m, dt, dist_init=0.02, dist_max=0.6):\n    \"\"\"\n    Check if distance between eastern or northern particle to center particle is bigger than `dist_max`\n    \"\"\"\n    nb_p = x.shape[0] // 3\n    delta = dist_max ** 2\n    for i in range(nb_p):\n        i0 = i * 3\n        i_n = i0 + 1\n        i_e = i0 + 2\n        # If particle already set, we skip\n        if m[i0] or m[i_n] or m[i_e]:\n            continue\n        # Distance with north\n        dxn, dyn = x[i0] - x[i_n], y[i0] - y[i_n]\n        dn = dxn ** 2 + dyn ** 2\n        # Distance with east\n        dxe, dye = x[i0] - x[i_e], y[i0] - y[i_e]\n        de = dxe ** 2 + dye ** 2\n\n        if dn >= delta or de >= delta:\n            s1 = dxe ** 2 + dxn ** 2 + dye ** 2 + dyn ** 2\n            s2 = ((dxn + dye) ** 2 + (dxe - dyn) ** 2) * (\n                (dxn - dye) ** 2 + (dxe + dyn) ** 2\n            )\n            g[i] = 1 / (2 * dt) * log2(1 / (2 * dist_init ** 2) * (s1 + s2 ** 0.5))\n            m[i0], m[i_n], m[i_e] = True, True, True\n\n\n@njit(cache=True)\ndef build_triplet(x, y, step=0.02):\n    \"\"\"\n    Triplet building for each position we add east and north point with defined step\n    \"\"\"\n    nb_x = x.shape[0]\n    x_ = empty(nb_x * 3, dtype=x.dtype)\n    y_ = empty(nb_x * 3, dtype=y.dtype)\n    for i in range(nb_x):\n        i0 = i * 3\n        i_n, i_e = i0 + 1, i0 + 2\n        x__, y__ = x[i], y[i]\n        x_[i0], y_[i0] = x__, y__\n        x_[i_n], y_[i_n] = x__, y__ + step\n        x_[i_e], y_[i_e] = x__ + step, y__\n    return x_, y_"
       ]
     },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Settings\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Step in degrees for ouput\nstep_grid_out = 1/25.\n# Initial separation in degrees\ndist_init = 1/50.\n# Final separation in degrees\ndist_max = 1/5.\n# Time of start\nt0 = 20268\n# Number of time step by days\ntime_step_by_days = 5\n# Maximal time of advection\n# Here we limit because our data cube cover only 3 month\nnb_days = 85\n# Backward or forward\nbackward=True"
+      ]
+    },
     {
       "cell_type": "markdown",
       "metadata": {},
@@ -80,7 +98,7 @@
       },
       "outputs": [],
       "source": [
-        "step = 0.02\nt0 = 20268\nx0_, y0_ = -5, 30\nlon_p, lat_p = arange(x0_, x0_ + 43, step), arange(y0_, y0_ + 16, step)\nx0, y0 = meshgrid(lon_p, lat_p)\ngrid_shape = x0.shape\nx0, y0 = x0.reshape(-1), y0.reshape(-1)\n# Identify all particle not on land\nm = ~isnan(c[t0].interp(\"adt\", x0, y0))\nx0, y0 = x0[m], y0[m]"
+        "x0_, y0_ = -5, 30\nlon_p, lat_p = arange(x0_, x0_ + 43, step_grid_out), arange(y0_, y0_ + 16, step_grid_out)\nx0, y0 = meshgrid(lon_p, lat_p)\ngrid_shape = x0.shape\nx0, y0 = x0.reshape(-1), y0.reshape(-1)\n# Identify all particle not on land\nm = ~isnan(c[t0].interp(\"adt\", x0, y0))\nx0, y0 = x0[m], y0[m]"
       ]
     },
     {
@@ -98,7 +116,7 @@
       },
       "outputs": [],
       "source": [
-        "time_step_by_days = 5\n# Array to compute fsle\nfsle = zeros(x0.shape[0], dtype=\"f4\")\nx, y = build_triplet(x0, y0)\nused = zeros(x.shape[0], dtype=\"bool\")\n\n# advection generator\nkw = dict(t_init=t0, nb_step=1, backward=True, mask_particule=used)\np = c.advect(x, y, \"u\", \"v\", time_step=86400 / time_step_by_days, **kw)\n\nnb_days = 85\n# We check at each step of advection if particle distance is over `dist_max`\nfor i in range(time_step_by_days * nb_days):\n    t, xt, yt = p.__next__()\n    dt = t / 86400.0 - t0\n    check_p(xt, yt, fsle, used, dt, dist_max=0.2, dist_init=step)\n\n# Get index with original_position\ni = ((x0 - x0_) / step).astype(\"i4\")\nj = ((y0 - y0_) / step).astype(\"i4\")\nfsle_ = empty(grid_shape, dtype=\"f4\")\nused_ = zeros(grid_shape, dtype=\"bool\")\nfsle_[j, i] = fsle\nused_[j, i] = used[::3]\n# Create a grid object\nfsle_custom = RegularGridDataset.with_array(\n    coordinates=(\"lon\", \"lat\"),\n    datas=dict(\n        fsle=ma.array(fsle_.T, mask=~used_.T),\n        lon=lon_p,\n        lat=lat_p,\n    ),\n    centered=True,\n)"
+        "# Array to compute fsle\nfsle = zeros(x0.shape[0], dtype=\"f4\")\nx, y = build_triplet(x0, y0, dist_init)\nused = zeros(x.shape[0], dtype=\"bool\")\n\n# advection generator\nkw = dict(t_init=t0, nb_step=1, backward=backward, mask_particule=used)\np = c.advect(x, y, \"u\", \"v\", time_step=86400 / time_step_by_days, **kw)\n\n# We check at each step of advection if particle distance is over `dist_max`\nfor i in range(time_step_by_days * nb_days):\n    t, xt, yt = p.__next__()\n    dt = t / 86400.0 - t0\n    check_p(xt, yt, fsle, used, dt, dist_max=dist_max, dist_init=dist_init)\n\n# Get index with original_position\ni = ((x0 - x0_) / step_grid_out).astype(\"i4\")\nj = ((y0 - y0_) / step_grid_out).astype(\"i4\")\nfsle_ = empty(grid_shape, dtype=\"f4\")\nused_ = zeros(grid_shape, dtype=\"bool\")\nfsle_[j, i] = fsle\nused_[j, i] = used[::3]\n# Create a grid object\nfsle_custom = RegularGridDataset.with_array(\n    coordinates=(\"lon\", \"lat\"),\n    datas=dict(\n        fsle=ma.array(fsle_.T, mask=~used_.T),\n        lon=lon_p,\n        lat=lat_p,\n    ),\n    centered=True,\n)"
       ]
     },
     {