Initial GeoIntel V1 foundation

backend/app/services/geojson_service.py

@@ -0,0 +1,134 @@
+from __future__ import annotations
+
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any
+
+import json
+from pyproj import Transformer, CRS
+from shapely.geometry import shape
+from shapely.geometry.base import BaseGeometry
+from shapely.ops import unary_union
+from shapely.ops import transform as _transform_geometry
+from shapely.validation import make_valid
+
+
+def parse_geojson_payload(raw_text: str | dict[str, Any]) -> dict[str, Any]:
+    if isinstance(raw_text, dict):
+        payload = raw_text
+    else:
+        try:
+            payload = json.loads(raw_text)
+        except Exception as exc:
+            raise ValueError("Uploaded dataset is not valid JSON") from exc
+
+    if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection":
+        raise ValueError("Upload must be a GeoJSON FeatureCollection")
+
+    features = payload.get("features") or []
+    if not isinstance(features, list):
+        raise ValueError("FeatureCollection features is invalid")
+
+    geometry_types: set[str] = set()
+    geometries = []
+    invalid_features = 0
+    polygon_area_m2: float | None = None
+    crs_assumed = None
+    for feature in features:
+        if not isinstance(feature, dict):
+            continue
+        geometry = feature.get("geometry")
+        if not geometry:
+            continue
+        try:
+            geom = shape(geometry)
+        except Exception as exc:
+            raise ValueError("Invalid feature geometry") from exc
+        if not geom.is_valid:
+            geom = make_valid(geom)
+            if not geom.is_valid:
+                invalid_features += 1
+                raise ValueError("Invalid geometry remains after repair")
+        geometry_types.add(str(geom.geom_type))
+        geometries.append(geom)
+
+    if geometries:
+        unioned = unary_union(geometries)
+        bounds = unioned.bounds
+        bounds_json = {
+            "min_x": float(bounds[0]),
+            "min_y": float(bounds[1]),
+            "max_x": float(bounds[2]),
+            "max_y": float(bounds[3]),
+        }
+    else:
+        bounds_json = None
+
+    crs = None
+    crs_assumed = False
+    raw_crs = payload.get("crs")
+    if isinstance(raw_crs, dict):
+        raw_name = raw_crs.get("properties", {}).get("name")
+        if isinstance(raw_name, str):
+            crs = raw_name
+    elif isinstance(raw_crs, str):
+        crs = raw_crs
+    if not crs:
+        crs = "EPSG:4326"
+        crs_assumed = True
+
+    polygon_area_m2 = _approximate_polygon_area_m2(geometries, crs)
+
+    return {
+        "feature_count": len(features),
+        "geometry_types": sorted(geometry_types),
+        "bounds_json": bounds_json,
+        "approximate_area_m2": polygon_area_m2,
+        "invalid_features": invalid_features,
+        "crs": crs,
+        "crs_assumed": crs_assumed,
+        "extracted_at": datetime.now(timezone.utc).isoformat(),
+        "feature_geometry_count": len(geometries),
+    }
+
+
+def load_dataset_text(file_path: str) -> str:
+    return Path(file_path).read_text(encoding="utf-8")
+
+
+def _approximate_polygon_area_m2(geometries: list[BaseGeometry], crs: str | None) -> float | None:
+    if not geometries:
+        return 0.0
+    try:
+        polygons = [geometry for geometry in geometries if geometry.geom_type.lower() in {"polygon", "multipolygon"}]
+        if not polygons:
+            return None
+        target_crs = CRS.from_epsg(31370)
+        source_crs = _crs_to_epsg(crs)
+        transformer = Transformer.from_crs(source_crs, target_crs, always_xy=True)
+        projected = [_transform_polygon_for_area(geometry, transformer) for geometry in polygons]
+        area = sum(item.area for item in projected)
+        if area < 0:
+            area = 0.0
+        return float(area)
+    except Exception:
+        return None
+
+
+def _crs_to_epsg(value: str | None) -> str:
+    if not value:
+        return "EPSG:4326"
+    normalized = value.upper().strip().replace(" ", "")
+    if normalized.startswith("EPSG:"):
+        return normalized
+    if normalized.replace("-", "").isdigit():
+        return f"EPSG:{normalized}"
+    return "EPSG:4326"
+
+
+def _transform_polygon_for_area(geometry: BaseGeometry, transformer: Transformer):
+    if geometry.is_empty:
+        return geometry
+    if geometry.geom_type.lower() in {"polygon", "multipolygon"}:
+        return _transform_geometry(transformer.transform, geometry)
+    return geometry