diff --git a/CHANGELOG.md b/CHANGELOG.md index 743f5bf29..7d8f40f36 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 ## [0.54.0] - 2026-10-28 +- Add public `group_geometries` helper in `openeo.extra.job_management` for geometry-grouping based job splitting, including optional tile-grid co-grouping. ### Changed diff --git a/docs/cookbook/job_manager.rst b/docs/cookbook/job_manager.rst index 20b76f912..b7de4f602 100644 --- a/docs/cookbook/job_manager.rst +++ b/docs/cookbook/job_manager.rst @@ -202,6 +202,7 @@ to resume where you left off. Parquet support requires the ``pyarrow`` package (see :ref:`optional dependencies `). + STAC API (experimental) ----------------------- @@ -247,7 +248,7 @@ Example: splitting a bounding box into 10x10 km tiles in EPSG:3857: # Split into 10x10 km tiles in EPSG:3857 gdf = split_area( - aoi=extent, + spatial_extent=extent, tile_size=10_000, # tile size in meters (units of the projection) projection="EPSG:3857", # projection for tiling ) @@ -300,7 +301,80 @@ or any other custom partitioning, you can pass a pre-computed tile_grid = gpd.read_parquet("my_tile_grid.parquet") # Only tiles intersecting the AOI are returned - gdf = split_area(aoi=extent , tile_grid=tile_grid) + gdf = split_area(spatial_extent=extent, tile_grid=tile_grid) + +Grouping geometries (e.g. point-based extractions) +-------------------------------------------------- + +When the work to split is a set of geometries (e.g. sample locations) +rather than a single large area, e.g. for point-based extraction jobs, +:py:func:`~openeo.extra.job_management._job_splitting.group_geometries` groups +the rows of a GeoDataFrame into batches, optionally of at most ``max_geometries`` rows. +Each row is treated as a single feature to group, based on its representative point, +regardless of geometry type (e.g. ``Point``, ``MultiPoint``, ``Polygon``). +Null or empty geometries are not allowed. +Geometries are not split themselves: use ``split_area`` for a large area of interest. +Like with ``split_area``, geometries without +CRS information are assumed to be in the CRS of the tile grid (if any). +Groups are formed by recursive bisection at the median along the longer axis, +so that batches stay geographically clustered and balanced. +The returned GeoDataFrames keep the original index, to link back to the input rows: + +.. code-block:: python + + import geopandas as gpd + from openeo.extra.job_management import group_geometries + + points = gpd.read_file("my_sample_locations.gpkg") + batches = group_geometries(points, max_geometries=500) + # batches is a list of GeoDataFrames, each with at most 500 rows + +Optionally, a ``tile_grid`` (e.g. a satellite tiling grid) can be passed so +that each feature is first assigned to the tile containing its representative point +and batches are formed per tile. +Features outside every tile are batched separately: + +.. code-block:: python + + tile_grid = gpd.read_parquet("my_tile_grid.parquet") + batches = group_geometries(points, max_geometries=500, tile_grid=tile_grid) + +Alternatively, a regular grid can be specified through ``tile_size`` and ``projection``. +This grid is anchored at the origin of the projection, +so the grouping of existing features does not change when other features are added: + +.. code-block:: python + + # group points per 50x50 km tile in EPSG:3857 + batches = group_geometries(points, max_geometries=500, tile_size=50_000, projection="EPSG:3857") + +Without ``max_geometries`` (default ``None``), there is no size limit and +geometries are simply grouped per tile: + +.. code-block:: python + + # one batch per tile + batches = group_geometries(points, tile_grid=tile_grid) + +Current limitations: + +- Geometries are never split: a polygon crossing tile edges, or larger than a tile, + ends up in a single batch and extends beyond its tile. + Use ``split_area`` for large areas of interest. +- Multi-part geometries (e.g. ``MultiPolygon``) are kept together as one feature, + placed by a single representative point, even if their parts are far apart. + To group parts individually, explode them into separate rows first + (the original index is kept, to link back to the input rows; + note that a single feature may then be spread over multiple batches): + + .. code-block:: python + + batches = group_geometries(points.explode(index_parts=False), max_geometries=500) + +- Batches are limited by number of geometries only, not by area or processing cost. + Use a tile grid to limit the spatial extent of batches. +- With overlapping tiles (e.g. Sentinel-2 MGRS), a feature is assigned to the first tile + (in grid order) containing its representative point, not to the tile with the largest overlap. Customizing Job Handling ======================== @@ -543,6 +617,7 @@ Job Database .. autoclass:: openeo.extra.job_management.ParquetJobDatabase + .. autofunction:: openeo.extra.job_management.create_job_db .. autofunction:: openeo.extra.job_management.get_job_db diff --git a/openeo/extra/job_management/__init__.py b/openeo/extra/job_management/__init__.py index 02806531e..6e2ac44e8 100644 --- a/openeo/extra/job_management/__init__.py +++ b/openeo/extra/job_management/__init__.py @@ -6,7 +6,7 @@ create_job_db, get_job_db, ) -from openeo.extra.job_management._job_splitting import split_area +from openeo.extra.job_management._job_splitting import group_geometries, split_area from openeo.extra.job_management._manager import MultiBackendJobManager from openeo.extra.job_management.process_based import ProcessBasedJobCreator @@ -19,5 +19,6 @@ "create_job_db", "get_job_db", "split_area", + "group_geometries", "MultiBackendJobManager", ] diff --git a/openeo/extra/job_management/_job_splitting.py b/openeo/extra/job_management/_job_splitting.py index 9cc854c17..56e31f94a 100644 --- a/openeo/extra/job_management/_job_splitting.py +++ b/openeo/extra/job_management/_job_splitting.py @@ -4,6 +4,8 @@ from typing import Dict, List, Optional, Tuple, Union import geopandas as gpd +import numpy as np +import pandas as pd import shapely import shapely.geometry.base from shapely.geometry import MultiPolygon, Polygon @@ -12,6 +14,7 @@ _log = logging.getLogger(__name__) + class JobSplittingFailure(Exception): pass @@ -72,25 +75,6 @@ def _parse_input_geometry( "Invalid bounding box: west must not equal east (zero-width bounding box)." ) - # if bbox["west"] > bbox["east"]: - # # In a geographic CRS this indicates an antimeridian crossing; - # # split into two polygons on each side of the ±180° meridian. - # if source_epsg is not None and _is_geographic_crs(source_epsg): - # _log.info( - # "Bounding box crosses the antimeridian (west=%s > east=%s); " - # "splitting into two polygons.", - # bbox["west"], - # bbox["east"], - # ) - # west_part = shapely.geometry.box(bbox["west"], bbox["south"], 180.0, bbox["north"]) - # east_part = shapely.geometry.box(-180.0, bbox["south"], bbox["east"], bbox["north"]) - # return MultiPolygon([west_part, east_part]), source_epsg - # raise JobSplittingFailure( - # "Invalid bounding box: west must be less than east. " - # "Antimeridian-crossing bounding boxes are only supported " - # "for geographic CRSs (e.g. EPSG:4326) with an explicit 'crs' field." - # ) - return bbox.as_geometry(), source_epsg elif isinstance(geometry, (Polygon, MultiPolygon)): return geometry, None @@ -374,3 +358,178 @@ def split_area( grid = _SizeBasedTileGrid(epsg=normalize_crs(projection), size=tile_size) return grid.get_tiles(aoi) + + +def _predefined_tile_keys(tile_grid: "_PredefinedTileGrid", points: np.ndarray) -> np.ndarray: + """For each point (in the grid CRS), the position of the first tile containing it (-1: none).""" + # TODO: polygons are assigned by their representative point only, so they may extend beyond their tile. + # Assigning them to e.g. the tile with the largest overlap could be a future addition. + n_tiles = len(tile_grid._gdf) + point_idx, tile_idx = tile_grid._gdf.sindex.query(points, predicate="intersects") + # Lowest tile position wins for points on shared tile edges. + keys = np.full(len(points), n_tiles) + np.minimum.at(keys, point_idx, tile_idx) + keys[keys == n_tiles] = -1 + return keys + + +def _bisect_group( + positions: np.ndarray, + *, + xs: np.ndarray, + ys: np.ndarray, + max_geometries: Optional[int], +) -> List[np.ndarray]: + """ + Recursively bisect a group of rows (given as positions) at the median of their locations, + along the longer axis of the locations' spread, until groups have at most *max_geometries* rows. + """ + result = [] + stack = [positions] + while stack: + pos = stack.pop() + if max_geometries is None or len(pos) <= max_geometries: + result.append(pos) + continue + px, py = xs[pos], ys[pos] + axis = px if np.ptp(px) >= np.ptp(py) else py + pos = pos[np.argsort(axis, kind="stable")] + cut = len(pos) // 2 + stack.append(pos[cut:]) + stack.append(pos[:cut]) + return result + + +def group_geometries( + geometries: gpd.GeoDataFrame, + *, + max_geometries: Optional[int] = None, + projection: Optional[str] = None, + tile_size: Optional[float] = None, + tile_grid: Optional[Union[_PredefinedTileGrid, _SizeBasedTileGrid, gpd.GeoDataFrame]] = None, +) -> List[gpd.GeoDataFrame]: + """ + Group geometries (e.g. sample locations for point-based extraction) + into spatially clustered groups of at most *max_geometries* each. + + Each row is treated as a single "feature" to group, regardless of its + geometry type (e.g. ``Point``, ``MultiPoint``, ``Polygon``): + *max_geometries* limits the number of rows (features) per group, + not the number of individual coordinates. To limit individual + coordinates instead, explode multi-part geometries into separate rows + first, e.g. with ``geometries.explode(index_parts=False)``. + + Features are located by their ``representative_point()`` (computed in the input CRS), + they are not split themselves. For a single large area of interest that has to be split, + use :func:`split_area` instead. + + Grouping is done by recursive bisection (k-d tree style): a group is cut at the + median of the feature locations along the longer axis of their spread, + until it has at most *max_geometries* features. + This keeps groups spatially clustered and balanced. + + If *tile_grid* is given (e.g. a Sentinel-2 MGRS/UTM tile grid), each feature is first + assigned to the tile containing its representative point, and every tile is bisected separately. + Features whose representative point is not in any tile are grouped (and bisected) separately. + Alternatively, *tile_size* and *projection* can be given to group + features in a regular grid of square tiles anchored at the origin of the projection + (so grouping does not change when other features are added). + + Current limitations: + + - Features are never split: a polygon crossing tile edges, or larger than a tile, + is assigned to a single tile and extends beyond it. Use :func:`split_area` for large areas. + - Groups are limited by feature count only, not by area or processing cost. + Use a tile grid to limit the spatial extent of groups. + - Multi-part geometries are placed by a single representative point, + so their parts can be far away from the rest of their group. + - With overlapping tiles, a feature is assigned to the first tile (in grid order) + containing its representative point, not to the tile with the largest overlap. + + :param geometries: GeoDataFrame with a ``geometry`` column. + Null or empty geometries are not allowed. + :param max_geometries: maximum number of geometries per group. + ``None`` (default) means no limit. + :param projection: EPSG string (e.g. ``"EPSG:3857"``) of the regular grid + to group by. Must be combined with *tile_size*. + :param tile_size: tile edge length in the unit of *projection*. + Must be combined with *projection*. + :param tile_grid: optional :class:`_PredefinedTileGrid` or + :class:`_SizeBasedTileGrid` instance, or GeoDataFrame of tile polygons + (with CRS set) to co-join geometries by before splitting, e.g. a + satellite tiling grid. + Mutually exclusive with *projection* / *tile_size*. + :return: list of GeoDataFrames, with the original index of *geometries* preserved. + :raises JobSplittingFailure: on invalid input (e.g. non-positive + *max_geometries*, null/empty geometries, contradictory grid arguments). + + Like with :func:`split_area`, geometries without CRS information + are assumed to be in the tile grid's CRS. + Without tile grid, no CRS is needed and the returned groups have no CRS either. + """ + if not isinstance(geometries, gpd.GeoDataFrame): + raise JobSplittingFailure(f"Expected a GeoDataFrame, got {type(geometries).__name__}.") + if "geometry" not in geometries.columns: + raise JobSplittingFailure("The GeoDataFrame must contain a 'geometry' column.") + if max_geometries is not None and max_geometries <= 0: + raise JobSplittingFailure(f"'max_geometries' must be positive, got {max_geometries!r}.") + + if tile_grid is not None: + if projection is not None or tile_size is not None: + raise JobSplittingFailure( + "Cannot combine 'tile_grid' with 'projection' or 'tile_size'. " + "Either pass a tile_grid, or pass projection + tile_size." + ) + if isinstance(tile_grid, gpd.GeoDataFrame): + tile_grid = _PredefinedTileGrid(tiles=tile_grid) + elif not isinstance(tile_grid, (_PredefinedTileGrid, _SizeBasedTileGrid)): + raise JobSplittingFailure( + "Unsupported tile_grid type: " + f"{type(tile_grid).__name__}. Expected _PredefinedTileGrid, _SizeBasedTileGrid, or GeoDataFrame." + ) + elif projection is not None or tile_size is not None: + if projection is None or tile_size is None: + raise JobSplittingFailure("'projection' and 'tile_size' must be provided together.") + tile_grid = _SizeBasedTileGrid.from_size_projection(size=tile_size, projection=projection) + + n_invalid = int((geometries.geometry.isna() | geometries.geometry.is_empty).sum()) + if n_invalid: + raise JobSplittingFailure(f"Found {n_invalid} null or empty geometries, which cannot be grouped.") + + if geometries.crs is None and tile_grid is not None: + _log.warning( + "Geometries have no CRS information; assuming they are already in the tile grid's CRS (EPSG:%d).", + tile_grid.crs, + ) + geometries = geometries.set_crs(f"EPSG:{tile_grid.crs}") + if geometries.crs is not None and geometries.crs.to_epsg() is not None: + _TileGridInterface._check_antimeridian_crossing(geometries.total_bounds, geometries.crs.to_epsg()) + + if geometries.empty: + return [] + + # TODO: limiting groups by area or processing cost (e.g. a weight per feature) is not supported yet. + # Work with positions (not index labels), so duplicated index labels are no problem. + locations = geometries.geometry.representative_point() + if tile_grid is not None: + locations = locations.to_crs(epsg=tile_grid.crs) + xs, ys = locations.x.to_numpy(), locations.y.to_numpy() + + all_positions = np.arange(len(geometries)) + if tile_grid is None: + position_groups = [all_positions] + else: + if isinstance(tile_grid, _SizeBasedTileGrid): + # Regular grid anchored at the origin of the CRS. + keys = pd.DataFrame({"x": np.floor(xs / tile_grid.size), "y": np.floor(ys / tile_grid.size)}) + elif isinstance(tile_grid, _PredefinedTileGrid): + keys = pd.DataFrame({"tile": _predefined_tile_keys(tile_grid, locations.to_numpy())}) + else: + raise JobSplittingFailure(f"Unsupported tile grid type: {type(tile_grid).__name__}.") + position_groups = list(keys.groupby(list(keys.columns), sort=True).indices.values()) + + result = [] + for positions in position_groups: + for group in _bisect_group(positions, xs=xs, ys=ys, max_geometries=max_geometries): + result.append(geometries.iloc[np.sort(group)]) + return result diff --git a/tests/extra/job_management/test_job_splitting.py b/tests/extra/job_management/test_job_splitting.py index f8f1ab227..ca334d733 100644 --- a/tests/extra/job_management/test_job_splitting.py +++ b/tests/extra/job_management/test_job_splitting.py @@ -1,3 +1,5 @@ +import warnings + import geopandas as gpd import pytest import shapely @@ -7,6 +9,7 @@ _PredefinedTileGrid, _SizeBasedTileGrid, _TileGridInterface, + group_geometries, split_area, ) from openeo.util import BBoxDict @@ -275,6 +278,12 @@ def test_returns_geodataframe(self): assert isinstance(result, gpd.GeoDataFrame) assert result.crs.to_epsg() == 4326 + def test_aoi_keyword_argument(self): + aoi = {"west": 0.0, "south": 0.0, "east": 10_000.0, "north": 10_000.0, "crs": "EPSG:3857"} + result = split_area(aoi=aoi, tile_size=10_000.0, projection="EPSG:3857") + assert len(result) == 1 + assert result.crs.to_epsg() == 3857 + def test_with_tile_size_and_projection(self): aoi = {"west": 0.0, "south": 0.0, "east": 1.0, "north": 1.0, "crs": "EPSG:4326"} result = split_area(aoi, projection="EPSG:4326", tile_size=1.0) @@ -386,3 +395,319 @@ def test_antimeridian_not_checked_for_projected_crs(self): aoi = {"west": 18_000_000.0, "south": 0.0, "east": 19_000_000.0, "north": 1_000_000.0, "crs": "EPSG:3857"} result = split_area(aoi, projection="EPSG:3857", tile_size=1_000_000.0) assert len(result) >= 1 + + +class TestGroupGeometries: + def _grid_points(self, n_side: int) -> gpd.GeoDataFrame: + points = [shapely.geometry.Point(x, y) for x in range(n_side) for y in range(n_side)] + return gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") + + def test_single_group_when_under_limit(self): + gdf = self._grid_points(5) # 25 points + result = group_geometries(gdf, max_geometries=100) + assert len(result) == 1 + assert len(result[0]) == 25 + + def test_splits_into_groups_of_max_geometries(self): + gdf = self._grid_points(10) # 100 points + result = group_geometries(gdf, max_geometries=30) + assert sum(len(g) for g in result) == 100 + assert all(len(g) <= 30 for g in result) + assert len(result) > 1 + + def test_groups_are_spatially_clustered(self): + """Points on the left half and right half should end up in different groups.""" + left = [shapely.geometry.Point(x, 0) for x in range(10)] + right = [shapely.geometry.Point(x + 100, 0) for x in range(10)] + gdf = gpd.GeoDataFrame(geometry=left + right, crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=10) + assert len(result) == 2 + for group in result: + xs = [geom.x for geom in group.geometry] + assert max(xs) - min(xs) < 50 + + def test_default_max_geometries_is_unlimited(self): + gdf = self._grid_points(30) # 900 points + result = group_geometries(gdf) + assert len(result) == 1 + assert len(result[0]) == 900 + + def test_default_max_geometries_groups_per_tile(self): + gdf = self._grid_points(30) + result = group_geometries(gdf, tile_size=15, projection="EPSG:4326") + assert sorted(len(g) for g in result) == [225, 225, 225, 225] + + def test_empty_geodataframe(self): + gdf = gpd.GeoDataFrame(geometry=[], crs="EPSG:4326") + assert group_geometries(gdf) == [] + + def test_coincident_points_are_split_by_row_order(self): + """Points that all share the same location are still split in groups of at most max_geometries.""" + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0, 0)] * 25, crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=10) + assert sum(len(g) for g in result) == 25 + assert all(len(g) <= 10 for g in result) + + def test_missing_geometry_column_raises(self): + gdf = gpd.GeoDataFrame({"foo": [1, 2]}) + with pytest.raises(JobSplittingFailure, match="geometry"): + group_geometries(gdf) + + def test_missing_crs_without_grid(self): + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0, 0)]) + result = group_geometries(gdf) + assert len(result) == 1 + assert result[0].crs is None + + def test_non_positive_max_geometries_raises(self): + gdf = self._grid_points(2) + with pytest.raises(JobSplittingFailure, match="max_geometries"): + group_geometries(gdf, max_geometries=0) + + def test_preserves_extra_columns(self): + gdf = gpd.GeoDataFrame( + {"id": [1, 2, 3], "geometry": [shapely.geometry.Point(x, 0) for x in range(3)]}, + crs="EPSG:4326", + ) + result = group_geometries(gdf, max_geometries=10) + assert len(result) == 1 + assert list(result[0]["id"]) == [1, 2, 3] + + def test_multipoint_rows_treated_as_single_feature(self): + """A row with a MultiPoint geometry counts as one feature towards max_geometries, not one per coordinate.""" + gdf = gpd.GeoDataFrame( + geometry=[ + shapely.geometry.MultiPoint([(0, 0), (0.1, 0.1), (0.2, 0.2)]), # 3 coordinates, 1 feature + shapely.geometry.Point(0.5, 0.5), + ], + crs="EPSG:4326", + ) + result = group_geometries(gdf, max_geometries=10) + assert len(result) == 1 + assert len(result[0]) == 2 # 2 rows/features, regardless of coordinate count + + def test_tile_grid_points_never_span_multiple_tiles(self): + tile_grid = gpd.GeoDataFrame( + {"tile": ["A", "B"]}, + geometry=[shapely.geometry.box(0, 0, 1, 1), shapely.geometry.box(1, 0, 2, 1)], + crs="EPSG:4326", + ) + points = [shapely.geometry.Point(0.1 * i, 0.5) for i in range(5)] # tile A + # Start at x=1.1 to avoid the x=1.0 shared-edge ambiguity. + points += [shapely.geometry.Point(1.1 + 0.1 * i, 0.5) for i in range(5)] # tile B + gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") + + result = group_geometries(gdf, max_geometries=100, tile_grid=tile_grid) + assert len(result) == 2 + for group in result: + xs = [geom.x for geom in group.geometry] + assert (max(xs) < 1) or (min(xs) >= 1) + + def test_tile_grid_oversized_tile_is_further_split(self): + tile_grid = gpd.GeoDataFrame(geometry=[shapely.geometry.box(0, 0, 10, 10)], crs="EPSG:4326") + points = [shapely.geometry.Point(x, y) for x in range(10) for y in range(10)] # 100 points, 1 tile + gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") + + result = group_geometries(gdf, max_geometries=30, tile_grid=tile_grid) + assert sum(len(g) for g in result) == 100 + assert all(len(g) <= 30 for g in result) + assert len(result) > 1 + + def test_tile_grid_unmatched_points_grouped_separately(self): + tile_grid = gpd.GeoDataFrame(geometry=[shapely.geometry.box(0, 0, 1, 1)], crs="EPSG:4326") + gdf = gpd.GeoDataFrame( + geometry=[shapely.geometry.Point(0.5, 0.5), shapely.geometry.Point(50, 50)], + crs="EPSG:4326", + ) + result = group_geometries(gdf, max_geometries=10, tile_grid=tile_grid) + assert len(result) == 2 + assert sum(len(g) for g in result) == 2 + + def test_tile_grid_reprojects_when_crs_differs(self): + tile_grid = gpd.GeoDataFrame( + geometry=[shapely.geometry.box(-111_320, -111_325, 111_320, 111_325)], crs="EPSG:3857" + ) + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0.1, 0.1)], crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=10, tile_grid=tile_grid) + assert len(result) == 1 + assert len(result[0]) == 1 + + def test_tile_grid_missing_crs_raises(self): + tile_grid = gpd.GeoDataFrame(geometry=[shapely.geometry.box(0, 0, 1, 1)]) + gdf = self._grid_points(2) + with pytest.raises(JobSplittingFailure, match="CRS"): + group_geometries(gdf, tile_grid=tile_grid) + + @pytest.mark.parametrize( + "geometries", + [ + [shapely.geometry.Point(x, 0) for x in range(5)], + shapely.geometry.MultiPoint([(x, 0) for x in range(5)]), + ], + ) + def test_non_geodataframe_raises(self, geometries): + with pytest.raises(JobSplittingFailure, match="Expected a GeoDataFrame"): + group_geometries(geometries, max_geometries=10) + + def test_points_without_crs_assume_grid_crs(self, caplog): + points = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(100_000 * x, 0) for x in range(3)]) + result = group_geometries(points, max_geometries=10, tile_size=150_000, projection="EPSG:3857") + assert [len(g) for g in result] == [2, 1] + assert all(g.crs.to_epsg() == 3857 for g in result) + assert "assuming they are already in the tile grid's CRS (EPSG:3857)" in caplog.text + + def test_tile_grid_interface_instance(self): + gdf = gpd.GeoDataFrame( + geometry=[shapely.geometry.Point(0.5, 0.5), shapely.geometry.Point(1.5, 0.5)], crs="EPSG:4326" + ) + tile_grid = _PredefinedTileGrid( + tiles=[shapely.geometry.box(0, 0, 1, 1), shapely.geometry.box(1, 0, 2, 1)], crs=4326 + ) + assert len(group_geometries(gdf, max_geometries=10, tile_grid=tile_grid)) == 2 + tile_grid = _SizeBasedTileGrid(epsg=4326, size=1.0) + assert len(group_geometries(gdf, max_geometries=10, tile_grid=tile_grid)) == 2 + + def test_unsupported_tile_grid_interface_raises(self): + class CustomTileGrid(_TileGridInterface): + @property + def crs(self) -> int: + return 4326 + + def get_tiles(self, geometry): + return gpd.GeoDataFrame(geometry=[shapely.geometry.box(0, 0, 1, 1)], crs="EPSG:4326") + + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0.5, 0.5)], crs="EPSG:4326") + with pytest.raises(JobSplittingFailure, match="Unsupported tile_grid type: CustomTileGrid"): + group_geometries(gdf, max_geometries=10, tile_grid=CustomTileGrid()) + + def test_null_geometry_raises(self): + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0, 0), None], crs="EPSG:4326") + with pytest.raises(JobSplittingFailure, match="1 null or empty"): + group_geometries(gdf) + + def test_empty_geometry_raises(self): + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(0, 0), shapely.geometry.Point()], crs="EPSG:4326") + with pytest.raises(JobSplittingFailure, match="1 null or empty"): + group_geometries(gdf) + + def test_duplicated_index_keeps_all_rows_once(self): + gdf = self._grid_points(4) + gdf.index = [0, 0, 1, 1, 2, 2, 3, 3] * 2 + result = group_geometries(gdf, max_geometries=3) + assert sum(len(g) for g in result) == 16 + coords = sorted((g.x, g.y) for group in result for g in group.geometry) + assert coords == sorted((g.x, g.y) for g in gdf.geometry) + assert sorted(i for group in result for i in group.index) == sorted(gdf.index) + + def test_original_index_preserved(self): + gdf = self._grid_points(3) + gdf.index = list("abcdefghi") + result = group_geometries(gdf, max_geometries=4) + assert sorted(i for group in result for i in group.index) == list("abcdefghi") + for group in result: + assert (gdf.loc[group.index].geometry.values == group.geometry.values).all() + + def test_polygons_in_geographic_crs_emit_no_warning(self): + gdf = gpd.GeoDataFrame( + geometry=[shapely.geometry.box(x, 0, x + 1, 1) for x in range(10)], + crs="EPSG:4326", + ) + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + result = group_geometries(gdf, max_geometries=4) + assert sum(len(g) for g in result) == 10 + + def test_concave_polygon_located_by_representative_point(self): + # U shape: its centroid (4.4, 3.7) lies in the notch, outside the polygon. + u_shape = shapely.geometry.box(0, 0, 10, 10).difference(shapely.geometry.box(2, 2, 9, 10)) + assert not u_shape.contains(u_shape.centroid) + gdf = gpd.GeoDataFrame( + geometry=[u_shape, shapely.geometry.Point(1, 6), shapely.geometry.Point(7, 0.5)], crs="EPSG:3857" + ) + result = group_geometries(gdf, tile_size=5, projection="EPSG:3857") + assert sorted(sorted(g.index) for g in result) == [[0, 1], [2]] + + def test_collinear_points_are_split_balanced(self): + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(3, y) for y in range(100)], crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=30) + assert [len(g) for g in result] == [25, 25, 25, 25] + + def test_balanced_bisection(self): + gdf = self._grid_points(10) + result = group_geometries(gdf, max_geometries=30) + assert sorted(len(g) for g in result) == [25, 25, 25, 25] + + def test_size_based_grid_is_globally_anchored(self): + def grouping(gdf): + result = group_geometries(gdf, tile_size=5, projection="EPSG:3857") + return {frozenset(g.index) for g in result} + + base = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(x, 1) for x in (1, 3, 6, 8, 11)], crs="EPSG:3857") + extended = gpd.GeoDataFrame( + geometry=list(base.geometry) + [shapely.geometry.Point(-1000, 1)], crs="EPSG:3857" + ) + assert grouping(extended) - {frozenset({5})} == grouping(base) + + def test_predefined_grid_polygon_assigned_by_representative_point(self): + tile_grid = gpd.GeoDataFrame( + geometry=[shapely.geometry.box(0, 0, 10, 10), shapely.geometry.box(10, 0, 20, 10)], crs="EPSG:3857" + ) + gdf = gpd.GeoDataFrame( + geometry=[ + shapely.geometry.Point(5, 5), + shapely.geometry.Point(15, 5), + shapely.geometry.box(8, 2, 13, 4), # representative point (10.5, 3) lies in tile B + ], + crs="EPSG:3857", + ) + result = group_geometries(gdf, tile_grid=tile_grid) + assert sorted(sorted(g.index) for g in result) == [[0], [1, 2]] + + def test_polygons_grouped_by_centroid(self): + polygons = [shapely.geometry.box(x, 0, x + 1, 1) for x in (0, 1, 100, 101)] + gdf = gpd.GeoDataFrame(geometry=polygons, crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=2) + assert len(result) == 2 + for group in result: + xs = [geom.centroid.x for geom in group.geometry] + assert max(xs) - min(xs) < 50 + + def test_tile_size_projection_groups_by_regular_grid(self): + points = [shapely.geometry.Point(0.1 * i, 0.5) for i in range(5)] + points += [shapely.geometry.Point(1.1 + 0.1 * i, 0.5) for i in range(5)] + gdf = gpd.GeoDataFrame(geometry=points, crs="EPSG:4326") + result = group_geometries(gdf, max_geometries=100, tile_size=1.0, projection="EPSG:4326") + assert len(result) == 2 + for group in result: + xs = [geom.x for geom in group.geometry] + assert (max(xs) < 1) or (min(xs) >= 1) + + def test_tile_size_projection_reprojects(self): + # ~0.1 degree apart near the equator is ~11 km in EPSG:3857: 50 km tiles keep them together + gdf = gpd.GeoDataFrame( + geometry=[shapely.geometry.Point(0.0, 0.0), shapely.geometry.Point(0.1, 0.0)], crs="EPSG:4326" + ) + result = group_geometries(gdf, max_geometries=10, tile_size=50_000, projection="EPSG:3857") + assert len(result) == 1 + result = group_geometries(gdf, max_geometries=10, tile_size=5_000, projection="EPSG:3857") + assert len(result) == 2 + + @pytest.mark.parametrize( + "kwargs", + [ + {"tile_size": 1.0}, + {"projection": "EPSG:4326"}, + {"tile_size": 1.0, "projection": "EPSG:4326", "tile_grid": "grid"}, + ], + ) + def test_invalid_grid_argument_combinations_raise(self, kwargs): + if kwargs.get("tile_grid") == "grid": + kwargs["tile_grid"] = gpd.GeoDataFrame(geometry=[shapely.geometry.box(0, 0, 1, 1)], crs="EPSG:4326") + with pytest.raises(JobSplittingFailure): + group_geometries(self._grid_points(2), **kwargs) + + def test_antimeridian_crossing_points_raises(self): + """Geographic points with longitude outside [-180, 180] are rejected to avoid antimeridian issues.""" + gdf = gpd.GeoDataFrame(geometry=[shapely.geometry.Point(190.0, 0.0)], crs="EPSG:4326") + with pytest.raises(JobSplittingFailure, match="antimeridian"): + group_geometries(gdf, max_geometries=10)