Bug with empty mask caught

smfm-project · Apr 9, 2018 · a723040 · a723040
1 parent 7922e8a
commit a723040
Show file tree

Hide file tree

Showing 2 changed files with 12 additions and 9 deletions.
diff --git a/biota/core.py b/biota/core.py
@@ -473,11 +473,13 @@ def getGamma0(self, polarisation = 'HV', units = 'natural', output = False, show
         if self.lee_filter:
             gamma0 = biota.filter.enhanced_lee_filter(gamma0, n_looks = self.nLooks)
 
+
         # Convert to natural units where specified
         if units == 'natural': gamma0 = 10 ** (gamma0 / 10.)
 
         # Keep masked values tidy
         gamma0.data[self.mask] = self.nodata
+        gamma0.mask = self.mask
 
         if output: self.__outputGeoTiff(gamma0, 'Gamma0')
 
@@ -548,6 +550,7 @@ def getWoodyCover(self, output = False, show = False):
                 WoodyCover.data[contiguous_area == False] = False
 
             WoodyCover.data[self.mask] = False
+            WoodyCover.mask = self.mask
 
             # Save output to class
             self.WoodyCover = WoodyCover
@@ -582,7 +585,7 @@ def getForestPatches(self, output = False, show = False):
             _, ForestPatches = biota.indices.getContiguousAreas(WoodyCover, True, min_pixels = min_pixels)
 
             # Tidy up masked pixels     
-            ForestPatches.data[ForestPatches.mask] = self.nodata
+            ForestPatches.data[np.ma.getmaskarray(ForestPatches)] = self.nodata
 
             # Save output to class
             self.ForestPatches = ForestPatches

diff --git a/biota/indices.py b/biota/indices.py
@@ -28,8 +28,8 @@ def getContiguousAreas(data, value, min_pixels = 1):
 
     # If any pixels are masked, we give them the value of the nearest valid pixel.
     if masked:
-        mask = data.mask
-        ind = ndimage.distance_transform_edt(data.mask, return_distances = False, return_indices = True)
+        mask = np.ma.getmaskarray(data)
+        ind = ndimage.distance_transform_edt(mask, return_distances = False, return_indices = True)
         data = data.data[tuple(ind)]
 
     # Extract area that meets condition
@@ -160,13 +160,13 @@ def _computeLDI(unique_ids):
         Calculates an index of landscape heterogeneity
         '''
 
-        # If masked array, separate mask
+        # If masked array, separate data from mask
         if np.ma.isMaskedArray(unique_ids):
-            mask = unique_ids.mask
-            unique_ids = unique_ids.data
+            mask = np.ma.getmask(unique_ids)
+            unique_ids = np.ma.getdata(unique_ids)
         else:
             mask = np.zeros_like(unique_ids, dtype=np.bool)
-
+                
         # Calculate LDI
         selection_1 = np.zeros_like(unique_ids,dtype = np.bool).ravel()
         selection_2 = np.zeros_like(unique_ids,dtype = np.bool).ravel()
@@ -228,7 +228,7 @@ def _computeLDI(unique_ids):
         #  If at least 50 % of data is present...
         if this_patch.mask.sum() <= ((patch_size ** 2) * 0.5):
 
-            # Compute LDI
+            ## Compute LDI
             LDI.data[n, m] = _computeLDI(unique_ids)
 
         else:
@@ -262,7 +262,7 @@ def calculateLDIChange(tile_change, patch_size = 'auto', output = False, show =
     from osgeo import gdal
 
     patch_size = _calculatePatchSize(tile_change, patch_size)
-
+    
     LDI_t1 = calculateLDI(tile_change.tile_t1, patch_size = patch_size)
     LDI_t2 = calculateLDI(tile_change.tile_t2, patch_size = patch_size)