Provide parameter for retrieving additional data around tile edges to…

… minimise edge effects (#104)

* initial commit of padding a tile

* add a test

* remove used print statements

* Add tests. Modify parameter name.

* Test if requested tile is aligned to internal tiles, if it is don't pad. Clarify param in docstring. Provide better default window.

* Add test for web-optimised COG

* Make black happy

* Add info to changes.txt

CHANGES.txt

@@ -1,3 +1,8 @@
+Next release 
+------------------
+
+- add tile_edge_padding option to be passed to rio_tiler.utils._tile_read to reduce sharp edges that occur due to resampling (#104)
+
 1.2.4 (2019-04-16)
 ------------------
 

rio_tiler/utils.py

@@ -17,11 +17,13 @@
 from rasterio.enums import Resampling, MaskFlags, ColorInterp
 from rasterio.io import DatasetReader, MemoryFile
 from rasterio import transform
+from rasterio import windows
 from rasterio.warp import calculate_default_transform, transform_bounds
 
 from rio_tiler.mercator import get_zooms
 
 from rio_tiler.errors import NoOverviewWarning
+from affine import Affine
 
 logger = logging.getLogger(__name__)
 
@@ -262,7 +264,13 @@ def has_alpha_band(src_dst):
 
 
 def _tile_read(
-    src_dst, bounds, tilesize, indexes=None, nodata=None, resampling_method="bilinear"
+    src_dst,
+    bounds,
+    tilesize,
+    indexes=None,
+    nodata=None,
+    resampling_method="bilinear",
+    tile_edge_padding=2,
 ):
     """
     Read data and mask.
@@ -281,6 +289,9 @@ def _tile_read(
     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.
 
     Returns
     -------
@@ -298,6 +309,27 @@ def _tile_read(
     )
 
     vrt_transform, vrt_width, vrt_height = get_vrt_transform(src_dst, bounds)
+    out_window = windows.Window(
+        col_off=0, row_off=0, width=vrt_width, height=vrt_height
+    )
+
+    if tile_edge_padding > 0 and not _requested_tile_aligned_with_internal_tile(
+        src_dst, bounds, tilesize
+    ):
+        vrt_transform = vrt_transform * Affine.translation(
+            -tile_edge_padding, -tile_edge_padding
+        )
+        orig__vrt_height = vrt_height
+        orig_vrt_width = vrt_width
+        vrt_height = vrt_height + 2 * tile_edge_padding
+        vrt_width = vrt_width + 2 * tile_edge_padding
+        out_window = windows.Window(
+            col_off=tile_edge_padding,
+            row_off=tile_edge_padding,
+            width=orig_vrt_width,
+            height=orig__vrt_height,
+        )
+
     vrt_params.update(dict(transform=vrt_transform, width=vrt_width, height=vrt_height))
 
     indexes = indexes if indexes is not None else src_dst.indexes
@@ -314,9 +346,11 @@ def _tile_read(
         data = vrt.read(
             out_shape=out_shape,
             indexes=indexes,
+            window=out_window,
             resampling=Resampling[resampling_method],
         )
-        mask = vrt.dataset_mask(out_shape=(tilesize, tilesize))
+
+        mask = vrt.dataset_mask(out_shape=(tilesize, tilesize), window=out_window)
 
         return data, mask
 
@@ -407,6 +441,27 @@ def tile_exists(bounds, tile_z, tile_x, tile_y):
     )
 
 
+def _requested_tile_aligned_with_internal_tile(src_dst, bounds, tilesize):
+    """Check if tile is aligned with internal tiles."""
+    if src_dst.crs != CRS.from_epsg(3857):
+        return False
+
+    col_off, row_off, w, h = windows.from_bounds(
+        *bounds, height=tilesize, transform=src_dst.transform, width=tilesize
+    ).flatten()
+
+    if round(w) % 64 and round(h) % 64:
+        return False
+
+    if (src_dst.width - round(col_off)) % 64:
+        return False
+
+    if (src_dst.height - round(row_off)) % 64:
+        return False
+
+    return True
+
+
 def _apply_discrete_colormap(arr, cmap):
     """
     Apply discrete colormap.

tests/test_utils.py

@@ -6,6 +6,7 @@
 from mock import patch
 
 import numpy as np
+import mercantile
 
 from rio_toa import toa_utils
 
@@ -39,6 +40,7 @@
 PIX4D_PATH = os.path.join(S3_LOCAL, KEY_PIX4D)
 
 COG_DST = os.path.join(os.path.dirname(__file__), "fixtures", "cog_name.tif")
+COG_WEB_TILED = os.path.join(os.path.dirname(__file__), "fixtures", "web.tif")
 
 
 with open("{}_MTL.txt".format(LANDSAT_PATH), "r") as f:
@@ -80,6 +82,65 @@ def test_tile_read_validResampling():
     assert mask.shape == (16, 16)
 
 
+def test_resampling_returns_different_results():
+    address = "{}_B2.TIF".format(LANDSAT_PATH)
+    bounds = (
+        -8844681.416934313,
+        3757032.814272982,
+        -8766409.899970293,
+        3835304.331237001,
+    )
+    tilesize = 16
+
+    arr, mask = utils.tile_read(address, bounds, tilesize)
+    arr2, mask2 = utils.tile_read(
+        address, bounds, tilesize, resampling_method="nearest"
+    )
+    assert not np.array_equal(arr, arr2)
+
+
+def test_resampling_with_diff_padding_returns_different_results():
+    address = "{}_B2.TIF".format(LANDSAT_PATH)
+    bounds = (
+        -8844681.416934313,
+        3757032.814272982,
+        -8766409.899970293,
+        3835304.331237001,
+    )
+    tilesize = 16
+
+    arr, mask = utils.tile_read(address, bounds, tilesize)
+    arr2, mask2 = utils.tile_read(address, bounds, tilesize, tile_edge_padding=0)
+    assert not np.array_equal(arr, arr2)
+
+
+def test_tile_padding_only_effects_edge_pixels():
+    """Adding tile padding should effect edge pixels only."""
+    address = "{}_B2.TIF".format(LANDSAT_PATH)
+    bounds = (
+        -8844681.416934313,
+        3757032.814272982,
+        -8766409.899970293,
+        3835304.331237001,
+    )
+    tilesize = 16
+
+    arr, mask = utils.tile_read(address, bounds, tilesize)
+    arr2, mask2 = utils.tile_read(address, bounds, tilesize, tile_edge_padding=0)
+    assert not np.array_equal(arr[0][0], arr2[0][0])
+    assert np.array_equal(arr[0][5:-5][5:-5], arr2[0][5:-5][5:-5])
+
+
+def test_that_tiling_ignores_padding_if_web_friendly_internal_tiles_exist():
+    address = COG_WEB_TILED
+    bounds = mercantile.bounds(147, 182, 9)
+    tilesize = 256
+
+    arr, mask = utils.tile_read(address, bounds, tilesize)
+    arr2, mask2 = utils.tile_read(address, bounds, tilesize, tile_edge_padding=0)
+    assert np.array_equal(arr, arr2)
+
+
 def test_tile_read_invalidResampling():
     """Should raise an error on invalid resampling method name."""
     address = "{}_B2.TIF".format(LANDSAT_PATH)