From b83378301706718de3fd21ee8ba4aaf04e61f32a Mon Sep 17 00:00:00 2001 From: Martino Boni Date: Tue, 29 Oct 2024 21:24:48 +0100 Subject: [PATCH] =?UTF-8?q?create=5Fcutline=20expose=20rasterio=20rowcol?= =?UTF-8?q?=20op=20parameter=20to=20fix=20pixel=20raste=E2=80=A6=20(#759)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * create_cutline expose rasterio rowcol op parameter to fix pixel raster x,y approximation * fix pydantic operator input type (got error with tox validation) * test update added op param to exists tests with default math.floor value. TODO: next need to test with specific feat and round or ceiling operators * UT create_cutline check callable operator passed rasies exception if not one of the accepted methods by rasterio.transform.rowcol * edit PR * update type * update changelog --------- Co-authored-by: Martino Co-authored-by: vincentsarago --- CHANGES.md | 4 ++- rio_tiler/utils.py | 24 +++++++++++++---- tests/test_io_rasterio.py | 1 + tests/test_utils.py | 55 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 78 insertions(+), 6 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index b9270087..5d2eb535 100644 --- a/CHANGES.md +++ b/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) diff --git a/rio_tiler/utils.py b/rio_tiler/utils.py index 098f0756..2c6b0abf 100644 --- a/rio_tiler/utils.py +++ b/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})" diff --git a/tests/test_io_rasterio.py b/tests/test_io_rasterio.py index 7bcba03d..7a8427f6 100644 --- a/tests/test_io_rasterio.py +++ b/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") diff --git a/tests/test_utils.py b/tests/test_utils.py index 37c356e1..abe65533 100644 --- a/tests/test_utils.py +++ b/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]