"""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,
    anchor: str = "sw",
) -> 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)``,
    optionally rotated by ``orientation_deg``, then projected back to WGS84
    lng/lat. Columns run east (+lng), rows run north (+lat) before rotation, and
    polygons are ordered row-major (r=0 c=0, r=0 c=1, ...) — the order plot refs
    are assigned in.

    ``anchor`` controls what the origin means:

    * ``"sw"`` (default) — the origin is the grid's south-west corner, so the
      whole grid grows north/east from it. Use when the caller picked an explicit
      corner to start from.
    * ``"center"`` — the origin is the CENTRE of the grid, so the grid straddles
      it symmetrically. Use when anchoring on a section's centroid, otherwise the
      grid spills off the north/east edges of the boundary (INDL-55 bug fix).
    """
    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

    # When centring, shift the whole grid (pre-rotation, in metre-space) so its
    # bounding box is centred on the origin instead of starting at it. The grid
    # spans from x=0 to the far edge of the last column ((cols-1)*step_x + width),
    # so half that extent is the offset; same for rows.
    if anchor == "center":
        off_x = ((cols - 1) * step_x + plot_width_m) / 2.0
        off_y = ((rows - 1) * step_y + plot_length_m) / 2.0
    else:
        off_x = off_y = 0.0

    polygons: list[Polygon] = []
    for r in range(rows):
        for c in range(cols):
            x0 = c * step_x - off_x
            y0 = r * step_y - off_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


def infer_lattice_anchor(plot_polygon: Polygon) -> dict:
    """Reverse ``build_plot_grid``: derive the ``origin_lat``/``origin_lng``,
    ``width_m``, ``length_m`` and ``orientation_deg`` that would reproduce
    ``plot_polygon`` as the ``(row=0, col=0)`` cell of an ``anchor="sw"`` grid.

    The returned origin is the polygon's own first ring vertex, which is *a*
    valid anchor for the infinite lattice this cell belongs to (any cell's own
    corner satisfies that), but ``anchor="sw"`` only ever extends outward in
    the +row/+col direction from whatever origin it's given — so anchoring on
    an arbitrary single cell can miss existing cells that sit "before" it in
    row/column order. Use :func:`resolve_lattice_sw_anchor` when re-anchoring
    onto a whole *set* of existing plots.
    """
    coords = list(plot_polygon.exterior.coords)
    (lng0, lat0), (lng1, lat1), (lng2, lat2) = coords[0], coords[1], coords[2]

    lat_clamped = max(min(lat0, 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

    dx1_m = (lng1 - lng0) * m_per_deg_lng
    dy1_m = (lat1 - lat0) * _METERS_PER_DEG_LAT
    dx2_m = (lng2 - lng1) * m_per_deg_lng
    dy2_m = (lat2 - lat1) * _METERS_PER_DEG_LAT

    return {
        "origin_lat": lat0,
        "origin_lng": lng0,
        "width_m": math.hypot(dx1_m, dy1_m),
        "length_m": math.hypot(dx2_m, dy2_m),
        "orientation_deg": math.degrees(math.atan2(dy1_m, dx1_m)) % 360,
    }


def dominant_orientation_deg(polygon: Polygon) -> float:
    """Angle (0-359°, measured the same way as ``build_plot_grid``'s
    ``orientation_deg`` — CCW from +x/east) of a boundary polygon's minimum
    rotated rectangle's longer edge.

    Used to pre-fill the grid generator's orientation so a new plot lattice
    aligns with the shape of the section boundary by default instead of 0°.
    Coordinates are projected into local metre-space first (same cos(lat)
    longitude scaling used elsewhere in this module) so the angle isn't
    skewed by longitude/latitude degrees having different real-world lengths.
    """
    centroid = polygon.centroid
    lat_clamped = max(min(centroid.y, 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

    local_coords = [
        ((lng - centroid.x) * m_per_deg_lng, (lat - centroid.y) * _METERS_PER_DEG_LAT)
        for lng, lat in polygon.exterior.coords
    ]
    rect = Polygon(local_coords).minimum_rotated_rectangle
    (x0, y0), (x1, y1), (x2, y2) = list(rect.exterior.coords)[:3]
    edge1 = math.hypot(x1 - x0, y1 - y0)
    edge2 = math.hypot(x2 - x1, y2 - y1)
    dx, dy = (x1 - x0, y1 - y0) if edge1 >= edge2 else (x2 - x1, y2 - y1)
    return math.degrees(math.atan2(dy, dx)) % 360


def resolve_lattice_sw_anchor(existing_plots: list[Polygon]) -> dict:
    """Given a section's existing plot polygons (assumed to share one lattice —
    same plot size/rotation, produced by one or more prior grid generations),
    return the ``origin_lat``/``origin_lng`` of that lattice's SW-most corner,
    plus ``width_m``/``length_m``/``orientation_deg`` taken from an arbitrary
    member.

    Anchoring a *new* ``build_plot_grid(..., anchor="sw")`` call there
    guarantees the new grid extends outward covering every existing cell
    (each one has a non-negative row/col offset from the SW corner by
    definition), not just the cells reachable from one arbitrarily-picked
    existing plot — picking a single interior cell as the anchor would only
    extend forward from it, silently missing existing cells that sit before
    it in row/column order.
    """
    anchors = [infer_lattice_anchor(p) for p in existing_plots]
    ref = anchors[0]

    theta = math.radians(ref["orientation_deg"])
    cos_t, sin_t = math.cos(theta), math.sin(theta)
    lat_clamped = max(min(ref["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

    # Project every existing cell's own origin into the grid's UNROTATED local
    # frame (relative to ref), then take the min on each axis — the true SW
    # corner of the lattice, even if no single existing plot sits exactly there
    # (e.g. it was clipped away by a section-boundary containment check).
    min_u = min_v = 0.0
    for a in anchors:
        dx_m = (a["origin_lng"] - ref["origin_lng"]) * m_per_deg_lng
        dy_m = (a["origin_lat"] - ref["origin_lat"]) * _METERS_PER_DEG_LAT
        u = dx_m * cos_t + dy_m * sin_t
        v = -dx_m * sin_t + dy_m * cos_t
        min_u = min(min_u, u)
        min_v = min(min_v, v)

    rx, ry = _rotate(min_u, min_v, cos_t, sin_t)
    return {
        "origin_lat": ref["origin_lat"] + ry / _METERS_PER_DEG_LAT,
        "origin_lng": ref["origin_lng"] + rx / m_per_deg_lng,
        "width_m": ref["width_m"],
        "length_m": ref["length_m"],
        "orientation_deg": ref["orientation_deg"],
    }