Mask resampling + binary mask option

.gitignore

@@ -104,4 +104,5 @@ ENV/
 .pytest_cache
 
 .benchmarks/
-tests/benchmarks/data/*
+tests/benchmarks/data/*
+tests/fixtures/mask*

rio_tiler/utils.py

@@ -346,41 +346,44 @@ def _tile_read(
     bounds_crs=None,
     minimum_tile_cover=None,
     warp_vrt_option={},
+    force_binary_mask=True,
 ):
     """
     Read data and mask.
 
     Attributes
     ----------
-    src_dst : rasterio.io.DatasetReader
-        rasterio.io.DatasetReader object
-    bounds : list
-        Output bounds (left, bottom, right, top) in target crs ("dst_crs").
-    tilesize : int
-        Output image size
-    indexes : list of ints or a single int, optional, (defaults: None)
-        If `indexes` is a list, the result is a 3D array, but is
-        a 2D array if it is a band index number.
-    nodata: int or float, optional (defaults: None)
-    resampling_method : str, optional (default: "bilinear")
-        Resampling algorithm.
-    tile_edge_padding : int, optional (default: 2)
-        Padding to apply to each edge of the tile when retrieving data
-        to assist in reducing resampling artefacts along edges.
-    dst_crs: CRS or str, optional
-        Target coordinate reference system (default "epsg:3857").
-    bounds_crs: CRS or str, optional
-        Overwrite bounds coordinate reference system (default None, equal to dst_crs).
-    minimum_tile_cover: float, optional (default: None)
-        Minimum % overlap for which to raise an error with dataset not
-        covering enought of the tile.
-    warp_vrt_option: dict, optional (default: {})
-        These will be passed to the rasterio.warp.WarpedVRT class.
+        src_dst : rasterio.io.DatasetReader
+            rasterio.io.DatasetReader object
+        bounds : list
+            Output bounds (left, bottom, right, top) in target crs ("dst_crs").
+        tilesize : int
+            Output image size
+        indexes : list of ints or a single int, optional, (defaults: None)
+            If `indexes` is a list, the result is a 3D array, but is
+            a 2D array if it is a band index number.
+        nodata: int or float, optional (defaults: None)
+        resampling_method : str, optional (default: "bilinear")
+            Resampling algorithm.
+        tile_edge_padding : int, optional (default: 2)
+            Padding to apply to each edge of the tile when retrieving data
+            to assist in reducing resampling artefacts along edges.
+        dst_crs: CRS or str, optional
+            Target coordinate reference system (default "epsg:3857").
+        bounds_crs: CRS or str, optional
+            Overwrite bounds coordinate reference system (default None, equal to dst_crs).
+        minimum_tile_cover: float, optional (default: None)
+            Minimum % overlap for which to raise an error with dataset not
+            covering enought of the tile.
+        warp_vrt_option: dict, optional (default: {})
+            These will be passed to the rasterio.warp.WarpedVRT class.
+        force_binary_mask, bool, optional (default: True)
+            If True, rio-tiler makes sure mask has only 0 or 255 values.
 
     Returns
     -------
-    data : numpy ndarray
-    mask: numpy array
+        data : numpy ndarray
+        mask: numpy array
 
     """
     if isinstance(indexes, int):
@@ -453,7 +456,26 @@ def _tile_read(
             window=out_window,
             resampling=Resampling[resampling_method],
         )
-        mask = vrt.dataset_mask(out_shape=(tilesize, tilesize), window=out_window)
+
+        # This is to overcome some problem with GDALGetMask()
+        if ColorInterp.alpha in vrt.colorinterp:
+            idx = vrt.colorinterp.index(ColorInterp.alpha) + 1
+            mask = vrt.read(
+                indexes=idx,
+                out_shape=(tilesize, tilesize),
+                window=out_window,
+                resampling=Resampling[resampling_method],
+                out_dtype="uint8",
+            )
+        else:
+            mask = vrt.dataset_mask(
+                out_shape=(tilesize, tilesize),
+                window=out_window,
+                resampling=Resampling[resampling_method],
+            )
+
+        if force_binary_mask:
+            mask = np.where(mask != 0, np.uint8(255), np.uint8(0))
 
         return data, mask
 

tests/benchmarks/test_benchmarks.py

@@ -18,7 +18,9 @@ def read_tile(src_path, tile):
     # We make sure to not store things in cache.
     with rasterio.Env(GDAL_CACHEMAX=0, NUM_THREADS="all"):
         with rasterio.open(src_path) as src_dst:
-            return utils._tile_read(src_dst, tile_bounds, 256)
+            return utils._tile_read(
+                src_dst, tile_bounds, 256, resampling_method="nearest"
+            )
 
 
 @pytest.mark.parametrize("tile_name", ["full", "boundless"])

tests/test_mask.py

@@ -0,0 +1,92 @@
+"""Benchmark."""
+
+import os
+import pytest
+
+import numpy
+import mercantile
+
+import rasterio
+from rasterio.crs import CRS
+from rasterio.coords import BoundingBox
+
+from rio_tiler import utils
+
+
+tiles = {
+    "masked": mercantile.Tile(x=535, y=498, z=10),
+    "boundless": mercantile.Tile(x=540, y=497, z=10),
+}
+equator = {
+    "name": "equator",
+    "bounds": BoundingBox(left=382792.5, bottom=362992.5, right=610507.5, top=595207.5),
+    "crs": CRS.from_epsg(32632),
+}
+
+dataset = [
+    dict(equator, dtype="uint8", nodata_type="alpha"),
+    dict(equator, dtype="uint8", nodata_type="nodata"),
+    dict(equator, dtype="uint8", nodata_type="mask"),
+    dict(equator, dtype="int8", nodata_type="alpha"),
+    dict(equator, dtype="int8", nodata_type="nodata"),
+    dict(equator, dtype="int8", nodata_type="mask"),
+    # dict(equator, dtype="uint16", nodata_type="alpha"), #fail
+    dict(equator, dtype="uint16", nodata_type="nodata"),
+    dict(equator, dtype="uint16", nodata_type="mask"),
+    # dict(equator, dtype="int16", nodata_type="alpha"), # Fail
+    dict(equator, dtype="int16", nodata_type="nodata"),
+    # dict(equator, dtype="int16", nodata_type="mask"), # Fail
+]
+
+cog_path = os.path.join(os.path.dirname(__file__), "fixtures", "mask")
+
+
+def test_mask_bilinear(cloudoptimized_geotiff):
+    """Test mask read with bilinear resampling"""
+    src_path = cloudoptimized_geotiff(
+        cog_path, **equator, dtype="uint8", nodata_type="mask"
+    )
+
+    tile = mercantile.Tile(x=535, y=498, z=10)
+    tile_bounds = mercantile.xy_bounds(tile)
+    with rasterio.open(src_path) as src_dst:
+        data, mask = utils._tile_read(
+            src_dst,
+            tile_bounds,
+            256,
+            resampling_method="bilinear",
+            force_binary_mask=True,
+        )
+        masknodata = (data[0] != 0).astype(numpy.uint8) * 255
+        numpy.testing.assert_array_equal(mask, masknodata)
+
+        data, mask = utils._tile_read(
+            src_dst,
+            tile_bounds,
+            256,
+            resampling_method="bilinear",
+            force_binary_mask=False,
+        )
+        masknodata = (data[0] != 0).astype(numpy.uint8) * 255
+        assert not numpy.array_equal(mask, masknodata)
+
+
+@pytest.mark.parametrize("resampling", ["bilinear", "nearest"])
+@pytest.mark.parametrize("tile_name", ["masked"])
+@pytest.mark.parametrize("dataset_info", dataset)
+def test_mask(dataset_info, tile_name, resampling, cloudoptimized_geotiff):
+    """Test tile read for multiple combination of datatype/mask/tile extent."""
+    src_path = cloudoptimized_geotiff(cog_path, **dataset_info)
+
+    tile = tiles[tile_name]
+    tile_bounds = mercantile.xy_bounds(tile)
+    with rasterio.open(src_path) as src_dst:
+        data, mask = utils._tile_read(
+            src_dst,
+            tile_bounds,
+            256,
+            resampling_method=resampling,
+            force_binary_mask=True,
+        )
+        masknodata = (data[0] != 0).astype(numpy.uint8) * 255
+        numpy.testing.assert_array_equal(mask, masknodata)

tests/test_utils.py

@@ -883,20 +883,20 @@ def test_aligned_with_internaltile():
 
 
 # See https://github.com/cogeotiff/rio-tiler/issues/105#issuecomment-492268836
-# def test_tile_read_validMask():
-#     """Dataset mask should be the same as the actual mask."""
-#     address = "{}_B2.TIF".format(LANDSAT_PATH)
-
-#     bounds = (
-#         -8844681.416934313,
-#         3757032.814272982,
-#         -8766409.899970293,
-#         3835304.331237001,
-#     )
-#     tilesize = 128
-#     arr, mask = utils.tile_read(address, bounds, tilesize, nodata=0)
-#     masknodata = (arr[0] != 0).astype(np.uint8) * 255
-#     np.testing.assert_array_equal(mask, masknodata)
+def test_tile_read_validMask():
+    """Dataset mask should be the same as the actual mask."""
+    address = "{}_B2.TIF".format(LANDSAT_PATH)
+
+    bounds = (
+        -8844681.416934313,
+        3757032.814272982,
+        -8766409.899970293,
+        3835304.331237001,
+    )
+    tilesize = 128
+    arr, mask = utils.tile_read(address, bounds, tilesize, nodata=0)
+    masknodata = (arr[0] != 0).astype(np.uint8) * 255
+    np.testing.assert_array_equal(mask, masknodata)
 
 
 def test_tile_read_crs():