create_cutline expose rasterio rowcol op parameter to fix pixel raste…

CHANGES.md

@@ -12,7 +12,9 @@
 
 * Cast Xarray `attrs` values in XarrayReader's `info()` response to avoid JSON encoding issues (https://github.com/cogeotiff/rio-tiler/pull/755)
 
-* refactor XarrayReader's `feature()` method to use the `part` method (https://github.com/cogeotiff/rio-tiler/pull/755)
+* Refactor XarrayReader's `feature()` method to use the `part` method (https://github.com/cogeotiff/rio-tiler/pull/755)
+
+* Allow `op` parameter for `create_cutline` and `_convert_to_raster_space` functions to better control rasterio's `rowcol` behaviour (author @Martenz, https://github.com/cogeotiff/rio-tiler/pull/759)
 
 # 7.0.1 (2024-10-22)
 

rio_tiler/utils.py

@@ -3,7 +3,17 @@
 import math
 import warnings
 from io import BytesIO
-from typing import Any, Dict, Generator, List, Optional, Sequence, Tuple, Union
+from typing import (
+    Any,
+    Callable,
+    Dict,
+    Generator,
+    List,
+    Optional,
+    Sequence,
+    Tuple,
+    Union,
+)
 
 import numpy
 import rasterio
@@ -642,12 +652,15 @@ def _calculateRatio(
 def _convert_to_raster_space(
     src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
     poly_coordinates: List,
+    op: Optional[Callable[[float], Any]] = None,
 ) -> List[str]:
+    # NOTE: we could remove this once we have rasterio >= 1.4.2
+    op = op or numpy.floor
     polygons = []
     for point in poly_coordinates:
         xs, ys = zip(*coords(point))
-        src_y, src_x = rowcol(src_dst.transform, xs, ys)
-        polygon = ", ".join([f"{x} {y}" for x, y in list(zip(src_x, src_y))])
+        src_y, src_x = rowcol(src_dst.transform, xs, ys, op=op)
+        polygon = ", ".join([f"{int(x)} {int(y)}" for x, y in list(zip(src_x, src_y))])
         polygons.append(f"({polygon})")
 
     return polygons
@@ -657,6 +670,7 @@ def create_cutline(
     src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT],
     geometry: Dict,
     geometry_crs: CRS = None,
+    op: Optional[Callable[[float], Any]] = None,
 ) -> str:
     """
     Create WKT Polygon Cutline for GDALWarpOptions.
@@ -678,13 +692,13 @@ def create_cutline(
         geometry = transform_geom(geometry_crs, src_dst.crs, geometry)
 
     if geom_type == "Polygon":
-        polys = ",".join(_convert_to_raster_space(src_dst, geometry["coordinates"]))
+        polys = ",".join(_convert_to_raster_space(src_dst, geometry["coordinates"], op))
         wkt = f"POLYGON ({polys})"
 
     elif geom_type == "MultiPolygon":
         multi_polys = []
         for poly in geometry["coordinates"]:
-            polys = ",".join(_convert_to_raster_space(src_dst, poly))
+            polys = ",".join(_convert_to_raster_space(src_dst, poly, op))
             multi_polys.append(f"({polys})")
         str_multipoly = ",".join(multi_polys)
         wkt = f"MULTIPOLYGON ({str_multipoly})"

tests/test_io_rasterio.py

@@ -762,6 +762,7 @@ def test_equality_part_feature():
         # I would assume this is due to rounding issue or reprojection of the cutline by GDAL
         # After some debugging locally I found out the rasterized mask is more precise
         # numpy.testing.assert_array_equal(img_part.mask, img_feat.mask)
+        # NOTE reply: can this rounding be fixed with a different operator passed to rasterio rowcol?
 
         # Re-Projection
         img_feat = src.feature(feat, dst_crs="epsg:3857")

tests/test_utils.py

@@ -389,6 +389,61 @@ def test_cutline():
         assert sum(mask[-1, :]) == 0  # last line
 
 
+def test_cutline_operator(dataset_fixture):
+    """Test rio_tiler.utils.create_cutline with operators."""
+    with MemoryFile(
+        dataset_fixture(
+            crs=CRS.from_epsg(4326),
+            bounds=(-175.0, -85, 175.0, 85.0),
+            dtype="uint8",
+            nodata_type="nodata",
+            width=720,
+            height=360,
+        )
+    ) as memfile:
+        with memfile.open() as src_dst:
+            feat = {
+                "type": "Polygon",
+                "coordinates": [
+                    [
+                        [-163.0, -83.0],
+                        [163.0, -83.0],
+                        [163.0, 83.0],
+                        [-163.0, 83.0],
+                        [-163.0, -83.0],
+                    ]
+                ],
+            }
+            cutline = utils.create_cutline(
+                src_dst,
+                feat,
+                geometry_crs="epsg:4326",
+            )
+            cutline_mathfloor = utils.create_cutline(
+                src_dst,
+                feat,
+                geometry_crs="epsg:4326",
+                op=math.floor,
+            )
+            assert cutline == cutline_mathfloor
+
+            cutline_npfloor = utils.create_cutline(
+                src_dst,
+                feat,
+                geometry_crs="epsg:4326",
+                op=np.floor,
+            )
+            assert cutline_npfloor == cutline_mathfloor
+
+            cutline_npceil = utils.create_cutline(
+                src_dst,
+                feat,
+                geometry_crs="epsg:4326",
+                op=np.ceil,
+            )
+            assert cutline_npceil != cutline_npfloor
+
+
 def test_parse_expression():
     """test parsing rio-tiler expression."""
     assert sorted(parse_expression("b1*b11+b3,b1*b2+B4")) == [1, 2, 3, 4, 11]