amplitude simplification

Author: AntSimi

src/py_eddy_tracker/dataset/grid.py

@@ -32,14 +32,6 @@ def raw_resample(datas, fixed_size):
     return interp(arange(fixed_size), arange(nb_value) * (fixed_size - 1) / (nb_value - 1) , datas)
 
 
-def contour_iter(self, anticyclonic_search):
-    for coll in self.collections[::1 if anticyclonic_search else -1]:
-        yield coll
-
-
-BaseQuadContourSet.iter_ = contour_iter
-
-
 @property
 def mean_coordinates(self):
     return self.vertices.mean(axis=0)
@@ -100,8 +92,8 @@ def uniform_resample(x_val, y_val, num_fac=2, fixed_size=None):
     dist = empty(x_val.shape)
     dist[0] = 0
     dist[1:] = distance(x_val[:-1], y_val[:-1], x_val[1:], y_val[1:])
-    # To be still monotonous
-    dist[1:][dist[1:]<1e-10] = 1e-10
+    # To be still monotonous (dist is store in m)
+    dist[1:][dist[1:]<1e-3] = 1e-3
     dist = dist.cumsum()
     # Get uniform distances
     if fixed_size is None:
@@ -121,13 +113,21 @@ def uniform_resample_stack(vertices, num_fac=2, fixed_size=None):
     data[:, 1] = y_new
     return data
 
+
 @njit(cache=True)
 def value_on_regular_contour(x_g, y_g, z_g, m_g, vertices, num_fac=2, fixed_size=None):
     x_val, y_val = vertices[:, 0], vertices[:, 1]
     x_new, y_new = uniform_resample(x_val, y_val, num_fac, fixed_size)
     return interp2d_geo(x_g, y_g, z_g, m_g, x_new[1:], y_new[1:])
 
 
+@njit(cache=True)
+def mean_on_regular_contour(x_g, y_g, z_g, m_g, vertices, num_fac=2, fixed_size=None):
+    x_val, y_val = vertices[:, 0], vertices[:, 1]
+    x_new, y_new = uniform_resample(x_val, y_val, num_fac, fixed_size)
+    return interp2d_geo(x_g, y_g, z_g, m_g, x_new[1:], y_new[1:]).mean()
+
+
 def fit_circle_path(self):
     if not hasattr(self, '_circle_params'):
         self._circle_params = _fit_circle_path(self.vertices)
@@ -715,7 +715,7 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     c_x, c_y = proj(speed_contour.lon, speed_contour.lat)
                     _, _, eddy_radius_s, aerr_s = fit_circle_c_numba(c_x, c_y)
 
-                    # Instantiate new EddyObservation object
+                    # Instantiate new EddyObservation object (high cost need to be review)
                     properties = EddiesObservations(size=1, track_extra_variables=track_extra_variables,
                                                     track_array_variables=array_sampling,
                                                     array_variables=array_variables)
@@ -780,7 +780,7 @@ def get_uavg(self, all_contours, centlon_e, centlat_e, original_contour, anticyc
         Calculate geostrophic speed around successive contours
         Returns the average
         """
-        max_average_speed = self.speed_coef(original_contour).mean()
+        max_average_speed = self.speed_coef_mean(original_contour)
         speed_array = [max_average_speed]
         pixel_min = 1
 
@@ -808,7 +808,7 @@ def get_uavg(self, all_contours, centlon_e, centlat_e, original_contour, anticyc
             if pixel_min > level_contour.nb_pixel:
                 break
             # Interpolate uspd to seglon, seglat, then get mean
-            level_average_speed = self.speed_coef(level_contour).mean()
+            level_average_speed = self.speed_coef_mean(level_contour)
             speed_array.append(level_average_speed)
             if level_average_speed >= max_average_speed:
                 max_average_speed = level_average_speed
@@ -944,12 +944,12 @@ def _low_filter(self, grid_name, x_cut, y_cut, factor=40.):
         z_interp = RectBivariateSpline(x_center, y_center, z_filtered, **opts_interpolation).ev(x, y)
         return ma.array(z_interp, mask=m_interp.ev(x, y) > 0.00001)
 
-    def speed_coef(self, contour):
+    def speed_coef_mean(self, contour):
         dist, idx = self.index_interp.query(uniform_resample_stack(contour.vertices)[1:], k=4)
         i_y = idx % self.x_c.shape[1]
         i_x = int_((idx - i_y) / self.x_c.shape[1])
         # A simplified solution to be change by a weight mean
-        return self._speed_norm[i_x, i_y].mean(axis=1)
+        return self._speed_norm[i_x, i_y].mean(axis=1).mean()
 
     def init_speed_coef(self, uname='u', vname='v'):
         self._speed_norm = (self.grid(uname) ** 2 + self.grid(vname) ** 2) ** .5
@@ -1348,11 +1348,11 @@ def add_uv(self, grid_height):
         d_x = self.EARTH_RADIUS * 2 * pi / 360 * d_x_degrees * cos(deg2rad(self.y_c))
         self.vars['v'] = d_hx / d_x * gof
 
-    def speed_coef(self, contour):
+    def speed_coef_mean(self, contour):
         """some nan can be compute over contour if we are near border,
         something to explore
         """
-        return value_on_regular_contour(
+        return mean_on_regular_contour(
             self.x_c, self.y_c,
             self._speed_ev, self._speed_ev.mask,
             contour.vertices)
@@ -1378,18 +1378,8 @@ def interp(self, grid_name, lons, lats):
         Returns:
             new z
         """
-        z = empty(lons.shape, dtype='f4').reshape(-1)
         g = self.grid(grid_name)
-        interp_numba(
-            self.x_c,
-            self.y_c,
-            g,
-            lons.reshape(-1),
-            lats.reshape(-1),
-            z,
-            g.fill_value
-        )
-        return z
+        return interp2d_geo(self.x_c, self.y_c, g, g.mask, lons, lats)
 
 
 @njit(cache=True, fastmath=True, parallel=True)

src/py_eddy_tracker/eddy_feature.py

@@ -46,16 +46,16 @@ class Amplitude(object):
         'pixel_mask',
         'sla',
         'contour',
-        'interval',
+        'interval_min',
         'amplitude',
         'mle',
         )
 
     def __init__(self, contour, contour_height, data, interval):
         # Height of the contour
         self.h_0 = contour_height
-        # Step between two level
-        self.interval = interval
+        # Step minimal to consider amplitude
+        self.interval_min = interval * 2
         # Indices of all pixels in contour
         self.contour = contour
         # Link on original grid (local view) or copy if it's on bound
@@ -71,10 +71,10 @@ def __init__(self, contour, contour_height, data, interval):
             self.grid_extract = data[x_start:x_stop, y_start:y_stop]
         # => maybe replace pixel out of contour by nan?
         self.pixel_mask = zeros(self.grid_extract.shape, dtype='bool')
+        i_x = contour.pixels_index[0] - x_start
         if on_bounds:
-            i_x = (contour.pixels_index[0] - x_start) % data.shape[0]
-        else:
-            i_x = contour.pixels_index[0] - x_start
+            i_x %= data.shape[0]
+
         self.pixel_mask[i_x, contour.pixels_index[1] - y_start] = True
 
         # Only pixel in contour
@@ -87,7 +87,7 @@ def __init__(self, contour, contour_height, data, interval):
     def within_amplitude_limits(self):
         """Need update
         """
-        return (self.interval * 2) <= self.amplitude
+        return self.interval_min <= self.amplitude
         return self.eddy.ampmin <= self.amplitude <= self.eddy.ampmax
 
     def all_pixels_below_h0(self, level):
@@ -96,26 +96,27 @@ def all_pixels_below_h0(self, level):
         are below a given SSH threshold for cyclonic eddies.
         """
         # In some case pixel value must be very near of contour bounds
-        if ((self.sla - self.h_0) > self.EPSILON).any() or self.sla.mask.any():
+        if self.sla.mask.any() or ((self.sla.data - self.h_0) > self.EPSILON).any():
             return False
         else:
             # All local extrema index on th box
-            lmi_i, lmi_j = detect_local_minima_(self.grid_extract, self.grid_extract.mask, self.pixel_mask, self.mle, 1)
+            lmi_i, lmi_j = detect_local_minima_(self.grid_extract.data, self.grid_extract.mask, self.pixel_mask, self.mle, 1)
+            # After we use grid.data because index are in contour and we check before than no pixel are hide
             nb = len(lmi_i)
-            (x_start, _), (y_start, _) = self.contour.bbox_slice
             if nb == 0:
                 logging.warning('No extrema found in contour in level %f', level)
                 return False
             elif nb == 1:
                 i, j = lmi_i[0], lmi_j[0]
             else:
                 # Verify if several extrema are seriously below contour
-                nb_real_extrema = ((level - self.grid_extract[lmi_i, lmi_j]) >= 2 * self.interval).sum()
+                nb_real_extrema = ((level - self.grid_extract.data[lmi_i, lmi_j]) >= self.interval_min).sum()
                 if nb_real_extrema > self.mle:
                     return False
-                index = self.grid_extract[lmi_i, lmi_j].argmin()
+                index = self.grid_extract.data[lmi_i, lmi_j].argmin()
                 i, j = lmi_i[index], lmi_j[index]
-            self.amplitude = abs(self.grid_extract[i, j] - self.h_0)
+            self.amplitude = abs(self.grid_extract.data[i, j] - self.h_0)
+            (x_start, _), (y_start, _) = self.contour.bbox_slice
             i += x_start
             j += y_start
             return i, j
@@ -126,26 +127,26 @@ def all_pixels_above_h0(self, level):
         are above a given SSH threshold for anticyclonic eddies.
         """
         # In some case pixel value must be very near of contour bounds
-        if ((self.sla - self.h_0) < - self.EPSILON).any() or self.sla.mask.any():
+        if self.sla.mask.any() or ((self.h_0 - self.sla.data) > self.EPSILON).any():
             return False
         else:
             # All local extrema index on th box
-            lmi_i, lmi_j = detect_local_minima_(self.grid_extract, self.grid_extract.mask, self.pixel_mask, self.mle, -1)
+            lmi_i, lmi_j = detect_local_minima_(self.grid_extract.data, self.grid_extract.mask, self.pixel_mask, self.mle, -1)
             nb = len(lmi_i)
-            (x_start, _), (y_start, _) = self.contour.bbox_slice
             if nb == 0:
                 logging.warning('No extrema found in contour in level %f', level)
                 return False
             elif nb == 1:
                 i, j = lmi_i[0], lmi_j[0]
             else:
                 # Verify if several extrema are seriously above contour
-                nb_real_extrema = ((self.grid_extract[lmi_i, lmi_j] - level) >= 2 * self.interval).sum()
+                nb_real_extrema = ((self.grid_extract.data[lmi_i, lmi_j] - level) >= self.interval_min).sum()
                 if nb_real_extrema > self.mle:
                     return False
-                index = self.grid_extract[lmi_i, lmi_j].argmax()
+                index = self.grid_extract.data[lmi_i, lmi_j].argmax()
                 i, j = lmi_i[index], lmi_j[index]
-            self.amplitude = abs(self.grid_extract[i, j] - self.h_0)
+            self.amplitude = abs(self.grid_extract.data[i, j] - self.h_0)
+            (x_start, _), (y_start, _) = self.contour.bbox_slice
             i += x_start
             j += y_start
             return i, j

src/py_eddy_tracker/observations.py

@@ -159,6 +159,14 @@ class EddiesObservations(object):
     *local maxima/minima* within a closed region of a sea level anomaly field.
 
     """
+    __slots__ = (
+        'track_extra_variables',
+        'track_array_variables',
+        'array_variables',
+        'observations',
+        'active',
+        'sign_type',
+    )
 
     DELTA_JJ_JE = 2448623 - (datetime(1992, 1, 1) - datetime(1950, 1, 1)).days 
 
@@ -851,6 +859,7 @@ def display(self, ax, ref=0, **kwargs):
 class VirtualEddiesObservations(EddiesObservations):
     """Class to work with virtual obs
     """
+    __slots__ = ()
 
     @property
     def elements(self):
@@ -862,6 +871,7 @@ def elements(self):
 class TrackEddiesObservations(EddiesObservations):
     """Class to practice Tracking on observations
     """
+    __slots__ = ()
 
     def filled_by_interpolation(self, mask):
         """Filled selected values by interpolation