60 lines
1.8 KiB
Python
60 lines
1.8 KiB
Python
from __future__ import annotations
|
|
|
|
from typing import Any
|
|
|
|
from pyproj import Transformer
|
|
from shapely import force_2d
|
|
from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, box, shape
|
|
from shapely.ops import transform
|
|
from shapely.validation import make_valid
|
|
|
|
|
|
def normalize_to_multipolygon(raw_geometry: dict[str, Any]) -> MultiPolygon:
|
|
geom = force_2d(shape(raw_geometry))
|
|
if geom.is_empty:
|
|
raise ValueError("Geometry is empty")
|
|
|
|
if not geom.is_valid:
|
|
geom = make_valid(geom)
|
|
|
|
if not geom.is_valid:
|
|
raise ValueError("Geometry is invalid and could not be repaired")
|
|
|
|
if geom.geom_type == "Polygon":
|
|
return MultiPolygon([geom])
|
|
if geom.geom_type == "MultiPolygon":
|
|
return MultiPolygon(geom.geoms)
|
|
if isinstance(geom, GeometryCollection):
|
|
polygons = [g for g in geom.geoms if isinstance(g, Polygon)]
|
|
multipolygons = [g for g in geom.geoms if g.geom_type == "MultiPolygon"]
|
|
if not polygons and not multipolygons:
|
|
raise ValueError("Only polygon geometries are supported for AOI")
|
|
normalized = []
|
|
normalized.extend(polygons)
|
|
for mp in multipolygons:
|
|
normalized.extend(mp.geoms)
|
|
return MultiPolygon(normalized)
|
|
|
|
raise ValueError("Only Polygon or MultiPolygon geometries are accepted")
|
|
|
|
|
|
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,
|
|
)
|
|
return float(projected.area)
|
|
|
|
|
|
def geometry_bbox_polygon(geom: MultiPolygon):
|
|
return box(*geom.bounds)
|