diff --git a/mapchete/formats/default/vector_file.py b/mapchete/formats/default/vector_file.py index f68737c1..203a1d12 100644 --- a/mapchete/formats/default/vector_file.py +++ b/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, diff --git a/mapchete/geometry/filter.py b/mapchete/geometry/filter.py index 11480c80..574c8845 100644 --- a/mapchete/geometry/filter.py +++ b/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,15 +39,19 @@ 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 @@ -55,7 +59,7 @@ def is_type( 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, ) diff --git a/mapchete/geometry/types.py b/mapchete/geometry/types.py index 91b900e2..6980f21f 100644 --- a/mapchete/geometry/types.py +++ b/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}") diff --git a/mapchete/io/vector.py b/mapchete/io/vector.py index 29391d0e..c55d633e 100644 --- a/mapchete/io/vector.py +++ b/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( diff --git a/mapchete/tile.py b/mapchete/tile.py index c13852e7..f01753cf 100644 --- a/mapchete/tile.py +++ b/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,7 +530,9 @@ 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 @@ -535,13 +540,36 @@ def snap_geometry_to_tiles(geometry=None, pyramid=None, zoom=None): 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, diff --git a/pyproject.toml b/pyproject.toml index e17c6450..856a6fe7 100644 --- a/pyproject.toml +++ b/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",