From c1c67150f88ba531c88534f71236819e168683da Mon Sep 17 00:00:00 2001 From: vschaffn Date: Thu, 13 Feb 2025 17:31:13 +0100 Subject: [PATCH 1/4] feat: add tiling function for raster --- geoutils/raster/georeferencing.py | 47 +++++++++++++++++++ geoutils/raster/raster.py | 75 ++++++++++++++++++++++++++++++- 2 files changed, 120 insertions(+), 2 deletions(-) diff --git a/geoutils/raster/georeferencing.py b/geoutils/raster/georeferencing.py index 604a718a..91859065 100644 --- a/geoutils/raster/georeferencing.py +++ b/geoutils/raster/georeferencing.py @@ -22,6 +22,7 @@ from __future__ import annotations +import math import warnings from typing import Iterable, Literal @@ -284,3 +285,49 @@ def _cast_nodata(out_dtype: DTypeLike, nodata: int | float | None) -> int | floa nodata = nodata return nodata + + +def _generate_tiling_grid( + row_min: int, + col_min: int, + row_max: int, + col_max: int, + row_split: int, + col_split: int, + overlap: int = 0, +) -> NDArrayNum: + """ + Generate a grid of positions by splitting [row_min, row_max] x + [col_min, col_max] into tiles of size row_split x col_split with optional overlap. + + :param row_min: Minimum row index of the bounding box to split. + :param col_min: Minimum column index of the bounding box to split. + :param row_max: Maximum row index of the bounding box to split. + :param col_max: Maximum column index of the bounding box to split. + :param row_split: Height of each tile. + :param col_split: Width of each tile. + :param overlap: size of overlapping between tiles (both vertically and horizontally). + :return: A numpy array grid with splits in two dimensions (0: row, 1: column), + where each cell contains [row_min, row_max, col_min, col_max]. + """ + # Calculate the number of splits considering overlap + nb_col_split = math.ceil((col_max - col_min) / (col_split - overlap)) + nb_row_split = math.ceil((row_max - row_min) / (row_split - overlap)) + + # Initialize the output grid + tiling_grid = np.zeros(shape=(nb_row_split, nb_col_split, 4), dtype=int) + + for row in range(nb_row_split): + for col in range(nb_col_split): + # Calculate the start of the tile + row_start = row_min + row * (row_split - overlap) + col_start = col_min + col * (col_split - overlap) + + # Calculate the end of the tile ensuring it doesn't exceed the bounds + row_end = min(row_max, row_start + row_split) + col_end = min(col_max, col_start + col_split) + + # Populate the grid with the tile boundaries + tiling_grid[row, col] = [row_start, row_end, col_start, col_end] + + return tiling_grid diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py index 78570f0a..65851268 100644 --- a/geoutils/raster/raster.py +++ b/geoutils/raster/raster.py @@ -41,6 +41,7 @@ import rioxarray import xarray as xr from affine import Affine +from matplotlib.patches import Rectangle from mpl_toolkits.axes_grid1 import make_axes_locatable from packaging.version import Version from rasterio.crs import CRS @@ -77,6 +78,7 @@ _cast_pixel_interpretation, _coords, _default_nodata, + _generate_tiling_grid, _ij2xy, _outside_image, _res, @@ -2655,6 +2657,44 @@ def translate( raster_copy.transform = translated_transform return raster_copy + def compute_tiling( + self, + tile_size: int, + raster_ref: RasterType, + overlap: int = 0, + ) -> NDArrayNum: + """ + Compute the raster tiling grid to coregister raster by block. + + :param tile_size: Size of each tile (square tiles) + :param raster_ref: The other raster to coregister, use to validate the shape + :param overlap: Size of overlap between tiles (optional) + :return: tiling_grid (array of tile boundaries), new_shape (shape of the tiled grid) + """ + if self.shape != raster_ref.shape: + raise Exception("Reference and secondary rasters do not have the same shape") + row_max, col_max = self.shape + + # Generate tiling + tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, tile_size, tile_size, overlap=overlap) + return tiling_grid + + def plot_tiling(self, tiling_grid: NDArrayNum) -> None: + """ + Plot raster with its tiling. + + :param tiling_grid: tiling given by Raster.compute_tiling. + """ + ax, caxes = self.plot(return_axes=True) + for tile in tiling_grid.reshape(-1, 4): + row_min, row_max, col_min, col_max = tile + x_min, y_min = self.transform * (col_min, row_min) # Bottom-left corner + x_max, y_max = self.transform * (col_max, row_max) # Top-right corne + rect = Rectangle( + (x_min, y_min), x_max - x_min, y_max - y_min, edgecolor="red", facecolor="none", linewidth=1.5 + ) + ax.add_patch(rect) + def save( self, filename: str | pathlib.Path | IO[bytes], @@ -2922,6 +2962,38 @@ def intersection(self, raster: str | Raster, match_ref: bool = True) -> tuple[fl # mypy raises a type issue, not sure how to address the fact that output of merge_bounds can be () return intersection # type: ignore + @overload + def plot( + self, + bands: int | tuple[int, ...] | None = None, + cmap: matplotlib.colors.Colormap | str | None = None, + vmin: float | int | None = None, + vmax: float | int | None = None, + alpha: float | int | None = None, + cbar_title: str | None = None, + add_cbar: bool = True, + ax: matplotlib.axes.Axes | Literal["new"] | None = None, + *, + return_axes: Literal[False] = False, + **kwargs: Any, + ) -> None: ... + + @overload + def plot( + self, + bands: int | tuple[int, ...] | None = None, + cmap: matplotlib.colors.Colormap | str | None = None, + vmin: float | int | None = None, + vmax: float | int | None = None, + alpha: float | int | None = None, + cbar_title: str | None = None, + add_cbar: bool = True, + ax: matplotlib.axes.Axes | Literal["new"] | None = None, + *, + return_axes: Literal[True], + **kwargs: Any, + ) -> tuple[matplotlib.axes.Axes, matplotlib.colors.Colormap]: ... + def plot( self, bands: int | tuple[int, ...] | None = None, @@ -3059,8 +3131,7 @@ def plot( # If returning axes if return_axes: return ax0, cax - else: - return None + return None def reduce_points( self, From b787b68f051255940dbd334a1c1b24f74f95fdc6 Mon Sep 17 00:00:00 2001 From: vschaffn Date: Thu, 13 Feb 2025 17:31:50 +0100 Subject: [PATCH 2/4] test: add test for the tiling function --- tests/test_raster/test_georeferencing.py | 59 +++++++++++++++++++++++- 1 file changed, 58 insertions(+), 1 deletion(-) diff --git a/tests/test_raster/test_georeferencing.py b/tests/test_raster/test_georeferencing.py index fe2fcba9..1a16a936 100644 --- a/tests/test_raster/test_georeferencing.py +++ b/tests/test_raster/test_georeferencing.py @@ -5,6 +5,7 @@ import geoutils as gu from geoutils import examples +from geoutils.raster.georeferencing import _generate_tiling_grid class TestGeoreferencing: @@ -173,7 +174,7 @@ def test_xy2ij(self) -> None: @pytest.mark.parametrize("example", [landsat_b4_path, aster_dem_path, landsat_rgb_path]) # type: ignore def test_coords(self, example: str) -> None: - img = gu.Raster(self.landsat_b4_path) + img = gu.Raster(example) # With lower left argument xx0, yy0 = img.coords(grid=False, force_offset="ll") @@ -209,3 +210,59 @@ def test_coords(self, example: str) -> None: xxgrid, yygrid = img.coords(grid=True, force_offset="ll") assert np.array_equal(xxgrid, np.repeat(xx0[np.newaxis, :], img.height, axis=0)) assert np.array_equal(yygrid, np.flipud(np.repeat(yy0[:, np.newaxis], img.width, axis=1))) + + @pytest.mark.parametrize("overlap", [0, 5]) # type: ignore + def test_tiling(self, overlap: int) -> None: + + # Test with mock data + tiling_grid_mock = _generate_tiling_grid(0, 0, 100, 100, 50, 50, overlap) + if overlap == 0: + expected_tiling = np.array([[[0, 50, 0, 50], [0, 50, 50, 100]], [[50, 100, 0, 50], [50, 100, 50, 100]]]) + assert np.array_equal(tiling_grid_mock, expected_tiling) + elif overlap == 5: + expected_tiling = np.array( + [ + [[0, 50, 0, 50], [0, 50, 45, 95], [0, 50, 90, 100]], + [[45, 95, 0, 50], [45, 95, 45, 95], [45, 95, 90, 100]], + [[90, 100, 0, 50], [90, 100, 45, 95], [90, 100, 90, 100]], + ] + ) + assert np.array_equal(tiling_grid_mock, expected_tiling) + + # Test with real data + img = gu.Raster(self.landsat_b4_path) + + # Define tiling parameters + row_split, col_split = 100, 100 + row_max, col_max = img.shape + + # Generate the tiling grid + tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, row_split, col_split, overlap) + + # Calculate expected number of tiles + nb_row_tiles = np.ceil(row_max / (row_split - overlap)).astype(int) + nb_col_tiles = np.ceil(col_max / (col_split - overlap)).astype(int) + + # Check that the tiling grid has the expected shape + assert tiling_grid.shape == (nb_row_tiles, nb_col_tiles, 4) + + # Check the boundaries of the first and last tile + assert np.array_equal(tiling_grid[0, 0], np.array([0, min(row_split, row_max), 0, min(col_split, col_max)])) + assert np.array_equal( + tiling_grid[-1, -1], + np.array( + [ + (nb_row_tiles - 1) * (row_split - overlap), + row_max, + (nb_col_tiles - 1) * (col_split - overlap), + col_max, + ] + ), + ) + + # Check if overlap is consistent between tiles + for row in range(nb_row_tiles - 1): + assert tiling_grid[row + 1, 0, 0] == tiling_grid[row, 0, 1] - overlap + + for col in range(nb_col_tiles - 1): # Skip last tile in column + assert tiling_grid[0, col + 1, 2] == tiling_grid[0, col, 3] - overlap From 7acb4e31e107810307b584b3ac454b75f29a9b4c Mon Sep 17 00:00:00 2001 From: vschaffn Date: Fri, 14 Feb 2025 16:55:30 +0100 Subject: [PATCH 3/4] fix: move and in , modify overlap handling --- geoutils/raster/__init__.py | 1 + geoutils/raster/georeferencing.py | 47 ------- geoutils/raster/raster.py | 2 +- geoutils/raster/sampling.py | 96 ------------- geoutils/raster/tiling.py | 168 +++++++++++++++++++++++ tests/test_raster/test_georeferencing.py | 57 -------- tests/test_raster/test_tiling.py | 94 +++++++++++++ 7 files changed, 264 insertions(+), 201 deletions(-) create mode 100644 geoutils/raster/tiling.py create mode 100644 tests/test_raster/test_tiling.py diff --git a/geoutils/raster/__init__.py b/geoutils/raster/__init__.py index 1ecfa46d..200bab9d 100644 --- a/geoutils/raster/__init__.py +++ b/geoutils/raster/__init__.py @@ -23,5 +23,6 @@ from geoutils.raster.multiraster import * # noqa from geoutils.raster.sampling import * # noqa from geoutils.raster.satimg import * # noqa +from geoutils.raster.tiling import * # noqa __all__ = ["RasterType", "Raster"] diff --git a/geoutils/raster/georeferencing.py b/geoutils/raster/georeferencing.py index 91859065..604a718a 100644 --- a/geoutils/raster/georeferencing.py +++ b/geoutils/raster/georeferencing.py @@ -22,7 +22,6 @@ from __future__ import annotations -import math import warnings from typing import Iterable, Literal @@ -285,49 +284,3 @@ def _cast_nodata(out_dtype: DTypeLike, nodata: int | float | None) -> int | floa nodata = nodata return nodata - - -def _generate_tiling_grid( - row_min: int, - col_min: int, - row_max: int, - col_max: int, - row_split: int, - col_split: int, - overlap: int = 0, -) -> NDArrayNum: - """ - Generate a grid of positions by splitting [row_min, row_max] x - [col_min, col_max] into tiles of size row_split x col_split with optional overlap. - - :param row_min: Minimum row index of the bounding box to split. - :param col_min: Minimum column index of the bounding box to split. - :param row_max: Maximum row index of the bounding box to split. - :param col_max: Maximum column index of the bounding box to split. - :param row_split: Height of each tile. - :param col_split: Width of each tile. - :param overlap: size of overlapping between tiles (both vertically and horizontally). - :return: A numpy array grid with splits in two dimensions (0: row, 1: column), - where each cell contains [row_min, row_max, col_min, col_max]. - """ - # Calculate the number of splits considering overlap - nb_col_split = math.ceil((col_max - col_min) / (col_split - overlap)) - nb_row_split = math.ceil((row_max - row_min) / (row_split - overlap)) - - # Initialize the output grid - tiling_grid = np.zeros(shape=(nb_row_split, nb_col_split, 4), dtype=int) - - for row in range(nb_row_split): - for col in range(nb_col_split): - # Calculate the start of the tile - row_start = row_min + row * (row_split - overlap) - col_start = col_min + col * (col_split - overlap) - - # Calculate the end of the tile ensuring it doesn't exceed the bounds - row_end = min(row_max, row_start + row_split) - col_end = min(col_max, col_start + col_split) - - # Populate the grid with the tile boundaries - tiling_grid[row, col] = [row_start, row_end, col_start, col_end] - - return tiling_grid diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py index 65851268..44e7754f 100644 --- a/geoutils/raster/raster.py +++ b/geoutils/raster/raster.py @@ -78,7 +78,6 @@ _cast_pixel_interpretation, _coords, _default_nodata, - _generate_tiling_grid, _ij2xy, _outside_image, _res, @@ -90,6 +89,7 @@ decode_sensor_metadata, parse_and_convert_metadata_from_filename, ) +from geoutils.raster.tiling import _generate_tiling_grid from geoutils.stats import nmad from geoutils.vector.vector import Vector diff --git a/geoutils/raster/sampling.py b/geoutils/raster/sampling.py index ae8085ae..c95584a5 100644 --- a/geoutils/raster/sampling.py +++ b/geoutils/raster/sampling.py @@ -109,99 +109,3 @@ def subsample_array( return unraveled_indices else: return array[unraveled_indices] - - -def _get_closest_rectangle(size: int) -> tuple[int, int]: - """ - Given a 1D array size, return a rectangular shape that is closest to a cube which the size fits in. - - If 'size' does not have an integer root, a rectangle is returned that is slightly larger than 'size'. - - :examples: - >>> _get_closest_rectangle(4) # size will be 4 - (2, 2) - >>> _get_closest_rectangle(9) # size will be 9 - (3, 3) - >>> _get_closest_rectangle(3) # size will be 4; needs padding afterward. - (2, 2) - >>> _get_closest_rectangle(55) # size will be 56; needs padding afterward. - (7, 8) - >>> _get_closest_rectangle(24) # size will be 25; needs padding afterward - (5, 5) - >>> _get_closest_rectangle(85620) # size will be 85849; needs padding afterward - (293, 293) - >>> _get_closest_rectangle(52011) # size will be 52212; needs padding afterward - (228, 229) - """ - close_cube = int(np.sqrt(size)) - - # If size has an integer root, return the respective cube. - if close_cube**2 == size: - return (close_cube, close_cube) - - # One of these rectangles/cubes will cover all cells, so return the first that does. - potential_rectangles = [(close_cube, close_cube + 1), (close_cube + 1, close_cube + 1)] - - for rectangle in potential_rectangles: - if np.prod(rectangle) >= size: - return rectangle - - raise NotImplementedError(f"Function criteria not met for rectangle of size: {size}") - - -def subdivide_array(shape: tuple[int, ...], count: int) -> NDArrayNum: - """ - Create indices for subdivison of an array in a number of blocks. - - If 'count' is divisible by the product of 'shape', the amount of cells in each block will be equal. - If 'count' is not divisible, the amount of cells in each block will be very close to equal. - - :param shape: The shape of a array to be subdivided. - :param count: The amount of subdivisions to make. - - :examples: - >>> subdivide_array((4, 4), 4) - array([[0, 0, 1, 1], - [0, 0, 1, 1], - [2, 2, 3, 3], - [2, 2, 3, 3]]) - - >>> subdivide_array((6, 4), 4) - array([[0, 0, 1, 1], - [0, 0, 1, 1], - [0, 0, 1, 1], - [2, 2, 3, 3], - [2, 2, 3, 3], - [2, 2, 3, 3]]) - - >>> subdivide_array((5, 4), 3) - array([[0, 0, 0, 0], - [0, 0, 0, 0], - [1, 1, 2, 2], - [1, 1, 2, 2], - [1, 1, 2, 2]]) - - :raises ValueError: If the 'shape' size (`np.prod(shape)`) is smallern than 'count' - If the shape is not a 2D shape. - - :returns: An array of shape 'shape' with 'count' unique indices. - """ - try: - import skimage.transform - except ImportError: - raise ImportError("Missing optional dependency, skimage.transform, required by this function.") - - if count > np.prod(shape): - raise ValueError(f"Shape '{shape}' size ({np.prod(shape)}) is smaller than 'count' ({count}).") - - if len(shape) != 2: - raise ValueError(f"Expected a 2D shape, got {len(shape)}D shape: {shape}") - - # Generate a small grid of indices, with the same unique count as 'count' - rect = _get_closest_rectangle(count) - small_indices = np.pad(np.arange(count), np.prod(rect) - count, mode="edge")[: int(np.prod(rect))].reshape(rect) - - # Upscale the grid to fit the output shape using nearest neighbour scaling. - indices = skimage.transform.resize(small_indices, shape, order=0, preserve_range=True).astype(int) - - return indices.reshape(shape) diff --git a/geoutils/raster/tiling.py b/geoutils/raster/tiling.py new file mode 100644 index 00000000..569f807c --- /dev/null +++ b/geoutils/raster/tiling.py @@ -0,0 +1,168 @@ +# Copyright (c) 2025 GeoUtils developers +# +# This file is part of the GeoUtils project: +# https://github.com/glaciohack/geoutils +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +"""Tiling tools for arrays and rasters.""" + +import math + +import numpy as np + +from geoutils._typing import NDArrayNum + + +def _get_closest_rectangle(size: int) -> tuple[int, int]: + """ + Given a 1D array size, return a rectangular shape that is closest to a cube which the size fits in. + + If 'size' does not have an integer root, a rectangle is returned that is slightly larger than 'size'. + + :examples: + >>> _get_closest_rectangle(4) # size will be 4 + (2, 2) + >>> _get_closest_rectangle(9) # size will be 9 + (3, 3) + >>> _get_closest_rectangle(3) # size will be 4; needs padding afterward. + (2, 2) + >>> _get_closest_rectangle(55) # size will be 56; needs padding afterward. + (7, 8) + >>> _get_closest_rectangle(24) # size will be 25; needs padding afterward + (5, 5) + >>> _get_closest_rectangle(85620) # size will be 85849; needs padding afterward + (293, 293) + >>> _get_closest_rectangle(52011) # size will be 52212; needs padding afterward + (228, 229) + """ + close_cube = int(np.sqrt(size)) + + # If size has an integer root, return the respective cube. + if close_cube**2 == size: + return close_cube, close_cube + + # One of these rectangles/cubes will cover all cells, so return the first that does. + potential_rectangles = [(close_cube, close_cube + 1), (close_cube + 1, close_cube + 1)] + + for rectangle in potential_rectangles: + if np.prod(rectangle) >= size: + return rectangle + + # Default return, should never reach here as one of the potential_rectangles will cover all cells. + return potential_rectangles[-1] + + +def subdivide_array(shape: tuple[int, ...], count: int) -> NDArrayNum: + """ + Create indices for subdivison of an array in a number of blocks. + + If 'count' is divisible by the product of 'shape', the amount of cells in each block will be equal. + If 'count' is not divisible, the amount of cells in each block will be very close to equal. + + :param shape: The shape of a array to be subdivided. + :param count: The amount of subdivisions to make. + + :examples: + >>> subdivide_array((4, 4), 4) + array([[0, 0, 1, 1], + [0, 0, 1, 1], + [2, 2, 3, 3], + [2, 2, 3, 3]]) + + >>> subdivide_array((6, 4), 4) + array([[0, 0, 1, 1], + [0, 0, 1, 1], + [0, 0, 1, 1], + [2, 2, 3, 3], + [2, 2, 3, 3], + [2, 2, 3, 3]]) + + >>> subdivide_array((5, 4), 3) + array([[0, 0, 0, 0], + [0, 0, 0, 0], + [1, 1, 2, 2], + [1, 1, 2, 2], + [1, 1, 2, 2]]) + + :raises ValueError: If the 'shape' size (`np.prod(shape)`) is smallern than 'count' + If the shape is not a 2D shape. + + :returns: An array of shape 'shape' with 'count' unique indices. + """ + try: + import skimage.transform + except ImportError: + raise ImportError("Missing optional dependency, skimage.transform, required by this function.") + + if count > np.prod(shape): + raise ValueError(f"Shape '{shape}' size ({np.prod(shape)}) is smaller than 'count' ({count}).") + + if len(shape) != 2: + raise ValueError(f"Expected a 2D shape, got {len(shape)}D shape: {shape}") + + # Generate a small grid of indices, with the same unique count as 'count' + rect = _get_closest_rectangle(count) + small_indices = np.pad(np.arange(count), np.prod(rect) - count, mode="edge")[: int(np.prod(rect))].reshape(rect) + + # Upscale the grid to fit the output shape using nearest neighbour scaling. + indices = skimage.transform.resize(small_indices, shape, order=0, preserve_range=True).astype(int) + + return indices.reshape(shape) + + +def _generate_tiling_grid( + row_min: int, + col_min: int, + row_max: int, + col_max: int, + row_split: int, + col_split: int, + overlap: int = 0, +) -> NDArrayNum: + """ + Generate a grid of positions by splitting [row_min, row_max] x + [col_min, col_max] into tiles of size row_split x col_split with optional overlap. + + :param row_min: Minimum row index of the bounding box to split. + :param col_min: Minimum column index of the bounding box to split. + :param row_max: Maximum row index of the bounding box to split. + :param col_max: Maximum column index of the bounding box to split. + :param row_split: Height of each tile. + :param col_split: Width of each tile. + :param overlap: size of overlapping between tiles (both vertically and horizontally). + :return: A numpy array grid with splits in two dimensions (0: row, 1: column), + where each cell contains [row_min, row_max, col_min, col_max]. + """ + # Calculate the number of splits considering overlap + nb_col_split = math.ceil((col_max - col_min) / col_split) + nb_row_split = math.ceil((row_max - row_min) / row_split) + + # Initialize the output grid + tiling_grid = np.zeros(shape=(nb_row_split, nb_col_split, 4), dtype=int) + + for row in range(nb_row_split): + for col in range(nb_col_split): + # Calculate the start of the tile + row_start = max(row_min + row * row_split - overlap, 0) + col_start = max(col_min + col * col_split - overlap, 0) + + # Calculate the end of the tile ensuring it doesn't exceed the bounds + row_end = min(row_max, (row + 1) * row_split + overlap) + col_end = min(col_max, (col + 1) * col_split + overlap) + + # Populate the grid with the tile boundaries + tiling_grid[row, col] = [row_start, row_end, col_start, col_end] + + return tiling_grid diff --git a/tests/test_raster/test_georeferencing.py b/tests/test_raster/test_georeferencing.py index 1a16a936..5b7514bb 100644 --- a/tests/test_raster/test_georeferencing.py +++ b/tests/test_raster/test_georeferencing.py @@ -5,7 +5,6 @@ import geoutils as gu from geoutils import examples -from geoutils.raster.georeferencing import _generate_tiling_grid class TestGeoreferencing: @@ -210,59 +209,3 @@ def test_coords(self, example: str) -> None: xxgrid, yygrid = img.coords(grid=True, force_offset="ll") assert np.array_equal(xxgrid, np.repeat(xx0[np.newaxis, :], img.height, axis=0)) assert np.array_equal(yygrid, np.flipud(np.repeat(yy0[:, np.newaxis], img.width, axis=1))) - - @pytest.mark.parametrize("overlap", [0, 5]) # type: ignore - def test_tiling(self, overlap: int) -> None: - - # Test with mock data - tiling_grid_mock = _generate_tiling_grid(0, 0, 100, 100, 50, 50, overlap) - if overlap == 0: - expected_tiling = np.array([[[0, 50, 0, 50], [0, 50, 50, 100]], [[50, 100, 0, 50], [50, 100, 50, 100]]]) - assert np.array_equal(tiling_grid_mock, expected_tiling) - elif overlap == 5: - expected_tiling = np.array( - [ - [[0, 50, 0, 50], [0, 50, 45, 95], [0, 50, 90, 100]], - [[45, 95, 0, 50], [45, 95, 45, 95], [45, 95, 90, 100]], - [[90, 100, 0, 50], [90, 100, 45, 95], [90, 100, 90, 100]], - ] - ) - assert np.array_equal(tiling_grid_mock, expected_tiling) - - # Test with real data - img = gu.Raster(self.landsat_b4_path) - - # Define tiling parameters - row_split, col_split = 100, 100 - row_max, col_max = img.shape - - # Generate the tiling grid - tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, row_split, col_split, overlap) - - # Calculate expected number of tiles - nb_row_tiles = np.ceil(row_max / (row_split - overlap)).astype(int) - nb_col_tiles = np.ceil(col_max / (col_split - overlap)).astype(int) - - # Check that the tiling grid has the expected shape - assert tiling_grid.shape == (nb_row_tiles, nb_col_tiles, 4) - - # Check the boundaries of the first and last tile - assert np.array_equal(tiling_grid[0, 0], np.array([0, min(row_split, row_max), 0, min(col_split, col_max)])) - assert np.array_equal( - tiling_grid[-1, -1], - np.array( - [ - (nb_row_tiles - 1) * (row_split - overlap), - row_max, - (nb_col_tiles - 1) * (col_split - overlap), - col_max, - ] - ), - ) - - # Check if overlap is consistent between tiles - for row in range(nb_row_tiles - 1): - assert tiling_grid[row + 1, 0, 0] == tiling_grid[row, 0, 1] - overlap - - for col in range(nb_col_tiles - 1): # Skip last tile in column - assert tiling_grid[0, col + 1, 2] == tiling_grid[0, col, 3] - overlap diff --git a/tests/test_raster/test_tiling.py b/tests/test_raster/test_tiling.py new file mode 100644 index 00000000..90e9db01 --- /dev/null +++ b/tests/test_raster/test_tiling.py @@ -0,0 +1,94 @@ +"""Test tiling tools for arrays and rasters.""" + +import numpy as np +import pytest + +import geoutils as gu +from geoutils import examples +from geoutils.raster.tiling import _generate_tiling_grid + + +class TestTiling: + + landsat_b4_path = examples.get_path("everest_landsat_b4") + + def test_subdivide_array(self) -> None: + test_shape = (6, 4) + test_count = 4 + subdivision_grid = gu.raster.subdivide_array(test_shape, test_count) + + assert subdivision_grid.shape == test_shape + assert np.unique(subdivision_grid).size == test_count + + assert np.unique(gu.raster.subdivide_array((3, 3), 3)).size == 3 + + with pytest.raises(ValueError, match=r"Expected a 2D shape, got 1D shape.*"): + gu.raster.subdivide_array((5,), 2) + + with pytest.raises(ValueError, match=r"Shape.*smaller than.*"): + gu.raster.subdivide_array((5, 2), 15) + + @pytest.mark.parametrize("overlap", [0, 5]) # type: ignore + def test_tiling(self, overlap: int) -> None: + + # Test with mock data + tiling_grid_mock = _generate_tiling_grid(0, 0, 100, 100, 50, 50, overlap) + if overlap == 0: + expected_tiling = np.array([[[0, 50, 0, 50], [0, 50, 50, 100]], [[50, 100, 0, 50], [50, 100, 50, 100]]]) + assert np.array_equal(tiling_grid_mock, expected_tiling) + elif overlap == 5: + expected_tiling = np.array( + [ + [[0, 55, 0, 55], [0, 55, 45, 100]], + [[45, 100, 0, 55], [45, 100, 45, 100]], + ] + ) + assert np.array_equal(tiling_grid_mock, expected_tiling) + + # Test with real data + img = gu.Raster(self.landsat_b4_path) + + # Define tiling parameters + row_split, col_split = 100, 100 + row_max, col_max = img.shape + + # Generate the tiling grid + tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, row_split, col_split, overlap) + + # Calculate expected number of tiles + nb_row_tiles = np.ceil(row_max / row_split).astype(int) + nb_col_tiles = np.ceil(col_max / col_split).astype(int) + + # Check that the tiling grid has the expected shape + assert tiling_grid.shape == (nb_row_tiles, nb_col_tiles, 4) + + # Check the boundaries of the first and last tile + assert np.array_equal( + tiling_grid[0, 0], + np.array( + [ + 0, + min(row_split + overlap, row_max), + 0, + min(col_split + overlap, col_max), + ] + ), + ) + assert np.array_equal( + tiling_grid[-1, -1], + np.array( + [ + (nb_row_tiles - 1) * row_split - overlap, + row_max, + (nb_col_tiles - 1) * col_split - overlap, + col_max, + ] + ), + ) + + # Check if overlap is consistent between tiles + for row in range(nb_row_tiles - 1): + assert tiling_grid[row + 1, 0, 0] == tiling_grid[row, 0, 1] - 2 * overlap + + for col in range(nb_col_tiles - 1): + assert tiling_grid[0, col + 1, 2] == tiling_grid[0, col, 3] - 2 * overlap From d861eb99f0ded14eb190d8849e379406e96e3f7f Mon Sep 17 00:00:00 2001 From: vschaffn Date: Tue, 18 Feb 2025 10:01:34 +0100 Subject: [PATCH 4/4] fix: move tiling functions from raster to tiling --- geoutils/raster/raster.py | 40 ------------------------------------- geoutils/raster/tiling.py | 42 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 40 deletions(-) diff --git a/geoutils/raster/raster.py b/geoutils/raster/raster.py index 44e7754f..99cf289d 100644 --- a/geoutils/raster/raster.py +++ b/geoutils/raster/raster.py @@ -41,7 +41,6 @@ import rioxarray import xarray as xr from affine import Affine -from matplotlib.patches import Rectangle from mpl_toolkits.axes_grid1 import make_axes_locatable from packaging.version import Version from rasterio.crs import CRS @@ -89,7 +88,6 @@ decode_sensor_metadata, parse_and_convert_metadata_from_filename, ) -from geoutils.raster.tiling import _generate_tiling_grid from geoutils.stats import nmad from geoutils.vector.vector import Vector @@ -2657,44 +2655,6 @@ def translate( raster_copy.transform = translated_transform return raster_copy - def compute_tiling( - self, - tile_size: int, - raster_ref: RasterType, - overlap: int = 0, - ) -> NDArrayNum: - """ - Compute the raster tiling grid to coregister raster by block. - - :param tile_size: Size of each tile (square tiles) - :param raster_ref: The other raster to coregister, use to validate the shape - :param overlap: Size of overlap between tiles (optional) - :return: tiling_grid (array of tile boundaries), new_shape (shape of the tiled grid) - """ - if self.shape != raster_ref.shape: - raise Exception("Reference and secondary rasters do not have the same shape") - row_max, col_max = self.shape - - # Generate tiling - tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, tile_size, tile_size, overlap=overlap) - return tiling_grid - - def plot_tiling(self, tiling_grid: NDArrayNum) -> None: - """ - Plot raster with its tiling. - - :param tiling_grid: tiling given by Raster.compute_tiling. - """ - ax, caxes = self.plot(return_axes=True) - for tile in tiling_grid.reshape(-1, 4): - row_min, row_max, col_min, col_max = tile - x_min, y_min = self.transform * (col_min, row_min) # Bottom-left corner - x_max, y_max = self.transform * (col_max, row_max) # Top-right corne - rect = Rectangle( - (x_min, y_min), x_max - x_min, y_max - y_min, edgecolor="red", facecolor="none", linewidth=1.5 - ) - ax.add_patch(rect) - def save( self, filename: str | pathlib.Path | IO[bytes], diff --git a/geoutils/raster/tiling.py b/geoutils/raster/tiling.py index 569f807c..2d034c9e 100644 --- a/geoutils/raster/tiling.py +++ b/geoutils/raster/tiling.py @@ -21,8 +21,10 @@ import math import numpy as np +from matplotlib.patches import Rectangle from geoutils._typing import NDArrayNum +from geoutils.raster import RasterType def _get_closest_rectangle(size: int) -> tuple[int, int]: @@ -166,3 +168,43 @@ def _generate_tiling_grid( tiling_grid[row, col] = [row_start, row_end, col_start, col_end] return tiling_grid + + +def compute_tiling( + tile_size: int, + raster_shape: tuple[int, int], + ref_shape: tuple[int, int], + overlap: int = 0, +) -> NDArrayNum: + """ + Compute the raster tiling grid to coregister raster by block. + + :param tile_size: Size of each tile (square tiles). + :param raster_shape: Shape of the raster to determine tiling parameters. + :param ref_shape: The shape of another raster to coregister, use to validate the shape. + :param overlap: Size of overlap between tiles (optional). + :return: tiling_grid (array of tile boundaries). + """ + if raster_shape != ref_shape: + raise Exception("Reference and secondary rasters do not have the same shape") + row_max, col_max = raster_shape + + # Generate tiling + tiling_grid = _generate_tiling_grid(0, 0, row_max, col_max, tile_size, tile_size, overlap=overlap) + return tiling_grid + + +def plot_tiling(raster: RasterType, tiling_grid: NDArrayNum) -> None: + """ + Plot raster with its tiling. + + :param raster: The raster to plot with its tiling. + :param tiling_grid: tiling given by compute_tiling. + """ + ax, caxes = raster.plot(return_axes=True) + for tile in tiling_grid.reshape(-1, 4): + row_min, row_max, col_min, col_max = tile + x_min, y_min = raster.transform * (col_min, row_min) # Bottom-left corner + x_max, y_max = raster.transform * (col_max, row_max) # Top-right corne + rect = Rectangle((x_min, y_min), x_max - x_min, y_max - y_min, edgecolor="red", facecolor="none", linewidth=1.5) + ax.add_patch(rect)