amplitude will not use pixel outside of contour

Author: AntSimi

src/py_eddy_tracker/dataset/grid.py

@@ -194,6 +194,7 @@ class GridDataset(object):
         'vars',
         'interpolators',
         'speed_coef',
+        'contours',
     )
 
     GRAVITY = 9.807
@@ -211,6 +212,7 @@ def __init__(self, filename, x_name, y_name, centered=None):
         self.x_dim = None
         self.y_dim = None
         self.centered = centered
+        self.contours = None
         self.xinterp = None
         self.yinterp = None
         self.filename = filename
@@ -413,16 +415,18 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
         x, y = self.x_c, self.y_c
 
         # Compute ssh contour
-        contours = Contours(x, y, data, levels, bbox_surface_min_degree=bbox_surface_min_degree, wrap_x=self.is_circular())
+        self.contours = Contours(x, y, data, levels, bbox_surface_min_degree=bbox_surface_min_degree, wrap_x=self.is_circular())
 
+        track_extra_variables = ['height_max_speed_contour', 'height_external_contour', 'height_inner_contour']
+        array_variables = ['contour_lon_e', 'contour_lat_e', 'contour_lon_s', 'contour_lat_s', 'uavg_profile']
         # Compute cyclonic and anticylonic research:
         a_and_c = list()
         for anticyclonic_search in [True, False]:
             eddies = list()
             iterator = 1 if anticyclonic_search else -1
 
             # Loop over each collection
-            for coll_ind, coll in enumerate(contours.iter(step=iterator)):
+            for coll_ind, coll in enumerate(self.contours.iter(step=iterator)):
                 corrected_coll_index = coll_ind
                 if iterator == -1:
                     corrected_coll_index = - coll_ind - 1
@@ -431,7 +435,7 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                 nb_paths = len(contour_paths)
                 if nb_paths == 0:
                     continue
-                cvalues = contours.cvalues[corrected_coll_index]
+                cvalues = self.contours.cvalues[corrected_coll_index]
                 logging.debug('doing collection %s, contour value %.4f, %d paths',
                               corrected_coll_index, cvalues, nb_paths)
 
@@ -440,6 +444,7 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     if current_contour.used:
                         continue
                     centlon_e, centlat_e, eddy_radius_e, aerr = current_contour.fit_circle()
+
                     # Filter for shape
                     if aerr < 0 or aerr > shape_error:
                         continue
@@ -465,13 +470,14 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     # Compute amplitude
                     reset_centroid, amp = self.get_amplitude(current_contour, cvalues, data,
                                                              anticyclonic_search=anticyclonic_search,
-                                                             level=contours.levels[corrected_coll_index], step=step)
-
+                                                             level=self.contours.levels[corrected_coll_index], step=step)
+                    print(reset_centroid, amp.amplitude)
                     # If we have a valid amplitude
                     if (not amp.within_amplitude_limits()) or (amp.amplitude == 0):
                         continue
 
                     if reset_centroid:
+
                         if self.is_circular():
                             centi = self.normalize_x_indice(reset_centroid[0])
                         else:
@@ -487,7 +493,8 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
 
                     # centlat_e and centlon_e must be index of maximum, we will loose some inner contour, if it's not
                     max_average_speed, speed_contour, inner_contour, speed_array, i_max_speed, i_inner = \
-                        self.get_uavg(contours, centlon_e, centlat_e, current_contour, anticyclonic_search, corrected_coll_index)
+                        self.get_uavg(self.contours, centlon_e, centlat_e, current_contour, anticyclonic_search,
+                                      corrected_coll_index)
 
                     # Use azimuth equal projection for radius
                     proj = Proj('+proj=aeqd +ellps=WGS84 +lat_0={1} +lon_0={0}'.format(*inner_contour.mean_coordinates))
@@ -502,16 +509,13 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     _, _, eddy_radius_s, aerr_s = fit_circle_c(c_x, c_y)
 
                     # Instantiate new EddyObservation object
-                    properties = EddiesObservations(
-                        size=1,
-                        track_extra_variables=[ 'height_max_speed_contour', 'height_external_contour', 'height_inner_contour'],
-                        track_array_variables=array_sampling,
-                        array_variables=['contour_lon_e', 'contour_lat_e', 'contour_lon_s', 'contour_lat_s', 'uavg_profile']
-                    )
-
-                    properties.obs['height_max_speed_contour'] = contours.cvalues[i_max_speed]
+                    properties = EddiesObservations(size=1, track_extra_variables=track_extra_variables,
+                                                    track_array_variables=array_sampling,
+                                                    array_variables=array_variables)
+
+                    properties.obs['height_max_speed_contour'] = self.contours.cvalues[i_max_speed]
                     properties.obs['height_external_contour'] = cvalues
-                    properties.obs['height_inner_contour'] = contours.cvalues[i_inner]
+                    properties.obs['height_inner_contour'] = self.contours.cvalues[i_inner]
                     array_size = speed_array.shape[0]
                     properties.obs['nb_contour_selected'] = array_size
                     if speed_array.shape[0] == 1:
@@ -537,7 +541,9 @@ def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape
                     # To reserve definitively the area
                     data.mask[i_x_in, i_y_in] = True
             if len(eddies) == 0:
-                eddies_collection = EddiesObservations()
+                eddies_collection = EddiesObservations(track_extra_variables=track_extra_variables,
+                                                       track_array_variables=array_sampling,
+                                                       array_variables=array_variables)
             else:
                 eddies_collection = EddiesObservations.concatenate(eddies)
             eddies_collection.sign_type = 1 if anticyclonic_search else -1
@@ -1168,6 +1174,11 @@ def init_speed_coef(self, uname='u', vname='v'):
         speed[speed.mask] = 0
         self._speed_ev = RectBivariateSpline(self.x_c, self.y_c, speed, kx=1, ky=1).ev
 
+    def display(self, ax, name, **kwargs):
+        if 'cmap' not in kwargs:
+            kwargs['cmap'] = 'coolwarm'
+        return ax.pcolormesh(self.x_bounds, self.y_bounds, self.grid(name).T, **kwargs)
+
     def interp(self, grid_name, lons, lats):
         """
         Compute z over lons, lats

src/py_eddy_tracker/eddy_feature.py

@@ -27,7 +27,7 @@
 """
 
 import logging
-from numpy import ones, where, empty, array, concatenate, ma
+from numpy import ones, where, empty, array, concatenate, ma, zeros
 from scipy.ndimage import minimum_filter
 from scipy.ndimage import binary_erosion
 from matplotlib.figure import Figure
@@ -43,6 +43,7 @@ class Amplitude(object):
     __slots__ = (
         'h_0',
         'grid_extract',
+        'mask',
         'sla',
         'contour',
         'interval',
@@ -65,6 +66,8 @@ def __init__(self, contour, contour_height, data, interval):
         else:
             self.grid_extract = data[slice_x, slice_y]
         # => maybe replace pixel out of contour by nan?
+        self.mask = zeros(self.grid_extract.shape, dtype='bool')
+        self.mask[contour.pixels_index[0] - slice_x.start, contour.pixels_index[1] - slice_y.start] = True
 
         # Only pixel in contour
         self.sla = data[contour.pixels_index]
@@ -83,12 +86,12 @@ def within_amplitude_limits(self):
     def _set_cyc_amplitude(self):
         """Get amplitude for cyclone
         """
-        self.amplitude = self.h_0 - self.grid_extract.min()
+        self.amplitude = self.h_0 - self.sla.min()
 
     def _set_acyc_amplitude(self):
         """Get amplitude for anticyclone
         """
-        self.amplitude = self.grid_extract.max() - self.h_0
+        self.amplitude = self.sla.max() - self.h_0
 
     def all_pixels_below_h0(self, level):
         """
@@ -100,6 +103,8 @@ def all_pixels_below_h0(self, level):
         else:
             self._set_local_extrema(1)
             lmi_i, lmi_j = where(self.local_extrema_inds)
+            m = self.mask[lmi_i, lmi_j]
+            lmi_i, lmi_j = lmi_i[m], lmi_j[m]
             nb_real_extrema = ((level - self.grid_extract[lmi_i, lmi_j]) >= 2 * self.interval).sum()
             if nb_real_extrema > self.mle:
                 return False
@@ -122,6 +127,8 @@ def all_pixels_above_h0(self, level):
         else:
             self._set_local_extrema(-1)
             lmi_i, lmi_j = where(self.local_extrema_inds)
+            m = self.mask[lmi_i, lmi_j]
+            lmi_i, lmi_j = lmi_i[m], lmi_j[m]
             nb_real_extrema = ((self.grid_extract[lmi_i, lmi_j] - level) >= 2 * self.interval).sum()
             if nb_real_extrema > self.mle:
                 return False
@@ -277,7 +284,7 @@ def __init__(self, x, y, z, levels, bbox_surface_min_degree, wrap_x=False, keep_
         logging.debug('Y shape : %s', y.shape)
         logging.debug('Z shape : %s', z.shape)
         logging.info('Start computing iso lines with %d levels from %f to %f ...', len(levels), levels[0], levels[-1])
-        self.contours = ax.contour(x, y, z.T, levels)
+        self.contours = ax.contour(x, y, z.T, levels, cmap='Spectral_r')
         if wrap_x:
             self.find_wrapcut_path_and_join(x[0], x[-1])
         logging.info('Finish computing iso lines')
@@ -410,8 +417,12 @@ def get_index_nearest_path_bbox_contain_pt(self, level, xpt, ypt):
             return self.contours.collections[level]._paths[index]
 
     def display(self, ax, **kwargs):
+        from matplotlib.collections import LineCollection
         for collection in self.contours.collections:
-            for path in collection.get_paths():
-                x, y= path.vertices[:, 0], path.vertices[:, 1]
-                ax.plot(x, y, **kwargs)
-
+            ax.add_collection(LineCollection(
+                (i.vertices for i in collection.get_paths()),
+                color=collection.get_color(),
+                **kwargs
+            ))
+        ax.update_datalim([self.contours._mins, self.contours._maxs])
+        ax.autoscale_view()

src/py_eddy_tracker/observations.py

@@ -843,6 +843,10 @@ def write_netcdf(self, path='./', filename='%(path)s/%(sign_type)s.nc'):
     def set_global_attr_netcdf(self, h_nc):
         pass
 
+    def display(self, ax, ref=0, **kwargs):
+        ax.plot((self.obs['contour_lon_s'].T - ref) % 360 + ref, self.obs['contour_lat_s'].T, ** kwargs)
+        ax.plot((self.obs['contour_lon_e'].T - ref) % 360 + ref, self.obs['contour_lat_e'].T, linestyle='-.', ** kwargs)
+
 
 class VirtualEddiesObservations(EddiesObservations):
     """Class to work with virtual obs