Repository: ungarj/mapchete (python)

Issue:
The extract function of a Grid can cause exception
I am using mapchete to read a data in a window from a TileDirectory.

I am using the extract function of the grid class for building a grid that I can use with mapchete.io.raster.read.read_raster_window.

However I have an exception when I call this function due to the fact that  the following line

`return Grid.from_bounds(bounds, Shape(height,width), self.crs)`

Is calling the Shape constructor with float width and height.

I fixed that with a cast.

Second point with this class is an exception raised when the width or height is 0. It might be interesting to have a something like a minimum grid size ?

Maybe I am using an unsupported path to do what I want though.



Pull request:
fix Grid.extract() behavior on small bounds
closes #681 

Patch:
diff --git a/mapchete/grid.py b/mapchete/grid.py
index d1277a7..8deb6d4 100644
--- a/mapchete/grid.py
+++ b/mapchete/grid.py
@@ -1,9 +1,13 @@
 from __future__ import annotations
 
+import math
 from typing import Tuple
 
 from affine import Affine
-from rasterio.transform import array_bounds, from_bounds, rowcol
+from rasterio.transform import array_bounds, from_bounds
+from rasterio.windows import from_bounds as window_from_bounds
+from rasterio.windows import bounds as window_bounds
+from rasterio.windows import Window
 from shapely.geometry import mapping, shape
 from tilematrix import Shape
 
@@ -30,15 +34,24 @@ class Grid:
 
     def extract(self, bounds: BoundsLike) -> Grid:
         bounds = Bounds.from_inp(bounds)
-        # I <3 axis orders!
-        (minrow, maxrow), (mincol, maxcol) = rowcol(
-            self.transform, [bounds.left, bounds.right], [bounds.top, bounds.bottom]
+
+        # use rasterio.window.Window to help with calculation
+        window = window_from_bounds(*bounds, transform=self.transform)
+
+        # now, properly round the window with the intention that every "touched"
+        # row or column is included and that height or width cannot be 0
+        window = Window(
+            col_off=math.floor(window.col_off),  # type: ignore
+            row_off=math.floor(window.row_off),  # type: ignore
+            width=math.ceil(window.width),  # type: ignore
+            height=math.ceil(window.height),  # type: ignore
+        )
+
+        return Grid.from_bounds(
+            window_bounds(window, self.transform),
+            Shape(window.height, window.width),
+            self.crs,
         )
-        width = maxcol - mincol
-        height = maxrow - minrow
-        if width < 0 or height < 0:  # pragma: no cover
-            raise ValueError(f"bounds {bounds} are outside of source grid")
-        return Grid.from_bounds(bounds, Shape(height, width), self.crs)
 
     @staticmethod
     def from_obj(obj):  # pragma: no cover


Test patch:
diff --git a/test/cli/default/__init__.py b/test/cli/default/__init__.py
index 9f81e5f..d70fb33 100644
--- a/test/cli/default/__init__.py
+++ b/test/cli/default/__init__.py
@@ -10,7 +10,7 @@ logger = logging.getLogger(__name__)
 def run_cli(
     args, expected_exit_code=0, output_contains=None, raise_exc=True, cli=mapchete_cli
 ):
-    result = CliRunner(env=dict(MAPCHETE_TEST="TRUE"), mix_stderr=True).invoke(
+    result = CliRunner(env=dict(MAPCHETE_TEST="TRUE")).invoke(
         cli, list(map(str, args)), catch_exceptions=True, standalone_mode=True
     )
     if output_contains:
diff --git a/test/test_grid.py b/test/test_grid.py
new file mode 100644
index 0000000..051d9d9
--- /dev/null
+++ b/test/test_grid.py
@@ -0,0 +1,32 @@
+import pytest
+
+from mapchete.bounds import Bounds
+from mapchete.grid import Grid
+
+
+@pytest.mark.parametrize(
+    "bounds",
+    [
+        (0, 0, 0.4, 0.6),
+        (0, 0, 0.5, 0.5),
+        (0, 0, 0.6, 0.4),
+        (0.4, 0.6, 1, 1),
+        (0.5, 0.5, 1, 1),
+        (0.6, 0.4, 1, 1),
+    ],
+)
+@pytest.mark.parametrize("grid_size", [1, 2, 3])
+def test_extract(bounds, grid_size):
+    # 3x3 grid
+    grid1 = Grid.from_bounds(
+        bounds=(0, 0, grid_size, grid_size),
+        shape=(grid_size, grid_size),
+        crs="EPSG:4326",
+    )
+
+    # extract from bounds smaller than the pixel size
+    grid2 = grid1.extract(bounds)
+
+    # make sure grid2 is snapped to grid1
+    assert grid2.shape == (1, 1)
+    assert grid2.bounds == Bounds.from_inp((0, 0, 1, 1))