Draft of tracking

Author: adelepoulle

share/eddy_identification.yaml

@@ -0,0 +1,95 @@
+# Two diagnostics, based on SLA (e.g., Chelton et al, 2011)
+DIAGNOSTIC_TYPE: 'SLA'
+
+# Specify domain
+DOMAIN:
+  #~ THE_DOMAIN: 'Global'
+  THE_DOMAIN: 'Regional'
+  LONMIN: -100
+  LONMAX: -5
+  LATMIN: 10
+  LATMAX: 45
+  DATE_STR: 2014-12-27
+  DATE_END: 2014-12-31
+
+DATASET:
+    # Path to data
+    #~ DATA_DIR: '/data/PVA/Externe/global/delayed-time/grids/msla/all-sat-merged/h/2014/'
+    DATA_DIR: '/home/adelepoulle/2014/'
+    # Model to select file
+    FILES_MODEL: dt_global_allsat_msla_h_????????_????????.nc
+    # Model to extract date (regular expression)
+    DATE_REGEXP: dt_global_allsat_msla_h_([0-9]*)_[0-9]*.nc
+    # Model to transfer date to python datetime
+    DATE_MODEL: '%Y%m%d'
+    # Grid Name to compute analysis
+    VAR_NAME: Grid_0001
+    # Longitude variable
+    LON_NAME: NbLongitudes
+    # Latitude variable
+    LAT_NAME: NbLatitudes
+    # Step to do SUBSAMPLING (1 no subsampling)
+    SUBSAMPLING: 1
+
+PATHS:
+  # Obtain file from:
+  # http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
+  RW_PATH: 'rossrad.dat'
+  # Path for saving of outputs
+  SAVE_DIR: 'test/'
+
+  # Reference Julian day (Julian date at Jan 1, 1992)
+JDAY_REFERENCE: 2448623.
+
+# Define contouring parameters
+CONTOUR_PARAMETER:
+  # Set SLA contour spacing (cm)
+  CONTOUR_PARAMETER_SLA:
+    MAX_SLA: 100.
+    INTERVAL: 1
+    #~ INTERVAL: 4
+
+# Shape error (see Kurian et al 2011)
+SHAPE_ERROR: 65. # % (Mason et al 2014 use 55%)
+    
+# Filtering
+SMOOTHING: Yes
+SMOOTHING_SLA:
+  TYPE: 'Gaussian'
+  ZWL: 20. # degrees, zonal wavelength (see Chelton etal 2011)
+  MWL: 10. # degrees, meridional wavelength
+
+# Min and max permitted eddy radii [degrees] and amplitudes (cm)
+# 0.4 for DT10, 0.35 for DT14
+RADMIN: 0.35
+RADMAX: 4.461
+AMPMIN: 0.
+AMPMAX: 150.
+
+# MPM14 use 1; CSS11 have no limit
+MAX_LOCAL_EXTREMA: 1
+
+# Typical parameters
+RAD0: 60000. # m
+AMP0: 2. # cm
+TEMP0: 15.
+SALT0: 35.
+
+# Parameters used by Chelton etal and Kurian etal (Sec. 3.2) to ensure 
+# the slow evolution of the eddies over time;
+# they use min and max values of 0.25 and 2.5
+EVOLVE_AMP_MIN: 0.0005 # min change in amplitude
+EVOLVE_AMP_MAX: 500 # max change in amplitude
+EVOLVE_AREA_MIN: 0.0005 # min change in area
+EVOLVE_AREA_MAX: 500 # max change in area
+
+# Define track_extra_variables to track and save:
+# - effective contour points
+# - speed-based contour points
+# - shape test values
+# - profiles of swirl velocity from effective contour inwards
+# Useful for working with ARGO data    
+TRACK_EXTRA_VARIABLES: No
+
+# 'RectBivariate' is faster
+INTERP_METHOD: 'RectBivariate'

share/tracking.yaml

@@ -0,0 +1,22 @@
+# Two diagnostics, based on SLA (e.g., Chelton et al, 2011)
+# or the Q-parameter (e.g., Kurian et al, 2011) are available
+DIAGNOSTIC_TYPE: 'SLA'
+
+PATHS:
+  FILES_PATTERN: /homelocal/adelepoulle/20151110_eddy_track/Anticyclonic_2014*.nc
+  #~ FILES_PATTERN: /homelocal/adelepoulle/20151110_eddy_track/Cyclonic_2014*.nc
+  # Path and filename of Chelton et al (1998) Rossby radius data
+  # Obtain file from:
+  # http://www-po.coas.oregonstate.edu/research/po/research/rossby_radius/
+  RW_PATH: 'rossrad.dat'
+  # Path for saving of outputs
+  SAVE_DIR: 'test/'
+
+TRACK_DURATION_MIN: 28 # Number of observations
+  
+SEPARATION_METHOD: 'ellipse'
+#SEPARATION_METHOD: 'sum_radii'
+
+
+# Typical parameters
+DIST0: 25000. # m separation distance after ~7 days (see CSS11 fig 22)

src/py_eddy_tracker/__init__.py

@@ -126,20 +126,20 @@ 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',
-            )
-        ),
+    #~ 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,
@@ -192,7 +192,7 @@ def parse_args(self, *args, **kwargs):
     track=dict(
         attr_name=None,
         nc_name='track',
-        nc_type='int32',
+        nc_type='uint32',
         nc_dims=('Nobs',),
         nc_attr=dict(
             long_name='track number',
@@ -203,7 +203,7 @@ def parse_args(self, *args, **kwargs):
     n=dict(
         attr_name=None,
         nc_name='n',
-        nc_type='int16',
+        nc_type='uint16',
         nc_dims=('Nobs',),
         nc_attr=dict(
             long_name='observation number',
@@ -254,3 +254,5 @@ def parse_args(self, *args, **kwargs):
             )
         ),
     )
+
+VAR_DESCR_inv = {VAR_DESCR[key]['nc_name'] : key for key in VAR_DESCR.keys()}

src/py_eddy_tracker/global_tracking.py

@@ -35,12 +35,10 @@ def create_variable(self, handler_nc, kwargs_variable,
         except ValueError:
             logging.warn('Data is empty')
 
-    def write_netcdf(self):
+    def write_netcdf(self, global_attr):
         """Write a netcdf with eddy
         """
         eddy_size = len(self.eddy.tmp_observations)
-        print dir(self)
-        exit()
         filename = '%s_%s.nc' % (self.sign_type, self.date.strftime('%Y%m%d'))
         with Dataset(filename, 'w', format='NETCDF4') as h_nc:
             logging.info('Create intermediary file %s', filename)
@@ -67,6 +65,10 @@ def write_netcdf(self):
                      dimensions=VAR_DESCR['type_cyc']['nc_dims']),
                 VAR_DESCR['type_cyc']['nc_attr'],
                 -1 if self.sign_type == 'Cyclonic' else 1)
+            # Global attr
+            h_nc.title = self.sign_type + ' eddy tracks'
+            for key, value in global_attr.iteritems():
+                h_nc.setncattr(key, str(value))
 
     def read_tracks(self):
         """

src/py_eddy_tracker/grid/__init__.py

@@ -61,20 +61,21 @@ def __init__(self):
         self._dy = None
         self._umask = None
         self._vmask = None
+        self.grid_filename = None
+        self.grid_date = None
 
-    @staticmethod
-    def read_nc(varfile, varname, indices=slice(None)):
+    def read_nc(self, varname, indices=slice(None)):
         """
         Read data from nectdf file
           varname : variable ('temp', 'mask_rho', etc) to read
-          indices : string of index ranges, eg. '[0,:,0]'
+          indices : slice
         """
-        with Dataset(varfile) as h_nc:
+        with Dataset(self.grid_filename) as h_nc:
             return h_nc.variables[varname][:][indices]
 
-    @staticmethod
-    def nc_variables(varfile):
-        with Dataset(varfile) as h_nc:
+    @property
+    def nc_variables(self):
+        with Dataset(self.grid_filename) as h_nc:
             return h_nc.variables.keys()
 
     @property
@@ -94,13 +95,13 @@ def view_unpad(self):
         return (self.slice_j_unpad,
                 self.slice_i_unpad)
 
-    def read_nc_att(self, varfile, varname, att):
+    def read_nc_att(self, varname, att):
         """
         Read data attribute from nectdf file
           varname : variable ('temp', 'mask_rho', etc) to read
           att : string of attribute, eg. 'valid_range'
         """
-        with Dataset(varfile) as h_nc:
+        with Dataset(self.grid_filename) as h_nc:
             return getattr(h_nc.variables[varname], att)
 
     @property
@@ -422,17 +423,10 @@ def set_interp_coeffs(self, sla, uspd):
         """
         Won't work for rotated grid
         """
-        if 'AVISO' in self.product:
-            self.sla_coeffs = interpolate.RectBivariateSpline(
-                self.lat[:, 0], self.lon[0], sla, kx=1, ky=1)
-            self.uspd_coeffs = interpolate.RectBivariateSpline(
-                self.lat[:, 0], self.lon[0], uspd, kx=1, ky=1)
-        elif 'ROMS' in self.product:
-            points = np.array([self.lon.ravel(), self.lat.ravel()]).T
-            self.sla_coeffs = interpolate.CloughTocher2DInterpolator(
-                points, sla.ravel())
-            self.uspd_coeffs = interpolate.CloughTocher2DInterpolator(
-                points, uspd.ravel())
+        self.sla_coeffs = interpolate.RectBivariateSpline(
+            self.lat[:, 0], self.lon[0], sla, kx=1, ky=1)
+        self.uspd_coeffs = interpolate.RectBivariateSpline(
+            self.lat[:, 0], self.lon[0], uspd, kx=1, ky=1)
         return self
 
     @staticmethod