From 86b524a5f4063f2f35db714d75b93ca1b89b0f19 Mon Sep 17 00:00:00 2001 From: Joachim Ungar Date: Mon, 4 Nov 2024 10:22:33 +0100 Subject: [PATCH] make sure CRS bounds have latlon CRS set; fix pyproj.CRS parsing issue --- mapchete/geometry/repair.py | 22 ++++++++++++---------- mapchete/geometry/reproject.py | 26 +++++++++++++++++++------- 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/mapchete/geometry/repair.py b/mapchete/geometry/repair.py index 7d1291aa..3e8ba48e 100644 --- a/mapchete/geometry/repair.py +++ b/mapchete/geometry/repair.py @@ -4,16 +4,18 @@ from mapchete.geometry.types import Geometry, MultiPolygon, Polygon -def repair(geometry: Geometry) -> Geometry: - repaired = ( - geometry.buffer(0) - if isinstance(geometry, (Polygon, MultiPolygon)) - else geometry - ) - if repaired.is_valid: - return repaired +def repair(geometry: Geometry, normalize: bool = True) -> Geometry: + if isinstance(geometry, (Polygon, MultiPolygon)): + out = geometry.buffer(0) + else: + out = geometry + + if normalize: + out = out.normalize() + + if out.is_valid: + return out else: raise TopologicalError( - "geometry is invalid (%s) and cannot be repaired" - % explain_validity(repaired) + f"geometry is invalid ({explain_validity(out)}) and cannot be repaired" ) diff --git a/mapchete/geometry/reproject.py b/mapchete/geometry/reproject.py index 6665564f..9dba36a2 100644 --- a/mapchete/geometry/reproject.py +++ b/mapchete/geometry/reproject.py @@ -3,6 +3,7 @@ import fiona import pyproj +from pyproj.exceptions import CRSError from fiona.transform import transform_geom from rasterio.crs import CRS from shapely.geometry import mapping, shape @@ -29,11 +30,11 @@ CRS_BOUNDS = { # http://spatialreference.org/ref/epsg/wgs-84/ - CRS.from_epsg(4326): Bounds(-180.0, -90.0, 180.0, 90.0), + CRS.from_epsg(4326): Bounds(-180.0, -90.0, 180.0, 90.0, crs=LATLON_CRS), # unknown source - CRS.from_epsg(3857): Bounds(-180.0, -85.0511, 180.0, 85.0511), + CRS.from_epsg(3857): Bounds(-180.0, -85.0511, 180.0, 85.0511, crs=LATLON_CRS), # http://spatialreference.org/ref/epsg/3035/ - CRS.from_epsg(3035): Bounds(-10.6700, 34.5000, 31.5500, 71.0500), + CRS.from_epsg(3035): Bounds(-10.6700, 34.5000, 31.5500, 71.0500, crs=LATLON_CRS), } @@ -42,15 +43,26 @@ def get_crs_bounds(crs: CRS) -> Bounds: # get bounds from known CRSes return CRS_BOUNDS[crs] except KeyError: + logger.debug("try to determine CRS bounds using pyproj ...") # try to get bounds using pyproj - area_of_use = pyproj.CRS(crs.to_epsg()).area_of_use - if area_of_use: - return Bounds.from_inp(area_of_use.bounds) + try: + # on UTM CRS, the area_of_use is None if pyproj.CRS is initialized with CRS.to_proj4(), thus + # prefer using CRS.to_epsg() and only use CRS.to_proj4() as backup + pyproj_crs = ( + pyproj.CRS(crs.to_epsg()) + if crs.is_epsg_code + else pyproj.CRS(crs.to_proj4()) + ) + if pyproj_crs.area_of_use: + return Bounds.from_inp(pyproj_crs.area_of_use.bounds, crs=LATLON_CRS) + except CRSError as exc: + logger.debug(exc) + pass raise ValueError(f"bounds of CRS {crs} could not be determined") def crs_is_epsg_4326(crs: CRS) -> bool: - return crs.is_epsg_code and crs.get("init") == "epsg:4326" + return crs == LATLON_CRS def reproject_geometry(