121 lines
4.8 KiB
Python
121 lines
4.8 KiB
Python
from __future__ import annotations
|
|
|
|
from math import isfinite
|
|
from numbers import Real
|
|
from typing import Any
|
|
|
|
from pyproj import CRS, Transformer
|
|
from shapely import force_2d, get_coordinates
|
|
from shapely.geometry import MultiPolygon, box, shape
|
|
from shapely.ops import transform
|
|
|
|
|
|
# This is a deliberately broad guard envelope around Belgium and the Belgian
|
|
# North Sea. Exact legal/regional clipping remains the responsibility of the
|
|
# persisted coverage Areas; this boundary prevents an AOI with valid-looking
|
|
# but globally misplaced coordinates from entering the workbench.
|
|
BELGIUM_AND_NORTH_SEA_GUARD_BOUNDS = (1.5, 48.5, 7.5, 52.5)
|
|
|
|
|
|
def _raw_coordinates_are_finite(value: Any) -> bool:
|
|
if isinstance(value, (list, tuple)):
|
|
return bool(value) and all(_raw_coordinates_are_finite(item) for item in value)
|
|
return isinstance(value, Real) and not isinstance(value, bool) and isfinite(float(value))
|
|
|
|
|
|
def normalize_to_multipolygon(raw_geometry: dict[str, Any]) -> MultiPolygon:
|
|
if isinstance(raw_geometry, dict) and "coordinates" in raw_geometry:
|
|
if not _raw_coordinates_are_finite(raw_geometry["coordinates"]):
|
|
raise ValueError("Geometry coordinates must be finite numbers")
|
|
try:
|
|
geom = force_2d(shape(raw_geometry))
|
|
except Exception as exc:
|
|
raise ValueError("Geometry is not valid GeoJSON") from exc
|
|
coordinates = get_coordinates(geom, include_z=False)
|
|
if coordinates.size == 0 or not all(isfinite(float(value)) for row in coordinates for value in row):
|
|
raise ValueError("Geometry coordinates must be finite numbers")
|
|
if geom.is_empty:
|
|
raise ValueError("Geometry is empty")
|
|
if not geom.is_valid:
|
|
raise ValueError("Geometry is invalid")
|
|
|
|
if geom.geom_type == "Polygon":
|
|
return MultiPolygon([geom])
|
|
if geom.geom_type == "MultiPolygon":
|
|
return MultiPolygon(geom.geoms)
|
|
|
|
raise ValueError("Only Polygon or MultiPolygon geometries are accepted")
|
|
|
|
|
|
def normalize_area_to_epsg4326(
|
|
raw_geometry: dict[str, Any],
|
|
source_crs: str,
|
|
) -> tuple[MultiPolygon, str]:
|
|
"""Validate an AOI and normalize its declared CRS to canonical WGS84.
|
|
|
|
The returned CRS string preserves the caller's declaration for provenance;
|
|
the returned geometry is always finite, polygonal and stored as EPSG:4326.
|
|
"""
|
|
|
|
declared_crs = str(source_crs or "").strip()
|
|
if not declared_crs:
|
|
raise ValueError("Area CRS is required")
|
|
try:
|
|
parsed_crs = CRS.from_user_input(declared_crs)
|
|
except Exception as exc:
|
|
raise ValueError("Area CRS is unknown or invalid") from exc
|
|
if not (parsed_crs.is_geographic or parsed_crs.is_projected):
|
|
raise ValueError("Area CRS must be a geographic or projected two-dimensional CRS")
|
|
if len(parsed_crs.axis_info) != 2:
|
|
raise ValueError("Area CRS must have exactly two spatial axes")
|
|
|
|
geometry = normalize_to_multipolygon(raw_geometry)
|
|
target_crs = CRS.from_epsg(4326)
|
|
if not parsed_crs.equals(target_crs):
|
|
try:
|
|
transformer = Transformer.from_crs(parsed_crs, target_crs, always_xy=True)
|
|
geometry = normalize_to_multipolygon(
|
|
transform(transformer.transform, geometry).__geo_interface__
|
|
)
|
|
except ValueError:
|
|
raise
|
|
except Exception as exc:
|
|
raise ValueError("Area geometry could not be transformed to EPSG:4326") from exc
|
|
|
|
min_x, min_y, max_x, max_y = geometry.bounds
|
|
if not all(isfinite(value) for value in (min_x, min_y, max_x, max_y)):
|
|
raise ValueError("Transformed area geometry contains non-finite coordinates")
|
|
world_bounds = (-180.0, -90.0, 180.0, 90.0)
|
|
if min_x < world_bounds[0] or min_y < world_bounds[1] or max_x > world_bounds[2] or max_y > world_bounds[3]:
|
|
raise ValueError("Transformed area geometry falls outside the EPSG:4326 coordinate domain")
|
|
guard = box(*BELGIUM_AND_NORTH_SEA_GUARD_BOUNDS)
|
|
if not guard.intersects(geometry):
|
|
raise ValueError("Area geometry falls outside Belgium and the Belgian North Sea workbench domain")
|
|
if not guard.covers(geometry):
|
|
raise ValueError("Area geometry must remain within the Belgium and Belgian North Sea workbench domain")
|
|
return geometry, declared_crs
|
|
|
|
|
|
def area_bounds_multipolygon(geom: MultiPolygon):
|
|
return {
|
|
"min_x": float(geom.bounds[0]),
|
|
"min_y": float(geom.bounds[1]),
|
|
"max_x": float(geom.bounds[2]),
|
|
"max_y": float(geom.bounds[3]),
|
|
}
|
|
|
|
|
|
def area_m2(geom: MultiPolygon) -> float:
|
|
projected = transform(
|
|
Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True).transform,
|
|
geom,
|
|
)
|
|
result = float(projected.area)
|
|
if not isfinite(result) or result <= 0:
|
|
raise ValueError("Area geometry must have a finite positive surface")
|
|
return result
|
|
|
|
|
|
def geometry_bbox_polygon(geom: MultiPolygon):
|
|
return box(*geom.bounds)
|