Files
geointel/backend/app/services/detection_georeferencing.py
T
Codex 0aff8e3b8c
GeoIntel release gates / Compile, test, contracts and builds (push) Canceled after 0s
GeoIntel release gates / Python and npm vulnerability policy (push) Canceled after 0s
GeoIntel release gates / GIS image, SBOM and container scan (push) Canceled after 0s
feat(scope): make Belgium and North Sea operational default
2026-07-22 02:11:48 +02:00

143 lines
6.1 KiB
Python

from __future__ import annotations
from typing import Any
from pyproj import Transformer
from shapely.geometry import Polygon
from app.core.errors import AppError
def pixel_bbox_to_epsg4326_polygon(bbox: list[float], tile: dict[str, Any], crs: str | None = None) -> Polygon:
if len(bbox) != 4:
raise AppError(code="DETECTION_INVALID_BBOX", message="YOLO detection bbox must contain four pixel coordinates", status_code=422)
x_min, y_min, x_max, y_max = [float(value) for value in bbox]
if x_max <= x_min or y_max <= y_min:
raise AppError(code="DETECTION_INVALID_BBOX", message="YOLO detection bbox must have positive width and height", status_code=422)
transform = tile.get("transform")
if isinstance(transform, list) and len(transform) >= 6:
corners = [
_apply_gdal_transform(transform, x_min, y_min),
_apply_gdal_transform(transform, x_max, y_min),
_apply_gdal_transform(transform, x_max, y_max),
_apply_gdal_transform(transform, x_min, y_max),
_apply_gdal_transform(transform, x_min, y_min),
]
else:
corners = _corners_from_bounds(bbox=[x_min, y_min, x_max, y_max], tile=tile)
source_crs = crs or tile.get("crs") or tile.get("source_crs") or "EPSG:4326"
if str(source_crs).upper() not in {"EPSG:4326", "4326"}:
transformer = Transformer.from_crs(source_crs, "EPSG:4326", always_xy=True)
corners = [transformer.transform(x, y) for x, y in corners]
polygon = Polygon(corners)
if polygon.is_empty or not polygon.is_valid:
raise AppError(code="DETECTION_INVALID_GEOMETRY", message="Georeferenced detection geometry is invalid", status_code=422)
return polygon
def pixel_points_to_epsg4326_polygon(points: list[list[float]], tile: dict[str, Any], crs: str | None = None) -> Polygon:
if not isinstance(points, list) or len(points) < 3:
raise AppError(
code="SEGMENTATION_INVALID_MASK",
message="Segmentation mask polygon must contain at least three pixel points",
status_code=422,
)
try:
pixel_points = [(float(point[0]), float(point[1])) for point in points]
except (TypeError, ValueError, IndexError) as exc:
raise AppError(
code="SEGMENTATION_INVALID_MASK",
message="Segmentation mask polygon points must be numeric [x, y] pairs",
status_code=422,
) from exc
transform = tile.get("transform")
if isinstance(transform, list) and len(transform) >= 6:
coordinates = [_apply_gdal_transform(transform, x, y) for x, y in pixel_points]
else:
coordinates = [_project_pixel_with_bounds(tile, x, y) for x, y in pixel_points]
source_crs = crs or tile.get("crs") or tile.get("source_crs") or "EPSG:4326"
if str(source_crs).upper() not in {"EPSG:4326", "4326"}:
transformer = Transformer.from_crs(source_crs, "EPSG:4326", always_xy=True)
coordinates = [transformer.transform(x, y) for x, y in coordinates]
if coordinates[0] != coordinates[-1]:
coordinates.append(coordinates[0])
polygon = Polygon(coordinates)
if not polygon.is_valid:
from shapely.validation import make_valid
repaired = make_valid(polygon)
polygon = _largest_polygon(repaired)
if polygon is None or polygon.is_empty or not polygon.is_valid or polygon.area <= 0:
raise AppError(
code="SEGMENTATION_INVALID_GEOMETRY",
message="Georeferenced segmentation geometry is invalid",
status_code=422,
)
return polygon
def _largest_polygon(geometry: Any) -> Polygon | None:
if isinstance(geometry, Polygon):
return geometry
candidates = [geom for geom in getattr(geometry, "geoms", []) if isinstance(geom, Polygon) and geom.area > 0]
if not candidates:
return None
return max(candidates, key=lambda geom: geom.area)
def _project_pixel_with_bounds(tile: dict[str, Any], px: float, py: float) -> tuple[float, float]:
bounds = tile.get("bounds")
pixel_window = tile.get("pixel_window")
if not (isinstance(bounds, list) and len(bounds) == 4 and isinstance(pixel_window, list) and len(pixel_window) == 4):
raise AppError(
code="DETECTION_TILE_MANIFEST_INVALID",
message="Tile manifest entries require transform or bounds plus pixel_window for georeferencing",
status_code=422,
)
left, bottom, right, top = [float(value) for value in bounds]
_, _, width, height = [float(value) for value in pixel_window]
if width <= 0 or height <= 0:
raise AppError(code="DETECTION_TILE_MANIFEST_INVALID", message="Tile pixel_window must have positive size", status_code=422)
return (left + (px / width) * (right - left), top - (py / height) * (top - bottom))
def _apply_gdal_transform(transform: list[float], x: float, y: float) -> tuple[float, float]:
c, a, b, f, d, e = [float(value) for value in transform[:6]]
return (a * x + b * y + c, d * x + e * y + f)
def _corners_from_bounds(bbox: list[float], tile: dict[str, Any]) -> list[tuple[float, float]]:
bounds = tile.get("bounds")
pixel_window = tile.get("pixel_window")
if not (isinstance(bounds, list) and len(bounds) == 4 and isinstance(pixel_window, list) and len(pixel_window) == 4):
raise AppError(
code="DETECTION_TILE_MANIFEST_INVALID",
message="Tile manifest entries require transform or bounds plus pixel_window for georeferencing",
status_code=422,
)
x_min, y_min, x_max, y_max = bbox
left, bottom, right, top = [float(value) for value in bounds]
_, _, width, height = [float(value) for value in pixel_window]
if width <= 0 or height <= 0:
raise AppError(code="DETECTION_TILE_MANIFEST_INVALID", message="Tile pixel_window must have positive size", status_code=422)
def project(px: float, py: float) -> tuple[float, float]:
x = left + (px / width) * (right - left)
y = top - (py / height) * (top - bottom)
return (x, y)
return [
project(x_min, y_min),
project(x_max, y_min),
project(x_max, y_max),
project(x_min, y_max),
project(x_min, y_min),
]