Updat of object oriented identification

Author: AntSimi

setup.py

@@ -16,6 +16,7 @@
     packages=find_packages('src'),
     package_dir={'': 'src'},
     scripts=[
+        'src/scripts/EddyId',
         'src/scripts/EddyIdentification',
         'src/scripts/EddyTracking',
         'src/scripts/EddyFinalTracking',
@@ -44,5 +45,6 @@
         'shapely',
         'pyyaml',
         'pyproj',
+        'pint'
         ],
     )

src/py_eddy_tracker/__init__.py

@@ -101,6 +101,20 @@ def parse_args(self, *args, **kwargs):
             reference_description="Julian date on Jan 1, 1992",
             )
         ),
+    time_jj=dict(
+        attr_name='time',
+        nc_name='time',
+        nc_type='int32',
+        nc_dims=('Nobs',),
+        nc_attr=dict(
+            standard_name='time',
+            units='days since 1950-01-01 00:00:00',
+            calendar='proleptic_gregorian',
+            axis='T',
+            long_name='time of gridded file',
+            description='date of this observation',
+            )
+        ),
     ocean_time=dict(
         attr_name=None,
         nc_name='ocean_time',
@@ -186,20 +200,6 @@ def parse_args(self, *args, **kwargs):
                         'around the contour defining the eddy perimeter',
             )
         ),
-    #~ radius=dict(
-        #~ attr_name=None,
-        #~ nc_name='L',
-        #~ nc_type='float32',
-        #~ nc_dims=('Nobs',),
-        #~ scale_factor=1e-3,
-        #~ nc_attr=dict(
-            #~ long_name='speed radius scale',
-            #~ units='km',
-            #~ description='radius of a circle whose area is equal to that '
-                        #~ 'enclosed by the contour of maximum circum-average'
-                        #~ ' speed',
-            #~ )
-        #~ ),
     speed_radius=dict(
         attr_name='speed_radius',
         scale_factor=100,

src/py_eddy_tracker/dataset/grid.py

@@ -4,20 +4,23 @@
 import logging
 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, nan, ceil, sinc, float64
+    linspace, errstate, int_, column_stack, interp, meshgrid, nan, ceil, sinc, float64, isnan
+from datetime import datetime
 from scipy.special import j1
-from scipy.signal import convolve2d
 from netCDF4 import Dataset
 from scipy.ndimage import gaussian_filter, convolve
-from scipy.interpolate import RectBivariateSpline
+from scipy.interpolate import RectBivariateSpline, interp1d
 from scipy.spatial import cKDTree
+from scipy.signal import welch, convolve2d
 from matplotlib.path import Path as BasePath
 from matplotlib.contour import QuadContourSet as BaseQuadContourSet
 from pyproj import Proj
+from pint import UnitRegistry
 from ..tools import fit_circle_c, distance_vector, winding_number_poly, poly_contain_poly, \
     distance, distance_point_vector
 from ..observations import EddiesObservations
 from ..eddy_feature import Amplitude, Contours
+from .. import VAR_DESCR
 
 
 def raw_resample(datas, fixed_size):
@@ -314,6 +317,12 @@ def is_circular(self):
         """
         return False
 
+    def units(self, varname):
+        with Dataset(self.filename) as h:
+            var = h.variables[varname]
+            if hasattr(var, 'units'):
+                return var.units
+
     def grid(self, varname):
         """give grid required
         """
@@ -351,8 +360,10 @@ def bounds(self):
         """
         return self.x_bounds.min(), self.x_bounds.max(), self.y_bounds.min(), self.y_bounds.max()
 
-    def eddy_identification(self, grid_height, uname, vname, step=0.005, shape_error=55,
+    def eddy_identification(self, grid_height, uname, vname, date, step=0.005, shape_error=55,
                             array_sampling=50, pixel_limit=None, bbox_surface_min_degree=.125**2):
+        if not isinstance(date, datetime):
+            raise Exception('Date argument be a datetime object')
         # The inf limit must be in pixel and  sup limit in surface
         if pixel_limit is None:
             pixel_limit = (4, 1000)
@@ -459,8 +470,7 @@ def eddy_identification(self, grid_height, uname, vname, step=0.005, shape_error
                     # Instantiate new EddyObservation object
                     properties = EddiesObservations(
                         size=1,
-                        track_extra_variables=['shape_error_e', 'shape_error_s', 'height_max_speed_contour',
-                                               'height_external_contour', 'height_inner_contour', 'nb_contour_selected'],
+                        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']
                     )
@@ -492,7 +502,23 @@ def eddy_identification(self, grid_height, uname, vname, step=0.005, shape_error
                     eddies.append(properties)
                     # To reserve definitively the area
                     data.mask[i_x_in, i_y_in] = True
-            a_and_c.append(EddiesObservations.concatenate(eddies))
+            if len(eddies) == 0:
+                eddies_collection = EddiesObservations()
+            else:
+                eddies_collection = EddiesObservations.concatenate(eddies)
+            eddies_collection.sign_type = 1 if anticyclonic_search else -1
+            eddies_collection.obs['time'] = (date - datetime(1950, 1, 1)).total_seconds() / 86400.
+
+            a_and_c.append(eddies_collection)
+        h_units = self.units(grid_height)
+        units = UnitRegistry()
+        in_unit = units.parse_expression(h_units)
+        if in_unit is not None:
+            for name in ['amplitude', 'height_max_speed_contour', 'height_external_contour', 'height_inner_contour']:
+                out_unit = units.parse_expression(VAR_DESCR[name]['nc_attr']['units'])
+                factor, _ = in_unit.to(out_unit).to_tuple()
+                a_and_c[0].obs[name] *= factor
+                a_and_c[1].obs[name] *= factor
         return a_and_c
 
     def get_uavg(self, all_contours, centlon_e, centlat_e, original_contour, anticyclonic_search, level_start,
@@ -528,7 +554,6 @@ def get_uavg(self, all_contours, centlon_e, centlat_e, original_contour, anticyc
             level_contour.pixels_in(self)
             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()
             speed_array.append(level_average_speed)
@@ -851,6 +876,10 @@ def kernel_bessel(self, lat, wave_length, order=1):
         # Estimate size of kernel
         step_y_km = self.ystep * distance(0, 0, 0, 1) / 1000
         step_x_km = self.xstep * distance(0, lat, 1, lat) / 1000
+        min_wave_length = max(step_x_km * 2, step_y_km * 2)
+        if wave_length < min_wave_length:
+            logging.error('Wave_length to short for resolution, must be > %d km', ceil(min_wave_length))
+            raise Exception()
         # half size will be multiply with by order
         half_x_pt, half_y_pt = ceil(wave_length / step_x_km).astype(int), ceil(wave_length / step_y_km).astype(int)
 
@@ -876,7 +905,7 @@ def kernel_bessel(self, lat, wave_length, order=1):
         kernel[dist_norm > order] = 0
         return kernel
 
-    def convolve_filter_with_dynamic_kernel(self, grid_name, kernel_func, lat_max, **kwargs_func):
+    def convolve_filter_with_dynamic_kernel(self, grid_name, kernel_func, lat_max=85, **kwargs_func):
         logging.warning('No filtering above %f degrees of latitude', lat_max)
         data = self.grid(grid_name).copy()
         # Matrix for result
@@ -943,6 +972,14 @@ def lanczos_low_filter(self, grid_name, wave_length, order=1, lat_max=85):
             grid_name, self.kernel_lanczos, lat_max=lat_max, wave_length=wave_length, order=order)
         self.vars[grid_name] = data_out
 
+    def bessel_band_filter(self, grid_name, wave_length_inf, wave_length_sup, **kwargs):
+        data_out = self.convolve_filter_with_dynamic_kernel(
+            grid_name, self.kernel_bessel, wave_length=wave_length_inf, **kwargs)
+        self.vars[grid_name] = data_out
+        data_out = self.convolve_filter_with_dynamic_kernel(
+            grid_name, self.kernel_bessel, wave_length=wave_length_sup, **kwargs)
+        self.vars[grid_name] -= data_out
+
     def bessel_high_filter(self, grid_name, wave_length, order=1, lat_max=85):
         data_out = self.convolve_filter_with_dynamic_kernel(
             grid_name, self.kernel_bessel, lat_max=lat_max, wave_length=wave_length, order=order)
@@ -953,6 +990,62 @@ def bessel_low_filter(self, grid_name, wave_length, order=1, lat_max=85):
             grid_name, self.kernel_bessel, lat_max=lat_max, wave_length=wave_length, order=order)
         self.vars[grid_name] = data_out
 
+    def spectrum_lonlat(self, grid_name, area=None, ref=None, **kwargs):
+        if area is None:
+            area = dict(llcrnrlon=190, urcrnrlon=280, llcrnrlat=-62, urcrnrlat=8)
+        scaling = kwargs.pop('scaling', 'density')
+        x0, y0 = self.nearest_grd_indice(area['llcrnrlon'], area['llcrnrlat'])
+        x1, y1 = self.nearest_grd_indice(area['urcrnrlon'], area['urcrnrlat'])
+
+        data = self.grid(grid_name)[x0:x1,y0:y1]
+
+        # Lat spectrum
+        pws = list()
+        step_y_km = self.ystep * distance(0, 0, 0, 1) / 1000
+        nb_invalid = 0
+        for i, _ in enumerate(self.x_c[x0:x1]):
+            f, pw = welch(data[i,:],  1 / step_y_km, scaling=scaling, **kwargs)
+            if isnan(pw).any():
+                nb_invalid += 1
+                continue
+            pws.append(pw)
+        if nb_invalid:
+            logging.warning('%d/%d columns invalid', nb_invalid, i + 1)
+        lat_content = 1 / f, array(pws).mean(axis=0)
+
+        # Lon spectrum
+        fs, pws = list(), list()
+        f_min, f_max = None, None
+        nb_invalid = 0
+        for i, lat in enumerate(self.y_c[y0:y1]):
+            step_x_km = self.xstep * distance(0, lat, 1, lat) / 1000
+            f, pw = welch(data[:,i], 1 / step_x_km, scaling=scaling, **kwargs)
+            if isnan(pw).any():
+                nb_invalid += 1
+                continue
+            if f_min is None:
+                f_min = f.min()
+                f_max = f.max()
+            else:
+                f_min = max(f_min, f.min())
+                f_max = min(f_max, f.max())
+            fs.append(f)
+            pws.append(pw)
+        if nb_invalid:
+            logging.warning('%d/%d lines invalid', nb_invalid, i + 1)
+        f_interp = linspace(f_min, f_max, f.shape[0])
+        pw_m = array(
+            [interp1d(f, pw, fill_value=0., bounds_error=False)(f_interp) for f, pw in zip(fs, pws)]).mean(axis=0)
+        lon_content = 1 / f_interp, pw_m
+        if ref is None:
+            return lon_content, lat_content
+        else:
+            ref_lon_content, ref_lat_content = ref.spectrum_lonlat(grid_name, area, **kwargs)
+            return (lon_content[0], lon_content[1] / ref_lon_content[1]), \
+                   (lat_content[0], lat_content[1] / ref_lat_content[1])
+
+
+
     def add_uv(self, grid_height):
         """Compute a u and v grid
         """

src/py_eddy_tracker/eddy_feature.py

@@ -49,7 +49,6 @@ class Amplitude(object):
         'iy',
         'interval',
         'amplitude',
-        'local_extrema',
         'local_extrema_inds',
         'mle',
         )
@@ -68,9 +67,10 @@ def __init__(self, i_contour_x, i_contour_y, contour_height, data, interval):
         self.iy -= self.slice_y.start
         # Link on original grid (local view)
         self.data = data[self.slice_x, self.slice_y]
+        # => maybe replace pixel out of contour by nan?
+
         # Amplitude which will be provide
         self.amplitude = 0
-        self.local_extrema = None
         self.local_extrema_inds = None
         # Maximum local extrema accepted
         self.mle = 1
@@ -101,23 +101,16 @@ def all_pixels_below_h0(self, level):
             return False
         else:
             self._set_local_extrema(1)
-            if 0 < self.local_extrema <= self.mle:
-                self._set_cyc_amplitude()
+            lmi_i, lmi_j = where(self.local_extrema_inds)
+            nb_real_extrema = ((level - self.data[lmi_i, lmi_j]) >= 2 * self.interval).sum()
+            if nb_real_extrema > self.mle:
                 return False
-            elif self.local_extrema > self.mle:
-                lmi_i, lmi_j = where(self.local_extrema_inds)
-                levnm2 = level - (2 * self.interval)
-                slamin = 1e5
-                for i, j in zip(lmi_i, lmi_j):
-                    if slamin >= self.data[i, j]:
-                        slamin = self.data[i, j]
-                        imin, jmin = i, j
-                    if self.data[i, j] >= levnm2:
-                        self._set_cyc_amplitude()
-                        # Prevent further calls to_set_cyc_amplitude
-                        levnm2 = 1e5
-                return imin + self.slice_x.start, jmin + self.slice_y.start
-        return False
+            amp_min = level - (2 * self.interval)
+            index = self.data[lmi_i, lmi_j].argmin()
+            jmin_, imin_ = lmi_j[index], lmi_i[index]
+            if self.data[imin_, jmin_] <= amp_min:
+                self._set_cyc_amplitude()
+            return imin_ + self.slice_x.start, jmin_ + self.slice_y.start
 
     def all_pixels_above_h0(self, level):
         """
@@ -130,39 +123,23 @@ def all_pixels_above_h0(self, level):
             return False
         else:
             self._set_local_extrema(-1)
-            # If we have a number of extrema avoid, we compute amplitude
-            if 0 < self.local_extrema <= self.mle:
-                self._set_acyc_amplitude()
+            lmi_i, lmi_j = where(self.local_extrema_inds)
+            nb_real_extrema = ((self.data[lmi_i, lmi_j] - level) >= 2 * self.interval).sum()
+            if nb_real_extrema > self.mle:
                 return False
-
-            # More than avoid
-            elif self.local_extrema > self.mle:
-                # index of extrema
-                lmi_i, lmi_j = where(self.local_extrema_inds)
-                levnp2 = level + (2 * self.interval)
-                slamax = -1e5
-                # Iteration on extrema
-                for i, j in zip(lmi_i, lmi_j):
-                    # We iterate on max and search the first sup of slamax
-                    if slamax <= self.data[i, j]:
-                        slamax = self.data[i, j]
-                        imax, jmax = i, j
-
-                    if self.data[i, j] <= levnp2:
-                        self._set_acyc_amplitude()
-                        # Prevent further calls to_set_acyc_amplitude
-                        levnp2 = -1e5
-                return imax + self.slice_x.start, jmax + self.slice_y.start
-        return False
+            amp_min = level + (2 * self.interval)
+            index = self.data[lmi_i, lmi_j].argmax()
+            jmin_, imin_ = lmi_j[index], lmi_i[index]
+            if self.data[imin_, jmin_] >= amp_min:
+                self._set_cyc_amplitude()
+            return imin_ + self.slice_x.start, jmin_ + self.slice_y.start
 
     def _set_local_extrema(self, sign):
         """
         Set count of local SLA maxima/minima within eddy
         """
         # mask of minima
         self.local_extrema_inds = self.detect_local_minima(self.data * sign)
-        # nb of minima
-        self.local_extrema = self.local_extrema_inds.sum()
 
     @staticmethod
     def detect_local_minima(grid):
@@ -174,8 +151,8 @@ def detect_local_minima(grid):
         """
         # Equivalent
         neighborhood = ones((3, 3), dtype='bool')
-        # ~ neighborhood = generate_binary_structure(grid.ndim, 2)
-
+        if hasattr(grid, 'mask'):
+            grid[grid.mask] = 2e10
         # Get local mimima
         detected_minima = minimum_filter(
             grid, footprint=neighborhood) == grid

src/py_eddy_tracker/grid/__init__.py

@@ -43,12 +43,17 @@
 
 
 def browse_dataset_in(data_dir, files_model, date_regexp, date_model,
-                      start_date=None, end_date=None, sub_sampling_step=1):
-    pattern_regexp = re_compile('.*/' + date_regexp)
-    full_path = join_path(data_dir, files_model)
-    logging.info('Search files : %s', full_path)
+                      start_date=None, end_date=None, sub_sampling_step=1,
+                      files=None):
+    if files is not None:
+        pattern_regexp = re_compile('.*/' + date_regexp)
+        filenames = bytes_(files)
+    else:
+        pattern_regexp = re_compile('.*/' + date_regexp)
+        full_path = join_path(data_dir, files_model)
+        logging.info('Search files : %s', full_path)
+        filenames = bytes_(glob(full_path))
 
-    filenames = bytes_(glob(full_path))
     dataset_list = empty(len(filenames),
                          dtype=[('filename', 'S500'),
                                 ('date', 'datetime64[D]'),