"""Geospatial helpers for the Cemetery Map feature (INDL-49 / shared with INDL-41).

Pure-Python geometry validation and conversion built on ``shapely``. PostGIS is
used for the authoritative ``ST_Area`` / ``ST_Centroid`` computation in the
service layer; this module handles GeoJSON <-> geometry conversion, validation,
and grid tiling for bulk plot generation.
"""
from __future__ import annotations

import math
from typing import Any, Optional

from shapely.geometry import Polygon, mapping, shape
from shapely.geometry.base import BaseGeometry

from src.core.exceptions import ValidationError

# One generation may not exceed this many plots (AC-10 / security A04).
MAX_GRID_PLOTS = 2000

# WGS84 metres-per-degree of latitude (near-constant). Longitude scales by cos(lat).
_METERS_PER_DEG_LAT = 111_320.0

_INVALID_BOUNDARY_MSG = "Invalid boundary shape — redraw and try again"


def validate_polygon_geojson(geojson: Any) -> Polygon:
    """Validate an incoming GeoJSON polygon and return a shapely ``Polygon``.

    Raises ``ValidationError`` (HTTP 422) for anything that is not a single,
    topologically valid, closed polygon ring — mirroring the PRD's server-side
    ring-closure / winding-order guard.
    """
    if not isinstance(geojson, dict):
        raise ValidationError(_INVALID_BOUNDARY_MSG)

    geom_type = geojson.get("type")
    if geom_type != "Polygon":
        raise ValidationError(
            f"Boundary must be a GeoJSON Polygon, got {geom_type!r}"
        )

    try:
        geom: BaseGeometry = shape(geojson)
    except Exception:  # malformed coordinates, etc.
        raise ValidationError(_INVALID_BOUNDARY_MSG)

    if geom.is_empty or not isinstance(geom, Polygon):
        raise ValidationError(_INVALID_BOUNDARY_MSG)

    # A valid polygon needs at least 3 distinct vertices.
    if len(geom.exterior.coords) < 4:
        raise ValidationError(_INVALID_BOUNDARY_MSG)

    if not geom.is_valid:
        from shapely.validation import explain_validity

        # A self-intersecting (crossing) ring is adversarial/malformed and must be
        # rejected outright — never silently "repaired" into a different shape
        # (security SAC-03 / A03). buffer(0) is only for benign self-touching rings.
        if "self-intersection" in explain_validity(geom).lower():
            raise ValidationError(_INVALID_BOUNDARY_MSG)

        # Attempt a zero-width buffer repair (fixes self-touching rings); if it
        # still isn't a clean polygon, reject.
        repaired = geom.buffer(0)
        if repaired.is_empty or not isinstance(repaired, Polygon) or not repaired.is_valid:
            raise ValidationError(_INVALID_BOUNDARY_MSG)
        geom = repaired

    # Normalise winding order (exterior counter-clockwise) for consistent storage.
    from shapely.geometry.polygon import orient

    return orient(geom, sign=1.0)


def polygon_to_wkt(geom: Polygon) -> str:
    """Serialise a shapely polygon to WKT for ``ST_GeogFromText``."""
    return geom.wkt


def geojson_to_wkt(geojson: Any) -> str:
    """Validate a GeoJSON polygon and return its WKT string."""
    return polygon_to_wkt(validate_polygon_geojson(geojson))


def geometry_to_geojson(value: Any) -> Optional[dict]:
    """Convert a stored geography value (WKBElement) to a GeoJSON dict, or None."""
    if value is None:
        return None
    try:
        from geoalchemy2.shape import to_shape

        geom = to_shape(value)
        return mapping(geom)
    except Exception:
        return None


def _rotate(x: float, y: float, cos_t: float, sin_t: float) -> tuple[float, float]:
    return (x * cos_t - y * sin_t, x * sin_t + y * cos_t)


def build_plot_grid(
    origin_lat: float,
    origin_lng: float,
    rows: int,
    cols: int,
    plot_width_m: float,
    plot_length_m: float,
    gap_m: float = 0.3,
    orientation_deg: float = 0.0,
) -> list[Polygon]:
    """Tile a rectangle into ``rows × cols`` non-overlapping plot polygons.

    Plots are laid out in metre-space anchored at ``(origin_lat, origin_lng)``
    (the grid's south-west corner), optionally rotated by ``orientation_deg``,
    then projected back to WGS84 lng/lat. Row 0/col 0 is the anchor; columns run
    east (+lng), rows run north (+lat) before rotation.

    Returned polygons are ordered row-major (r=0 c=0, r=0 c=1, ...), which is the
    order plot refs are assigned in.
    """
    theta = math.radians(orientation_deg or 0.0)
    cos_t, sin_t = math.cos(theta), math.sin(theta)

    # Guard against degenerate latitude (poles) for the longitude scale factor.
    lat_clamped = max(min(origin_lat, 89.9), -89.9)
    m_per_deg_lng = _METERS_PER_DEG_LAT * math.cos(math.radians(lat_clamped))
    if m_per_deg_lng < 1e-6:
        m_per_deg_lng = 1e-6

    step_x = plot_width_m + gap_m
    step_y = plot_length_m + gap_m

    polygons: list[Polygon] = []
    for r in range(rows):
        for c in range(cols):
            x0 = c * step_x
            y0 = r * step_y
            # Local rectangle corners (metres), CCW.
            local_corners = [
                (x0, y0),
                (x0 + plot_width_m, y0),
                (x0 + plot_width_m, y0 + plot_length_m),
                (x0, y0 + plot_length_m),
            ]
            ring = []
            for lx, ly in local_corners:
                rx, ry = _rotate(lx, ly, cos_t, sin_t)
                dlng = rx / m_per_deg_lng
                dlat = ry / _METERS_PER_DEG_LAT
                ring.append((origin_lng + dlng, origin_lat + dlat))
            polygons.append(Polygon(ring))
    return polygons