Merge pull request AntSimi#77 from CoriPegliasco/master

add ellipse fit

Author: CoriPegliasco

examples/01_general_things/pet_storage.py

@@ -32,7 +32,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

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
 
 import py_eddy_tracker.gui
 from py_eddy_tracker.appli.gui import Anim
 from py_eddy_tracker.data import get_demo_path
 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"])

src/py_eddy_tracker/dataset/grid.py

@@ -609,17 +609,17 @@ def eddy_identification(
         :param str vname: Grid name of v speed component
         :param datetime.datetime date: Date which will be stored in object to date data
         :param float,int step: Height between two layers in m
-        :param float,int shape_error: Maximal error allowed for outter contour in %
+        :param float,int shape_error: Maximal error allowed for outermost contour in %
         :param int sampling: Number of points to store contours and speed profile
         :param str sampling_method: Method to resample 'uniform' or 'visvalingam'
         :param (int,int),None pixel_limit:
-            Min and max number of pixels inside the inner and the outer contour to be considered as an eddy
+            Min and max number of pixels inside the inner and the outermost contour to be considered as an eddy
         :param float,None precision: Truncate values at the defined precision in m
         :param str force_height_unit: Unit used for height unit
         :param str force_speed_unit: Unit used for speed unit
         :param dict kwargs: Argument given to amplitude
 
-        :return: Return a list of 2 elements: Anticyclone and Cyclone
+        :return: Return a list of 2 elements: Anticyclones and Cyclones
         :rtype: py_eddy_tracker.observations.observation.EddiesObservations
 
         .. minigallery:: py_eddy_tracker.GridDataset.eddy_identification
@@ -729,7 +729,8 @@ def eddy_identification(
                 for contour in contour_paths:
                     if contour.used:
                         continue
-                    # FIXME : center could be not in contour and fit on raw sampling
+                    # FIXME : center could be outside the contour due to the fit
+                    # FIXME : warning : the fit is made on raw sampling
                     _, _, _, aerr = contour.fit_circle()
 
                     # Filter for shape
@@ -752,6 +753,7 @@ def eddy_identification(
                     ):
                         continue
 
+                    # Test the number of pixels within the outermost contour
                     # FIXME : Maybe limit max must be replace with a maximum of surface
                     if (
                         contour.nb_pixel < pixel_limit[0]
@@ -760,6 +762,9 @@ def eddy_identification(
                         contour.reject = 3
                         continue
 
+                    # Here the considered contour passed shape_error test, masked_pixels test,
+                    # values strictly above (AEs) or below (CEs) the contour, number_pixels test)
+
                     # Compute amplitude
                     reset_centroid, amp = self.get_amplitude(
                         contour,

src/py_eddy_tracker/eddy_feature.py

@@ -31,7 +31,7 @@
 class Amplitude(object):
     """
     Class to calculate *amplitude* and counts of *local maxima/minima*
-    within a closed region of a sea level anomaly field.
+    within a closed region of a sea surface height field.
     """
 
     EPSILON = 1e-8
@@ -66,8 +66,8 @@ def __init__(
         :param array data:
         :param float interval:
         :param int mle: maximum number of local maxima in contour
-        :param int nb_step_min: number of interval to consider like an eddy
-        :param int nb_step_to_be_mle: number of interval to be consider like another maxima
+        :param int nb_step_min: number of intervals to consider an eddy
+        :param int nb_step_to_be_mle: number of intervals to be considered as an another maxima
         """
 
         # Height of the contour
@@ -102,8 +102,9 @@ def __init__(
         self.nb_pixel = i_x.shape[0]
 
         # Only pixel in contour
+        # FIXME : change sla by ssh as the grid can be adt?
         self.sla = data[contour.pixels_index]
-        # Amplitude which will be provide
+        # Amplitude which will be provided
         self.amplitude = 0
         # Maximum local extrema accepted
         self.mle = mle
@@ -115,9 +116,10 @@ def within_amplitude_limits(self):
     def all_pixels_below_h0(self, level):
         """
         Check CSS11 criterion 1: The SSH values of all of the pixels
-        are below a given SSH threshold for cyclonic eddies.
+        are below (above) a given SSH threshold for cyclonic (anticyclonic)
+        eddies.
         """
-        # In some case pixel value must be very near of contour bounds
+        # In some cases pixel value may be very close to the contour bounds
         if self.sla.mask.any() or ((self.sla.data - self.h_0) > self.EPSILON).any():
             return False
         else:
@@ -293,10 +295,10 @@ class Contours(object):
 
     Attributes:
       contour:
-        A matplotlib contour object of high-pass filtered SLA
+        A matplotlib contour object of high-pass filtered SSH
 
       eddy:
-        A tracklist object holding the SLA data
+        A tracklist object holding the SSH data
 
       grd:
         A grid object
@@ -406,7 +408,7 @@ def __init__(self, x, y, z, levels, wrap_x=False, keep_unclose=False):
         fig = Figure()
         ax = fig.add_subplot(111)
         if wrap_x:
-            logger.debug("wrapping activate to compute contour")
+            logger.debug("wrapping activated to compute contour")
             x = concatenate((x, x[:1] + 360))
             z = ma.concatenate((z, z[:1]))
         logger.debug("X shape : %s", x.shape)
@@ -602,8 +604,8 @@ def display(
             Must be 'shape_error', 'x', 'y' or 'radius'.
             If define display_criterion is not use.
             bins argument must be define
-        :param array bins: bins use to colorize contour
-        :param str cmap: Name of cmap to use for field display
+        :param array bins: bins used to colorize contour
+        :param str cmap: Name of cmap for field display
         :param dict kwargs: look at :py:meth:`matplotlib.collections.LineCollection`
 
         .. minigallery:: py_eddy_tracker.Contours.display
@@ -688,7 +690,7 @@ def display(
             ax.autoscale_view()
 
     def label_contour_unused_which_contain_eddies(self, eddies):
-        """Select contour which contain several eddies"""
+        """Select contour containing several eddies"""
         if eddies.sign_type == 1:
             # anticyclonic
             sl = slice(None, -1)

src/py_eddy_tracker/featured_tracking/old_tracker_reference.py

@@ -21,13 +21,13 @@ def cost_function(records_in, records_out, distance):
         return distance
 
     def mask_function(self, other, distance):
-        """We mask link with ellips and ratio"""
-        # Compute Parameter of ellips
+        """We mask link with ellipse and ratio"""
+        # Compute Parameter of ellipse
         minor, major = 1.05, 1.5
-        y = self.basic_formula_ellips_major_axis(
+        y = self.basic_formula_ellipse_major_axis(
             self.lat, degrees=True, c0=minor, cmin=minor, cmax=major, lat1=23, lat2=5
         )
-        # mask from ellips
+        # mask from ellipse
         mask = self.shifted_ellipsoid_degrees_mask(
             other, minor=minor, major=y  # Minor can be bigger than major??
         )

src/py_eddy_tracker/observations/network.py

@@ -513,7 +513,6 @@ def relatives(self, obs, order=2):
             else:
                 segments_connexion[seg][0] = i_slice
 
-
             if i_p != -1:
 
                 if p_seg not in segments_connexion:
@@ -531,12 +530,7 @@ def relatives(self, obs, order=2):
                 segments_connexion[seg][1].append(n_seg)
                 segments_connexion[n_seg][1].append(seg)
 
-
-        i_obs = (
-            [obs]
-            if not hasattr(obs, "__iter__")
-            else obs
-        )
+        i_obs = [obs] if not hasattr(obs, "__iter__") else obs
         import numpy as np
 
         distance = zeros(segment.size, dtype=np.uint16) - 1