Merge pull request #645 from ungarj/singlepart_multipart_filtering

Singlepart multipart filtering

mapchete/formats/default/vector_file.py

@@ -289,9 +289,11 @@ def _read_from_cache(self, validity_check=True, clip_to_crs_bounds=False):
         if checked not in self._cache:
             self._cache[checked] = list(
                 read_vector_window(
-                    self._in_memory_features
-                    if self._memory_cache_active
-                    else self.path,
+                    (
+                        self._in_memory_features
+                        if self._memory_cache_active
+                        else self.path
+                    ),
                     self.tile,
                     validity_check=validity_check,
                     clip_to_crs_bounds=clip_to_crs_bounds,

mapchete/geometry/filter.py

@@ -1,4 +1,4 @@
-from typing import Generator, Union
+from typing import Generator, Type, Union
 
 from mapchete.errors import GeometryTypeError
 from mapchete.geometry.types import (
@@ -39,23 +39,27 @@ def multipart_to_singleparts(
 
 def is_type(
     geometry: Geometry,
-    target_type: Union[str, Geometry],
+    target_type: Union[str, Type[Geometry]],
     allow_multipart: bool = True,
 ) -> bool:
+    """
+    Checks whether geometry type is in alignment with target type.
+    """
     target_type = get_geometry_type(target_type)
     if isinstance(geometry, target_type):
         return True
     elif isinstance(geometry, GeometryCollection):
         return False
 
+    # SinglePart is MultiPart or MultiPart is SinglePart
     if allow_multipart:
         return isinstance(geometry, get_multipart_type(target_type))
     return False
 
 
 def filter_by_geometry_type(
     geometry: Geometry,
-    target_type: Union[str, Geometry],
+    target_type: Union[str, Type[Geometry]],
     allow_multipart: bool = True,
 ):
     """Yields geometries only if they match the target type.
@@ -64,12 +68,16 @@ def filter_by_geometry_type(
     If allow_multipart is set to False, multipart geometries are broken down into their
     subgeometries. If set to True, a MultiPoint will be yielded if the
     """
-    target_type = get_geometry_type(target_type)
-    if is_type(geometry, target_type=target_type, allow_multipart=allow_multipart):
+    target_geometry_type = get_geometry_type(target_type)
+    if is_type(
+        geometry, target_type=target_geometry_type, allow_multipart=allow_multipart
+    ):
         yield geometry
 
     elif isinstance(geometry, MultipartGeometry):
         for subgeometry in multipart_to_singleparts(geometry):
             yield from filter_by_geometry_type(
-                subgeometry, target_type=target_type, allow_multipart=allow_multipart
+                subgeometry,
+                target_type=target_geometry_type,
+                allow_multipart=allow_multipart,
             )

mapchete/geometry/types.py

@@ -53,7 +53,7 @@ class GeoInterface(Protocol):
 CoordArrays = Tuple[Iterable[float], Iterable[float]]
 
 
-def get_multipart_type(geometry: Geometry) -> MultipartGeometry:
+def get_multipart_type(geometry_type: Union[Type[Geometry], str]) -> MultipartGeometry:
     try:
         return {
             Point: MultiPoint,
@@ -62,14 +62,38 @@ def get_multipart_type(geometry: Geometry) -> MultipartGeometry:
             MultiPoint: MultiPoint,
             MultiLineString: MultiLineString,
             MultiPolygon: MultiPolygon,
-        }[geometry]
+        }[
+            get_geometry_type(geometry_type)
+        ]  # type: ignore
     except KeyError:
         raise GeometryTypeError(
-            f"geometry type {geometry.type} has no corresponding multipart type"
+            f"geometry type {geometry_type} has no corresponding multipart type"
         )
 
 
-def get_geometry_type(geometry_type: Union[Type[Geometry], str]) -> Type[Geometry]:
+def get_singlepart_type(
+    geometry_type: Union[Type[Geometry], str]
+) -> SinglepartGeometry:
+    try:
+        return {
+            Point: Point,
+            LineString: LineString,
+            Polygon: Polygon,
+            MultiPoint: Point,
+            MultiLineString: LineString,
+            MultiPolygon: Polygon,
+        }[
+            get_geometry_type(geometry_type)
+        ]  # type: ignore
+    except KeyError:  # pragma: no cover
+        raise GeometryTypeError(
+            f"geometry type {geometry_type} has no corresponding multipart type"
+        )
+
+
+def get_geometry_type(
+    geometry_type: Union[Type[Geometry], str, dict, Geometry]
+) -> Geometry:
     if isinstance(geometry_type, str):
         try:
             return {
@@ -87,5 +111,5 @@ def get_geometry_type(geometry_type: Union[Type[Geometry], str]) -> Type[Geometr
                 f"geometry type cannot be determined from {geometry_type}"
             )
     elif issubclass(geometry_type, Geometry):
-        return geometry_type
+        return geometry_type  # type: ignore
     raise GeometryTypeError(f"geometry type cannot be determined from {geometry_type}")

mapchete/io/vector.py

@@ -25,7 +25,7 @@
     segmentize_geometry,
     to_shape,
 )
-from mapchete.geometry.types import get_geometry_type
+from mapchete.geometry.types import get_geometry_type, get_singlepart_type
 from mapchete.io import copy
 from mapchete.path import MPath, fs_from_path
 from mapchete.settings import IORetrySettings
@@ -323,7 +323,7 @@ def _get_reprojected_features(
                 clipped_geom = original_geom.intersection(dst_bbox)
                 for checked_geom in filter_by_geometry_type(
                     clipped_geom,
-                    original_geom.geom_type,
+                    get_singlepart_type(original_geom.geom_type),
                 ):
                     # reproject each feature to tile CRS
                     reprojected_geom = reproject_geometry(

mapchete/tile.py

@@ -1,4 +1,5 @@
 """Mapchtete handling tiles."""
+
 from __future__ import annotations
 
 import logging
@@ -9,9 +10,11 @@
 from affine import Affine
 from rasterio.enums import Resampling
 from rasterio.features import rasterize, shapes
+from rasterio.transform import from_bounds
 from rasterio.warp import reproject
 from shapely import clip_by_rect
 from shapely.geometry import box, shape
+from shapely.geometry.base import BaseGeometry
 from shapely.ops import unary_union
 from tilematrix import Tile, TilePyramid
 from tilematrix._conf import ROUND
@@ -527,21 +530,46 @@ def _count_cells(pyramid, geometry, minzoom, maxzoom):
     return int(count)
 
 
-def snap_geometry_to_tiles(geometry=None, pyramid=None, zoom=None):
+def snap_geometry_to_tiles(
+    geometry: BaseGeometry, pyramid: BufferedTilePyramid, zoom: int
+) -> BaseGeometry:
     if geometry.is_empty:
         return geometry
     # calculate everything using an unbuffered pyramid, because otherwise the Affine
     # object cannot be calculated
     unbuffered_pyramid = BufferedTilePyramid.from_dict(
         dict(pyramid.to_dict(), pixelbuffer=0)
     )
-    transform = unbuffered_pyramid.matrix_affine(zoom)
+
+    # use subset because otherwise memory usage can crash the program
+    left, bottom, right, top = geometry.bounds
+    # clip geometry bounds to pyramid bounds
+    left = max([left, unbuffered_pyramid.bounds.left])
+    bottom = max([bottom, unbuffered_pyramid.bounds.bottom])
+    right = min([right, unbuffered_pyramid.bounds.right])
+    top = min([top, unbuffered_pyramid.bounds.top])
+    lb_tile = unbuffered_pyramid.tile_from_xy(left, bottom, zoom, on_edge_use="rt")
+    rt_tile = unbuffered_pyramid.tile_from_xy(right, top, zoom, on_edge_use="lb")
+    snapped_left, snapped_south, snapped_east, snapped_north = (
+        lb_tile.bounds.left,
+        lb_tile.bounds.bottom,
+        rt_tile.bounds.right,
+        rt_tile.bounds.top,
+    )
+    width = abs(rt_tile.col - lb_tile.col) + 1
+    height = abs(rt_tile.row - lb_tile.row) + 1
+    out_shape = (height, width)
+    transform = from_bounds(
+        west=snapped_left,
+        south=snapped_south,
+        east=snapped_east,
+        north=snapped_north,
+        width=width,
+        height=height,
+    )
     raster = rasterize(
         [(geometry, 1)],
-        out_shape=(
-            unbuffered_pyramid.matrix_height(zoom),
-            unbuffered_pyramid.matrix_width(zoom),
-        ),
+        out_shape=out_shape,
         fill=0,
         transform=transform,
         dtype=np.uint8,

pyproject.toml

@@ -32,7 +32,7 @@ dependencies = [
     "geojson-pydantic",    
     "importlib-metadata",
     "importlib-resources",
-    "numpy>=1.16",
+    "numpy>=1.16,<=2.0.1",
     "oyaml",
     "pydantic>=2.3.0",
     "pydantic_settings>=2.0.0",