Skip to content

Commit

Permalink
create_cutline expose rasterio rowcol op parameter to fix pixel raste… (
Browse files Browse the repository at this point in the history
#759)

* 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 <[email protected]>
Co-authored-by: vincentsarago <[email protected]>
  • Loading branch information
3 people authored Oct 29, 2024
1 parent d3a552e commit b833783
Show file tree
Hide file tree
Showing 4 changed files with 78 additions and 6 deletions.
4 changes: 3 additions & 1 deletion CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
24 changes: 19 additions & 5 deletions rio_tiler/utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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.
Expand All @@ -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})"
Expand Down
1 change: 1 addition & 0 deletions tests/test_io_rasterio.py
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
55 changes: 55 additions & 0 deletions tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down

0 comments on commit b833783

Please sign in to comment.