Last Jan CLS commit

Author: evanmason

make_eddy_track_AVISO.py

@@ -35,7 +35,7 @@
 from py_eddy_tracker_classes import plt, np, dt, Dataset, ndimage, time, \
                                     datestr2datetime, gaussian_resolution, \
                                     get_cax, collection_loop, track_eddies, \
-                                    anim_figure
+                                    anim_figure, interpolate
 import make_eddy_tracker_list_obj as eddy_tracker
 from dateutil import parser
 from mpl_toolkits.basemap import Basemap
@@ -505,8 +505,8 @@ def __init__(self, AVISO_FILE, LONMIN, LONMAX, LATMIN, LATMAX,
         self.set_u_v_eke()
         pad2 = 2 * self.pad
         self.shape = (self.f().shape[0] - pad2, self.f().shape[1] - pad2)
+    
 
-
     def get_AVISO_data(self, AVISO_FILE):
         """
         Read nc data from AVISO file
@@ -517,7 +517,7 @@ def get_AVISO_data(self, AVISO_FILE):
                 ssh1 = self.read_nc(AVISO_FILE, 'SLA',
                        indices='[:, self.jp0:self.jp1, :self.ip0]')
                 ssh0 = self.read_nc(AVISO_FILE, 'SLA',
-                       indices='[:,self.jp0:self.jp1, self.ip1:]')
+                       indices='[:, self.jp0:self.jp1, self.ip1:]')
                 ssh0, ssh1 = ssh0.squeeze(), ssh1.squeeze()
                 ssh0 *= 100. # m to cm
                 ssh1 *= 100. # m to cm
@@ -526,7 +526,7 @@ def get_AVISO_data(self, AVISO_FILE):
                 ssh1 = self.read_nc(AVISO_FILE, 'Grid_0001',
                        indices='[:self.ip0, self.jp0:self.jp1]').T
                 ssh0 = self.read_nc(AVISO_FILE, 'Grid_0001',
-                       indices='[self.ip1:,self.jp0:self.jp1]').T
+                       indices='[self.ip1:, self.jp0:self.jp1]').T
 
             zeta = np.ma.concatenate((ssh0, ssh1), axis=1)
 
@@ -591,19 +591,19 @@ def fillmask(self, x, mask):
 
         # Multiply the queried good points by the weight, selecting only the
         # nearest points. Divide by the number of nearest points to get average
-        xfill = weight * x[igood[:,0][iquery], igood[:,1][iquery]]
-        xfill = (xfill / weight.sum(axis=1)[:,np.newaxis]).sum(axis=1)
+        xfill = weight * x[igood[:, 0][iquery], igood[:,1][iquery]]
+        xfill = (xfill / weight.sum(axis=1)[:, np.newaxis]).sum(axis=1)
 
         # Place average of nearest good points, xfill, into bad point locations
-        x[ibad[:,0], ibad[:,1]] = xfill
+        x[ibad[:, 0], ibad[:, 1]] = xfill
         return x
 
 
     def lon(self):
         if self.ZERO_CROSSING:
             # TO DO: These concatenations are possibly expensive, they
             # shouldn't need to happen with every call to self.lon()
-            lon0 = self._lon[self.j0:self.j1,  self.i1:]
+            lon0 = self._lon[self.j0:self.j1, self.i1:]
             lon1 = self._lon[self.j0:self.j1, :self.i0]
             print 'fix this so called only once'
             return np.concatenate((lon0 - 360., lon1), axis=1)
@@ -612,23 +612,23 @@ def lon(self):
 
     def lat(self):
         if self.ZERO_CROSSING:
-            lat0 = self._lat[self.j0:self.j1,  self.i1:]
+            lat0 = self._lat[self.j0:self.j1, self.i1:]
             lat1 = self._lat[self.j0:self.j1, :self.i0]
             return np.concatenate((lat0, lat1), axis=1)
         else:            
             return self._lat[self.j0:self.j1, self.i0:self.i1]
 
     def lonpad(self):
         if self.ZERO_CROSSING:
-            lon0 = self._lon[self.jp0:self.jp1,  self.ip1:]
+            lon0 = self._lon[self.jp0:self.jp1, self.ip1:]
             lon1 = self._lon[self.jp0:self.jp1, :self.ip0]
             return np.concatenate((lon0 - 360., lon1), axis=1)
         else:
             return self._lon[self.jp0:self.jp1, self.ip0:self.ip1]
 
     def latpad(self):
         if self.ZERO_CROSSING:
-            lat0 = self._lat[self.jp0:self.jp1,  self.ip1:]
+            lat0 = self._lat[self.jp0:self.jp1, self.ip1:]
             lat1 = self._lat[self.jp0:self.jp1, :self.ip0]
             return np.concatenate((lat0, lat1), axis=1)
         else:            
@@ -665,7 +665,7 @@ def pn(self): # Reciprocal of dy
 
     def get_resolution(self):
         return np.sqrt(np.diff(self.lon()[1:], axis=1) *
-                       np.diff(self.lat()[:,1:], axis=0)).mean()
+                       np.diff(self.lat()[:, 1:], axis=0)).mean()
 
 
     def boundary(self):
@@ -676,10 +676,10 @@ def boundary(self):
         Returns:
           lon/lat boundary points
         """
-        lon = np.r_[(self.lon()[:,0],     self.lon()[-1],
-                     self.lon()[::-1,-1], self.lon()[0,::-1])]
-        lat = np.r_[(self.lat()[:,0],     self.lat()[-1],
-                     self.lat()[::-1,-1], self.lat()[0,::-1])]
+        lon = np.r_[(self.lon()[:, 0], self.lon()[-1],
+                     self.lon()[::-1, -1], self.lon()[0, ::-1])]
+        lat = np.r_[(self.lat()[:, 0], self.lat()[-1],
+                     self.lat()[::-1, -1], self.lat()[0, ::-1])]
         return lon, lat
 
 
@@ -693,32 +693,15 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
         return path.Path(brypath)
 
 
-    #def pcol_2dxy(self, x, y):
-        #"""
-        #Function to shift x, y for subsequent use with pcolor
-        #by Jeroen Molemaker UCLA 2008
-        #"""
-        #Mp, Lp = x.shape
-        #M = Mp - 1
-        #L = Lp - 1
-        #x_pcol = np.zeros((Mp, Lp))
-        #y_pcol = np.zeros((Mp, Lp))
-        #x_tmp = self.half_interp(x[:,:L], x[:,1:Lp])
-        #x_pcol[1:Mp,1:Lp] = self.half_interp(x_tmp[0:M,:], x_tmp[1:Mp,:])
-        #x_pcol[0,:] = 2. * x_pcol[1,:] - x_pcol[2,:]
-        #x_pcol[:,0] = 2. * x_pcol[:,1] - x_pcol[:,2]
-        #y_tmp = self.half_interp(y[:,0:L], y[:,1:Lp]    )
-        #y_pcol[1:Mp,1:Lp] = self.half_interp(y_tmp[0:M,:], y_tmp[1:Mp,:])
-        #y_pcol[0,:] = 2. * y_pcol[1,:] - y_pcol[2,:]
-        #y_pcol[:,0] = 2. * y_pcol[:,1] - y_pcol[:,2]
-        #return x_pcol, y_pcol
-
-
-
-
-
-
-
+        
+    def set_interp_coeffs(self, sla, uspd):
+        """
+        """
+        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
 
 
 #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@@ -840,16 +823,16 @@ def brypath(self, imin=0, imax=-1, jmin=0, jmax=-1):
     #--------------------------------------------------------------------------
 
     assert DATE_STR < DATE_END, 'DATE_END must be larger than DATE_STR'
-    assert DIAGNOSTIC_TYPE in ('Q', 'SLA'), 'DIAGNOSTIC_TYPE not properly defined'
+    assert DIAGNOSTIC_TYPE in ('Q', 'SLA'), 'Undefined DIAGNOSTIC_TYPE'
 
     thestartdate = dt.date2num(datestr2datetime(str(DATE_STR)))
     theenddate = dt.date2num(datestr2datetime(str(DATE_END)))
 
     # Get complete AVISO file list
     AVISO_FILES = sorted(glob.glob(DATA_DIR + AVISO_FILES))
 
-    # Use this for subsampling to get identical list as old_AVISO
-    #AVISO_FILES = AVISO_FILES[5:-5:7]
+    # For subsampling to get identical list as old_AVISO use:
+    # AVISO_FILES = AVISO_FILES[5:-5:7]
     if AVISO_DT14 and AVISO_DT14_SUBSAMP:
         AVISO_FILES = AVISO_FILES[5:-5:np.int(DAYS_BTWN_RECORDS)]
 
@@ -1118,14 +1101,20 @@ 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.ma.masked_where(
+            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] == False,
-                                      Uspd)
-            A_eddy.Uspd = Uspd.copy()
-            C_eddy.Uspd = Uspd.copy()
+                                      uspd)
+            A_eddy.uspd = uspd.copy()
+            C_eddy.uspd = uspd.copy()
 
+            # Set interpolation coefficients
+            sla_grd.set_interp_coeffs(sla, uspd)
+            A_eddy.sla_coeffs = sla_grd.sla_coeffs
+            A_eddy.uspd_coeffs = sla_grd.uspd_coeffs
+            C_eddy.sla_coeffs = sla_grd.sla_coeffs
+            C_eddy.uspd_coeffs = sla_grd.uspd_coeffs
 
             # Get contours of Q/sla parameter
             if 'first_record' not in locals():

py_eddy_tracker_amplitude.py

@@ -111,10 +111,7 @@ def __init__(self, contlon, contlat, eddy, grd):
         self.jmin, self.jmax = self.eddy.jmin, self.eddy.jmax
         self.sla = self.eddy.sla[self.jmin:self.jmax,
                                  self.imin:self.imax].copy()
-        self.rbspline = interpolate.RectBivariateSpline(
-                          self.eddy.grd.lat()[self.jmin:self.jmax, 0],
-                          self.eddy.grd.lon()[0, self.imin:self.imax],
-                          self.sla, kx=1, ky=1)
+        self.rbspline = grd.sla_coeffs
         h0 = self.rbspline.ev(self.contlat, self.contlon)
         self.h0 = h0[np.isfinite(h0)].mean()
         self.amplitude = np.atleast_1d(0.)

py_eddy_tracker_classes.py

@@ -473,10 +473,10 @@ def get_uavg(Eddy, CS, collind, centlon_e, centlat_e, poly_eff,
     
     If save_all_uavg == True we want uavg for every contour
     """
-    def calc_uavg(rbspline, lon, lat):
-        """
-        """
-        return rbspline.ev(lon, lat).mean()
+    #def calc_uavg(rbspline, lon, lat):
+        #"""
+        #"""
+        #return rbspline.ev(lon, lat).mean()
 
     # True for debug figures
     #debug_U = False
@@ -495,20 +495,25 @@ def calc_uavg(rbspline, lon, lat):
     points = np.array([grd.lon()[jmin:jmax,imin:imax].ravel(),
                        grd.lat()[jmin:jmax,imin:imax].ravel()]).T
 
-    rbspline = interpolate.RectBivariateSpline(
-                           grd.lat()[jmin:jmax,0],
-                           grd.lon()[0,imin:imax],
-                           Eddy.Uspd[jmin:jmax,imin:imax], kx=1, ky=1)
+    #rbspline = interpolate.RectBivariateSpline(
+                           #grd.lat()[jmin:jmax,0],
+                           #grd.lon()[0,imin:imax],
+                           #Eddy.uspd[jmin:jmax,imin:imax], kx=1, ky=1)
 
     # First contour is the outer one (effective)
     theseglon, theseglat = poly_eff.vertices[:,0].copy(), \
                            poly_eff.vertices[:,1].copy()
     theseglon, theseglat = eddy_tracker.uniform_resample(
                                 theseglon, theseglat, method='akima')
-    uavg = calc_uavg(rbspline, theseglon[:-1], theseglat[:-1])
+    uavg = Eddy.uspd_coeffs.ev(theseglon[:-1], theseglat[:-1]).mean()
 
     if save_all_uavg:
         all_uavg = [uavg]
+        pixel_min = 1 # iterate until 1 pixel
+    
+    else:
+        # Iterate until PIXEL_THRESHOLD[0] number of pixels
+        pixel_min = Eddy.PIXEL_THRESHOLD[0]
 
     start = True
     citer = np.nditer(CS.cvalues, flags=['f_index'])
@@ -526,16 +531,10 @@ def calc_uavg(rbspline, lon, lat):
             # 1. Ensure polygon_i is within polygon_e
             # 2. Ensure polygon_i contains point centlon_e, centlat_e
             # 3. Respect size range
-            if not save_all_uavg:
-                # Iterate until PIXEL_THRESHOLD[0] number of pixels
-                pixel_min = Eddy.PIXEL_THRESHOLD[0]
-            else:
-                pixel_min = 1 # iterate until 1 pixel
-            
             if np.all([poly_eff.contains_path(poly_i),
                        poly_i.contains_point([centlon_e, centlat_e]),
-                       np.logical_and(mask_i_sum >= pixel_min,
-                                      mask_i_sum <= Eddy.PIXEL_THRESHOLD[1])]):
+                       (mask_i_sum >= pixel_min and
+                        mask_i_sum <= Eddy.PIXEL_THRESHOLD[1])]):
                 proceed = True
             else: 
                 proceed = False
@@ -547,9 +546,8 @@ def calc_uavg(rbspline, lon, lat):
                 seglon, seglat = eddy_tracker.uniform_resample(
                                       seglon, seglat, method='akima')
 
-                # Interpolate Uspd to seglon, seglat, then get mean
-                #uavgseg = calc_uavg(points, Eddy.Uspd[jmin:jmax,imin:imax], seglon[:-1], seglat[:-1])
-                uavgseg = calc_uavg(rbspline, seglon[:-1], seglat[:-1])
+                # Interpolate uspd to seglon, seglat, then get mean
+                uavgseg = Eddy.uspd_coeffs.ev(seglon[:-1], seglat[:-1]).mean()
 
                 if save_all_uavg:
                     all_uavg.append(uavgseg)
@@ -765,7 +763,7 @@ def collection_loop(CS, grd, rtime, A_list_obj, C_list_obj,
 
                                 if 'Q' in Eddy.DIAGNOSTIC_TYPE:
                                     #print 'change to rectbispline'
-                                    #uavg = interpolate.griddata(points, Eddy.Uspd[jmin:jmax,imin:imax].ravel(),
+                                    #uavg = interpolate.griddata(points, Eddy.uspd[jmin:jmax,imin:imax].ravel(),
                                                                 #(contlon_e, contlat_e), 'linear')
                                     #uavg = np.nan_to_num(uavg).max()
                                     ##uavg = 0; print 'fix me'