From 551871d94db75ef75f0c512f44e36f5aee92f291 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Tue, 25 Feb 2025 02:55:58 +0000 Subject: [PATCH 01/19] Added Cellular Automata-based Enclosed Tessellation --- momepy/functional/_elements.py | 445 +++++++++++++++++++++++++++++++-- 1 file changed, 431 insertions(+), 14 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 22a17a0a..5598b6ac 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -2,15 +2,19 @@ import geopandas as gpd import libpysal +import math import numpy as np import pandas as pd import shapely +from collections import deque, Counter +from enum import Enum from geopandas import GeoDataFrame, GeoSeries from joblib import Parallel, delayed from libpysal.cg import voronoi_frames from libpysal.graph import Graph from packaging.version import Version from pandas import MultiIndex, Series +from typing import Tuple, List GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev") @@ -129,6 +133,11 @@ def enclosed_tessellation( segment: float = 0.5, threshold: float = 0.05, n_jobs: int = -1, + use_ca: bool = False, + cell_size: float = 1.0, + neighbor_mode: str = "moore", + fulfill: bool = False, + barriers_for_inner: GeoSeries | GeoDataFrame = None, ) -> GeoDataFrame: """Generate enclosed tessellation @@ -265,7 +274,19 @@ def enclosed_tessellation( # generate tessellation in parallel new = Parallel(n_jobs=n_jobs)( - delayed(_tess)(*t, threshold, shrink, segment, index_name) for t in tuples + delayed(_tess)( + *t, + threshold, + shrink, + segment, + index_name, + use_ca, + cell_size, + neighbor_mode, + fulfill, + barriers_for_inner, + ) + for t in tuples ) new_df = pd.concat(new, axis=0) @@ -297,24 +318,57 @@ def enclosed_tessellation( return pd.concat([new_df, singles.drop(columns="position"), clean_blocks]) -def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id): +def _tess( + ix, + poly, + blg, + threshold, + shrink, + segment, + enclosure_id, + use_ca, + cell_size, + neighbor_mode, + fulfill, + barriers_for_inner, +): """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" # check if threshold is set and filter buildings based on the threshold if threshold: - blg = blg[ - shapely.area(shapely.intersection(blg.geometry.array, poly)) - > (shapely.area(blg.geometry.array) * threshold) - ] + if isinstance(poly, shapely.geometry.Polygon): + blg = blg[ + shapely.area(shapely.intersection(blg.geometry.array, poly)) + > (shapely.area(blg.geometry.array) * threshold) + ] + else: + blg = blg[ + shapely.area( + shapely.intersection( + blg.geometry.array, shapely.geometry.Polygon(poly.boundary) + ) + ) + > (shapely.area(blg.geometry.array) * threshold) + ] if len(blg) > 1: - tess = voronoi_frames( - blg, - clip=poly, - shrink=shrink, - segment=segment, - return_input=False, - as_gdf=True, - ) + if not use_ca: + tess = voronoi_frames( + blg, + clip=poly, + shrink=shrink, + segment=segment, + return_input=False, + as_gdf=True, + ) + else: + tess = _voronoi_by_ca( + seed_geoms=blg, + barrier_geoms=poly, + cell_size=cell_size, + neighbor_mode=neighbor_mode, + fulfill=fulfill, + barriers_for_inner=barriers_for_inner, + ) tess[enclosure_id] = ix return tess @@ -332,6 +386,369 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id): ) +def _voronoi_by_ca( + seed_geoms: GeoSeries | GeoDataFrame, + barrier_geoms: GeoSeries | GeoDataFrame, + cell_size: float = 1.0, + neighbor_mode: str = "moore", + fulfill: bool = True, + barriers_for_inner: GeoSeries | GeoDataFrame = None, +) -> GeoDataFrame: + """ + Generate an aggregated Voronoi tessellation as a GeoDataFrame via a cellular automata. + + This unified function performs the following: + - Ensures that the CRS of seed and barrier geometries are aligned. + - Combines inner barriers with the enclosure, . + - Computes grid bounds covering both seed and barrier geometries. + - Marks barrier cells using a prepared geometry for fast intersection. + - Seeds the grid with seed geometries. + - Propagates seed values via a BFS expansion that respects barriers. + - Uses a voting mechanism to finalize boundary cell assignments. + - Converts grid cells into a GeoDataFrame and dissolves adjacent cells + with the same seed id. + - Clips the output cells by the barrier boundaries. + + Parameters: + seed_geoms: GeoDataFrame containing seed features. + barrier_geoms: GeoDataFrame containing barrier features or a shapely Polygon. + cell_size: Grid cell size. By default it is 1.0. + neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. + fulfill: Whether to assign adjacent seed cells to the same seed id. By default it is True. + barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None. + + Returns: + A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by barriers. + """ + + # If there is barriers_for_inner, add the intersected or contained barriers to the barrier_geoms + if barriers_for_inner is not None: + # get inner barriers + inner_barriers = _get_inner_barriers( + barrier_geoms, barriers_for_inner.to_crs(seed_geoms.crs) + ) + + if barrier_geoms.geom_type == "Polygon": + # Wrap a single barrier geometry (Polygon, LineString, or GeometryCollection) + barrier_geoms = GeoSeries([barrier_geoms.boundary], crs=seed_geoms.crs) + + if len(inner_barriers) > 0: + # add inner barriers as a list of geometry to the barrier_geoms + barrier_geoms = pd.concat( + [barrier_geoms, inner_barriers], ignore_index=True + ) + else: + if barrier_geoms.geom_type == "Polygon": + # Wrap a single barrier geometry (Polygon, LineString, or GeometryCollection) + barrier_geoms = GeoSeries([barrier_geoms.boundary], crs=seed_geoms.crs) + + # Compute grid bounds + origin, grid_width, grid_height = _get_grid_bounds( + seed_geoms, barrier_geoms, cell_size + ) + + # Prepare barrier geometries if provided. + barrier_union = None + prep_barrier = None + if not barrier_geoms.empty: + barrier_union = shapely.ops.unary_union(barrier_geoms.geometry) + prep_barrier = shapely.prepared.prep(barrier_union) + + # Initialize grid states with UNKNOWN values. + states = np.full((grid_height, grid_width), CellState.UNKNOWN.value, dtype=int) + + # Identify barrier cells in the grid + if prep_barrier is not None: + xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) + cell_polys = GeoSeries( + [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xs.flatten(), ys.flatten()) + ] + ) + barrier_mask = cell_polys.intersects(barrier_union).values.reshape( + (grid_height, grid_width) + ) + states[barrier_mask] = CellState.BARRIER.value + + # Seed the grid with seed geometries. + for site_id, geom in enumerate(seed_geoms.geometry): + if not geom.is_empty: + cells = _geom_to_cells(geom, origin, cell_size, grid_width, grid_height) + valid_cells = [ + (x, y) for x, y in cells if states[y, x] == CellState.UNKNOWN.value + ] + if valid_cells: + indices = np.array(valid_cells) + states[indices[:, 1], indices[:, 0]] = site_id + + # Initialize the BFS queue with all seeded cells’ neighbors. + queue = deque() + seed_indices = np.argwhere(states >= 0) + for y, x in seed_indices: + _enqueue_neighbors(x, y, states, grid_width, grid_height, neighbor_mode, queue) + + # Process BFS to propagate seed values. + while queue: + # Dequeue the current cell and skip if it is not a frontier cell. + x_current, y_current = queue.popleft() + if states[y_current, x_current] != CellState.FRONTIER.value: + continue + + # Get neighbor cells that were already assigned a seed id or still unknown (state >= 0). + # Note that boundary or barrier cells are skipped (state < 0). + neighbor_seeds = [ + states[ny, nx] + for nx, ny in _get_neighbors( + x_current, y_current, grid_width, grid_height, mode=neighbor_mode + ) + if states[ny, nx] >= 0 + ] + if not neighbor_seeds: + continue + + # Assign as a boundary if multiple seed ids are found. + if len(set(neighbor_seeds)) > 1: + states[y_current, x_current] = CellState.BOUNDARY.value + # EIf not, equeue neighbor cells for further propagation. + else: + assigned_seed = set(neighbor_seeds).pop() + states[y_current, x_current] = assigned_seed + _enqueue_neighbors( + x_current, + y_current, + states, + grid_width, + grid_height, + neighbor_mode, + queue, + ) + + # Post-process barrier and boundary cells using a voting mechanism. + if fulfill: + states = _assign_adjacent_seed_cells(states, neighbor_mode) + + # Create grid cell polygons and build a GeoDataFrame. + xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) + grid_polys = [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xs.flatten(), ys.flatten()) + ] + grid_gdf = GeoDataFrame( + {"site_id": states.flatten()}, geometry=grid_polys, crs=seed_geoms.crs + ) + + # Include only cells with valid seed assignments and dissolve contiguous regions. + grid_gdf = ( + grid_gdf[grid_gdf["site_id"] >= 0] + .dissolve(by="site_id") + .reset_index() + .drop("site_id", axis=1) + ) + + # Clip by barriers + if barrier_geoms is not None and (not barrier_geoms.empty): + # Create a union of the barrier geometries. + barrier_union = shapely.ops.unary_union(barrier_geoms.geometry) + # If the barrier union is not a polygon (e.g., it's a MultiLineString), polygonize it. + if not isinstance( + barrier_union, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon) + ): + barrier_polys = list(shapely.ops.polygonize(barrier_union)) + if barrier_polys: + barrier_union = shapely.ops.unary_union(barrier_polys) + # Clip each polygon in the grid using the barrier boundary. + grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) + + return grid_gdf + + +def _get_inner_barriers(enclosure, barriers): + """ + Get inner barriers that intersect or are contained within an enclosure. + + Args: + enclosure (GeoSeries or GeoDataFrame): The enclosure geometry. + barriers (GeoDataFrame): The barriers GeoDataFrame. + + Returns: + shapely.geometry.Polygon: A single Polygon combining the enclosure + and any intersecting barriers. + """ + # Find barriers intersecting or contained in the enclosure + inner_barriers = barriers[ + barriers.intersects(enclosure) | barriers.contains(enclosure) + ] + + # Clip those segments to stay within the enclosure + inner_barriers = gpd.clip(inner_barriers, enclosure) + + return GeoSeries(inner_barriers.geometry, crs=barriers.crs) + + +class CellState(Enum): + """ + Enumeration of cell states for grid processing for improving the readability, instead of integers. + + Attributes: + UNKNOWN: Cell has not been processed. + BOUNDARY: Cell is at a junction between different seed regions. + BARRIER: Cell originally designated as a barrier. + FRONTIER: Cell queued for BFS expansion. + """ + + UNKNOWN = -1 + BOUNDARY = -2 + BARRIER = -3 + FRONTIER = -4 + + +def _get_cell_polygon( + x_idx: int, y_idx: int, cell_size: float = 1.0, origin: Tuple[float, float] = (0, 0) +) -> shapely.geometry.Polygon: + """ + Generate a grid cell polygon based on the given indices, cell size, and origin. + """ + ox, oy = origin + return shapely.geometry.Polygon( + [ + (ox + x_idx * cell_size, oy + y_idx * cell_size), + (ox + (x_idx + 1) * cell_size, oy + y_idx * cell_size), + (ox + (x_idx + 1) * cell_size, oy + (y_idx + 1) * cell_size), + (ox + x_idx * cell_size, oy + (y_idx + 1) * cell_size), + ] + ) + + +def _get_neighbors( + x: int, y: int, max_x: int, max_y: int, mode: str = "moore" +) -> List[Tuple[int, int]]: + """ + Retrieve valid neighboring cell indices based on connectivity. + + Parameters: + x, y: The current cell indices. + max_x, max_y: Dimensions of the grid. + mode: "moore" for 8-connected or "neumann" for 4-connected neighbors. + + Returns: + A list of (x, y) tuples for valid neighbor indices. + """ + neighbor_dirs = { + "moore": [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)], + "neumann": [(0, -1), (-1, 0), (1, 0), (0, 1)], + } + directions = neighbor_dirs.get(mode) + if directions is None: + raise ValueError("Invalid neighbor_mode: choose 'moore' or 'neumann'") + return [ + (x + dx, y + dy) + for dx, dy in directions + if 0 <= x + dx < max_x and 0 <= y + dy < max_y + ] + + +def _get_grid_bounds( + seed_geoms: GeoSeries | GeoDataFrame, + barrier_geoms: GeoSeries | GeoDataFrame, + cell_size: float, +) -> Tuple[Tuple[float, float], int, int]: + """ + Compute the grid bounds required to cover both seed and barrier geometries. + """ + seed_bounds = seed_geoms.total_bounds # [xmin, ymin, xmax, ymax] + barrier_bounds = barrier_geoms.total_bounds + + xmin = min(seed_bounds[0], barrier_bounds[0]) + ymin = min(seed_bounds[1], barrier_bounds[1]) + xmax = max(seed_bounds[2], barrier_bounds[2]) + ymax = max(seed_bounds[3], barrier_bounds[3]) + grid_width = math.ceil((xmax - xmin) / cell_size) + grid_height = math.ceil((ymax - ymin) / cell_size) + return (xmin, ymin), grid_width, grid_height + + +def _geom_to_cells( + geom: shapely.geometry, + origin: Tuple[float, float], + cell_size: float, + grid_width: int, + grid_height: int, +) -> List[Tuple[int, int]]: + """ + Determine grid cell indices that intersect the given geometry. + """ + if isinstance(geom, shapely.geometry.Point): + sx = int((geom.x - origin[0]) // cell_size) + sy = int((geom.y - origin[1]) // cell_size) + return [(sx, sy)] if 0 <= sx < grid_width and 0 <= sy < grid_height else [] + + else: + minx, miny, maxx, maxy = geom.bounds + start_x = max(0, int((minx - origin[0]) // cell_size)) + start_y = max(0, int((miny - origin[1]) // cell_size)) + end_x = min(grid_width, int(math.ceil((maxx - origin[0]) / cell_size))) + end_y = min(grid_height, int(math.ceil((maxy - origin[1]) / cell_size))) + + x_range = np.arange(start_x, end_x) + y_range = np.arange(start_y, end_y) + xx, yy = np.meshgrid(x_range, y_range) + candidate_polys = GeoSeries( + [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xx.flatten(), yy.flatten()) + ] + ) + mask = candidate_polys.intersects(geom) + return list(zip(xx.flatten()[mask], yy.flatten()[mask])) + + +def _enqueue_neighbors( + x: int, + y: int, + states: np.ndarray, + grid_width: int, + grid_height: int, + neighbor_mode: str, + queue: deque, +) -> None: + """ + Enqueue valid neighboring cells for BFS expansion. + """ + for nx, ny in _get_neighbors(x, y, grid_width, grid_height, mode=neighbor_mode): + if states[ny, nx] == CellState.UNKNOWN.value: + states[ny, nx] = CellState.FRONTIER.value + queue.append((nx, ny)) + + +def _assign_adjacent_seed_cells( + states: np.ndarray, neighbor_mode: str = "moore" +) -> np.ndarray: + """ + Reassign border and barrier cells to the proximate seed areas using a voting mechanism. + """ + new_states = states.copy() + indices = np.argwhere( + np.isin(states, [CellState.BARRIER.value, CellState.BOUNDARY.value]) + ) + grid_height, grid_width = states.shape + + for y, x in indices: + neighbor_seeds = [ + states[ny, nx] + for nx, ny in _get_neighbors( + x, y, grid_width, grid_height, mode=neighbor_mode + ) + if states[ny, nx] >= 0 + ] + if neighbor_seeds: + cnt = Counter(neighbor_seeds) + # In case of ties, choose the smaller seed id. + chosen_seed = min(cnt.items(), key=lambda item: (-item[1], item[0]))[0] + new_states[y, x] = chosen_seed + return new_states + + def verify_tessellation(tessellation, geometry): """Check whether result matches buildings and contains only Polygons. From 5b4a06ee51f1516e7c06ef31cc4979513563ab80 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Tue, 4 Mar 2025 21:56:54 +0000 Subject: [PATCH 02/19] delete fulfill and make it as a default behaviour --- momepy/functional/_elements.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 5598b6ac..7d7894f4 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -136,7 +136,6 @@ def enclosed_tessellation( use_ca: bool = False, cell_size: float = 1.0, neighbor_mode: str = "moore", - fulfill: bool = False, barriers_for_inner: GeoSeries | GeoDataFrame = None, ) -> GeoDataFrame: """Generate enclosed tessellation @@ -283,7 +282,6 @@ def enclosed_tessellation( use_ca, cell_size, neighbor_mode, - fulfill, barriers_for_inner, ) for t in tuples @@ -329,7 +327,6 @@ def _tess( use_ca, cell_size, neighbor_mode, - fulfill, barriers_for_inner, ): """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" @@ -366,7 +363,6 @@ def _tess( barrier_geoms=poly, cell_size=cell_size, neighbor_mode=neighbor_mode, - fulfill=fulfill, barriers_for_inner=barriers_for_inner, ) tess[enclosure_id] = ix @@ -391,7 +387,6 @@ def _voronoi_by_ca( barrier_geoms: GeoSeries | GeoDataFrame, cell_size: float = 1.0, neighbor_mode: str = "moore", - fulfill: bool = True, barriers_for_inner: GeoSeries | GeoDataFrame = None, ) -> GeoDataFrame: """ @@ -414,7 +409,6 @@ def _voronoi_by_ca( barrier_geoms: GeoDataFrame containing barrier features or a shapely Polygon. cell_size: Grid cell size. By default it is 1.0. neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. - fulfill: Whether to assign adjacent seed cells to the same seed id. By default it is True. barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None. Returns: @@ -525,8 +519,7 @@ def _voronoi_by_ca( ) # Post-process barrier and boundary cells using a voting mechanism. - if fulfill: - states = _assign_adjacent_seed_cells(states, neighbor_mode) + states = _assign_adjacent_seed_cells(states, neighbor_mode) # Create grid cell polygons and build a GeoDataFrame. xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) From 0b3a5197d65c70d0b9ba0802e4d10166cd42ae3c Mon Sep 17 00:00:00 2001 From: yu-ta-sato Date: Wed, 19 Mar 2025 14:14:35 +0000 Subject: [PATCH 03/19] Add coverage_simplify with modification of coverage in _voronoi_by_ca --- momepy/functional/_elements.py | 119 +++++++++++++++++++++++---------- 1 file changed, 83 insertions(+), 36 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 7d7894f4..7274c05e 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -19,6 +19,7 @@ GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev") LPS_GE_411 = Version(libpysal.__version__) >= Version("4.11.dev") +SPL_GE_210 = Version(shapely.__version__) >= Version("2.1.0rc1") __all__ = [ "morphological_tessellation", @@ -137,6 +138,7 @@ def enclosed_tessellation( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, + cell_tolerance: float = 2.0 ) -> GeoDataFrame: """Generate enclosed tessellation @@ -283,6 +285,7 @@ def enclosed_tessellation( cell_size, neighbor_mode, barriers_for_inner, + cell_tolerance ) for t in tuples ) @@ -328,6 +331,7 @@ def _tess( cell_size, neighbor_mode, barriers_for_inner, + cell_tolerance ): """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" # check if threshold is set and filter buildings based on the threshold @@ -364,6 +368,7 @@ def _tess( cell_size=cell_size, neighbor_mode=neighbor_mode, barriers_for_inner=barriers_for_inner, + cell_tolerance=cell_tolerance * cell_size, ) tess[enclosure_id] = ix return tess @@ -388,6 +393,7 @@ def _voronoi_by_ca( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, + cell_tolerance: float = 2.0, ) -> GeoDataFrame: """ Generate an aggregated Voronoi tessellation as a GeoDataFrame via a cellular automata. @@ -409,7 +415,9 @@ def _voronoi_by_ca( barrier_geoms: GeoDataFrame containing barrier features or a shapely Polygon. cell_size: Grid cell size. By default it is 1.0. neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. - barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None. + barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None + cell_tolerance: Tolerance for simplifying the grid cells. By default it is 2.0. + Returns: A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by barriers. @@ -422,48 +430,71 @@ def _voronoi_by_ca( barrier_geoms, barriers_for_inner.to_crs(seed_geoms.crs) ) - if barrier_geoms.geom_type == "Polygon": - # Wrap a single barrier geometry (Polygon, LineString, or GeometryCollection) - barrier_geoms = GeoSeries([barrier_geoms.boundary], crs=seed_geoms.crs) + # Handle barrier_geoms if it is a Polygon or MultiPolygon + if barrier_geoms.geom_type == "Polygon": + # Take buffer of polygon and extract its exterior boundary + barrier_geoms_buffered = GeoSeries( + [barrier_geoms.buffer(10).exterior], crs=seed_geoms.crs + ) + barrier_geoms = GeoSeries([barrier_geoms], crs=seed_geoms.crs) + + elif barrier_geoms.geom_type == "MultiPolygon": + # Process each polygon: take buffer then exterior boundary (to ensure there's no gap between enclosures) + barrier_geoms_buffered = GeoSeries( + [poly.buffer(10).exterior for poly in barrier_geoms.geoms], + crs=seed_geoms.crs, + ) + barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) - if len(inner_barriers) > 0: - # add inner barriers as a list of geometry to the barrier_geoms - barrier_geoms = pd.concat( - [barrier_geoms, inner_barriers], ignore_index=True - ) else: - if barrier_geoms.geom_type == "Polygon": - # Wrap a single barrier geometry (Polygon, LineString, or GeometryCollection) - barrier_geoms = GeoSeries([barrier_geoms.boundary], crs=seed_geoms.crs) + raise ValueError("barrier_geoms must be a Polygon or MultiPolygon") + + outer_union = shapely.ops.unary_union(barrier_geoms_buffered) + + # Compute inner barriers union if available + if ( + "inner_barriers" in locals() + and inner_barriers is not None + and not inner_barriers.empty + ): + inner_union = shapely.ops.unary_union(inner_barriers.geometry) + else: + inner_union = None + + # Combine outer barrier with inner barriers + if outer_union and inner_union: + prep_barrier = shapely.ops.unary_union([outer_union, inner_union]) + elif outer_union: + prep_barrier = outer_union + elif inner_union: + prep_barrier = inner_union + else: + prep_barrier = None # Compute grid bounds origin, grid_width, grid_height = _get_grid_bounds( - seed_geoms, barrier_geoms, cell_size + seed_geoms, barrier_geoms_buffered, cell_size ) - # Prepare barrier geometries if provided. - barrier_union = None - prep_barrier = None - if not barrier_geoms.empty: - barrier_union = shapely.ops.unary_union(barrier_geoms.geometry) - prep_barrier = shapely.prepared.prep(barrier_union) - # Initialize grid states with UNKNOWN values. states = np.full((grid_height, grid_width), CellState.UNKNOWN.value, dtype=int) + xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) + cell_polys = GeoSeries( + [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xs.flatten(), ys.flatten()) + ] + ) + # Identify barrier cells in the grid if prep_barrier is not None: - xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) - cell_polys = GeoSeries( - [ - _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xs.flatten(), ys.flatten()) - ] + barrier_mask = cell_polys.intersects(prep_barrier).values.reshape( + grid_height, grid_width ) - barrier_mask = cell_polys.intersects(barrier_union).values.reshape( - (grid_height, grid_width) - ) - states[barrier_mask] = CellState.BARRIER.value + else: + barrier_mask = np.zeros((grid_height, grid_width), dtype=bool) + states[barrier_mask] = CellState.BARRIER.value # Seed the grid with seed geometries. for site_id, geom in enumerate(seed_geoms.geometry): @@ -553,6 +584,17 @@ def _voronoi_by_ca( # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) + if cell_tolerance is not None: + if SPL_GE_210: + # Simplify coverages with the parameter cell_tolerance as the area threshold of Visvalingam-Whyatt algorithm. + grid_gdf["geometry"] = shapely.coverage_simplify( + grid_gdf["geometry"].array, tolerance=cell_tolerance + ) + else: + warnings.warn( + "Shapely 2.1.0 or higher is required for coverage_simplify. Skipping." + ) + return grid_gdf @@ -576,6 +618,9 @@ def _get_inner_barriers(enclosure, barriers): # Clip those segments to stay within the enclosure inner_barriers = gpd.clip(inner_barriers, enclosure) + # Only keep the geometry which is within the enclosure + inner_barriers = inner_barriers[inner_barriers.within(enclosure)] + return GeoSeries(inner_barriers.geometry, crs=barriers.crs) @@ -646,9 +691,6 @@ def _get_grid_bounds( barrier_geoms: GeoSeries | GeoDataFrame, cell_size: float, ) -> Tuple[Tuple[float, float], int, int]: - """ - Compute the grid bounds required to cover both seed and barrier geometries. - """ seed_bounds = seed_geoms.total_bounds # [xmin, ymin, xmax, ymax] barrier_bounds = barrier_geoms.total_bounds @@ -656,9 +698,14 @@ def _get_grid_bounds( ymin = min(seed_bounds[1], barrier_bounds[1]) xmax = max(seed_bounds[2], barrier_bounds[2]) ymax = max(seed_bounds[3], barrier_bounds[3]) - grid_width = math.ceil((xmax - xmin) / cell_size) - grid_height = math.ceil((ymax - ymin) / cell_size) - return (xmin, ymin), grid_width, grid_height + # expand bounds by 1 cell in each direction + new_xmin = xmin - cell_size + new_ymin = ymin - cell_size + new_xmax = xmax + cell_size + new_ymax = ymax + cell_size + grid_width = math.ceil((new_xmax - new_xmin) / cell_size) + grid_height = math.ceil((new_ymax - new_ymin) / cell_size) + return (new_xmin, new_ymin), grid_width, grid_height def _geom_to_cells( From d754ce4743dec95580fcd329e71c9a8a3027e5c5 Mon Sep 17 00:00:00 2001 From: yu-ta-sato Date: Wed, 19 Mar 2025 15:38:36 +0000 Subject: [PATCH 04/19] Derive torelance for coverage_simplify from the cell_size and let shapely < 2.1.0 raise an error --- momepy/functional/_elements.py | 27 +++++++++++---------------- 1 file changed, 11 insertions(+), 16 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 7274c05e..932ab690 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -138,7 +138,6 @@ def enclosed_tessellation( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, - cell_tolerance: float = 2.0 ) -> GeoDataFrame: """Generate enclosed tessellation @@ -285,7 +284,6 @@ def enclosed_tessellation( cell_size, neighbor_mode, barriers_for_inner, - cell_tolerance ) for t in tuples ) @@ -331,7 +329,6 @@ def _tess( cell_size, neighbor_mode, barriers_for_inner, - cell_tolerance ): """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" # check if threshold is set and filter buildings based on the threshold @@ -368,7 +365,6 @@ def _tess( cell_size=cell_size, neighbor_mode=neighbor_mode, barriers_for_inner=barriers_for_inner, - cell_tolerance=cell_tolerance * cell_size, ) tess[enclosure_id] = ix return tess @@ -393,7 +389,6 @@ def _voronoi_by_ca( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, - cell_tolerance: float = 2.0, ) -> GeoDataFrame: """ Generate an aggregated Voronoi tessellation as a GeoDataFrame via a cellular automata. @@ -416,7 +411,6 @@ def _voronoi_by_ca( cell_size: Grid cell size. By default it is 1.0. neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None - cell_tolerance: Tolerance for simplifying the grid cells. By default it is 2.0. Returns: @@ -584,16 +578,17 @@ def _voronoi_by_ca( # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) - if cell_tolerance is not None: - if SPL_GE_210: - # Simplify coverages with the parameter cell_tolerance as the area threshold of Visvalingam-Whyatt algorithm. - grid_gdf["geometry"] = shapely.coverage_simplify( - grid_gdf["geometry"].array, tolerance=cell_tolerance - ) - else: - warnings.warn( - "Shapely 2.1.0 or higher is required for coverage_simplify. Skipping." - ) + if SPL_GE_210: + # Simplify coverages with coverage_simplify. + # torelance set as the square root of the isosceles right triangle with 2 cells_size edges. + # For the behavior of coverage_simplify see https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html + grid_gdf["geometry"] = shapely.coverage_simplify( + grid_gdf["geometry"].array, tolerance= ((2 * cell_size) ** 2 / 2) ** 0.5 + ) + else: + raise ImportError( + "Shapely 2.1.0 or higher is required for coverage_simplify." + ) return grid_gdf From 37cd109a53fdada91b36069772d05010c0f88ee6 Mon Sep 17 00:00:00 2001 From: Yuta Sato <65852668+yu-ta-sato@users.noreply.github.com> Date: Wed, 26 Mar 2025 00:03:31 +0000 Subject: [PATCH 05/19] Apply suggestions from code review Co-authored-by: Martin Fleischmann --- momepy/functional/_elements.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 932ab690..7d6e7ac8 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -134,10 +134,8 @@ def enclosed_tessellation( segment: float = 0.5, threshold: float = 0.05, n_jobs: int = -1, - use_ca: bool = False, - cell_size: float = 1.0, - neighbor_mode: str = "moore", - barriers_for_inner: GeoSeries | GeoDataFrame = None, + inner_barriers: GeoSeries | GeoDataFrame = None, + **kwargs ) -> GeoDataFrame: """Generate enclosed tessellation @@ -435,7 +433,7 @@ def _voronoi_by_ca( elif barrier_geoms.geom_type == "MultiPolygon": # Process each polygon: take buffer then exterior boundary (to ensure there's no gap between enclosures) barrier_geoms_buffered = GeoSeries( - [poly.buffer(10).exterior for poly in barrier_geoms.geoms], + shapely.buffer(shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10) crs=seed_geoms.crs, ) barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) @@ -443,7 +441,7 @@ def _voronoi_by_ca( else: raise ValueError("barrier_geoms must be a Polygon or MultiPolygon") - outer_union = shapely.ops.unary_union(barrier_geoms_buffered) + outer_union = barrier_geoms_buffered.union_all() # Compute inner barriers union if available if ( @@ -451,13 +449,13 @@ def _voronoi_by_ca( and inner_barriers is not None and not inner_barriers.empty ): - inner_union = shapely.ops.unary_union(inner_barriers.geometry) + inner_union = inner_barriers.union_all() else: inner_union = None # Combine outer barrier with inner barriers if outer_union and inner_union: - prep_barrier = shapely.ops.unary_union([outer_union, inner_union]) + prep_barrier = shapely.union(outer_union, inner_union) elif outer_union: prep_barrier = outer_union elif inner_union: @@ -567,14 +565,14 @@ def _voronoi_by_ca( # Clip by barriers if barrier_geoms is not None and (not barrier_geoms.empty): # Create a union of the barrier geometries. - barrier_union = shapely.ops.unary_union(barrier_geoms.geometry) + barrier_union = barrier_geoms.union_all() # If the barrier union is not a polygon (e.g., it's a MultiLineString), polygonize it. if not isinstance( barrier_union, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon) ): - barrier_polys = list(shapely.ops.polygonize(barrier_union)) - if barrier_polys: - barrier_union = shapely.ops.unary_union(barrier_polys) + barrier_polys = shapely.get_parts(shapely.polygonize(barrier_union)) + if not barrier_polys.empty: + barrier_union = shapely.union_all(barrier_polys) # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) @@ -643,7 +641,7 @@ def _get_cell_polygon( Generate a grid cell polygon based on the given indices, cell size, and origin. """ ox, oy = origin - return shapely.geometry.Polygon( + return shapely.Polygon( [ (ox + x_idx * cell_size, oy + y_idx * cell_size), (ox + (x_idx + 1) * cell_size, oy + y_idx * cell_size), From 3cab43b861c778866f6bb134edbb41c084511ae9 Mon Sep 17 00:00:00 2001 From: yu-ta-sato Date: Wed, 26 Mar 2025 01:05:51 +0000 Subject: [PATCH 06/19] Modified API for enclosed_tessellation with kwargs and fixed bugs --- momepy/functional/_elements.py | 91 +++++++++++++++------------------- 1 file changed, 39 insertions(+), 52 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 7d6e7ac8..5c8d6219 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -135,7 +135,7 @@ def enclosed_tessellation( threshold: float = 0.05, n_jobs: int = -1, inner_barriers: GeoSeries | GeoDataFrame = None, - **kwargs + **kwargs, ) -> GeoDataFrame: """Generate enclosed tessellation @@ -185,6 +185,10 @@ def enclosed_tessellation( n_jobs : int, optional The number of jobs to run in parallel. -1 means using all available cores. By default -1 + inner_barriers: GeoSeries | GeoDataFrame, optional + Barriers that should be included in the tessellation process. By default None. + **kwargs + Additional keyword arguments passed to the Cellular Automata-based Voronoi Diagram. Warnings -------- @@ -273,15 +277,7 @@ def enclosed_tessellation( # generate tessellation in parallel new = Parallel(n_jobs=n_jobs)( delayed(_tess)( - *t, - threshold, - shrink, - segment, - index_name, - use_ca, - cell_size, - neighbor_mode, - barriers_for_inner, + *t, threshold, shrink, segment, index_name, inner_barriers, **kwargs ) for t in tuples ) @@ -316,38 +312,24 @@ def enclosed_tessellation( def _tess( - ix, - poly, - blg, - threshold, - shrink, - segment, - enclosure_id, - use_ca, - cell_size, - neighbor_mode, - barriers_for_inner, + ix, poly, blg, threshold, shrink, segment, enclosure_id, inner_barriers, **kwargs ): """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" + + # Extract CA parameters from kwargs + cell_size = kwargs.pop("cell_size", 1.0) + neighbor_mode = kwargs.pop("neighbor_mode", "moore") + simplify = kwargs.pop("simplify", True) + # check if threshold is set and filter buildings based on the threshold if threshold: - if isinstance(poly, shapely.geometry.Polygon): - blg = blg[ - shapely.area(shapely.intersection(blg.geometry.array, poly)) - > (shapely.area(blg.geometry.array) * threshold) - ] - else: - blg = blg[ - shapely.area( - shapely.intersection( - blg.geometry.array, shapely.geometry.Polygon(poly.boundary) - ) - ) - > (shapely.area(blg.geometry.array) * threshold) - ] + blg = blg[ + shapely.area(shapely.intersection(blg.geometry.array, poly)) + > (shapely.area(blg.geometry.array) * threshold) + ] if len(blg) > 1: - if not use_ca: + if inner_barriers is None or inner_barriers.empty: tess = voronoi_frames( blg, clip=poly, @@ -362,7 +344,8 @@ def _tess( barrier_geoms=poly, cell_size=cell_size, neighbor_mode=neighbor_mode, - barriers_for_inner=barriers_for_inner, + barriers_for_inner=inner_barriers, + simplify=simplify, ) tess[enclosure_id] = ix return tess @@ -387,6 +370,7 @@ def _voronoi_by_ca( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, + simplify: bool = True, ) -> GeoDataFrame: """ Generate an aggregated Voronoi tessellation as a GeoDataFrame via a cellular automata. @@ -409,31 +393,38 @@ def _voronoi_by_ca( cell_size: Grid cell size. By default it is 1.0. neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None + simplify: Whether to simplify the coverage using shapely's coverage_simplify. By default it is True. Returns: A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by barriers. """ - # If there is barriers_for_inner, add the intersected or contained barriers to the barrier_geoms - if barriers_for_inner is not None: - # get inner barriers - inner_barriers = _get_inner_barriers( - barrier_geoms, barriers_for_inner.to_crs(seed_geoms.crs) + # Check if shapely version is 2.1.0 or higher for coverage_simplify + if not SPL_GE_210: + raise ImportError( + "Shapely 2.1.0 or higher is required for tessellation with inner barriers provided." ) + # Get inner barriers as intersection or containment of the barrier_geoms + inner_barriers = _get_inner_barriers( + barrier_geoms, barriers_for_inner.to_crs(seed_geoms.crs) + ) + # Handle barrier_geoms if it is a Polygon or MultiPolygon if barrier_geoms.geom_type == "Polygon": - # Take buffer of polygon and extract its exterior boundary + # Take buffer of polygon and extract its exterior boundary (10 cells) barrier_geoms_buffered = GeoSeries( - [barrier_geoms.buffer(10).exterior], crs=seed_geoms.crs + [barrier_geoms.buffer(10 * cell_size).exterior], crs=seed_geoms.crs ) barrier_geoms = GeoSeries([barrier_geoms], crs=seed_geoms.crs) elif barrier_geoms.geom_type == "MultiPolygon": # Process each polygon: take buffer then exterior boundary (to ensure there's no gap between enclosures) barrier_geoms_buffered = GeoSeries( - shapely.buffer(shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10) + shapely.buffer( + shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10 * cell_size + ), crs=seed_geoms.crs, ) barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) @@ -576,16 +567,12 @@ def _voronoi_by_ca( # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) - if SPL_GE_210: + if simplify: # Simplify coverages with coverage_simplify. # torelance set as the square root of the isosceles right triangle with 2 cells_size edges. # For the behavior of coverage_simplify see https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html grid_gdf["geometry"] = shapely.coverage_simplify( - grid_gdf["geometry"].array, tolerance= ((2 * cell_size) ** 2 / 2) ** 0.5 - ) - else: - raise ImportError( - "Shapely 2.1.0 or higher is required for coverage_simplify." + grid_gdf["geometry"].array, tolerance=((2 * cell_size) ** 2 / 2) ** 0.5 ) return grid_gdf @@ -614,7 +601,7 @@ def _get_inner_barriers(enclosure, barriers): # Only keep the geometry which is within the enclosure inner_barriers = inner_barriers[inner_barriers.within(enclosure)] - return GeoSeries(inner_barriers.geometry, crs=barriers.crs) + return inner_barriers class CellState(Enum): From 5f77453a77a52e20e2e1c3a9d5aae222c814b56f Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 25 Sep 2025 10:34:52 +0200 Subject: [PATCH 07/19] various changes --- momepy/functional/_elements.py | 257 +++++++++++++++++++++------------ 1 file changed, 163 insertions(+), 94 deletions(-) diff --git a/momepy/functional/_elements.py b/momepy/functional/_elements.py index 5c8d6219..0e61a2b2 100644 --- a/momepy/functional/_elements.py +++ b/momepy/functional/_elements.py @@ -1,20 +1,19 @@ +import math import warnings +from collections import Counter, deque +from enum import Enum import geopandas as gpd import libpysal -import math import numpy as np import pandas as pd import shapely -from collections import deque, Counter -from enum import Enum from geopandas import GeoDataFrame, GeoSeries from joblib import Parallel, delayed from libpysal.cg import voronoi_frames from libpysal.graph import Graph from packaging.version import Version from pandas import MultiIndex, Series -from typing import Tuple, List GPD_GE_013 = Version(gpd.__version__) >= Version("0.13.0") GPD_GE_10 = Version(gpd.__version__) >= Version("1.0dev") @@ -135,7 +134,8 @@ def enclosed_tessellation( threshold: float = 0.05, n_jobs: int = -1, inner_barriers: GeoSeries | GeoDataFrame = None, - **kwargs, + cell_size: float = 1, + neighbor_mode: str = "moore", ) -> GeoDataFrame: """Generate enclosed tessellation @@ -186,9 +186,16 @@ def enclosed_tessellation( The number of jobs to run in parallel. -1 means using all available cores. By default -1 inner_barriers: GeoSeries | GeoDataFrame, optional - Barriers that should be included in the tessellation process. By default None. - **kwargs - Additional keyword arguments passed to the Cellular Automata-based Voronoi Diagram. + Barriers that should be included in the tessellation process. By passing + inner barriers, tessellation will be derived using less performant + Cellular Automata implememtation but it will recognise dangling barries as valid + limits of the cell growth. See the user guide for details. By default None. + cell_size : float, optional + Grid cell size when ``inner_barriers`` is not None. Otherwise ignored. + By default 1.0 + neighbor_mode : str, optional + Choice of neighbor connectivity ('moore' or 'neumann') when ``inner_barriers`` + is not None. Otherwise ignored.. By default 'moore'. Warnings -------- @@ -277,7 +284,14 @@ def enclosed_tessellation( # generate tessellation in parallel new = Parallel(n_jobs=n_jobs)( delayed(_tess)( - *t, threshold, shrink, segment, index_name, inner_barriers, **kwargs + *t, + threshold, + shrink, + segment, + index_name, + inner_barriers, + cell_size, + neighbor_mode, ) for t in tuples ) @@ -312,14 +326,48 @@ def enclosed_tessellation( def _tess( - ix, poly, blg, threshold, shrink, segment, enclosure_id, inner_barriers, **kwargs + ix, + poly, + blg, + threshold, + shrink, + segment, + enclosure_id, + inner_barriers, + cell_size, + neighbor_mode, ): - """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" + """Generate tessellation for a single enclosure. Helper for enclosed_tessellation. - # Extract CA parameters from kwargs - cell_size = kwargs.pop("cell_size", 1.0) - neighbor_mode = kwargs.pop("neighbor_mode", "moore") - simplify = kwargs.pop("simplify", True) + Parameters + ---------- + ix : int + Enclosure index. + poly : shapely.Geometry + Enclosure geometry. + blg : GeoSeries | GeoDataFrame + Buildings within the enclosure. + threshold : float + Threshold for building inclusion. + shrink : float + Shrink distance for tessellation. + segment : float + Segmentation distance. + enclosure_id : str + Column name for enclosure ID. + inner_barriers : GeoSeries | GeoDataFrame + Inner barriers for tessellation. + cell_size : float + Grid cell size when inner_barriers is not None. + neighbor_mode : str + Choice of neighbor connectivity ('moore' or 'neumann') + when inner_barriers is not None. + + Returns + ------- + GeoDataFrame + Tessellation for the enclosure. + """ # check if threshold is set and filter buildings based on the threshold if threshold: @@ -345,7 +393,6 @@ def _tess( cell_size=cell_size, neighbor_mode=neighbor_mode, barriers_for_inner=inner_barriers, - simplify=simplify, ) tess[enclosure_id] = ix return tess @@ -370,46 +417,53 @@ def _voronoi_by_ca( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, - simplify: bool = True, ) -> GeoDataFrame: """ - Generate an aggregated Voronoi tessellation as a GeoDataFrame via a cellular automata. + Generate an aggregated Voronoi tessellation as a GeoDataFrame via cellular automata. This unified function performs the following: - - Ensures that the CRS of seed and barrier geometries are aligned. - - Combines inner barriers with the enclosure, . - - Computes grid bounds covering both seed and barrier geometries. - - Marks barrier cells using a prepared geometry for fast intersection. - - Seeds the grid with seed geometries. - - Propagates seed values via a BFS expansion that respects barriers. - - Uses a voting mechanism to finalize boundary cell assignments. - - Converts grid cells into a GeoDataFrame and dissolves adjacent cells + + - Ensures that the CRS of seed and barrier geometries are aligned. + - Combines inner barriers with the enclosure. + - Computes grid bounds covering both seed and barrier geometries. + - Marks barrier cells using a prepared geometry for fast intersection. + - Seeds the grid with seed geometries. + - Propagates seed values via a BFS expansion that respects barriers. + - Uses a voting mechanism to finalize boundary cell assignments. + - Converts grid cells into a GeoDataFrame and dissolves adjacent cells with the same seed id. - - Clips the output cells by the barrier boundaries. + - Clips the output cells by the barrier boundaries. - Parameters: - seed_geoms: GeoDataFrame containing seed features. - barrier_geoms: GeoDataFrame containing barrier features or a shapely Polygon. - cell_size: Grid cell size. By default it is 1.0. - neighbor_mode: Choice of neighbor connectivity ('moore' or 'neumann'). By default it is 'moore'. - barriers_for_inner: GeoDataFrame containing inner barriers to be included. By default it is None - simplify: Whether to simplify the coverage using shapely's coverage_simplify. By default it is True. + Parameters + ---------- + seed_geoms : GeoSeries | GeoDataFrame + GeoDataFrame containing seed features. + barrier_geoms : GeoSeries | GeoDataFrame + GeoDataFrame containing barrier features or a shapely Polygon. + cell_size : float, optional + Grid cell size. By default 1.0 + neighbor_mode : str, optional + Choice of neighbor connectivity ('moore' or 'neumann'). By default 'moore'. + barriers_for_inner : GeoSeries | GeoDataFrame, optional + GeoDataFrame containing inner barriers to be included. By default None - Returns: - A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by barriers. + Returns + ------- + GeoDataFrame + A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by + barriers. """ # Check if shapely version is 2.1.0 or higher for coverage_simplify if not SPL_GE_210: raise ImportError( - "Shapely 2.1.0 or higher is required for tessellation with inner barriers provided." + "Shapely 2.1.0 or higher is required for tessellation with inner " + "barriers provided." ) # Get inner barriers as intersection or containment of the barrier_geoms - inner_barriers = _get_inner_barriers( - barrier_geoms, barriers_for_inner.to_crs(seed_geoms.crs) - ) + inner_barriers = _get_inner_barriers(barrier_geoms, barriers_for_inner) # Handle barrier_geoms if it is a Polygon or MultiPolygon if barrier_geoms.geom_type == "Polygon": @@ -420,7 +474,8 @@ def _voronoi_by_ca( barrier_geoms = GeoSeries([barrier_geoms], crs=seed_geoms.crs) elif barrier_geoms.geom_type == "MultiPolygon": - # Process each polygon: take buffer then exterior boundary (to ensure there's no gap between enclosures) + # Process each polygon: take buffer then exterior boundary + # (to ensure there's no gap between enclosures) barrier_geoms_buffered = GeoSeries( shapely.buffer( shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10 * cell_size @@ -430,16 +485,12 @@ def _voronoi_by_ca( barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) else: - raise ValueError("barrier_geoms must be a Polygon or MultiPolygon") + raise ValueError("Enclosure must be a Polygon or MultiPolygon") outer_union = barrier_geoms_buffered.union_all() # Compute inner barriers union if available - if ( - "inner_barriers" in locals() - and inner_barriers is not None - and not inner_barriers.empty - ): + if inner_barriers is not None and not inner_barriers.empty: inner_union = inner_barriers.union_all() else: inner_union = None @@ -466,7 +517,7 @@ def _voronoi_by_ca( cell_polys = GeoSeries( [ _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xs.flatten(), ys.flatten()) + for x, y in zip(xs.flatten(), ys.flatten(), strict=True) ] ) @@ -503,7 +554,8 @@ def _voronoi_by_ca( if states[y_current, x_current] != CellState.FRONTIER.value: continue - # Get neighbor cells that were already assigned a seed id or still unknown (state >= 0). + # Get neighbor cells that were already assigned a seed id or still unknown + # (state >= 0). # Note that boundary or barrier cells are skipped (state < 0). neighbor_seeds = [ states[ny, nx] @@ -518,7 +570,7 @@ def _voronoi_by_ca( # Assign as a boundary if multiple seed ids are found. if len(set(neighbor_seeds)) > 1: states[y_current, x_current] = CellState.BOUNDARY.value - # EIf not, equeue neighbor cells for further propagation. + # If not, equeue neighbor cells for further propagation. else: assigned_seed = set(neighbor_seeds).pop() states[y_current, x_current] = assigned_seed @@ -539,7 +591,7 @@ def _voronoi_by_ca( xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) grid_polys = [ _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xs.flatten(), ys.flatten()) + for x, y in zip(xs.flatten(), ys.flatten(), strict=True) ] grid_gdf = GeoDataFrame( {"site_id": states.flatten()}, geometry=grid_polys, crs=seed_geoms.crs @@ -557,9 +609,10 @@ def _voronoi_by_ca( if barrier_geoms is not None and (not barrier_geoms.empty): # Create a union of the barrier geometries. barrier_union = barrier_geoms.union_all() - # If the barrier union is not a polygon (e.g., it's a MultiLineString), polygonize it. + # If the barrier union is not a polygon (e.g., it's a MultiLineString), + # polygonize it. if not isinstance( - barrier_union, (shapely.geometry.Polygon, shapely.geometry.MultiPolygon) + barrier_union, shapely.geometry.Polygon | shapely.geometry.MultiPolygon ): barrier_polys = shapely.get_parts(shapely.polygonize(barrier_union)) if not barrier_polys.empty: @@ -567,36 +620,34 @@ def _voronoi_by_ca( # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) - if simplify: - # Simplify coverages with coverage_simplify. - # torelance set as the square root of the isosceles right triangle with 2 cells_size edges. - # For the behavior of coverage_simplify see https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html - grid_gdf["geometry"] = shapely.coverage_simplify( - grid_gdf["geometry"].array, tolerance=((2 * cell_size) ** 2 / 2) ** 0.5 - ) + # Simplify coverages with coverage_simplify. + # torelance set as the square root of the isosceles right triangle with 2 + # cells_size edges. For the behavior of coverage_simplify see + # https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html + grid_gdf["geometry"] = shapely.coverage_simplify( + grid_gdf["geometry"].array, tolerance=((2 * cell_size) ** 2 / 2) ** 0.5 + ) return grid_gdf def _get_inner_barriers(enclosure, barriers): - """ - Get inner barriers that intersect or are contained within an enclosure. + """Get inner barriers that intersect or are contained within an enclosure. - Args: - enclosure (GeoSeries or GeoDataFrame): The enclosure geometry. - barriers (GeoDataFrame): The barriers GeoDataFrame. + Parameters + ---------- + enclosure : GeoSeries | GeoDataFrame + The enclosure geometry. + barriers : GeoDataFrame + The barriers GeoDataFrame. - Returns: - shapely.geometry.Polygon: A single Polygon combining the enclosure - and any intersecting barriers. + Returns + ------- + shapely.Polygon + A single Polygon combining the enclosure and any intersecting barriers. """ - # Find barriers intersecting or contained in the enclosure - inner_barriers = barriers[ - barriers.intersects(enclosure) | barriers.contains(enclosure) - ] - # Clip those segments to stay within the enclosure - inner_barriers = gpd.clip(inner_barriers, enclosure) + inner_barriers = gpd.clip(barriers, enclosure) # Only keep the geometry which is within the enclosure inner_barriers = inner_barriers[inner_barriers.within(enclosure)] @@ -606,13 +657,19 @@ def _get_inner_barriers(enclosure, barriers): class CellState(Enum): """ - Enumeration of cell states for grid processing for improving the readability, instead of integers. + Enumeration of cell states for grid processing for improving the readability, + instead of integers. - Attributes: - UNKNOWN: Cell has not been processed. - BOUNDARY: Cell is at a junction between different seed regions. - BARRIER: Cell originally designated as a barrier. - FRONTIER: Cell queued for BFS expansion. + Attributes + ---------- + UNKNOWN : int + Cell has not been processed. + BOUNDARY : int + Cell is at a junction between different seed regions. + BARRIER : int + Cell originally designated as a barrier. + FRONTIER : int + Cell queued for BFS expansion. """ UNKNOWN = -1 @@ -622,7 +679,7 @@ class CellState(Enum): def _get_cell_polygon( - x_idx: int, y_idx: int, cell_size: float = 1.0, origin: Tuple[float, float] = (0, 0) + x_idx: int, y_idx: int, cell_size: float = 1.0, origin: tuple[float, float] = (0, 0) ) -> shapely.geometry.Polygon: """ Generate a grid cell polygon based on the given indices, cell size, and origin. @@ -640,16 +697,27 @@ def _get_cell_polygon( def _get_neighbors( x: int, y: int, max_x: int, max_y: int, mode: str = "moore" -) -> List[Tuple[int, int]]: +) -> list[tuple[int, int]]: """ Retrieve valid neighboring cell indices based on connectivity. - Parameters: - x, y: The current cell indices. - max_x, max_y: Dimensions of the grid. - mode: "moore" for 8-connected or "neumann" for 4-connected neighbors. + Parameters + ---------- + x : int + The current cell x index. + y : int + The current cell y index. + max_x : int + The maximum x dimension of the grid. + max_y : int + The maximum y dimension of the grid. + mode : str, optional + The connectivity mode, "moore" for 8-connected or "neumann" for 4-connected + neighbors. By default "moore". - Returns: + Returns + ------- + list[tuple[int, int]] A list of (x, y) tuples for valid neighbor indices. """ neighbor_dirs = { @@ -670,7 +738,7 @@ def _get_grid_bounds( seed_geoms: GeoSeries | GeoDataFrame, barrier_geoms: GeoSeries | GeoDataFrame, cell_size: float, -) -> Tuple[Tuple[float, float], int, int]: +) -> tuple[tuple[float, float], int, int]: seed_bounds = seed_geoms.total_bounds # [xmin, ymin, xmax, ymax] barrier_bounds = barrier_geoms.total_bounds @@ -690,11 +758,11 @@ def _get_grid_bounds( def _geom_to_cells( geom: shapely.geometry, - origin: Tuple[float, float], + origin: tuple[float, float], cell_size: float, grid_width: int, grid_height: int, -) -> List[Tuple[int, int]]: +) -> list[tuple[int, int]]: """ Determine grid cell indices that intersect the given geometry. """ @@ -716,11 +784,11 @@ def _geom_to_cells( candidate_polys = GeoSeries( [ _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xx.flatten(), yy.flatten()) + for x, y in zip(xx.flatten(), yy.flatten(), strict=True) ] ) mask = candidate_polys.intersects(geom) - return list(zip(xx.flatten()[mask], yy.flatten()[mask])) + return list(zip(xx.flatten()[mask], yy.flatten()[mask], strict=True)) def _enqueue_neighbors( @@ -745,7 +813,8 @@ def _assign_adjacent_seed_cells( states: np.ndarray, neighbor_mode: str = "moore" ) -> np.ndarray: """ - Reassign border and barrier cells to the proximate seed areas using a voting mechanism. + Reassign border and barrier cells to the proximate seed areas using a voting + mechanism. """ new_states = states.copy() indices = np.argwhere( @@ -789,7 +858,7 @@ def verify_tessellation(tessellation, geometry): Returns ------- tuple(excluded, multipolygons) - Tuple of indices of building IDs not present in tessellations and MultiPolygons. + tuple of indices of building IDs not present in tessellations and MultiPolygons. Examples -------- From 7da1bc87c0f58c1fae39f76ce855f23dff35bb5d Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 25 Sep 2025 10:41:17 +0200 Subject: [PATCH 08/19] move stuff --- momepy/elements.py | 533 +++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 518 insertions(+), 15 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index c0abe3da..3c21cb60 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -1,4 +1,7 @@ +import math import warnings +from collections import Counter, deque +from enum import Enum import geopandas as gpd import numpy as np @@ -153,6 +156,9 @@ def enclosed_tessellation( threshold: float = 0.05, simplify: bool = True, n_jobs: int = -1, + inner_barriers: GeoSeries | GeoDataFrame = None, + cell_size: float = 1, + neighbor_mode: str = "moore", **kwargs, ) -> GeoDataFrame: """Generate enclosed tessellation @@ -206,6 +212,17 @@ def enclosed_tessellation( n_jobs : int, optional The number of jobs to run in parallel. -1 means using all available cores. By default -1 + inner_barriers: GeoSeries | GeoDataFrame, optional + Barriers that should be included in the tessellation process. By passing + inner barriers, tessellation will be derived using less performant + Cellular Automata implememtation but it will recognise dangling barries as valid + limits of the cell growth. See the user guide for details. By default None. + cell_size : float, optional + Grid cell size when ``inner_barriers`` is not None. Otherwise ignored. + By default 1.0 + neighbor_mode : str, optional + Choice of neighbor connectivity ('moore' or 'neumann') when ``inner_barriers`` + is not None. Otherwise ignored.. By default 'moore'. **kwargs Additional keyword arguments pased to libpysal.cg.voronoi_frames, such as ``grid_size``. @@ -299,7 +316,18 @@ def enclosed_tessellation( # generate tessellation in parallel new = Parallel(n_jobs=n_jobs)( - delayed(_tess)(*t, threshold, shrink, segment, index_name, simplify, kwargs) + delayed(_tess)( + *t, + threshold, + shrink, + segment, + index_name, + simplify, + inner_barriers, + cell_size, + neighbor_mode, + kwargs, + ) for t in tuples ) @@ -332,8 +360,51 @@ def enclosed_tessellation( return pd.concat([new_df, singles.drop(columns="position"), clean_blocks]) -def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify, kwargs): - """Generate tessellation for a single enclosure. Helper for enclosed_tessellation""" +def _tess( + ix, + poly, + blg, + threshold, + shrink, + segment, + enclosure_id, + to_simplify, + inner_barriers, + cell_size, + neighbor_mode, + kwargs, +): + """Generate tessellation for a single enclosure. Helper for enclosed_tessellation. + + Parameters + ---------- + ix : int + Enclosure index. + poly : shapely.Geometry + Enclosure geometry. + blg : GeoSeries | GeoDataFrame + Buildings within the enclosure. + threshold : float + Threshold for building inclusion. + shrink : float + Shrink distance for tessellation. + segment : float + Segmentation distance. + enclosure_id : str + Column name for enclosure ID. + inner_barriers : GeoSeries | GeoDataFrame + Inner barriers for tessellation. + cell_size : float + Grid cell size when inner_barriers is not None. + neighbor_mode : str + Choice of neighbor connectivity ('moore' or 'neumann') + when inner_barriers is not None. + + Returns + ------- + GeoDataFrame + Tessellation for the enclosure. + """ # check if threshold is set and filter buildings based on the threshold if threshold: blg = blg[ @@ -342,19 +413,30 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify, ] if len(blg) > 1: - tess = voronoi_frames( - blg, - clip=poly, - shrink=shrink, - segment=segment, - return_input=False, - as_gdf=True, - **kwargs, - ) - if to_simplify: - tess.geometry = shapely.coverage_simplify( - tess.geometry, tolerance=segment / 2, simplify_boundary=False + if inner_barriers is None or inner_barriers.empty: + tess = voronoi_frames( + blg, + clip=poly, + shrink=shrink, + segment=segment, + return_input=False, + as_gdf=True, + **kwargs, ) + if to_simplify: + tess.geometry = shapely.coverage_simplify( + tess.geometry, tolerance=segment / 2, simplify_boundary=False + ) + else: + tess = _voronoi_by_ca( + seed_geoms=blg, + barrier_geoms=poly, + cell_size=cell_size, + neighbor_mode=neighbor_mode, + barriers_for_inner=inner_barriers, + to_simplify=to_simplify, + ) + tess[enclosure_id] = ix return tess @@ -372,6 +454,427 @@ def _tess(ix, poly, blg, threshold, shrink, segment, enclosure_id, to_simplify, ) +def _voronoi_by_ca( + seed_geoms: GeoSeries | GeoDataFrame, + barrier_geoms: GeoSeries | GeoDataFrame, + cell_size: float = 1.0, + neighbor_mode: str = "moore", + barriers_for_inner: GeoSeries | GeoDataFrame = None, + to_simplify: bool = True, +) -> GeoDataFrame: + """ + Generate an aggregated Voronoi tessellation as a GeoDataFrame via cellular automata. + + This unified function performs the following: + + - Ensures that the CRS of seed and barrier geometries are aligned. + - Combines inner barriers with the enclosure. + - Computes grid bounds covering both seed and barrier geometries. + - Marks barrier cells using a prepared geometry for fast intersection. + - Seeds the grid with seed geometries. + - Propagates seed values via a BFS expansion that respects barriers. + - Uses a voting mechanism to finalize boundary cell assignments. + - Converts grid cells into a GeoDataFrame and dissolves adjacent cells + with the same seed id. + - Clips the output cells by the barrier boundaries. + + Parameters + ---------- + seed_geoms : GeoSeries | GeoDataFrame + GeoDataFrame containing seed features. + barrier_geoms : GeoSeries | GeoDataFrame + GeoDataFrame containing barrier features or a shapely Polygon. + cell_size : float, optional + Grid cell size. By default 1.0 + neighbor_mode : str, optional + Choice of neighbor connectivity ('moore' or 'neumann'). By default 'moore'. + barriers_for_inner : GeoSeries | GeoDataFrame, optional + GeoDataFrame containing inner barriers to be included. By default None + + + Returns + ------- + GeoDataFrame + A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by + barriers. + """ + # Get inner barriers as intersection or containment of the barrier_geoms + inner_barriers = _get_inner_barriers(barrier_geoms, barriers_for_inner) + + # Handle barrier_geoms if it is a Polygon or MultiPolygon + if barrier_geoms.geom_type == "Polygon": + # Take buffer of polygon and extract its exterior boundary (10 cells) + barrier_geoms_buffered = GeoSeries( + [barrier_geoms.buffer(10 * cell_size).exterior], crs=seed_geoms.crs + ) + barrier_geoms = GeoSeries([barrier_geoms], crs=seed_geoms.crs) + + elif barrier_geoms.geom_type == "MultiPolygon": + # Process each polygon: take buffer then exterior boundary + # (to ensure there's no gap between enclosures) + barrier_geoms_buffered = GeoSeries( + shapely.buffer( + shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10 * cell_size + ), + crs=seed_geoms.crs, + ) + barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) + + else: + raise ValueError("Enclosure must be a Polygon or MultiPolygon") + + outer_union = barrier_geoms_buffered.union_all() + + # Compute inner barriers union if available + if inner_barriers is not None and not inner_barriers.empty: + inner_union = inner_barriers.union_all() + else: + inner_union = None + + # Combine outer barrier with inner barriers + if outer_union and inner_union: + prep_barrier = shapely.union(outer_union, inner_union) + elif outer_union: + prep_barrier = outer_union + elif inner_union: + prep_barrier = inner_union + else: + prep_barrier = None + + # Compute grid bounds + origin, grid_width, grid_height = _get_grid_bounds( + seed_geoms, barrier_geoms_buffered, cell_size + ) + + # Initialize grid states with UNKNOWN values. + states = np.full((grid_height, grid_width), CellState.UNKNOWN.value, dtype=int) + + xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) + cell_polys = GeoSeries( + [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xs.flatten(), ys.flatten(), strict=True) + ] + ) + + # Identify barrier cells in the grid + if prep_barrier is not None: + barrier_mask = cell_polys.intersects(prep_barrier).values.reshape( + grid_height, grid_width + ) + else: + barrier_mask = np.zeros((grid_height, grid_width), dtype=bool) + states[barrier_mask] = CellState.BARRIER.value + + # Seed the grid with seed geometries. + for site_id, geom in enumerate(seed_geoms.geometry): + if not geom.is_empty: + cells = _geom_to_cells(geom, origin, cell_size, grid_width, grid_height) + valid_cells = [ + (x, y) for x, y in cells if states[y, x] == CellState.UNKNOWN.value + ] + if valid_cells: + indices = np.array(valid_cells) + states[indices[:, 1], indices[:, 0]] = site_id + + # Initialize the BFS queue with all seeded cells’ neighbors. + queue = deque() + seed_indices = np.argwhere(states >= 0) + for y, x in seed_indices: + _enqueue_neighbors(x, y, states, grid_width, grid_height, neighbor_mode, queue) + + # Process BFS to propagate seed values. + while queue: + # Dequeue the current cell and skip if it is not a frontier cell. + x_current, y_current = queue.popleft() + if states[y_current, x_current] != CellState.FRONTIER.value: + continue + + # Get neighbor cells that were already assigned a seed id or still unknown + # (state >= 0). + # Note that boundary or barrier cells are skipped (state < 0). + neighbor_seeds = [ + states[ny, nx] + for nx, ny in _get_neighbors( + x_current, y_current, grid_width, grid_height, mode=neighbor_mode + ) + if states[ny, nx] >= 0 + ] + if not neighbor_seeds: + continue + + # Assign as a boundary if multiple seed ids are found. + if len(set(neighbor_seeds)) > 1: + states[y_current, x_current] = CellState.BOUNDARY.value + # If not, equeue neighbor cells for further propagation. + else: + assigned_seed = set(neighbor_seeds).pop() + states[y_current, x_current] = assigned_seed + _enqueue_neighbors( + x_current, + y_current, + states, + grid_width, + grid_height, + neighbor_mode, + queue, + ) + + # Post-process barrier and boundary cells using a voting mechanism. + states = _assign_adjacent_seed_cells(states, neighbor_mode) + + # Create grid cell polygons and build a GeoDataFrame. + xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) + grid_polys = [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xs.flatten(), ys.flatten(), strict=True) + ] + grid_gdf = GeoDataFrame( + {"site_id": states.flatten()}, geometry=grid_polys, crs=seed_geoms.crs + ) + + # Include only cells with valid seed assignments and dissolve contiguous regions. + grid_gdf = ( + grid_gdf[grid_gdf["site_id"] >= 0] + .dissolve(by="site_id") + .reset_index() + .drop("site_id", axis=1) + ) + + # Clip by barriers + if barrier_geoms is not None and (not barrier_geoms.empty): + # Create a union of the barrier geometries. + barrier_union = barrier_geoms.union_all() + # If the barrier union is not a polygon (e.g., it's a MultiLineString), + # polygonize it. + if not isinstance( + barrier_union, shapely.geometry.Polygon | shapely.geometry.MultiPolygon + ): + barrier_polys = shapely.get_parts(shapely.polygonize(barrier_union)) + if not barrier_polys.empty: + barrier_union = shapely.union_all(barrier_polys) + # Clip each polygon in the grid using the barrier boundary. + grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) + + # Simplify coverages with coverage_simplify. + # torelance set as the square root of the isosceles right triangle with 2 + # cells_size edges. For the behavior of coverage_simplify see + # https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html + if to_simplify: + grid_gdf["geometry"] = shapely.coverage_simplify( + grid_gdf["geometry"].array, tolerance=((2 * cell_size) ** 2 / 2) ** 0.5 + ) + + return grid_gdf + + +def _get_inner_barriers(enclosure, barriers): + """Get inner barriers that intersect or are contained within an enclosure. + + Parameters + ---------- + enclosure : GeoSeries | GeoDataFrame + The enclosure geometry. + barriers : GeoDataFrame + The barriers GeoDataFrame. + + Returns + ------- + shapely.Polygon + A single Polygon combining the enclosure and any intersecting barriers. + """ + # Clip those segments to stay within the enclosure + inner_barriers = gpd.clip(barriers, enclosure) + + # Only keep the geometry which is within the enclosure + inner_barriers = inner_barriers[inner_barriers.within(enclosure)] + + return inner_barriers + + +class CellState(Enum): + """ + Enumeration of cell states for grid processing for improving the readability, + instead of integers. + + Attributes + ---------- + UNKNOWN : int + Cell has not been processed. + BOUNDARY : int + Cell is at a junction between different seed regions. + BARRIER : int + Cell originally designated as a barrier. + FRONTIER : int + Cell queued for BFS expansion. + """ + + UNKNOWN = -1 + BOUNDARY = -2 + BARRIER = -3 + FRONTIER = -4 + + +def _get_cell_polygon( + x_idx: int, y_idx: int, cell_size: float = 1.0, origin: tuple[float, float] = (0, 0) +) -> shapely.geometry.Polygon: + """ + Generate a grid cell polygon based on the given indices, cell size, and origin. + """ + ox, oy = origin + return shapely.Polygon( + [ + (ox + x_idx * cell_size, oy + y_idx * cell_size), + (ox + (x_idx + 1) * cell_size, oy + y_idx * cell_size), + (ox + (x_idx + 1) * cell_size, oy + (y_idx + 1) * cell_size), + (ox + x_idx * cell_size, oy + (y_idx + 1) * cell_size), + ] + ) + + +def _get_neighbors( + x: int, y: int, max_x: int, max_y: int, mode: str = "moore" +) -> list[tuple[int, int]]: + """ + Retrieve valid neighboring cell indices based on connectivity. + + Parameters + ---------- + x : int + The current cell x index. + y : int + The current cell y index. + max_x : int + The maximum x dimension of the grid. + max_y : int + The maximum y dimension of the grid. + mode : str, optional + The connectivity mode, "moore" for 8-connected or "neumann" for 4-connected + neighbors. By default "moore". + + Returns + ------- + list[tuple[int, int]] + A list of (x, y) tuples for valid neighbor indices. + """ + neighbor_dirs = { + "moore": [(-1, -1), (-1, 0), (-1, 1), (0, -1), (0, 1), (1, -1), (1, 0), (1, 1)], + "neumann": [(0, -1), (-1, 0), (1, 0), (0, 1)], + } + directions = neighbor_dirs.get(mode) + if directions is None: + raise ValueError("Invalid neighbor_mode: choose 'moore' or 'neumann'") + return [ + (x + dx, y + dy) + for dx, dy in directions + if 0 <= x + dx < max_x and 0 <= y + dy < max_y + ] + + +def _get_grid_bounds( + seed_geoms: GeoSeries | GeoDataFrame, + barrier_geoms: GeoSeries | GeoDataFrame, + cell_size: float, +) -> tuple[tuple[float, float], int, int]: + seed_bounds = seed_geoms.total_bounds # [xmin, ymin, xmax, ymax] + barrier_bounds = barrier_geoms.total_bounds + + xmin = min(seed_bounds[0], barrier_bounds[0]) + ymin = min(seed_bounds[1], barrier_bounds[1]) + xmax = max(seed_bounds[2], barrier_bounds[2]) + ymax = max(seed_bounds[3], barrier_bounds[3]) + # expand bounds by 1 cell in each direction + new_xmin = xmin - cell_size + new_ymin = ymin - cell_size + new_xmax = xmax + cell_size + new_ymax = ymax + cell_size + grid_width = math.ceil((new_xmax - new_xmin) / cell_size) + grid_height = math.ceil((new_ymax - new_ymin) / cell_size) + return (new_xmin, new_ymin), grid_width, grid_height + + +def _geom_to_cells( + geom: shapely.geometry, + origin: tuple[float, float], + cell_size: float, + grid_width: int, + grid_height: int, +) -> list[tuple[int, int]]: + """ + Determine grid cell indices that intersect the given geometry. + """ + if isinstance(geom, shapely.geometry.Point): + sx = int((geom.x - origin[0]) // cell_size) + sy = int((geom.y - origin[1]) // cell_size) + return [(sx, sy)] if 0 <= sx < grid_width and 0 <= sy < grid_height else [] + + else: + minx, miny, maxx, maxy = geom.bounds + start_x = max(0, int((minx - origin[0]) // cell_size)) + start_y = max(0, int((miny - origin[1]) // cell_size)) + end_x = min(grid_width, int(math.ceil((maxx - origin[0]) / cell_size))) + end_y = min(grid_height, int(math.ceil((maxy - origin[1]) / cell_size))) + + x_range = np.arange(start_x, end_x) + y_range = np.arange(start_y, end_y) + xx, yy = np.meshgrid(x_range, y_range) + candidate_polys = GeoSeries( + [ + _get_cell_polygon(x, y, cell_size, origin) + for x, y in zip(xx.flatten(), yy.flatten(), strict=True) + ] + ) + mask = candidate_polys.intersects(geom) + return list(zip(xx.flatten()[mask], yy.flatten()[mask], strict=True)) + + +def _enqueue_neighbors( + x: int, + y: int, + states: np.ndarray, + grid_width: int, + grid_height: int, + neighbor_mode: str, + queue: deque, +) -> None: + """ + Enqueue valid neighboring cells for BFS expansion. + """ + for nx, ny in _get_neighbors(x, y, grid_width, grid_height, mode=neighbor_mode): + if states[ny, nx] == CellState.UNKNOWN.value: + states[ny, nx] = CellState.FRONTIER.value + queue.append((nx, ny)) + + +def _assign_adjacent_seed_cells( + states: np.ndarray, neighbor_mode: str = "moore" +) -> np.ndarray: + """ + Reassign border and barrier cells to the proximate seed areas using a voting + mechanism. + """ + new_states = states.copy() + indices = np.argwhere( + np.isin(states, [CellState.BARRIER.value, CellState.BOUNDARY.value]) + ) + grid_height, grid_width = states.shape + + for y, x in indices: + neighbor_seeds = [ + states[ny, nx] + for nx, ny in _get_neighbors( + x, y, grid_width, grid_height, mode=neighbor_mode + ) + if states[ny, nx] >= 0 + ] + if neighbor_seeds: + cnt = Counter(neighbor_seeds) + # In case of ties, choose the smaller seed id. + chosen_seed = min(cnt.items(), key=lambda item: (-item[1], item[0]))[0] + new_states[y, x] = chosen_seed + return new_states + + def verify_tessellation(tessellation, geometry): """Check whether result matches buildings and contains only Polygons. From 89c7c554a7871d9e857b974c4ecc1857720e5045 Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 25 Sep 2025 10:59:28 +0200 Subject: [PATCH 09/19] cleanup and simplify --- momepy/elements.py | 44 ++++++++++++++++++++++++++++++-------------- 1 file changed, 30 insertions(+), 14 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index 3c21cb60..39f64747 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -423,10 +423,8 @@ def _tess( as_gdf=True, **kwargs, ) - if to_simplify: - tess.geometry = shapely.coverage_simplify( - tess.geometry, tolerance=segment / 2, simplify_boundary=False - ) + tolerance = segment / 2 + else: tess = _voronoi_by_ca( seed_geoms=blg, @@ -434,7 +432,14 @@ def _tess( cell_size=cell_size, neighbor_mode=neighbor_mode, barriers_for_inner=inner_barriers, - to_simplify=to_simplify, + ) + # torelance set as the square root of the isosceles right triangle with 2 + # cells_size edges + tolerance = ((2 * cell_size) ** 2 / 2) ** 0.5 + + if to_simplify: + tess.geometry = shapely.coverage_simplify( + tess.geometry, tolerance=tolerance, simplify_boundary=False ) tess[enclosure_id] = ix @@ -460,7 +465,6 @@ def _voronoi_by_ca( cell_size: float = 1.0, neighbor_mode: str = "moore", barriers_for_inner: GeoSeries | GeoDataFrame = None, - to_simplify: bool = True, ) -> GeoDataFrame: """ Generate an aggregated Voronoi tessellation as a GeoDataFrame via cellular automata. @@ -656,15 +660,27 @@ def _voronoi_by_ca( # Clip each polygon in the grid using the barrier boundary. grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) - # Simplify coverages with coverage_simplify. - # torelance set as the square root of the isosceles right triangle with 2 - # cells_size edges. For the behavior of coverage_simplify see - # https://shapely.readthedocs.io/en/latest/reference/shapely.coverage_simplify.html - if to_simplify: - grid_gdf["geometry"] = shapely.coverage_simplify( - grid_gdf["geometry"].array, tolerance=((2 * cell_size) ** 2 / 2) ** 0.5 - ) + # cleanup possible collections + def sanitize_geometry(geom): + if geom.geom_type in ["Polygon", "MultiPolygon"]: + return geom + elif geom.geom_type == "GeometryCollection": + parts = shapely.get_parts(geom) + valid_parts = [ + p for p in parts if p.geom_type in ["Polygon", "MultiPolygon"] + ] + if len(valid_parts) == 0: + return None + elif len(valid_parts) == 1: + return valid_parts[0] + else: + return shapely.MultiPolygon(valid_parts) + else: + # Drop points, lines, etc. + return None + grid_gdf["geometry"] = grid_gdf["geometry"].apply(sanitize_geometry) + grid_gdf = grid_gdf.dropna(subset=["geometry"]) return grid_gdf From bc4d83ce748f212478c953a5d3a0df6e898cc4ff Mon Sep 17 00:00:00 2001 From: Martin Fleischmann Date: Thu, 25 Sep 2025 11:24:40 +0200 Subject: [PATCH 10/19] make mypy happy --- momepy/elements.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/momepy/elements.py b/momepy/elements.py index 39f64747..91863c1e 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -582,7 +582,9 @@ def _voronoi_by_ca( states[indices[:, 1], indices[:, 0]] = site_id # Initialize the BFS queue with all seeded cells’ neighbors. - queue = deque() + # ... existing code ... + + queue: deque[tuple[int, int]] = deque() seed_indices = np.argwhere(states >= 0) for y, x in seed_indices: _enqueue_neighbors(x, y, states, grid_width, grid_height, neighbor_mode, queue) From 552cead76dcb1016a62a6342daaf857290b6c8dd Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Thu, 25 Sep 2025 16:44:37 +0100 Subject: [PATCH 11/19] Replace per-cell polygon loops with vectorisations --- momepy/elements.py | 36 ++++++++++++++++++------------------ test_tess.py | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 18 deletions(-) create mode 100644 test_tess.py diff --git a/momepy/elements.py b/momepy/elements.py index 91863c1e..d5f1e782 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -554,16 +554,18 @@ def _voronoi_by_ca( states = np.full((grid_height, grid_width), CellState.UNKNOWN.value, dtype=int) xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) - cell_polys = GeoSeries( - [ - _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xs.flatten(), ys.flatten(), strict=True) - ] + xs_flat = xs.ravel() + ys_flat = ys.ravel() + cell_polys_array = shapely.box( + origin[0] + xs_flat * cell_size, + origin[1] + ys_flat * cell_size, + origin[0] + (xs_flat + 1) * cell_size, + origin[1] + (ys_flat + 1) * cell_size, ) # Identify barrier cells in the grid if prep_barrier is not None: - barrier_mask = cell_polys.intersects(prep_barrier).values.reshape( + barrier_mask = shapely.intersects(cell_polys_array, prep_barrier).reshape( grid_height, grid_width ) else: @@ -630,11 +632,7 @@ def _voronoi_by_ca( states = _assign_adjacent_seed_cells(states, neighbor_mode) # Create grid cell polygons and build a GeoDataFrame. - xs, ys = np.meshgrid(np.arange(grid_width), np.arange(grid_height)) - grid_polys = [ - _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xs.flatten(), ys.flatten(), strict=True) - ] + grid_polys = GeoSeries(cell_polys_array, crs=seed_geoms.crs) grid_gdf = GeoDataFrame( {"site_id": states.flatten()}, geometry=grid_polys, crs=seed_geoms.crs ) @@ -836,14 +834,16 @@ def _geom_to_cells( x_range = np.arange(start_x, end_x) y_range = np.arange(start_y, end_y) xx, yy = np.meshgrid(x_range, y_range) - candidate_polys = GeoSeries( - [ - _get_cell_polygon(x, y, cell_size, origin) - for x, y in zip(xx.flatten(), yy.flatten(), strict=True) - ] + xx_flat = xx.ravel() + yy_flat = yy.ravel() + candidate_polys = shapely.box( + origin[0] + xx_flat * cell_size, + origin[1] + yy_flat * cell_size, + origin[0] + (xx_flat + 1) * cell_size, + origin[1] + (yy_flat + 1) * cell_size, ) - mask = candidate_polys.intersects(geom) - return list(zip(xx.flatten()[mask], yy.flatten()[mask], strict=True)) + mask = shapely.intersects(candidate_polys, geom) + return list(zip(xx_flat[mask], yy_flat[mask], strict=True)) def _enqueue_neighbors( diff --git a/test_tess.py b/test_tess.py new file mode 100644 index 00000000..f56a2a84 --- /dev/null +++ b/test_tess.py @@ -0,0 +1,34 @@ +import time + +import geopandas as gpd +import momepy + +# Load data +path = momepy.datasets.get_path("bubenec") +buildings = gpd.read_file(path, layer="buildings") +streets = gpd.read_file(path, layer="streets") + +# Generate enclosures +enclosures = momepy.enclosures(streets) + +# Use streets as inner_barriers for testing +inner_barriers = streets + +# Time the enclosed_tessellation without inner_barriers +start = time.time() +tess_no_inner = momepy.enclosed_tessellation(buildings, enclosures) +end = time.time() +time_no_inner = end - start + +# Time the enclosed_tessellation with inner_barriers +start = time.time() +tess_with_inner = momepy.enclosed_tessellation( + buildings, enclosures, inner_barriers=inner_barriers +) +end = time.time() +time_with_inner = end - start + +print(f"Time without inner_barriers: {time_no_inner} seconds") +print(f"Time with inner_barriers: {time_with_inner} seconds") +print(f"Number of tessellation cells without: {len(tess_no_inner)}") +print(f"Number of tessellation cells with: {len(tess_with_inner)}") From 9879473b79333abc6f909d2b2dbd9ed180efbfaa Mon Sep 17 00:00:00 2001 From: Yuta Sato <65852668+yu-ta-sato@users.noreply.github.com> Date: Sat, 27 Sep 2025 23:35:46 +0100 Subject: [PATCH 12/19] Delete test_tess.py --- test_tess.py | 34 ---------------------------------- 1 file changed, 34 deletions(-) delete mode 100644 test_tess.py diff --git a/test_tess.py b/test_tess.py deleted file mode 100644 index f56a2a84..00000000 --- a/test_tess.py +++ /dev/null @@ -1,34 +0,0 @@ -import time - -import geopandas as gpd -import momepy - -# Load data -path = momepy.datasets.get_path("bubenec") -buildings = gpd.read_file(path, layer="buildings") -streets = gpd.read_file(path, layer="streets") - -# Generate enclosures -enclosures = momepy.enclosures(streets) - -# Use streets as inner_barriers for testing -inner_barriers = streets - -# Time the enclosed_tessellation without inner_barriers -start = time.time() -tess_no_inner = momepy.enclosed_tessellation(buildings, enclosures) -end = time.time() -time_no_inner = end - start - -# Time the enclosed_tessellation with inner_barriers -start = time.time() -tess_with_inner = momepy.enclosed_tessellation( - buildings, enclosures, inner_barriers=inner_barriers -) -end = time.time() -time_with_inner = end - start - -print(f"Time without inner_barriers: {time_no_inner} seconds") -print(f"Time with inner_barriers: {time_with_inner} seconds") -print(f"Number of tessellation cells without: {len(tess_no_inner)}") -print(f"Number of tessellation cells with: {len(tess_with_inner)}") From 9e8527917096bdfddbca7cf495306a7eb7173d5c Mon Sep 17 00:00:00 2001 From: Yuta Sato <65852668+yu-ta-sato@users.noreply.github.com> Date: Sat, 27 Sep 2025 23:39:55 +0100 Subject: [PATCH 13/19] Remove typo Co-authored-by: James Gaboardi --- momepy/elements.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/momepy/elements.py b/momepy/elements.py index d5f1e782..9a9b4d3e 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -222,7 +222,7 @@ def enclosed_tessellation( By default 1.0 neighbor_mode : str, optional Choice of neighbor connectivity ('moore' or 'neumann') when ``inner_barriers`` - is not None. Otherwise ignored.. By default 'moore'. + is not None. Otherwise ignored. By default 'moore'. **kwargs Additional keyword arguments pased to libpysal.cg.voronoi_frames, such as ``grid_size``. From df212451be3f548442c39442b9739f727ee97eef Mon Sep 17 00:00:00 2001 From: Yuta Sato <65852668+yu-ta-sato@users.noreply.github.com> Date: Sat, 27 Sep 2025 23:40:26 +0100 Subject: [PATCH 14/19] Remove unnecessary empty line Co-authored-by: James Gaboardi --- momepy/elements.py | 1 - 1 file changed, 1 deletion(-) diff --git a/momepy/elements.py b/momepy/elements.py index 9a9b4d3e..bc44f698 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -495,7 +495,6 @@ def _voronoi_by_ca( barriers_for_inner : GeoSeries | GeoDataFrame, optional GeoDataFrame containing inner barriers to be included. By default None - Returns ------- GeoDataFrame From c9a79e3a4a39b1b721d6156341b1d8f2852f1c21 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Wed, 8 Oct 2025 16:38:14 +0200 Subject: [PATCH 15/19] Add intersection for inner barrier generation --- momepy/elements.py | 3 +- momepy/tests/test_elements.py | 131 +++++++++++++++++++++++++++++++++- 2 files changed, 132 insertions(+), 2 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index bc44f698..1ad07970 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -702,7 +702,8 @@ def _get_inner_barriers(enclosure, barriers): inner_barriers = gpd.clip(barriers, enclosure) # Only keep the geometry which is within the enclosure - inner_barriers = inner_barriers[inner_barriers.within(enclosure)] + inner_barriers = inner_barriers[inner_barriers.intersects(enclosure) | + inner_barriers.within(enclosure)] return inner_barriers diff --git a/momepy/tests/test_elements.py b/momepy/tests/test_elements.py index 44d20dda..5ec0fb12 100644 --- a/momepy/tests/test_elements.py +++ b/momepy/tests/test_elements.py @@ -10,9 +10,10 @@ from packaging.version import Version from pandas.testing import assert_index_equal from shapely import LineString, affinity -from shapely.geometry import MultiPoint, Polygon, box +from shapely.geometry import MultiPoint, MultiPolygon, Polygon, box import momepy as mm +import momepy.elements as elements SHPLY_GE_210 = Version(shapely.__version__) >= Version("2.1.0") LPS_G_4_13_0 = Version(libpysal.__version__) > Version("4.13.0") @@ -175,6 +176,134 @@ def test_enclosed_tessellation(self): assert (self.df_buildings.index.isin(tessellation_custom_index.index)).all() assert tessellation_custom_index.enclosure_index.isin(custom_index.index).all() + tessellation_inner_barrier = mm.enclosed_tessellation( + self.df_buildings, + self.enclosures.geometry, + simplify=False, + inner_barriers=self.df_streets, + cell_size=5, + neighbor_mode="neumann", + n_jobs=1, + ) + assert set(tessellation_inner_barrier.geom_type.unique()) <= {"Polygon", "MultiPolygon"} + assert tessellation_inner_barrier.crs == self.df_buildings.crs + assert len(tessellation_inner_barrier) == len(tessellation) + assert tessellation_inner_barrier.enclosure_index.isin( + self.enclosures.index + ).all() + + def test_enclosed_tessellation_inner_barrier_cellular(self): + crs = self.df_buildings.crs + buildings = gpd.GeoDataFrame( + {"uID": [0, 1, 2]}, + geometry=[ + box(10, 10, 30, 30), + box(70, 10, 90, 30), + shapely.Point(15, 65), + ], + crs=crs, + ) + multipolygon = shapely.MultiPolygon([box(0, 0, 40, 80), box(60, 0, 100, 80)]) + enclosures = gpd.GeoDataFrame({"geometry": [multipolygon]}, crs=crs) + inner_barriers = gpd.GeoDataFrame( + { + "geometry": [ + LineString([(50, -10), (50, 90)]), + LineString([(0, 40), (100, 40)]), + ] + }, + crs=crs, + ) + + tess = mm.enclosed_tessellation( + buildings, + enclosures, + simplify=False, + inner_barriers=inner_barriers, + cell_size=1, + neighbor_mode="neumann", + threshold=None, + n_jobs=1, + ) + + assert set(tess.geom_type.unique()) <= {"Polygon", "MultiPolygon"} + assert tess.enclosure_index.unique().tolist() == [0] + assert set(buildings.index).issubset(tess.index) + + point_barriers = gpd.GeoDataFrame( + {"geometry": [shapely.Point(20, 40), shapely.Point(80, 40)]}, crs=crs + ) + tess_point = mm.enclosed_tessellation( + buildings, + enclosures, + simplify=False, + inner_barriers=point_barriers, + cell_size=1, + neighbor_mode="moore", + threshold=None, + n_jobs=1, + ) + assert set(tess_point.geom_type.unique()) <= {"Polygon", "MultiPolygon"} + assert set(buildings.index).issubset(tess_point.index) + + def test_enclosed_tessellation_invalid_enclosure_geometry(self): + crs = self.df_buildings.crs + buildings = gpd.GeoDataFrame( + {"uID": [1, 2]}, + geometry=[box(10, 10, 30, 30), box(70, 10, 90, 30)], + crs=crs, + ) + invalid_enclosures = gpd.GeoDataFrame( + {"geometry": [LineString([(0, 20), (100, 20)])]}, crs=crs + ) + with pytest.raises(ValueError, match="Enclosure must be a Polygon"): + mm.enclosed_tessellation( + buildings, + invalid_enclosures, + simplify=False, + inner_barriers=self.df_streets.head(1), + cell_size=1, + threshold=None, + n_jobs=1, + ) + + def test_enclosed_tessellation_invalid_neighbor_mode(self): + crs = self.df_buildings.crs + buildings = gpd.GeoDataFrame( + {"uID": [1, 2]}, + geometry=[box(10, 10, 30, 30), box(70, 10, 90, 30)], + crs=crs, + ) + enclosures = gpd.GeoDataFrame( + { + "geometry": [ + shapely.MultiPolygon([box(0, 0, 40, 80), box(60, 0, 100, 80)]) + ] + }, + crs=crs, + ) + with pytest.raises(ValueError, match="Invalid neighbor_mode"): + mm.enclosed_tessellation( + buildings, + enclosures, + simplify=False, + inner_barriers=self.df_streets.head(1), + cell_size=1, + neighbor_mode="invalid", + threshold=None, + n_jobs=1, + ) + + def test_morphological_and_enclosed_tessellation_require_shapely(self, monkeypatch): + monkeypatch.setattr(elements, "SHPLY_GE_210", False) + with pytest.raises(ImportError): + mm.morphological_tessellation(self.df_buildings.head(1)) + with pytest.raises(ImportError): + mm.enclosed_tessellation( + self.df_buildings.head(1), + self.enclosures.geometry, + ) + def test_verify_tessellation(self): df = self.df_buildings b = df.total_bounds From 374e226f0139fbc0e60d2a3d7bbefdeb44806242 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Wed, 8 Oct 2025 17:43:18 +0200 Subject: [PATCH 16/19] Refactor _voronoi_by_ca to improve code coverage --- momepy/elements.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index 1ad07970..5405f2fc 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -501,30 +501,31 @@ def _voronoi_by_ca( A GeoDataFrame representing the aggregated Voronoi tessellation, clipped by barriers. """ + geom_type = barrier_geoms.geom_type + if geom_type not in {"Polygon", "MultiPolygon"}: + raise ValueError("Enclosure must be a Polygon or MultiPolygon") + # Get inner barriers as intersection or containment of the barrier_geoms inner_barriers = _get_inner_barriers(barrier_geoms, barriers_for_inner) # Handle barrier_geoms if it is a Polygon or MultiPolygon - if barrier_geoms.geom_type == "Polygon": + if geom_type == "Polygon": # Take buffer of polygon and extract its exterior boundary (10 cells) barrier_geoms_buffered = GeoSeries( [barrier_geoms.buffer(10 * cell_size).exterior], crs=seed_geoms.crs ) barrier_geoms = GeoSeries([barrier_geoms], crs=seed_geoms.crs) - elif barrier_geoms.geom_type == "MultiPolygon": + elif geom_type == "MultiPolygon": # Process each polygon: take buffer then exterior boundary # (to ensure there's no gap between enclosures) + parts = list(shapely.get_parts(barrier_geoms)) + exteriors = shapely.get_exterior_ring(parts) + buffered_exteriors = shapely.buffer(exteriors, 10 * cell_size) barrier_geoms_buffered = GeoSeries( - shapely.buffer( - shapely.get_exterior(shapely.get_parts(barrier_geoms)), 10 * cell_size - ), - crs=seed_geoms.crs, + list(buffered_exteriors), crs=seed_geoms.crs ) - barrier_geoms = GeoSeries(barrier_geoms, crs=seed_geoms.crs) - - else: - raise ValueError("Enclosure must be a Polygon or MultiPolygon") + barrier_geoms = GeoSeries(parts, crs=seed_geoms.crs) outer_union = barrier_geoms_buffered.union_all() From 63f11a3f1ab87246f739b046d46cd9b3a93d2d63 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Wed, 8 Oct 2025 18:29:40 +0200 Subject: [PATCH 17/19] Refactor masking / barrier generation to improve code coverage --- momepy/elements.py | 94 ++++++++++++++++------------------------------ 1 file changed, 33 insertions(+), 61 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index 5405f2fc..be450e54 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -528,22 +528,18 @@ def _voronoi_by_ca( barrier_geoms = GeoSeries(parts, crs=seed_geoms.crs) outer_union = barrier_geoms_buffered.union_all() + inner_union = ( + inner_barriers.union_all() + if inner_barriers is not None and not inner_barriers.empty + else None + ) - # Compute inner barriers union if available - if inner_barriers is not None and not inner_barriers.empty: - inner_union = inner_barriers.union_all() - else: - inner_union = None - - # Combine outer barrier with inner barriers - if outer_union and inner_union: - prep_barrier = shapely.union(outer_union, inner_union) - elif outer_union: - prep_barrier = outer_union - elif inner_union: - prep_barrier = inner_union - else: - prep_barrier = None + geoms_to_union = [ + geom + for geom in (outer_union, inner_union) + if geom is not None and not geom.is_empty + ] + prep_barrier = shapely.union_all(geoms_to_union) if geoms_to_union else None # Compute grid bounds origin, grid_width, grid_height = _get_grid_bounds( @@ -564,12 +560,11 @@ def _voronoi_by_ca( ) # Identify barrier cells in the grid - if prep_barrier is not None: + barrier_mask = np.zeros((grid_height, grid_width), dtype=bool) + if prep_barrier is not None and not prep_barrier.is_empty: barrier_mask = shapely.intersects(cell_polys_array, prep_barrier).reshape( grid_height, grid_width ) - else: - barrier_mask = np.zeros((grid_height, grid_width), dtype=bool) states[barrier_mask] = CellState.BARRIER.value # Seed the grid with seed geometries. @@ -647,37 +642,31 @@ def _voronoi_by_ca( # Clip by barriers if barrier_geoms is not None and (not barrier_geoms.empty): - # Create a union of the barrier geometries. barrier_union = barrier_geoms.union_all() - # If the barrier union is not a polygon (e.g., it's a MultiLineString), - # polygonize it. - if not isinstance( - barrier_union, shapely.geometry.Polygon | shapely.geometry.MultiPolygon - ): - barrier_polys = shapely.get_parts(shapely.polygonize(barrier_union)) - if not barrier_polys.empty: - barrier_union = shapely.union_all(barrier_polys) - # Clip each polygon in the grid using the barrier boundary. - grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) + polygonized = [] + if barrier_union is not None and not barrier_union.is_empty: + polygonized = [ + geom + for geom in shapely.get_parts(shapely.polygonize([barrier_union])) + if geom.geom_type in {"Polygon", "MultiPolygon"} + ] + if polygonized: + barrier_union = shapely.union_all(polygonized) + if barrier_union is not None and not barrier_union.is_empty: + grid_gdf["geometry"] = grid_gdf["geometry"].intersection(barrier_union) # cleanup possible collections def sanitize_geometry(geom): - if geom.geom_type in ["Polygon", "MultiPolygon"]: - return geom - elif geom.geom_type == "GeometryCollection": - parts = shapely.get_parts(geom) - valid_parts = [ - p for p in parts if p.geom_type in ["Polygon", "MultiPolygon"] + parts = ( + [ + part + for part in shapely.get_parts(geom) + if part.geom_type in {"Polygon", "MultiPolygon"} ] - if len(valid_parts) == 0: - return None - elif len(valid_parts) == 1: - return valid_parts[0] - else: - return shapely.MultiPolygon(valid_parts) - else: - # Drop points, lines, etc. - return None + if geom is not None + else [] + ) + return shapely.union_all(parts) if parts else None grid_gdf["geometry"] = grid_gdf["geometry"].apply(sanitize_geometry) grid_gdf = grid_gdf.dropna(subset=["geometry"]) @@ -732,23 +721,6 @@ class CellState(Enum): FRONTIER = -4 -def _get_cell_polygon( - x_idx: int, y_idx: int, cell_size: float = 1.0, origin: tuple[float, float] = (0, 0) -) -> shapely.geometry.Polygon: - """ - Generate a grid cell polygon based on the given indices, cell size, and origin. - """ - ox, oy = origin - return shapely.Polygon( - [ - (ox + x_idx * cell_size, oy + y_idx * cell_size), - (ox + (x_idx + 1) * cell_size, oy + y_idx * cell_size), - (ox + (x_idx + 1) * cell_size, oy + (y_idx + 1) * cell_size), - (ox + x_idx * cell_size, oy + (y_idx + 1) * cell_size), - ] - ) - - def _get_neighbors( x: int, y: int, max_x: int, max_y: int, mode: str = "moore" ) -> list[tuple[int, int]]: From 7c588b02844b6f9cbe7c923f3ab891f4e882f6e9 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Wed, 8 Oct 2025 18:29:51 +0200 Subject: [PATCH 18/19] Add cul-de-sac test cases --- momepy/tests/test_elements.py | 78 ++++++++++++++++++++++++++++------- 1 file changed, 62 insertions(+), 16 deletions(-) diff --git a/momepy/tests/test_elements.py b/momepy/tests/test_elements.py index 5ec0fb12..fcf66c57 100644 --- a/momepy/tests/test_elements.py +++ b/momepy/tests/test_elements.py @@ -193,26 +193,57 @@ def test_enclosed_tessellation(self): ).all() def test_enclosed_tessellation_inner_barrier_cellular(self): - crs = self.df_buildings.crs + # Define sample building geometries (including Points to cover Point handling). + blg_polygons = [ + shapely.geometry.Polygon([(15, 32), (35, 32), (35, 38), (15, 38)]), + shapely.geometry.Polygon([(15, 22), (35, 22), (35, 28), (15, 28)]), + shapely.geometry.Polygon([(15, 92), (35, 92), (35, 98), (15, 98)]), + shapely.geometry.Polygon([(15, 82), (35, 82), (35, 88), (15, 88)]), + shapely.geometry.Polygon([(45, 62), (65, 62), (65, 68), (45, 68)]), + shapely.geometry.Polygon([(45, 52), (65, 52), (65, 58), (45, 58)]), + shapely.geometry.Point(25, 50), # Point building to test Point handling + shapely.geometry.Point(55, 80), # Another Point building + ] buildings = gpd.GeoDataFrame( - {"uID": [0, 1, 2]}, - geometry=[ - box(10, 10, 30, 30), - box(70, 10, 90, 30), - shapely.Point(15, 65), - ], - crs=crs, + { + "building_id": list(range(1, len(blg_polygons) + 1)), + "geometry": blg_polygons, + }, + crs="EPSG:3857", ) - multipolygon = shapely.MultiPolygon([box(0, 0, 40, 80), box(60, 0, 100, 80)]) - enclosures = gpd.GeoDataFrame({"geometry": [multipolygon]}, crs=crs) + + # Define sample barrier geometries. + barrier_geoms = [ + shapely.geometry.LineString([(0, 0), (80, 0)]), + shapely.geometry.LineString([(80, 0), (80, 120)]), + shapely.geometry.LineString([(80, 120), (0, 120)]), + shapely.geometry.LineString([(0, 120), (0, 0)]), + shapely.geometry.LineString([(40, 0), (40, 110)]), + shapely.geometry.LineString([(10, 30), (40, 30)]), + shapely.geometry.LineString([(10, 90), (40, 90)]), + shapely.geometry.LineString([(40, 60), (70, 60)]), + ] + inner_barriers = gpd.GeoDataFrame( { - "geometry": [ - LineString([(50, -10), (50, 90)]), - LineString([(0, 40), (100, 40)]), - ] + "name": [ + "Bottom Edge", + "Right Edge", + "Top Edge", + "Left Edge", + "Main Vertical", + "Left Cul-de-Sac (Bottom)", + "Left Cul-de-Sac (Top)", + "Right Cul-de-Sac (Middle)", + ], + "geometry": barrier_geoms, }, - crs=crs, + crs="EPSG:3857", + ) + + # Create enclosure from the barrier boundaries + enclosures = gpd.GeoDataFrame( + {"geometry": [box(0, 0, 80, 120)]}, crs="EPSG:3857" ) tess = mm.enclosed_tessellation( @@ -231,7 +262,8 @@ def test_enclosed_tessellation_inner_barrier_cellular(self): assert set(buildings.index).issubset(tess.index) point_barriers = gpd.GeoDataFrame( - {"geometry": [shapely.Point(20, 40), shapely.Point(80, 40)]}, crs=crs + {"geometry": [shapely.Point(20, 40), shapely.Point(80, 40)]}, + crs="EPSG:3857", ) tess_point = mm.enclosed_tessellation( buildings, @@ -246,6 +278,20 @@ def test_enclosed_tessellation_inner_barrier_cellular(self): assert set(tess_point.geom_type.unique()) <= {"Polygon", "MultiPolygon"} assert set(buildings.index).issubset(tess_point.index) + empty_barriers = gpd.GeoDataFrame(geometry=[], crs="EPSG:3857") + tess_empty = mm.enclosed_tessellation( + buildings, + enclosures, + simplify=False, + inner_barriers=empty_barriers, + cell_size=1, + neighbor_mode="moore", + threshold=None, + n_jobs=1, + ) + assert set(tess_empty.geom_type.unique()) <= {"Polygon", "MultiPolygon"} + assert set(buildings.index).issubset(tess_empty.index) + def test_enclosed_tessellation_invalid_enclosure_geometry(self): crs = self.df_buildings.crs buildings = gpd.GeoDataFrame( From 7d534c7158a1eff83f6c2a86d6797156595fa6d8 Mon Sep 17 00:00:00 2001 From: Yuta Sato Date: Wed, 8 Oct 2025 18:38:47 +0200 Subject: [PATCH 19/19] Update docstrings --- momepy/elements.py | 32 ++++++++++++++++++++++---------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/momepy/elements.py b/momepy/elements.py index be450e54..6eeced5f 100644 --- a/momepy/elements.py +++ b/momepy/elements.py @@ -169,6 +169,11 @@ def enclosed_tessellation( street network, railway). Original morphological tessellation is used under the hood to partition each enclosure. + When ``inner_barriers`` are provided, the tessellation is derived using a Cellular + Automata implementation that recognizes dangling barriers (such as dead-end streets) + as valid limits of cell growth. This is more computationally intensive but handles + complex barrier configurations more accurately. + Tessellation requires data of relatively high level of precision and there are three particular patterns causing issues: @@ -213,19 +218,26 @@ def enclosed_tessellation( The number of jobs to run in parallel. -1 means using all available cores. By default -1 inner_barriers: GeoSeries | GeoDataFrame, optional - Barriers that should be included in the tessellation process. By passing - inner barriers, tessellation will be derived using less performant - Cellular Automata implememtation but it will recognise dangling barries as valid - limits of the cell growth. See the user guide for details. By default None. + Barriers that should be included in the tessellation process. When provided, + tessellation will be derived using a Cellular Automata implementation that + recognizes dangling barriers (such as dead-end streets or cul-de-sacs) as valid + limits of cell growth. This is more computationally intensive than the default + Voronoi-based approach but can handle inner barriers. By default None. cell_size : float, optional - Grid cell size when ``inner_barriers`` is not None. Otherwise ignored. - By default 1.0 + Grid cell size for the Cellular Automata implementation when ``inner_barriers`` + is not None. Smaller values provide higher resolution but increase computational + cost. This parameter controls the spatial granularity of the tessellation grid. + When ``inner_barriers`` is None, this parameter is ignored. By default 1.0 neighbor_mode : str, optional - Choice of neighbor connectivity ('moore' or 'neumann') when ``inner_barriers`` - is not None. Otherwise ignored. By default 'moore'. + Choice of neighbor connectivity for the Cellular Automata implementation when + ``inner_barriers`` is not None. Options are 'moore' (8-connected, including + diagonal neighbors) or 'neumann' (4-connected, only orthogonal neighbors). + When ``inner_barriers`` is None, this parameter is ignored. By default 'moore'. **kwargs - Additional keyword arguments pased to libpysal.cg.voronoi_frames, such as - ``grid_size``. + Additional keyword arguments passed to :func:`libpysal.cg.voronoi_frames` when + ``inner_barriers`` is None, such as ``grid_size``. These arguments are ignored + when ``inner_barriers`` is provided and the Cellular Automata implementation is + used. Warnings --------