Merge branch 'master' into visvalingam

Author: AntSimi

examples/01_general_things/pet_storage.py

@@ -16,8 +16,10 @@
 """
 
 import py_eddy_tracker_sample
+from matplotlib import pyplot as plt
+from numpy import arange, outer
 
-from py_eddy_tracker.data import get_demo_path, get_remote_demo_sample
+from py_eddy_tracker.data import get_demo_path
 from py_eddy_tracker.observations.network import NetworkObservations
 from py_eddy_tracker.observations.observation import EddiesObservations, Table
 from py_eddy_tracker.observations.tracking import TrackEddiesObservations
@@ -32,7 +34,7 @@
 # array field like contour/profile are 2D column.
 
 # %%
-# Eddies files (zarr or netcdf) could be loaded with ```load_file``` method:
+# Eddies files (zarr or netcdf) can be loaded with ```load_file``` method:
 eddies_collections = EddiesObservations.load_file(get_demo_path("Cyclonic_20160515.nc"))
 eddies_collections.field_table()
 # offset and scale_factor are used only when data is stored in zarr or netCDF4
@@ -58,6 +60,53 @@
 # ---------------
 # All contours are stored on the same number of points, and are resampled if needed with an algorithm to be stored as objects
 
+# %%
+# Speed profile storage
+# ---------------------
+# Speed profile is an interpolation of speed mean along each contour.
+# For each contour included in eddy, we compute mean of speed along the contour,
+# and after we interpolate speed mean array on a fixed size array.
+#
+# Several field are available to understand "uavg_profile" :
+#  0. - num_contours : Number of contour in eddies, must be equal to amplitude divide by isoline step
+#  1. - height_inner_contour : height of inner contour used
+#  2. - height_max_speed_contour : height of max speed contour used
+#  3. - height_external_contour : height of outter contour used
+#
+# Last value of "uavg_profile" is for inner contour and first value for outter contour.
+
+# Observations selection of "uavg_profile" with high number of contour(Eddy with high amplitude)
+e = eddies_collections.extract_with_mask(eddies_collections.num_contours > 15)
+
+# %%
+
+# Raw display of profiles with more than 15 contours
+ax = plt.subplot(111)
+_ = ax.plot(e.uavg_profile.T, lw=0.5)
+
+# %%
+
+# Profile from inner to outter
+ax = plt.subplot(111)
+ax.plot(e.uavg_profile[:, ::-1].T, lw=0.5)
+_ = ax.set_xlabel("From inner to outter contour"), ax.set_ylabel("Speed (m/s)")
+
+# %%
+
+# If we normalize indice of contour to set speed contour to 1 and inner contour to 0
+ax = plt.subplot(111)
+h_in = e.height_inner_contour
+h_s = e.height_max_speed_contour
+h_e = e.height_external_contour
+r = (h_e - h_in) / (h_s - h_in)
+nb_pt = e.uavg_profile.shape[1]
+# Create an x array for each profile
+x = outer(arange(nb_pt) / nb_pt, r)
+
+ax.plot(x, e.uavg_profile[:, ::-1].T, lw=0.5)
+_ = ax.set_xlabel("From inner to outter contour"), ax.set_ylabel("Speed (m/s)")
+
+
 # %%
 # Trajectories
 # ------------
@@ -86,11 +135,7 @@
 # - next_obs : Index of the next observation in the full dataset, if -1 there are no next observation (the segment ends)
 # - previous_cost : Result of the cost function (1 is a good association, 0 is bad) with previous observation
 # - next_cost : Result of the cost function (1 is a good association, 0 is bad) with next observation
-eddies_network = NetworkObservations.load_file(
-    get_remote_demo_sample(
-        "eddies_med_adt_allsat_dt2018_err70_filt500_order1/Anticyclonic_network.nc"
-    )
-)
+eddies_network = NetworkObservations.load_file(get_demo_path("network_med.nc"))
 eddies_network.field_table()
 
 # %%

examples/14_generic_tools/pet_fit_contour.py

@@ -15,30 +15,30 @@
 from py_eddy_tracker import data
 from py_eddy_tracker.generic import coordinates_to_local, local_to_coordinates
 from py_eddy_tracker.observations.observation import EddiesObservations
-from py_eddy_tracker.poly import fit_circle_, fit_ellips
+from py_eddy_tracker.poly import fit_circle_, fit_ellipse
 
 # %%
 # Load example identification file
 a = EddiesObservations.load_file(data.get_demo_path("Anticyclonic_20190223.nc"))
 
 
 # %%
-# Function to draw circle or ellips from parameter
+# Function to draw circle or ellipse from parameter
 def build_circle(x0, y0, r):
     angle = radians(linspace(0, 360, 50))
     x_norm, y_norm = cos(angle), sin(angle)
     return local_to_coordinates(x_norm * r, y_norm * r, x0, y0)
 
 
-def build_ellips(x0, y0, a, b, theta):
+def build_ellipse(x0, y0, a, b, theta):
     angle = radians(linspace(0, 360, 50))
     x = a * cos(theta) * cos(angle) - b * sin(theta) * sin(angle)
     y = a * sin(theta) * cos(angle) + b * cos(theta) * sin(angle)
     return local_to_coordinates(x, y, x0, y0)
 
 
 # %%
-# Plot fitted circle or ellips on stored contour
+# Plot fitted circle or ellipse on stored contour
 xs, ys = a.contour_lon_s, a.contour_lat_s
 
 fig = plt.figure(figsize=(15, 15))
@@ -51,9 +51,9 @@ def build_ellips(x0, y0, a, b, theta):
     ax = fig.add_subplot(4, 4, j)
     ax.grid(), ax.set_aspect("equal")
     ax.plot(x, y, label="store", color="black")
-    x0, y0, a, b, theta = fit_ellips(x_, y_)
+    x0, y0, a, b, theta = fit_ellipse(x_, y_)
     x0, y0 = local_to_coordinates(x0, y0, x0_, y0_)
-    ax.plot(*build_ellips(x0, y0, a, b, theta), label="ellips", color="green")
+    ax.plot(*build_ellipse(x0, y0, a, b, theta), label="ellipse", color="green")
     x0, y0, radius, shape_error = fit_circle_(x_, y_)
     x0, y0 = local_to_coordinates(x0, y0, x0_, y0_)
     ax.plot(*build_circle(x0, y0, radius), label="circle", color="red", lw=0.5)

examples/16_network/pet_follow_particle.py

@@ -33,7 +33,7 @@ def _repr_html_(self, *args, **kwargs):
 
     def save(self, *args, **kwargs):
         if args[0].endswith("gif"):
-            # In this case gif is use to create thumbnail which are not use but consume same time than video
+            # In this case gif is used to create thumbnail which are not used but consumes same time than video
             # So we create an empty file, to save time
             with open(args[0], "w") as _:
                 pass
@@ -147,7 +147,7 @@ def particle_candidate(x, y, c, eddies, t_start, i_target, pct, **kwargs):
     e = eddies.extract_with_mask(m_start)
     # to be able to get global index
     translate_start = where(m_start)[0]
-    # Identify particle in eddies(only in core)
+    # Identify particle in eddies (only in core)
     i_start = e.contains(x, y, intern=True)
     m = i_start != -1
     x, y, i_start = x[m], y[m], i_start[m]
@@ -158,9 +158,9 @@ def particle_candidate(x, y, c, eddies, t_start, i_target, pct, **kwargs):
     e_end = eddies.extract_with_mask(m_end)
     # to be able to get global index
     translate_end = where(m_end)[0]
-    # Id eddies for each alive particle(in core and extern)
+    # Id eddies for each alive particle (in core and extern)
     i_end = e_end.contains(x, y)
-    # compute matrix and filled target array
+    # compute matrix and fill target array
     get_matrix(i_start, i_end, translate_start, translate_end, i_target, pct)
 
 
@@ -169,7 +169,7 @@ def get_matrix(i_start, i_end, translate_start, translate_end, i_target, pct):
     nb_start, nb_end = translate_start.size, translate_end.size
     # Matrix which will store count for every couple
     count = zeros((nb_start, nb_end), dtype=nb_types.int32)
-    # Number of particle in each origin observation
+    # Number of particles in each origin observation
     ref = zeros(nb_start, dtype=nb_types.int32)
     # For each particle
     for i in range(i_start.size):
@@ -181,7 +181,7 @@ def get_matrix(i_start, i_end, translate_start, translate_end, i_target, pct):
     for i in range(nb_start):
         for j in range(nb_end):
             pct_ = count[i, j]
-            # If there are particle from i to j
+            # If there are particles from i to j
             if pct_ != 0:
                 # Get percent
                 pct_ = pct_ / ref[i] * 100.0

examples/16_network/pet_ioannou_2017_case.py

@@ -14,15 +14,19 @@
 from matplotlib import pyplot as plt
 from matplotlib.animation import FuncAnimation
 from matplotlib.ticker import FuncFormatter
-from numpy import arange, where
+from numpy import arange, where, array, pi
 
 from py_eddy_tracker.appli.gui import Anim
 from py_eddy_tracker.data import get_demo_path
 from py_eddy_tracker.gui import GUI_AXES
 from py_eddy_tracker.observations.network import NetworkObservations
 
+from py_eddy_tracker.generic import coordinates_to_local, local_to_coordinates
+from py_eddy_tracker.poly import fit_ellipse
 
 # %%
+
+
 class VideoAnimation(FuncAnimation):
     def _repr_html_(self, *args, **kwargs):
         """To get video in html and have a player"""
@@ -188,3 +192,47 @@ def update_axes(ax, mappable=None):
 m = close_to_i3.scatter_timeline(ax, "shape_error_e", vmin=14, vmax=70, **kw)
 cb = update_axes(ax, m["scatter"])
 cb.set_label("Effective shape error")
+
+# %%
+# Rotation angle
+# --------------
+# For each obs, fit an ellipse to the contour, with theta the angle from the x-axis,
+# a the semi ax in x direction and b the semi ax in y dimension
+
+theta_ = list()
+a_ = list()
+b_ = list()
+for obs in close_to_i3:
+    x, y = obs["contour_lon_s"], obs["contour_lat_s"]
+    x0_, y0_ = x.mean(), y.mean()
+    x_, y_ = coordinates_to_local(x, y, x0_, y0_)
+    x0, y0, a, b, theta = fit_ellipse(x_, y_)
+    theta_.append(theta)
+    a_.append(a)
+    b_.append(b)
+a_ = array(a_)
+b_ = array(b_)
+
+# %%
+# Theta
+ax = timeline_axes()
+m = close_to_i3.scatter_timeline(ax, theta_, vmin=-pi / 2, vmax=pi / 2, cmap="hsv")
+cb = update_axes(ax, m["scatter"])
+
+# %%
+# a
+ax = timeline_axes()
+m = close_to_i3.scatter_timeline(ax, a_ * 1e-3, vmin=0, vmax=80, cmap="Spectral_r")
+cb = update_axes(ax, m["scatter"])
+
+# %%
+# b
+ax = timeline_axes()
+m = close_to_i3.scatter_timeline(ax, b_ * 1e-3, vmin=0, vmax=80, cmap="Spectral_r")
+cb = update_axes(ax, m["scatter"])
+
+# %%
+# a/b
+ax = timeline_axes()
+m = close_to_i3.scatter_timeline(ax, a_ / b_, vmin=1, vmax=2, cmap="Spectral_r")
+cb = update_axes(ax, m["scatter"])

notebooks/python_module/01_general_things/pet_storage.ipynb

@@ -26,7 +26,7 @@
       },
       "outputs": [],
       "source": [
-        "import py_eddy_tracker_sample\n\nfrom py_eddy_tracker.data import get_demo_path, get_remote_demo_sample\nfrom py_eddy_tracker.observations.network import NetworkObservations\nfrom py_eddy_tracker.observations.observation import EddiesObservations, Table\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations"
+        "import py_eddy_tracker_sample\nfrom matplotlib import pyplot as plt\nfrom numpy import arange, outer\n\nfrom py_eddy_tracker.data import get_demo_path\nfrom py_eddy_tracker.observations.network import NetworkObservations\nfrom py_eddy_tracker.observations.observation import EddiesObservations, Table\nfrom py_eddy_tracker.observations.tracking import TrackEddiesObservations"
       ]
     },
     {
@@ -108,6 +108,57 @@
         "## Contour storage\nAll contours are stored on the same number of points, and are resampled if needed with an algorithm to be stored as objects\n\n"
       ]
     },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## Speed profile storage\nSpeed profile is an interpolation of speed mean along each contour.\nFor each contour included in eddy, we compute mean of speed along the contour,\nand after we interpolate speed mean array on a fixed size array.\n\nSeveral field are available to understand \"uavg_profile\" :\n 0. - num_contours : Number of contour in eddies, must be equal to amplitude divide by isoline step\n 1. - height_inner_contour : height of inner contour used\n 2. - height_max_speed_contour : height of max speed contour used\n 3. - height_external_contour : height of outter contour used\n\nLast value of \"uavg_profile\" is for inner contour and first value for outter contour.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Observations selection of \"uavg_profile\" with high number of contour(Eddy with high amplitude)\ne = eddies_collections.extract_with_mask(eddies_collections.num_contours > 15)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Raw display of profiles with more than 15 contours\nax = plt.subplot(111)\n_ = ax.plot(e.uavg_profile.T, lw=0.5)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Profile from inner to outter\nax = plt.subplot(111)\nax.plot(e.uavg_profile[:, ::-1].T, lw=0.5)\n_ = ax.set_xlabel(\"From inner to outter contour\"), ax.set_ylabel(\"Speed (m/s)\")"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# If we normalize indice of contour to set speed contour to 1 and inner contour to 0\nax = plt.subplot(111)\nh_in = e.height_inner_contour\nh_s = e.height_max_speed_contour\nh_e = e.height_external_contour\nr = (h_e - h_in) / (h_s - h_in)\nnb_pt = e.uavg_profile.shape[1]\n# Create an x array for each profile\nx = outer(arange(nb_pt) / nb_pt, r)\n\nax.plot(x, e.uavg_profile[:, ::-1].T, lw=0.5)\n_ = ax.set_xlabel(\"From inner to outter contour\"), ax.set_ylabel(\"Speed (m/s)\")"
+      ]
+    },
     {
       "cell_type": "markdown",
       "metadata": {},
@@ -141,7 +192,7 @@
       },
       "outputs": [],
       "source": [
-        "eddies_network = NetworkObservations.load_file(\n    get_remote_demo_sample(\n        \"eddies_med_adt_allsat_dt2018_err70_filt500_order1/Anticyclonic_network.nc\"\n    )\n)\neddies_network.field_table()"
+        "eddies_network = NetworkObservations.load_file(get_demo_path(\"network_med.nc\"))\neddies_network.field_table()"
       ]
     },
     {
@@ -179,7 +230,7 @@
       "name": "python",
       "nbconvert_exporter": "python",
       "pygments_lexer": "ipython3",
-      "version": "3.7.7"
+      "version": "3.9.2"
     }
   },
   "nbformat": 4,

src/py_eddy_tracker/appli/network.py

@@ -37,7 +37,7 @@ def build_network():
 
 
 def divide_network():
-    parser = EddyParser("Separate path for a same group(network)")
+    parser = EddyParser("Separate path for a same group (network)")
     parser.add_argument("input", help="input network file")
     parser.add_argument("out", help="output file")
     parser.contour_intern_arg()
@@ -66,7 +66,7 @@ def subset_network():
         "--length",
         nargs=2,
         type=int,
-        help="Nb of day which must be cover by network, first minimum number of day and last maximum number of day,"
+        help="Nb of days that must be covered by the network, first minimum number of day and last maximum number of day,"
         "if value is negative, this bound won't be used",
     )
     parser.add_argument(
@@ -85,8 +85,8 @@ def subset_network():
         "--period",
         nargs=2,
         type=int,
-        help="Start day and end day, if it's negative value we will add to day min and add to day max,"
-        "if 0 it s not use",
+        help="Start day and end day, if it's a negative value we will add to day min and add to day max,"
+        "if 0 it is not used",
     )
     args = parser.parse_args()
     n = NetworkObservations.load_file(args.input, raw_data=True)