New masked_value parameter for cropping

telluric/georaster.py

@@ -85,7 +85,8 @@ class PixelStrategy(Enum):
 
 def merge_all(rasters, roi=None, dest_resolution=None, merge_strategy=MergeStrategy.UNION,
               shape=None, ul_corner=None, crs=None,
-              pixel_strategy=PixelStrategy.FIRST):
+              pixel_strategy=PixelStrategy.FIRST,
+              crop_masked_value=None):
     """Merge a list of rasters, cropping by a region of interest.
        There are cases that the roi is not precise enough for this cases one can use,
        the upper left corner the shape and crs to precisely define the roi.
@@ -100,7 +101,7 @@ def merge_all(rasters, roi=None, dest_resolution=None, merge_strategy=MergeStrat
         dtype=rasters[0].dtype, shape=shape, ul_corner=ul_corner, crs=crs)
 
     # Create a list of single band rasters
-    all_band_names, projected_rasters = _prepare_rasters(rasters, merge_strategy, empty)
+    all_band_names, projected_rasters = _prepare_rasters(rasters, merge_strategy, empty, crop_masked_value)
     assert len(projected_rasters) == len(rasters)
 
     prepared_rasters = _apply_pixel_strategy(projected_rasters, pixel_strategy)
@@ -172,8 +173,12 @@ def key(rs):
     return rasters_final
 
 
-def _prepare_rasters(rasters, merge_strategy, first):
-    # type: (List[GeoRaster2], MergeStrategy, GeoRaster2) -> Tuple[IndexedSet[str], List[Optional[_Raster]]]
+def _prepare_rasters(rasters,                # type: List[GeoRaster2]
+                     merge_strategy,         # type: MergeStrategy
+                     first,                  # type: GeoRaster2
+                     crop_masked_value=None  # type: Union[int, float]
+                     ):
+    # type: (...) -> Tuple[IndexedSet[str], List[Optional[_Raster]]]
     """Prepares the rasters according to the baseline (first) raster and the merge strategy.
 
     The baseline (first) raster is used to crop and reproject the other rasters,
@@ -185,7 +190,7 @@ def _prepare_rasters(rasters, merge_strategy, first):
     all_band_names = IndexedSet(first.band_names)
     projected_rasters = []
     for raster in rasters:
-        projected_raster = _prepare_other_raster(first, raster)
+        projected_raster = _prepare_other_raster(first, raster, crop_masked_value)
 
         # Modify the bands only if an intersecting raster was returned
         if projected_raster:
@@ -216,12 +221,14 @@ def _explode_raster(raster, band_names=[]):
     return [_Raster(image=raster.bands_data([band_name]), band_names=[band_name]) for band_name in band_names]
 
 
-def _prepare_other_raster(one, other):
-    # type: (GeoRaster2, GeoRaster2) -> Union[_Raster, None]
+def _prepare_other_raster(one, other, crop_masked_value=None):
+    # type: (GeoRaster2, GeoRaster2, Union[int, float]) -> Union[_Raster, None]
     # Crop and reproject the second raster, if necessary
     if not (one.crs == other.crs and one.affine.almost_equals(other.affine) and one.shape == other.shape):
         if one.footprint().intersects(other.footprint()):
-            other = other.crop(one.footprint(), resolution=one.resolution())
+            other = other.crop(one.footprint(),
+                               resolution=one.resolution(),
+                               masked_value=crop_masked_value)
             other = other._reproject(new_width=one.width, new_height=one.height,
                                      dest_affine=one.affine, dst_crs=one.crs,
                                      resampling=Resampling.nearest)
@@ -931,11 +938,14 @@ def type_min(dtype):
         dst_array = dst_array.astype(dst_type)
         return self.copy_with(image=dst_array)
 
-    def crop(self, vector, resolution=None, masked=True, bands=None):
+    def crop(self, vector, resolution=None, masked=True, masked_value=None, bands=None):
         """
         crops raster outside vector (convex hull)
         :param vector: GeoVector
         :param resolution: output resolution, None for full resolution
+        :param masked: boolean, if `True` the return value will be a masked array. Default is True
+        :param masked_value: only this pixel value will be masked, optional
+        :param bands: list of indices of requested bands, default None which returns all bands
         :return: GeoRaster
         """
         bounds, window = self._vector_to_raster_bounds(vector, boundless=self._image is None)
@@ -944,7 +954,9 @@ def crop(self, vector, resolution=None, masked=True, bands=None):
         else:
             xsize, ysize = (None, None)
 
-        return self.pixel_crop(bounds, xsize, ysize, window=window, masked=masked, bands=bands)
+        return self.pixel_crop(bounds, xsize, ysize,
+                               window=window, masked=masked,
+                               masked_value=masked_value, bands=bands)
 
     def _window(self, bounds, to_round=True):
         # self.window expects to receive the arguments west, south, east, north,
@@ -994,13 +1006,15 @@ def _resolution_to_output_shape(self, bounds, resolution):
         ysize = round(height / yscale)
         return xsize, ysize
 
-    def pixel_crop(self, bounds, xsize=None, ysize=None, window=None, masked=True, bands=None):
+    def pixel_crop(self, bounds, xsize=None, ysize=None, window=None,
+                   masked=True, masked_value=None, bands=None):
         """Crop raster outside vector (convex hull).
 
         :param bounds: bounds of requester portion of the image in image pixels
         :param xsize: output raster width, None for full resolution
         :param ysize: output raster height, None for full resolution
         :param windows: the bounds representation window on image in image pixels, Optional
+        :param masked_value: only this pixel value will be masked, optional
         :param bands: list of indices of requested bands, default None which returns all bands
         :return: GeoRaster
         """
@@ -1015,7 +1029,8 @@ def pixel_crop(self, bounds, xsize=None, ysize=None, window=None, masked=True, b
                                                        bounds[1],
                                                        bounds[2] - bounds[0] + 1,
                                                        bounds[3] - bounds[1] + 1)
-            return self.get_window(window, xsize=xsize, ysize=ysize, masked=masked, bands=bands)
+            return self.get_window(window, xsize=xsize, ysize=ysize,
+                                   masked=masked, masked_value=masked_value, bands=bands)
 
     def _crop(self, bounds, xsize=None, ysize=None):
         """Crop raster outside vector (convex hull).
@@ -1579,8 +1594,9 @@ def _get_window_out_shape(self, bands, window, xsize, ysize):
 
     def get_window(self, window, bands=None,
                    xsize=None, ysize=None,
-                   resampling=Resampling.cubic, masked=True, affine=None
-                   ):
+                   resampling=Resampling.cubic,
+                   masked=True, masked_value=None,
+                   affine=None):
         """Get window from raster.
 
         :param window: requested window
@@ -1589,6 +1605,7 @@ def get_window(self, window, bands=None,
         :param ysize: tile y size default None, for full resolution pass None
         :param resampling: which Resampling to use on reading, default Resampling.cubic
         :param masked: boolean, if `True` the return value will be a masked array. Default is True
+        :masked_value: only this pixel value will be masked, optional
         :return: GeoRaster2 of tile
         """
         bands = bands or list(range(1, self.num_bands + 1))
@@ -1605,6 +1622,8 @@ def get_window(self, window, bands=None,
             }
             with self._raster_opener(self._filename) as raster:  # type: rasterio.io.DatasetReader
                 array = raster.read(bands, **read_params)
+                if masked and (masked_value is not None):
+                    array = np.ma.masked_where(array == masked_value, array)
             affine = affine or self._calculate_new_affine(window, out_shape[2], out_shape[1])
 
             raster = self.copy_with(image=array, affine=affine)

tests/test_merge_all.py

@@ -5,6 +5,7 @@
 import numpy as np
 from numpy.testing import assert_array_equal
 from affine import Affine
+from shapely import geometry
 
 from telluric import FeatureCollection, GeoFeature
 from telluric.constants import WEB_MERCATOR_CRS, WGS84_CRS
@@ -140,6 +141,25 @@ def test_crop_for_merging(main_r, cropping_r):
     assert(rr.affine.almost_equals(cropping_r.affine))
 
 
+def test_crop_masked_value():
+    rr = GeoRaster2.open("tests/data/raster/nomask.jp2")
+    bounds = [-6530057, -3574558, -6492196, -3639376]
+    roi = GeoVector(geometry.Polygon.from_bounds(*bounds), crs=WEB_MERCATOR_CRS)
+
+    # https://github.com/mapbox/rasterio/issues/1449
+    out = rr.crop(roi, masked=True)
+    assert(out.image.mask is np.ma.nomask)
+
+    out = rr.crop(roi, masked=True, masked_value=0)
+    assert(out.image.mask[0, -1, -1])
+
+    out = rr.crop(roi, masked=False)
+    # It is very confused behaviour that crop returns masked array even
+    # if masked=False. This is due to default value of nodata
+    # (nodata=0) in GeoRaster2 constructor.
+    assert(out.image.mask[0, -1, -1])
+
+
 def test_pixel_crop():
     rr = black_and_white_raster([1, 2, 3], height=1000, width=1000)
     out = rr.pixel_crop((0, 0, 1000, 1000))