Updated soil moisture interpolation routines

biota/SM.py

@@ -85,7 +85,7 @@ def _interpolateTime(sm):
 
     return sm_interp
 
-def _resampleSM(sm_interp, tile, geo_t):
+def _resampleSM(sm_interp, tile, geo_t, interpolation = 'avearge'):
     '''
     '''
 
@@ -101,16 +101,22 @@ def _resampleSM(sm_interp, tile, geo_t):
     ds_dest.SetGeoTransform(tile.geo_t)
     ds_dest.SetProjection(tile.proj)
 
-    # Reproject input GeoTiff to match the ALOS tile
-    gdal.ReprojectImage(ds_source, ds_dest, tile.proj, tile.proj, gdal.GRA_NearestNeighbour)
+    # Select interpolation type
+    if interpolation == 'average' or interpolation == 'nearest':
+        interp_type = gdal.GRA_NearestNeighbour
+    elif interpolation == 'cubic':
+        interp_type = gdal.GRA_CubicSpline
 
+    # Reproject input GeoTiff to match the ALOS tile    
+    gdal.ReprojectImage(ds_source, ds_dest, tile.proj, tile.proj, interp_type)
+
     # Load resampled image into memory
     sm_resampled = ds_dest.GetRasterBand(1).ReadAsArray()
 
     return np.ma.array(sm_resampled, mask = sm_resampled == -9999.)
 
 
-def _loadSM(tile, date, search_days = 7):
+def _loadSM(tile, date, search_days = 7, interpolation = 'average'):
     '''
     Load and reproject soil moisture data
     '''
@@ -122,12 +128,12 @@ def _loadSM(tile, date, search_days = 7):
     sm_interp = _interpolateTime(sm)
 
     # Reproject to match ALOS tile (nearest neighbor)
-    sm_resampled = _resampleSM(sm_interp, tile, geo_t)
+    sm_resampled = _resampleSM(sm_interp, tile, geo_t, interpolation = interpolation)
 
     return sm_resampled
 
 
-def getSM(tile, search_days = 7):
+def getSM(tile, search_days = 7, interpolation = 'average'):
     """
     Function to load a resampled soil moisture image from the ESA CCI soil moisture product. The function returns a masked array with estimated volumetric soil moisture content (m^2/m^2), averaged for each satellite overpass in a tile.
     
@@ -137,7 +143,9 @@ def getSM(tile, search_days = 7):
     Returns:
         A masked array of estimated soil moisture (m^2/m^2)
     """
-
+
+    assert interpolation in ['average', 'nearest', 'cubic'], "Soil moisture interpolation type must be 'average', 'nearest', or 'cubic'."
+
     # Build an array of dates in string format, which is quicker to search
     dates = tile.getDate().astype(np.int)
 
@@ -151,12 +159,17 @@ def getSM(tile, search_days = 7):
             continue
 
         # Interpolate through time and draw out the central value, and zoom to scale of tile
-        sm_resampled = _loadSM(tile, dt.datetime(1970,1,1).date() + dt.timedelta(date), search_days = search_days)
+        sm_resampled = _loadSM(tile, dt.datetime(1970,1,1).date() + dt.timedelta(date), search_days = search_days, interpolation = interpolation)
 
         # Only output soil moisture where > 25 % data exists        
         if float((sm_resampled.mask[dates == date]).sum()) / (sm_resampled.mask[dates == date]).shape[0] <= 0.25:
-            out[dates == date] = np.ma.mean(sm_resampled[dates == date])        
-
+
+            if interpolation == 'average':
+                # Calculate the mean average of values in overpass
+                out[dates == date] = np.ma.mean(sm_resampled[dates == date])
+            elif interpolation == 'nearest' or interpolation == 'cubic':
+                out[dates == date] = sm_resampled[dates == date]
+
     # Add the mask
     out = np.ma.array(out, mask = np.logical_or(tile.mask, out == 999999.))
 

biota/core.py

@@ -38,7 +38,7 @@ class LoadTile(object):
         mask:
     """
 
-    def __init__(self, data_dir, lat, lon, year, forest_threshold = 10., area_threshold = 0., downsample_factor = 1, lee_filter = False, window_size = 5, sm_dir = os.getcwd(), output_dir = os.getcwd()):
+    def __init__(self, data_dir, lat, lon, year, forest_threshold = 10., area_threshold = 0., downsample_factor = 1, lee_filter = False, window_size = 5, sm_dir = os.getcwd(), sm_interpolation = 'average', output_dir = os.getcwd()):
         """
         Loads data and metadata for an ALOS mosaic tile.
         """
@@ -56,6 +56,7 @@ def __init__(self, data_dir, lat, lon, year, forest_threshold = 10., area_thresh
         assert window_size % 2 == 1, "Option window_size must be an odd integer."
         assert os.path.isdir(os.path.expanduser(data_dir)), "Specified data directory (%s) does not exist"%str(data_dir)
         assert os.path.isdir(os.path.expanduser(sm_dir)), "Specified soil moisture directory (%s) does not exist"%str(sm_dir)
+        assert sm_interpolation in ['nearest', 'average', 'cubic'], "Soil moisture interpolation type must be one of 'nearest', 'average' or 'cubic'."
         assert os.path.isdir(os.path.expanduser(output_dir)), "Specified output directory (%s) does not exist"%str(output_dir)
         assert type(forest_threshold) == float or type(forest_threshold) == int, "Forest threshold must be numeric."
         assert type(area_threshold) == float or type(area_threshold) == int, "Area threshold must be numeric."
@@ -86,6 +87,7 @@ def __init__(self, data_dir, lat, lon, year, forest_threshold = 10., area_thresh
 
         # Set up soil moisture data location
         self.SM_dir = os.path.expanduser(sm_dir.rstrip('/'))
+        self.SM_interpolation = sm_interpolation
 
         # Set up locations for file output
         self.output_pattern = self.__getOutputPattern()
@@ -531,15 +533,36 @@ def getDOY(self, output = False, show = False):
 
         return self.DOY       
 
-    def getSM(self, output = False, show = False):
+    def getYearArray(self, output = False, show = False):
+        """
+        Loads year of overpass into a numpy array
+        """
+
+        # Don't rerun processing if already present in memory
+        if not hasattr(self, 'YearArray'):
+
+            dates = self.getDate()       
+
+            YearArray = dates.astype('datetime64[Y]').astype(int) + 1970
+
+            # Save output to class
+            self.YearArray = YearArray
+
+        if output: self.__outputGeoTiff(self.YearArray, 'YearArray', dtype = gdal.GDT_Int32)
+
+        if show: self.__showArray(self.YearArray, title = 'YearArray', cbartitle = 'Year', vmin = self.year-1, vmax = self.year+1, cmap = 'Spectral')
+
+        return self.YearArray
+
+    def getSM(self, output = False, show = False, search_days = 7):
         """
         Loads a soil moisture map using the ESA CCI soil moisture product.
         """
 
         # Don't rerun processing if already present in memory
         if not hasattr(self, 'SM'):
 
-            SM = biota.SM.getSM(self)
+            SM = biota.SM.getSM(self, search_days = search_days, interpolation = self.SM_interpolation)
 
             # Keep masked values tidy
             SM.data[self.mask] = self.nodata

biota/indices.py

@@ -388,7 +388,7 @@ def calculateProportionalChange(tile_change, change_type, patch_size = 'auto', o
 
             # Calculate proportion of patch subject to change_type
             change_downsampled.data[n, m] = int(round((float(C.sum()) / ((patch_size ** 2) - C.mask.sum())) * 100))
-
+            
         else:
             change_downsampled.mask[n, m] = True