Add function to extract u profil

Author: adelepoulle

src/py_eddy_tracker/__init__.py

@@ -213,6 +213,26 @@ def parse_args(self, *args, **kwargs):
                         'scale L',
             )
         ),
+    i=dict(
+        attr_name='i',
+        nc_name='i',
+        nc_type='uint16',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            long_name='longitude index in the grid of the detection',
+            description='longitude index in the grid of the detection',
+            )
+        ),
+    j=dict(
+        attr_name='j',
+        nc_name='j',
+        nc_type='uint16',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            long_name='latitude index in the grid of the detection',
+            description='latitude index in the grid of the detection',
+            )
+        ),
     eke=dict(
         attr_name='eke',
         nc_name='Teke',
@@ -327,7 +347,7 @@ def parse_args(self, *args, **kwargs):
         attr_name=None,
         nc_name='uavg_profile',
         nc_type='f4',
-        nc_dims=('uavg_contour_count', 'Nobs',),
+        nc_dims=('Nobs', 'NbSample'),
         nc_attr=dict(
             long_name='radial profile of uavg',
             description='all uavg values from effective contour inwards to '
@@ -356,6 +376,42 @@ def parse_args(self, *args, **kwargs):
             units='%',
             )
         ),
+    height_max_speed_contour=dict(
+        attr_name=None,
+        nc_name='height_max_speed_contour',
+        nc_type='f4',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            units='m',
+            )
+        ),
+    height_external_contour=dict(
+        attr_name=None,
+        nc_name='height_external_contour',
+        nc_type='f4',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            units='m',
+            )
+        ),
+    height_inner_contour=dict(
+        attr_name=None,
+        nc_name='height_inner_contour',
+        nc_type='f4',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            units='m',
+            )
+        ),
+    nb_contour_selected=dict(
+        attr_name=None,
+        nc_name='nb_contour_selected',
+        nc_type='u2',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            units='1',
+            )
+        ),
     )
 
 for key in VAR_DESCR.keys():

src/py_eddy_tracker/dataset/grid.py

@@ -2,7 +2,6 @@
 """
 """
 import logging
-from scipy.interpolate import interp1d
 from numpy import concatenate, int32, empty, maximum, where, array, \
     sin, deg2rad, pi, ones, cos, ma, int8, histogram2d, arange, float_, \
     linspace, errstate, int_, column_stack, interp, meshgrid, unique, nan
@@ -18,6 +17,11 @@
 from ..eddy_feature import Amplitude, get_uavg, Contours
 
 
+def raw_resample(datas, fixed_size):
+    nb_value = datas.shape[0]
+    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
@@ -441,7 +445,7 @@ def eddy_identification(self, grid_height, uname, vname,
                             centlat_e = y[centi, centj]
 
                     # centlat_e and centlon_e must be index of maximum, we will loose some inner contour, if it's not
-                    average_speed, speed_contour, inner_contour, speed_array = \
+                    max_average_speed, speed_contour, inner_contour, speed_array, i_max_speed, i_inner = \
                         get_uavg(self, contours, centlon_e, centlat_e, current_contour, anticyclonic_search, corrected_coll_index)
 
                     # Use azimuth equal projection for radius
@@ -459,13 +463,31 @@ def eddy_identification(self, grid_height, uname, vname,
                     # Instantiate new EddyObservation object
                     properties = EddiesObservations(
                         size=1,
-                        track_extra_variables=['shape_error_e', 'shape_error_s'],
+                        track_extra_variables=['shape_error_e', 'shape_error_s', 'height_max_speed_contour',
+                                               'height_external_contour', 'height_inner_contour', 'nb_contour_selected'],
                         track_array_variables=array_sampling,
-                        array_variables=['contour_lon_e', 'contour_lat_e', 'contour_lon_s', 'contour_lat_s']
+                        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.obs['height_external_contour'] = cvalues
+                    properties.obs['height_inner_contour'] = contours.cvalues[i_inner]
+                    array_size = speed_array.shape[0]
+                    properties.obs['nb_contour_selected'] = array_size
+                    properties.obs['uavg_profile'] = raw_resample(speed_array, array_sampling)
+                    # from matplotlib import pyplot as plt
+                    # if array_size > 10:
+                    #     plt.figure()
+                    #     plt.plot(linspace(properties.obs['height_external_contour'],properties.obs['height_inner_contour'], speed_array.shape[0]), speed_array, 'b')
+                    #     plt.axvline(properties.obs['height_inner_contour'], color='g')
+                    #     plt.axvline(properties.obs['height_max_speed_contour'], color='r')
+                    #     plt.axvline(properties.obs['height_external_contour'], color='k')
+                    #     plt.title('%d' % array_size)
+                    #     plt.ylim(0,None)
+                    #     plt.show()
                     properties.obs['amplitude'] = amp.amplitude
                     properties.obs['radius_s'] = eddy_radius_s / 1000
-                    properties.obs['speed_radius'] = average_speed
+                    properties.obs['speed_radius'] = max_average_speed
                     properties.obs['radius_e'] = eddy_radius_e / 1000
                     properties.obs['shape_error_e'] = aerr
                     properties.obs['shape_error_s'] = aerr_s
@@ -820,5 +842,7 @@ def speed_coef(self, contour):
     def init_speed_coef(self, uname='u', vname='v'):
         """Draft
         """
-        uspd = (self.grid(uname) ** 2 + self.grid(vname) ** 2) ** .5
-        self._speed_ev = RectBivariateSpline(self.x_c, self.y_c, uspd, kx=1, ky=1).ev
+        speed = (self.grid(uname) ** 2 + self.grid(vname) ** 2) ** .5
+        # Evaluation near masked value will be smoothed to 0 !!!, not perfect
+        speed[speed.mask] = 0
+        self._speed_ev = RectBivariateSpline(self.x_c, self.y_c, speed, kx=1, ky=1).ev

src/py_eddy_tracker/eddy_feature.py

@@ -98,7 +98,7 @@ def all_pixels_below_h0(self, level):
         are below a given SSH threshold for cyclonic eddies.
         """
         sla = self.data[self.ix, self.iy]
-        if (sla > self.h_0).any():
+        if (sla > self.h_0).any() or (hasattr(sla, 'mask') and sla.mask.any()):
             return False
         else:
             self._set_local_extrema(1)
@@ -126,7 +126,7 @@ def all_pixels_above_h0(self, level):
         are above a given SSH threshold for anticyclonic eddies.
         """
         sla = self.data[self.ix, self.iy]
-        if (sla < self.h_0).any():
+        if (sla < self.h_0).any() or (hasattr(sla, 'mask') and sla.mask.any()):
             # i.e.,with self.amplitude == 0
             return False
         else:
@@ -335,21 +335,25 @@ def get_uavg(grid, all_contours, centlon_e, centlat_e, original_contour, anticyc
     Calculate geostrophic speed around successive contours
     Returns the average
     """
-    average_speed = grid.speed_coef(original_contour).mean()
-    speed_array = [average_speed]
+    max_average_speed = grid.speed_coef(original_contour).mean()
+    speed_array = [max_average_speed]
     pixel_min = 1
 
     eddy_contours = [original_contour]
     inner_contour = selected_contour = original_contour
     # Must start only on upper or lower contour, no need to test the two part
-    for coll in all_contours.iter(start=level_start + 1, step=1 if anticyclonic_search else -1):
+    step = 1 if anticyclonic_search else -1
+    i_inner = i_max_speed = -1
+
+    for i, coll in enumerate(all_contours.iter(start=level_start + step, step=step)):
         level_contour = coll.get_nearest_path_bbox_contain_pt(centlon_e, centlat_e)
         # Leave loop if no contours at level
         if level_contour is None:
             break
-        # 1. Ensure polygon_i contains point centlon_e, centlat_e
-        if winding_number_poly(centlon_e, centlat_e, level_contour.vertices) == 0:
-            break
+        # 1. Ensure polygon_i contains point centlon_e, centlat_e (Maybe we loose some inner contour if eddy
+        #        core are not centered)
+        # if winding_number_poly(centlon_e, centlat_e, level_contour.vertices) == 0:
+        #     break
         # 2. Ensure polygon_i is within polygon_e
         if not poly_contain_poly(original_contour.vertices, level_contour.vertices):
             break
@@ -360,11 +364,15 @@ def get_uavg(grid, all_contours, centlon_e, centlat_e, original_contour, anticyc
         # Interpolate uspd to seglon, seglat, then get mean
         level_average_speed = grid.speed_coef(level_contour).mean()
         speed_array.append(level_average_speed)
-        if level_average_speed >= average_speed:
-            average_speed = level_average_speed
+        if level_average_speed >= max_average_speed:
+            max_average_speed = level_average_speed
+            i_max_speed = i
             selected_contour = level_contour
         inner_contour = level_contour
         eddy_contours.append(level_contour)
+        i_inner = i
     for contour in eddy_contours:
         contour.used = True
-    return average_speed, selected_contour, inner_contour, speed_array
+    i_max_speed = level_start + step + step * i_max_speed
+    i_inner = level_start + step + step * i_inner
+    return max_average_speed, selected_contour, inner_contour, array(speed_array), i_max_speed, i_inner

src/py_eddy_tracker/observations.py

@@ -356,7 +356,7 @@ def load_from_netcdf(cls, filename):
                 if variable == 'cyc':
                     continue
                 var_inv = VAR_DESCR_inv[variable]
-                if var_inv not in cls.ELEMENTS and var_inv not in kwargs['array_variables']:
+                if var_inv not in cls.ELEMENTS and var_inv not in kwargs.get('array_variables', list()):
                     kwargs['track_extra_variables'].append(var_inv)
 
             eddies = cls(size=nb_obs, ** kwargs)
@@ -366,6 +366,9 @@ def load_from_netcdf(cls, filename):
                 eddies.obs[VAR_DESCR_inv[variable]
                            ] = h_nc.variables[variable][:]
             eddies.sign_type = h_nc.variables['cyc'][0]
+            if eddies.sign_type == 0:
+                logging.info('File come from another algorithm')
+                eddies.sign_type = -1
         return eddies
 
     @classmethod