From 1562a1efa8b2560116f5d7cbf961a672f2b0f4dc Mon Sep 17 00:00:00 2001 From: Mike Taves Date: Thu, 15 Aug 2024 15:07:59 +1200 Subject: [PATCH] Feature: allow different snap modes used for grid size and offset (#38) --- gridit/__main__.py | 1 + gridit/classmethods.py | 154 ++++++++++++++++++++++++----- gridit/cli.py | 21 ++++ tests/test_classmethods.py | 192 ++++++++++++++++++++++++++++++++++--- tests/test_cli.py | 18 ++++ 5 files changed, 348 insertions(+), 38 deletions(-) diff --git a/gridit/__main__.py b/gridit/__main__.py index 70dea40..dbb8946 100644 --- a/gridit/__main__.py +++ b/gridit/__main__.py @@ -366,6 +366,7 @@ def write_output_with_parts(ar): error(err, show_usage=False) logger.info("%s", grid) + logger.info("bounds: %s", grid.bounds) logger.info("has mask: %s", mask is not None) # Process array from * options diff --git a/gridit/classmethods.py b/gridit/classmethods.py index b85956a..430e9ff 100644 --- a/gridit/classmethods.py +++ b/gridit/classmethods.py @@ -1,13 +1,57 @@ """Grid from_* classmethods.""" +from decimal import Decimal +from itertools import product from math import ceil, floor -from typing import Optional +from typing import Optional, Union from gridit.logger import get_logger +_lr = ["left", "right"] +_tb = ["top", "bottom"] +snap_modes = ["full", "half"] + list("-".join(two) for two in product(_tb, _lr)) -def get_shape_top_left(bounds, resolution, buffer=0.0): - minx, miny, maxx, maxy = bounds + +def get_shape_top_left( + bounds: tuple, + resolution: Union[float, Decimal], + buffer: Union[float, Decimal] = Decimal("0.0"), + snap: Union[str, tuple] = "full", +): + """Get shape and top-left coordinate to definea grid from bounds. + + Parameters + ---------- + bounds : tuple of float or Decimal + Bounding box, ordered (minx, miny, maxx, maxy). + resolution : float or Decimal + Grid resolution for x- and y-directions. + buffer : float or Decmial, default 0.0 + Optional buffer distance to expand bounds. + snap : {full, half, top-left, top-right, bottom-left, bottom-right} or tuple + Snap mode used to evaluate grid size and offset. Default 'full' will + snap bounds to a multiple of resolution, and 'half' will snap to + half-resolution. Corner specifications, e.g. 'bottom-left' will + snap the grid to align with this coordinate. Alternatively, + a coordinate tuple (snapx, snapy) can be provided to snap the grid, + although the grid does not necessarily include the coordinate. + + Returns + ------- + shape: tuple of int + Dimensions of grid (ny, nx). + top_left: tuple of float + Snapped top-left corner of grid (minx, maxy). + + """ + if all(isinstance(x, Decimal) for x in bounds): + minx, miny, maxx, maxy = bounds + else: + minx, miny, maxx, maxy = map(lambda x: Decimal(str(x)), bounds) + if not isinstance(resolution, Decimal): + resolution = Decimal(str(resolution)) + if not isinstance(buffer, Decimal): + buffer = Decimal(str(buffer)) if not (minx <= maxx): raise ValueError("'minx' must be less than 'maxx'") elif not (miny <= maxy): @@ -22,24 +66,55 @@ def get_shape_top_left(bounds, resolution, buffer=0.0): maxx += buffer maxy += buffer dx = dy = resolution - if buffer > 0.0: - minx = dx * round(minx / dx) - miny = dy * round(miny / dy) - maxx = dx * round(maxx / dx) - maxy = dy * round(maxy / dy) + if snap == "full": + snapx = snapy = Decimal("0.0") + elif snap == "half": + snapx = dx / Decimal("2.0") + snapy = dy / Decimal("2.0") + elif isinstance(snap, tuple): + if len(snap) != 2: + raise TypeError("'snap' tuple must have 2 items: (snapx, snapy)") + if all(isinstance(x, Decimal) for x in snap): + snapx, snapy = snap + else: + snapx, snapy = map(lambda x: Decimal(str(x)), snap) + elif snap in snap_modes: + if "top" in snap: + snapy = maxy + else: + assert "bottom" in snap, snap + leny = maxy - miny + ny = ceil(leny / dy) or 1 + snapy = miny + ny * dy + if leny == 0.0: + miny += dy + maxy += dy + if "left" in snap: + snapx = minx + else: + assert "right" in snap, snap + lenx = maxx - minx + nx = ceil(lenx / dx) or 1 + snapx = maxx - nx * dx + if lenx == 0.0: + minx -= dx + maxx -= dx else: - minx = dx * floor(minx / dx) - miny = dy * floor(miny / dy) - maxx = dx * ceil(maxx / dx) - maxy = dy * ceil(maxy / dy) - lenx = maxx - minx - leny = maxy - miny - assert lenx % dx == 0.0 - assert leny % dy == 0.0 - nx = int(lenx / dx) - ny = int(leny / dy) + raise ValueError(f"'snap' must be one of {snap_modes} or tuple (snapx, snapy)") + snapx %= dx + snapy %= dy + minx = dx * floor((minx - snapx) / dx) + snapx + maxx = dx * ceil((maxx - snapx) / dx) + snapx + miny = dy * floor((miny - snapy) / dy) + snapy + maxy = dy * ceil((maxy - snapy) / dy) + snapy + nx = int((maxx - minx) / dx) or 1 + ny = int((maxy - miny) / dy) or 1 shape = ny, nx - top_left = (minx, maxy) + top_left = (float(minx), float(maxy)) + # print( + # f"POLYGON (({minx} {maxy}, {minx} {miny}, " + # f"{maxx} {miny}, {maxx} {maxy}, {minx} {maxy}))" + # ) return shape, top_left @@ -51,7 +126,9 @@ def from_bbox( maxx: float, maxy: float, resolution: float, + *, buffer: float = 0.0, + snap: Union[str, tuple] = "full", projection: Optional[str] = None, logger=None, ): @@ -69,6 +146,13 @@ def from_bbox( A grid resolution, e.g. 250.0 for 250m x 250m buffer : float, default 0.0 Add buffer to extents of bounding box. + snap : {full, half, top-left, top-right, bottom-left, bottom-right} or tuple + Snap mode used to evaluate grid size and offset. Default 'full' will + snap bounds to a multiple of resolution, and 'half' will snap to + half-resolution. Corner specifications, e.g. 'bottom-left' will + snap the grid to align with this coordinate. Alternatively, + a coordinate tuple (snapx, snapy) can be provided to snap the grid, + although the grid does not necessarily include the coordinate. projection : optional str, default None Coordinate reference system described as a string either as (e.g.) EPSG:2193 or a WKT string. @@ -97,7 +181,7 @@ def from_bbox( logger = get_logger(cls.__name__) logger.info("creating from a bounding box") bounds = minx, miny, maxx, maxy - shape, top_left = get_shape_top_left(bounds, resolution, buffer) + shape, top_left = get_shape_top_left(bounds, resolution, buffer, snap) return cls( resolution=resolution, shape=shape, @@ -109,7 +193,13 @@ def from_bbox( @classmethod def from_raster( - cls, fname: str, resolution: float = None, buffer: float = 0.0, logger=None + cls, + fname: str, + resolution: Optional[float] = None, + *, + buffer: float = 0.0, + snap: Union[str, tuple] = "full", + logger=None, ): """Fetch grid information from a raster. @@ -123,6 +213,13 @@ def from_raster( expanded and "snapped" to a multiple of the resolution. buffer : float, default 0.0. Add buffer to extents of raster. + snap : {full, half, top-left, top-right, bottom-left, bottom-right} or tuple + Snap mode used to evaluate grid size and offset. Default 'full' will + snap bounds to a multiple of resolution, and 'half' will snap to + half-resolution. Corner specifications, e.g. 'bottom-left' will + snap the grid to align with this coordinate. Alternatively, + a coordinate tuple (snapx, snapy) can be provided to snap the grid, + although the grid does not necessarily include the coordinate. logger : logging.Logger, optional Logger to show messages. @@ -154,7 +251,7 @@ def from_raster( resolution = t.a ny, nx = shape bounds = t.c, t.f + ny * t.e, t.c + nx * t.a, t.f - shape, top_left = get_shape_top_left(bounds, resolution, buffer) + shape, top_left = get_shape_top_left(bounds, resolution, buffer, snap) else: resolution = t.a top_left = t.c, t.f @@ -172,8 +269,10 @@ def from_vector( cls, fname: str, resolution: float, - filter: dict = None, + *, + filter: Union[dict, str, None] = None, buffer: float = 0.0, + snap: Union[str, tuple] = "full", layer=None, logger=None, ): @@ -194,6 +293,13 @@ def from_vector( used if Fiona 1.9 or later is installed. buffer : float, default 0.0 Add buffer to extents of vector data. + snap : {full, half, top-left, top-right, bottom-left, bottom-right} or tuple + Snap mode used to evaluate grid size and offset. Default 'full' will + snap bounds to a multiple of resolution, and 'half' will snap to + half-resolution. Corner specifications, e.g. 'bottom-left' will + snap the grid to align with this coordinate. Alternatively, + a coordinate tuple (snapx, snapy) can be provided to snap the grid, + although the grid does not necessarily include the coordinate. layer : int or str, default None The integer index or name of a layer in a multi-layer dataset. logger : logging.Logger, optional @@ -237,7 +343,7 @@ def from_vector( flt.close() else: # full shapefile bounds bounds = ds.bounds - shape, top_left = get_shape_top_left(bounds, resolution, buffer) + shape, top_left = get_shape_top_left(bounds, resolution, buffer, snap) return cls( resolution=resolution, shape=shape, diff --git a/gridit/cli.py b/gridit/cli.py index 5823156..86e5734 100644 --- a/gridit/cli.py +++ b/gridit/cli.py @@ -2,6 +2,8 @@ __all__ = [] +import re +from decimal import Decimal from importlib.util import find_spec from gridit.grid import Grid @@ -63,6 +65,18 @@ def add_grid_parser_arguments(parser): default=0.0, help="Add buffer to extents of grid, default 0.", ) + grid_group.add_argument( + "--snap", + metavar="SNAP", + type=str, + default="full", + help="Snap mode used to evaluate grid size and offset. Default 'full' will " + "snap bounds to a multiple of resolution, and 'half' will snap to " + "half-resolution. Corner specifications, e.g. 'bottom-left' will " + "snap the grid to align with this coordinate. Alternatively, " + "a coordinate tuple (snapx, snapy) can be provided to snap the grid, " + "although the grid does not necessarily include the coordinate.", + ) grid_group.add_argument( "--projection", metavar="STR", @@ -101,6 +115,13 @@ def error_msg(msg: str, name: str = ""): grid_args["resolution"] = args.resolution if args.buffer: grid_args["buffer"] = args.buffer + if args.snap: + if snapxy := re.findall(r"([\-\+]?[\d\.]+)", args.snap): + if len(snapxy) != 2: + raise ValueError("snap tuple must have two floats") + grid_args["snap"] = tuple(map(Decimal, snapxy)) + else: + grid_args["snap"] = args.snap if args.projection: grid_args["projection"] = args.projection from_grid_methods = ["bbox", "raster", "vector"] diff --git a/tests/test_classmethods.py b/tests/test_classmethods.py index 6887607..ec27827 100644 --- a/tests/test_classmethods.py +++ b/tests/test_classmethods.py @@ -1,9 +1,138 @@ import pytest from gridit import Grid +from gridit.classmethods import get_shape_top_left from .conftest import datadir, requires_pkg +# low-level method + + +@pytest.mark.parametrize( + "resolution, snap, top_left", + [ + (0.2, "full", (42.0, 73.0)), + (0.2, "half", (41.9, 73.1)), + (0.2, "top-left", (42.0, 73.0)), + (0.2, "top-right", (41.8, 73.0)), + (0.2, "bottom-left", (42.0, 73.2)), + (0.2, "bottom-right", (41.8, 73.2)), + (0.2, (1.0, 1.0), (42.0, 73.0)), + (0.2, (2.0, 2.0), (42.0, 73.0)), + (1.0, "full", (42.0, 73.0)), + (1.0, "half", (41.5, 73.5)), + (1.0, "top-left", (42.0, 73.0)), + (1.0, "top-right", (41.0, 73.0)), + (1.0, "bottom-left", (42.0, 74.0)), + (1.0, "bottom-right", (41.0, 74.0)), + (1.0, (1.0, 1.0), (42.0, 73.0)), + (1.0, (2.0, 2.0), (42.0, 73.0)), + (2.0, "full", (42.0, 74.0)), + (2.0, "half", (41.0, 73.0)), + (2.0, "top-left", (42.0, 73.0)), + (2.0, "top-right", (40.0, 73.0)), + (2.0, "bottom-left", (42.0, 75.0)), + (2.0, "bottom-right", (40.0, 75.0)), + (2.0, (1.0, 1.0), (41.0, 73.0)), + (2.0, (2.0, 2.0), (42.0, 74.0)), + (5.0, "full", (40.0, 75.0)), + (5.0, "half", (37.5, 77.5)), + (5.0, (1.0, 2.0), (41.0, 77.0)), + (5.0, (2.8, 2.7), (37.8, 77.7)), + (10.0, "full", (40.0, 80.0)), + (20.0, "full", (40.0, 80.0)), + (33.3, "full", (33.3, 99.9)), + ], +) +def test_get_shape_top_left_point(resolution, snap, top_left): + x, y = 42.0, 73.0 # point + buffer = 0.0 + shape, (minx, maxy) = get_shape_top_left((x, y, x, y), resolution, buffer, snap) + maxx = minx + resolution + miny = maxy - resolution + assert top_left == (minx, maxy) + assert shape == (1, 1) + # check if point is within bounds + assert minx <= x <= maxx + assert miny <= y <= maxy + # check snapped edges touch + if isinstance(snap, str): + if "top" in snap: + assert maxy == y + if "bottom" in snap: + assert miny == y + if "left" in snap: + assert minx == x + if "right" in snap: + assert maxx == x + else: + # tuple snaps to top left + snapx, snapy = snap + assert (snapx % resolution) == pytest.approx(minx % resolution) + assert (snapy % resolution) == pytest.approx(maxy % resolution) + + +@pytest.mark.parametrize( + "resolution, buffer, snap, top_left, shape", + [ + (0.2, 0.01, "full", (41.8, 73.2), (2, 2)), + (0.2, 0.01, "half", (41.9, 73.1), (1, 1)), + (0.2, 0.01, "top-left", (41.99, 73.01), (1, 1)), + (0.2, 0.01, "top-right", (41.81, 73.01), (1, 1)), + (0.2, 0.01, "bottom-left", (41.99, 73.19), (1, 1)), + (0.2, 0.01, "bottom-right", (41.81, 73.19), (1, 1)), + (0.2, 0.01, (1.0, 1.0), (41.8, 73.2), (2, 2)), + (0.2, 0.01, (12.9, 12.71), (41.9, 73.11), (1, 1)), + (1.0, 0.5, "half", (41.5, 73.5), (1, 1)), + (1.0, 1.0, "full", (41.0, 74.0), (2, 2)), + (1.0, 1.0, "half", (40.5, 74.5), (3, 3)), + (1.0, 1.0, "top-left", (41.0, 74.0), (2, 2)), + (1.0, 1.0, "top-right", (41.0, 74.0), (2, 2)), + (1.0, 1.0, "bottom-left", (41.0, 74.0), (2, 2)), + (1.0, 1.0, "bottom-right", (41.0, 74.0), (2, 2)), + (1.0, 1.0, (-1.0, -1.0), (41.0, 74.0), (2, 2)), + (1.0, 1.0, (2.0, 2.0), (41.0, 74.0), (2, 2)), + (2.0, 2.0, "full", (40.0, 76.0), (3, 2)), + (2.0, 2.0, "half", (39.0, 75.0), (2, 3)), + (2.0, 2.0, "top-left", (40.0, 75.0), (2, 2)), + (2.0, 2.0, "top-right", (40.0, 75.0), (2, 2)), + (2.0, 2.0, "bottom-left", (40.0, 75.0), (2, 2)), + (2.0, 2.0, "bottom-right", (40.0, 75.0), (2, 2)), + (2.0, 2.0, (-1.0, -1.0), (39.0, 75.0), (2, 3)), + (2.0, 2.0, (2022.0, 2222.0), (40.0, 76.0), (3, 2)), + (33.3, 0.01, "full", (33.3, 99.9), (1, 1)), + (33.3, 6.7, "full", (33.3, 99.9), (2, 1)), + (33.3, 9.0, "full", (0.0, 99.9), (2, 2)), + (33.3, 9.0, "full", (0.0, 99.9), (2, 2)), + ], +) +def test_get_shape_top_left_buffer(resolution, buffer, snap, top_left, shape): + x, y = 42.0, 73.0 # point + (ny, nx), (minx, maxy) = get_shape_top_left((x, y, x, y), resolution, buffer, snap) + assert top_left == (minx, maxy) + assert shape == (ny, nx) + # check if point is within bounds + maxx = minx + resolution * nx + miny = maxy - resolution * ny + assert minx <= x <= maxx + assert miny <= y <= maxy + # check snapped edges touch + if isinstance(snap, str): + if "top" in snap: + assert maxy == pytest.approx(y + buffer) + if "bottom" in snap: + assert miny == pytest.approx(y - buffer) + if "left" in snap: + assert minx == pytest.approx(x - buffer) + if "right" in snap: + assert maxx == pytest.approx(x + buffer) + else: + # tuple snaps to top left + snapx, snapy = snap + assert (snapx % resolution) == pytest.approx(minx % resolution) + assert (snapy % resolution) == pytest.approx(maxy % resolution) + + mana_dem_path = datadir / "Mana.tif" mana_polygons_path = datadir / "Mana_polygons.shp" lines_path = datadir / "waitaku2_lines.shp" @@ -27,8 +156,10 @@ def test_grid_from_bbox_point(): def test_grid_from_bbox_buffer(): - grid = Grid.from_bbox(1748762.8, 5448908.9, 1749509, 5449749, 25, 20, "EPSG:2193") - expected = Grid(25.0, (35, 31), (1748750.0, 5449775.0), "EPSG:2193") + grid = Grid.from_bbox( + 1748762.8, 5448908.9, 1749509, 5449749, 25, buffer=20, projection="EPSG:2193" + ) + expected = Grid(25.0, (36, 33), (1748725.0, 5449775.0), "EPSG:2193") assert grid == expected @@ -57,9 +188,27 @@ def test_grid_from_raster_buffer(): expected = Grid(8.0, (282, 213), (1748672.0, 5451112.0), grid.projection) assert grid == expected + # with snap mode + grid = Grid.from_raster(mana_dem_path, buffer=16.0, snap="bottom-right") + assert grid == expected + + # with snap half mode + grid = Grid.from_raster(mana_dem_path, buffer=16.0, snap="half") + expected = Grid(8.0, (283, 214), (1748668.0, 5451116.0), grid.projection) + assert grid == expected + # with buffer - grid = Grid.from_raster(mana_dem_path, 10.0, 20.0) - expected = Grid(10.0, (227, 171), (1748670.0, 5451120.0), grid.projection) + grid = Grid.from_raster(mana_dem_path, 10.0, buffer=20.0) + expected = Grid(10.0, (227, 172), (1748660.0, 5451120.0), grid.projection) + assert grid == expected + + # with buffer + snap modes + grid = Grid.from_raster(mana_dem_path, 10.0, buffer=20.0, snap="top-left") + expected = Grid(10.0, (227, 172), (1748668.0, 5451116.0), grid.projection) + assert grid == expected + + grid = Grid.from_raster(mana_dem_path, 10.0, buffer=20.0, snap="bottom-left") + expected = Grid(10.0, (227, 172), (1748668.0, 5451122.0), grid.projection) assert grid == expected @@ -78,18 +227,23 @@ def test_grid_from_vector_point(): expected = Grid(250.0, (50, 28), (1810500.0, 5877750.0), grid.projection) assert grid == expected + # all + half snap mode + grid = Grid.from_vector(points_path, 250, snap="half") + expected = Grid(250.0, (49, 27), (1810625.0, 5877625.0), grid.projection) + assert grid == expected + # filter - grid = Grid.from_vector(points_path, 250, {"id": [5, 9]}) + grid = Grid.from_vector(points_path, 250, filter={"id": [5, 9]}) expected = Grid(250.0, (19, 7), (1810500.0, 5873750.0), grid.projection) assert grid == expected - grid = Grid.from_vector(points_path, 250, {"id": 5}) + grid = Grid.from_vector(points_path, 250, filter={"id": 5}) expected = Grid(250.0, (1, 1), (1812000.0, 5869250.0), grid.projection) assert grid == expected # filter + buffer - grid = Grid.from_vector(points_path, 250, {"id": 5}, buffer=240) - expected = Grid(250.0, (2, 2), (1811750.0, 5869250.0), grid.projection) + grid = Grid.from_vector(points_path, 250, filter={"id": 5}, buffer=240) + expected = Grid(250.0, (3, 3), (1811750.0, 5869500.0), grid.projection) assert grid == expected @@ -113,14 +267,19 @@ def test_grid_from_vector_polygon(): grid = Grid.from_vector(datadir, 100, layer="mana_polygons") assert grid == expected + # layer + half snap mode + grid = Grid.from_vector(datadir, 100, layer="mana_polygons", snap="half") + expected = Grid(100.0, (23, 18), (1748650.0, 5451150.0), grid.projection) + assert grid == expected + # filter - grid = Grid.from_vector(mana_polygons_path, 100, {"name": "South-east"}) + grid = Grid.from_vector(mana_polygons_path, 100, filter={"name": "South-east"}) expected = Grid(100.0, (14, 13), (1749100.0, 5450400.0), grid.projection) assert grid == expected # buffer grid = Grid.from_vector(mana_polygons_path, 100, buffer=500) - expected = Grid(100.0, (32, 27), (1748200.0, 5451600.0), grid.projection) + expected = Grid(100.0, (34, 28), (1748100.0, 5451700.0), grid.projection) assert grid == expected @@ -132,17 +291,22 @@ def test_grid_from_vector_line(): assert grid == expected # filter - grid = Grid.from_vector(lines_path, 250, {"StreamOrde": 5}) + grid = Grid.from_vector(lines_path, 250, filter={"StreamOrde": 5}) expected = Grid(250.0, (19, 14), (1808750.0, 5877000.0), grid.projection) assert grid == expected - grid = Grid.from_vector(lines_path, 250, {"StreamOrde": [4, 5]}) + grid = Grid.from_vector(lines_path, 250, filter={"StreamOrde": [4, 5]}) expected = Grid(250.0, (28, 41), (1804750.0, 5877000.0), grid.projection) assert grid == expected # buffer grid = Grid.from_vector(lines_path, 250, buffer=500) - expected = Grid(250.0, (70, 66), (1803000.0, 5878750.0), grid.projection) + expected = Grid(250.0, (71, 68), (1802750.0, 5879000.0), grid.projection) + assert grid == expected + + # buffer + half snap mode + grid = Grid.from_vector(lines_path, 250, buffer=500, snap="half") + expected = Grid(250.0, (71, 67), (1802875.0, 5878875.0), grid.projection) assert grid == expected @@ -154,6 +318,6 @@ def test_grid_from_vector_filter_sql_where(): pytest.skip("Fiona 1.9 or later required to use SQL WHERE") # filter - grid = Grid.from_vector(lines_path, 250, "StreamOrde>=5") + grid = Grid.from_vector(lines_path, 250, filter="StreamOrde>=5") expected = Grid(250.0, (19, 14), (1808750.0, 5877000.0), grid.projection) assert grid == expected diff --git a/tests/test_cli.py b/tests/test_cli.py index 9699d89..b3fde99 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -35,6 +35,24 @@ def test_grid_from_bbox(grid_from_bbox_args): assert returncode == 0 +def test_grid_from_bbox_more_args(grid_from_bbox_args): + stdout, stderr, returncode = run_cli(grid_from_bbox_args + ["--snap", "half"]) + assert len(stderr) == 0 + assert len(stdout) > 0 + assert returncode == 0 + + stdout, stderr, returncode = run_cli(grid_from_bbox_args + ["--snap", "(1.1 -2.9)"]) + assert len(stderr) == 0 + assert len(stdout) > 0 + assert returncode == 0 + + # bad arg + stdout, stderr, returncode = run_cli(grid_from_bbox_args + ["--snap", "HALF"]) + assert len(stderr) > 0 + assert len(stdout) > 0 + assert returncode != 0 + + @requires_pkg("rasterio") def test_grid_from_bbox_array_from_raster(tmp_path, grid_from_bbox_args): out_path = tmp_path / "out.tif"