Version 2.0.3

evanmason · evanmason · commit 7662ea736f0b · 2015-06-02T18:55:18.000+02:00
diff --git a/haversine.f90 b/haversine.f90
@@ -27,7 +27,7 @@ module haversine
 !
 !    If you have ifort on your system, change 'gfortran' to 'intelem'!
 !
-!   Version 2.0.1
+!   Version 2.0.3
 !
 !===========================================================================
 implicit none
diff --git a/make_eddy_track_AVISO.py b/make_eddy_track_AVISO.py
@@ -24,7 +24,7 @@
 
 make_eddy_track_AVISO.py
 
-Version 2.0.1
+Version 2.0.3
 
 
 Scroll down to line ~640 to get started
@@ -87,7 +87,6 @@ class PyEddyTracker (object):
     def __init__(self):
         """
         Set some constants
-
           *ZERO_CROSSING*:
             Boolean, *True* if THE_DOMAIN crosses 0 degree meridian
         """
@@ -318,15 +317,15 @@ def get_AVISO_f_pm_pn(self):
         # Get pm and pn
         pm = np.zeros_like(self.lonpad())
         pm[:, 1:-1] = self.haversine_dist(lonu[:, :-1], latu[:, :-1],
-                                          lonu[:, 1:],  latu[:, 1:])
+                                          lonu[:, 1:], latu[:, 1:])
         pm[:, 0] = pm[:, 1]
         pm[:, -1] = pm[:, -2]
         self._dx = pm
         self._pm = np.reciprocal(pm)
 
         pn = np.zeros_like(self.lonpad())
         pn[1:-1] = self.haversine_dist(lonv[:-1], latv[:-1],
-                                       lonv[1:],  latv[1:])
+                                       lonv[1:], latv[1:])
         pn[0] = pn[1]
         pn[-1] = pn[-2]
         self._dy = pn
@@ -499,7 +498,7 @@ def __init__(self, AVISO_FILE, THE_DOMAIN, PRODUCT,
         Initialise the grid object
         """
         super(AvisoGrid, self).__init__()
-        print '\nInitialising the *AVISO_grid*'
+        print('\nInitialising the *AVISO_grid*')
         self.THE_DOMAIN = THE_DOMAIN
         if 'AVISO' in PRODUCT:
             self.PRODUCT = 'AVISO'
@@ -601,7 +600,6 @@ def get_AVISO_data(self, AVISO_FILE):
 
         try:  # Extrapolate over land points with fillmask
             zeta = fillmask(zeta, self.mask == 1)
-            #zeta = fillmask(zeta, 1 + (-1 * zeta.mask))
         except Exception:  # In case no landpoints
             zeta = np.ma.masked_array(zeta)
         return zeta.astype(np.float64)
@@ -847,7 +845,7 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
     else:
         PRODUCT = 'AVISO_DT10'
         DAYS_BTWN_RECORDS = 7.  # old seven day AVISO
-
+    config['DAYS_BTWN_RECORDS'] = DAYS_BTWN_RECORDS
     AVISO_FILES = config['AVISO']['AVISO_FILES']
 
     #TRACK_DURATION_MIN = config['TRACK_DURATION_MIN']
@@ -972,10 +970,10 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
 
     # See Chelton section B2 (0.4 degree radius)
     # These should give 8 and 1000 for 0.25 deg resolution
-    PIXMIN = np.round((np.pi * config['RADMIN']**2) /
-                      sla_grd.get_resolution()**2)
-    PIXMAX = np.round((np.pi * config['RADMAX']**2) /
-                      sla_grd.get_resolution()**2)
+    PIXMIN = np.round((np.pi * config['RADMIN'] ** 2) /
+                      sla_grd.get_resolution() ** 2)
+    PIXMAX = np.round((np.pi * config['RADMAX'] ** 2) /
+                      sla_grd.get_resolution() ** 2)
     print '--- Pixel range = %s-%s' % (np.int(PIXMIN),
                                        np.int(PIXMAX))
 
@@ -1093,7 +1091,7 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
                                    xi / sla_grd.f()[sla_grd.jup0:sla_grd.jup1,
                                                     sla_grd.iup0:sla_grd.iup1])
 
-                xicopy = np.ma.copy(xi)
+                xicopy = xi.copy()
 
             elif 'SLA' in DIAGNOSTIC_TYPE:
 
@@ -1103,11 +1101,11 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
                 C_eddy.slacopy = sla.copy()
 
             # Get scalar speed
-            uspd = np.sqrt(sla_grd.u**2 + sla_grd.v**2)
+            uspd = np.sqrt(sla_grd.u ** 2 + sla_grd.v ** 2)
             uspd = np.ma.masked_where(
                 sla_grd.mask[sla_grd.jup0:sla_grd.jup1,
                              sla_grd.iup0:sla_grd.iup1] == 0,
-                uspd)
+                             uspd)
             A_eddy.uspd = uspd.copy()
             C_eddy.uspd = uspd.copy()
 
@@ -1121,7 +1119,7 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
             # Get contours of Q/sla parameter
             if 'first_record' not in locals():
 
-                print '------ processing SLA contours for eddies'
+                print('------ processing SLA contours for eddies')
                 contfig = plt.figure(99)
                 ax = contfig.add_subplot(111)
 
@@ -1252,7 +1250,7 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
                 A_eddy.set_old_variables()
                 C_eddy.set_old_variables()
                 START = False
-                print '------ tracking eddies'
+                print('------ tracking eddies')
             else:
                 first_record = False
 
@@ -1278,9 +1276,9 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
             if not first_record:
 
                 if A_eddy.VERBOSE:
-                    print '--- saving to nc', A_eddy.SAVE_DIR
-                    print '--- saving to nc', C_eddy.SAVE_DIR
-                    print '+++'
+                    print('--- saving to nc', A_eddy.SAVE_DIR)
+                    print('--- saving to nc', C_eddy.SAVE_DIR)
+                    print('+++')
 
                 A_eddy.write2netcdf(rtime)
                 C_eddy.write2netcdf(rtime)
diff --git a/make_eddy_track_ROMS.py b/make_eddy_track_ROMS.py
@@ -24,7 +24,7 @@
 
 make_eddy_track_ROMS.py
 
-Version 2.0.0
+Version 2.0.3
 
 
 
diff --git a/make_eddy_tracker_list_obj.py b/make_eddy_tracker_list_obj.py
@@ -25,7 +25,7 @@
 
 make_eddy_tracker_list_obj.py
 
-Version 2.0.1
+Version 2.0.3
 
 
 ===========================================================================
@@ -483,7 +483,7 @@ def get_active_tracks(self, rtime):
                 active_tracks.append(i)
         return active_tracks
 
-    def get_inactive_tracks(self, rtime):
+    def get_inactive_tracks(self, rtime, stopper=0):
         """
         Return list of indices to inactive tracks.
         This call also identifies and removes
@@ -495,6 +495,15 @@ def get_inactive_tracks(self, rtime):
                 inactive_tracks.append(i)
         return inactive_tracks
 
+    def kill_all_tracks(self):
+        """
+        Mark all tracks as not alive
+        """
+        for track in self.tracklist:
+            track.alive = False
+        print('------ all %s tracks killed for final saving'
+               % self.SIGN_TYPE.replace('one', 'onic').lower())
+
     def create_netcdf(self, directory, savedir,
                       grd=None, YMIN=None, YMAX=None,
                       MMIN=None, MMAX=None, MODEL=None,
@@ -742,13 +751,17 @@ def _remove_inactive_tracks(self):
         self.tracklist = tracklist.tolist()
         return alive_inds
 
-    def write2netcdf(self, rtime):
+    def write2netcdf(self, rtime, stopper=0):
         """
         Write inactive tracks to netcdf file.
         'ncind' is important because prevents writing of
         already written tracks.
         Each inactive track is 'emptied' after saving
+        
+        rtime - current timestamp
+        stopper - dummy value (either 0 or 1)
         """
+        rtime += stopper
         tracks2save = np.array([self.get_inactive_tracks(rtime)])
         DBR = self.DAYS_BTWN_RECORDS
 
@@ -833,7 +846,12 @@ def write2netcdf(self, rtime):
                             self.ncind += tsize
                             self.ch_index += 1
                             nc.sync()
-        
+
+        # Print final message and return
+        if stopper:
+            print('All %ss saved' % self.SIGN_TYPE.replace('one', 'onic').lower())
+            return
+
         # Get index to first currently active track
         #try:
             #lasti = self.get_active_tracks(rtime)[0]
diff --git a/py_eddy_tracker_classes.py b/py_eddy_tracker_classes.py
@@ -24,7 +24,7 @@
 
 py_eddy_tracker_classes.py
 
-Version 2.0.1
+Version 2.0.3
 ===========================================================================
 
 
@@ -376,7 +376,7 @@ def fit_circle(xvec, yvec):
 
     # Indices of polygon points outside circle
     # p_inon_? : polygon x or y points inside & on the circle
-    pout_id = np.nonzero(dist_poly > r)
+    pout_id = (dist_poly > r).nonzero()
 
     p_inon_x = xvec  # init
     p_inon_y = yvec  # init
@@ -954,9 +954,10 @@ def track_eddies(Eddy, first_record):
         new_cntr_e, new_cntr_s, new_uavg_prf = [], [], []
         new_shp_err = np.array([])
 
-        new_inds, backup_ind, backup_dist = np.array([], dtype=np.int16), \
-                                            np.array([], dtype=np.int16), \
-                                            np.array([])
+        #new_inds, backup_ind, backup_dist = np.array([], dtype=np.int16), \
+                                            #np.array([], dtype=np.int16), \
+                                            #np.array([])
+        backup_ind = np.array([], dtype=np.int16)
 
         '''
         # See http://stackoverflow.com/questions/11528078/determining-duplicate-values-in-an-array
@@ -983,6 +984,7 @@ def track_eddies(Eddy, first_record):
             plt.xlabel('new')
             plt.ylabel('old')
             plt.title('original')
+            plt.show()
 
         # Make an ellipse at current old_eddy location
         # (See CSS11 sec. B4, pg. 208)
@@ -997,7 +999,7 @@ def track_eddies(Eddy, first_record):
         for new_ind in new_sparse_inds:
             
             new_dist = dist_mat[old_ind, new_ind]
-            
+
             within_range = False
             #print far_away, Eddy.search_ellipse.rw_c_mod, new_dist
             if new_dist < Eddy.search_ellipse.rw_c_mod:#far_away:
@@ -1018,54 +1020,38 @@ def track_eddies(Eddy, first_record):
 
             # Pass only the eddies within ellipse or sep_dist
             if within_range:
+                
+                dist_arr = np.r_[dist_arr, new_dist]
+                new_ln.append(Eddy.new_lon_tmp[new_ind])
+                new_lt.append(Eddy.new_lat_tmp[new_ind])
+                new_rd_s.append(Eddy.new_radii_tmp_s[new_ind])
+                new_rd_e.append(Eddy.new_radii_tmp_e[new_ind])
+                new_am.append(Eddy.new_amp_tmp[new_ind])
+                new_Ua.append(Eddy.new_uavg_tmp[new_ind])
+                new_ek.append(Eddy.new_teke_tmp[new_ind])
+                new_tm.append(Eddy.new_time_tmp[new_ind])
 
-                old_amplitude[:] = Eddy.old_amp[old_ind]
-                new_amplitude[:] = Eddy.new_amp_tmp[new_ind]
-
-                # Following CSS11, we use effective radius rather than speed based...
-                old_area[:] = Eddy.old_radii_s[old_ind]**2
-                old_area *= np.pi
-                new_area[:] = Eddy.new_radii_tmp_s[new_ind]**2
-                new_area *= np.pi
-
-                # Pass only the eddies within min and max times old amplitude
-                # and area (KCCMC11 and CSS11 use 0.25 and 2.5, respectively)
-                if (new_amplitude >= (Eddy.EVOLVE_AMP_MIN * old_amplitude) and
-                    new_amplitude <= (Eddy.EVOLVE_AMP_MAX * old_amplitude)) and \
-                   (new_area >= (Eddy.EVOLVE_AREA_MIN * old_area) and
-                    new_area <= (Eddy.EVOLVE_AREA_MAX * old_area)):
-
-                    dist_arr = np.r_[dist_arr, new_dist]
-                    new_ln.append(Eddy.new_lon_tmp[new_ind])
-                    new_lt.append(Eddy.new_lat_tmp[new_ind])
-                    new_rd_s.append(Eddy.new_radii_tmp_s[new_ind])
-                    new_rd_e.append(Eddy.new_radii_tmp_e[new_ind])
-                    new_am.append(Eddy.new_amp_tmp[new_ind])
-                    new_Ua.append(Eddy.new_uavg_tmp[new_ind])
-                    new_ek.append(Eddy.new_teke_tmp[new_ind])
-                    new_tm.append(Eddy.new_time_tmp[new_ind])
-
-                    if 'ROMS' in Eddy.PRODUCT:
-                        #new_tp = np.r_[new_tp, Eddy.new_temp_tmp[new_ind]]
-                        #new_st = np.r_[new_st, Eddy.new_salt_tmp[new_ind]]
-                        pass
-
-                    if Eddy.TRACK_EXTRA_VARIABLES:
-                        new_cntr_e.append(Eddy.new_contour_e_tmp[new_ind])
-                        new_cntr_s.append(Eddy.new_contour_s_tmp[new_ind])
-                        new_uavg_prf.append(Eddy.new_uavg_profile_tmp[new_ind])
-                        new_shp_err = np.r_[new_shp_err,
-                                            Eddy.new_shape_error_tmp[new_ind]]
-
-                    new_inds = np.r_[new_inds, new_ind]
-                    backup_ind = np.r_[backup_ind, new_ind]
-
-                    # An old (active) eddy has been detected, so
-                    # corresponding new_eddy_inds set to False
-                    new_eddy_inds[np.nonzero(Eddy.new_lon_tmp ==
-                                             Eddy.new_lon_tmp[new_ind])] = False
-
-                    dist_mat[:, new_ind] = far_away # km
+
+                if 'ROMS' in Eddy.PRODUCT:
+                    #new_tp = np.r_[new_tp, Eddy.new_temp_tmp[new_ind]]
+                    #new_st = np.r_[new_st, Eddy.new_salt_tmp[new_ind]]
+                    pass
+
+                if Eddy.TRACK_EXTRA_VARIABLES:
+                    new_cntr_e.append(Eddy.new_contour_e_tmp[new_ind])
+                    new_cntr_s.append(Eddy.new_contour_s_tmp[new_ind])
+                    new_uavg_prf.append(Eddy.new_uavg_profile_tmp[new_ind])
+                    new_shp_err = np.r_[new_shp_err,
+                                        Eddy.new_shape_error_tmp[new_ind]]
+
+                backup_ind = np.r_[backup_ind, new_ind]
+
+                # An old (active) eddy has been detected, so
+                # corresponding new_eddy_inds set to False
+                new_eddy_inds[np.nonzero(Eddy.new_lon_tmp ==
+                                            Eddy.new_lon_tmp[new_ind])] = False
+
+                dist_mat[:, new_ind] = far_away # km
 
         if Eddy.TRACK_EXTRA_VARIABLES:
             kwargs = {'contour_e': new_cntr_e, 'contour_s': new_cntr_s,
@@ -1075,7 +1061,6 @@ def track_eddies(Eddy, first_record):
 
         # Only one eddy within range
         if dist_arr.size == 1:  # then update the eddy track only
-
             # Use index 0 here because type is list and we want the scalar
             args = (Eddy, old_ind, new_ln[0], new_lt[0], new_rd_s[0],
                     new_rd_e[0], new_am[0], new_Ua[0], new_ek[0], new_tm[0],
@@ -1103,8 +1088,8 @@ def track_eddies(Eddy, first_record):
 
                 # Choice of using effective or speed-based...
                 delta_area_tmp = np.array([np.pi *
-                                           (Eddy.old_radii_e[old_ind]**2),
-                                           np.pi * (new_rd_e[i]**2)]).ptp()
+                                           (Eddy.old_radii_e[old_ind] ** 2),
+                                           np.pi * (new_rd_e[i] ** 2)]).ptp()
                 delta_amp_tmp = np.array([Eddy.old_amp[old_ind],
                                           new_am[i]]).ptp()
                 delta_area = np.r_[delta_area, delta_area_tmp]
@@ -1120,9 +1105,9 @@ def track_eddies(Eddy, first_record):
                         #np.abs(np.diff([Eddy.old_salt[old_ind], new_st[i]]))]
 
             # This from Penven etal (2005)
-            deltaX = np.sqrt((delta_area / AREA0)**2 +
-                             (delta_amp / AMP0)**2 +
-                             (dist_arr / DIST0)**2)
+            deltaX = np.sqrt((delta_area / AREA0) ** 2 +
+                             (delta_amp / AMP0) ** 2 +
+                             (dist_arr / DIST0) ** 2)
 
             dx = deltaX.argsort()
             dx0 = dx[0]  # index to the nearest eddy
diff --git a/py_eddy_tracker_property_classes.py b/py_eddy_tracker_property_classes.py
@@ -25,7 +25,7 @@
 
 py_eddy_tracker_amplitude.py
 
-Version 2.0.1
+Version 2.0.3
 
 
 ===========================================================================
@@ -297,6 +297,8 @@ def debug_figure(self, grd):
             #cmin, cmax = (self.sla.min(), self.sla.max())
             cm = plt.cm.gist_ncar
             #cm = plt.cm.hsv
+            print self.h0_check.shape
+            print self.contlon.shape
             
             plt.title('Local max/min count: %s, Amp: %s' % (
                              self.local_extrema, self.amplitude))
diff --git a/roms_grid.py b/roms_grid.py

Original file line number	Diff line number	Diff line change
`@@ -27,7 +27,7 @@ module haversine`
`27`	`27`	`!`
`28`	`28`	`! If you have ifort on your system, change 'gfortran' to 'intelem'!`
`29`	`29`	`!`
`30`		`-! Version 2.0.1`
	`30`	`+! Version 2.0.3`
`31`	`31`	`!`
`32`	`32`	`!===========================================================================`
`33`	`33`	`implicit none`
Original file line number	Diff line number	Diff line change
`@@ -24,7 +24,7 @@`
`24`	`24`
`25`	`25`	`make_eddy_track_ROMS.py`
`26`	`26`
`27`		`-Version 2.0.0`
	`27`	`+Version 2.0.3`
`28`	`28`
`29`	`29`
`30`	`30`