#!/usr/bin/env python3
"""
Site Layers-overlay  -  ELVTR AI for Architects, Assignment 3 student tool.

Adds the layers archiMap does not have (statutory planning + climate/sun) to an
archiMap DXF export, aligned into the same drawing, then renders styled
diagrams and a verification register.

    pip install ezdxf shapely pyproj matplotlib requests
    python site_layers_overlay.py --lat -33.9107 --lon 151.1432 --address "14 Kays Ave West, Dulwich Hill" \
        --jurisdiction NSW --dxf my_archimap.dxf --style style_card.toml --radius 200

Rules this tool keeps:
  * keyless public endpoints only
  * one layer failing never stops the run - it becomes a register row
  * every pulled or computed fact becomes a register row with an EMPTY verdict
  * the style card, not the software, decides how the diagrams look
"""
from __future__ import annotations

VERSION = "0.1.0"

import argparse
import csv
import datetime as dt
import hashlib
import json
import math
import re
import sys
import time
import traceback
import zipfile
from pathlib import Path

try:
    import tomllib  # Python 3.11+
except ModuleNotFoundError:  # pragma: no cover
    tomllib = None

import numpy as np
import requests
from shapely.geometry import (LineString, MultiPolygon, Point, Polygon, box,
                              mapping, shape)
from shapely.ops import unary_union
from shapely import affinity
from pyproj import Transformer
import ezdxf
from ezdxf import bbox as ezbbox
from ezdxf import path as ezpath

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection, PatchCollection
from matplotlib.patches import Polygon as MplPolygon, PathPatch
from matplotlib.path import Path as MplPath

plt.rcParams["svg.fonttype"] = "none"      # keep text as text in the SVG

UA = f"ELVTR-SiteLayers/{VERSION} (architecture coursework; keyless public endpoints)"
TODAY = dt.date.today().isoformat()

# ---------------------------------------------------------------------------
# 1. Jurisdiction adapter table.  NSW filled.  Others are explicit, empty slots:
#    each produces a register row saying so, and the run carries on.
# ---------------------------------------------------------------------------
_PP = "https://mapprod3.environment.nsw.gov.au/arcgis/rest/services/ePlanning/"
_PRIN = _PP + "Planning_Portal_Principal_Planning/MapServer"
_HAZ = _PP + "Planning_Portal_Hazard/MapServer"
_PROT = _PP + "Planning_Portal_Protection/MapServer"

ADAPTERS: dict[str, dict | None] = {
    "NSW": {
        "name": "New South Wales (NSW Planning Portal + NSW Spatial Services)",
        "geocoder": "nsw",
        "cadastre": {
            # Primary: Spatial Services portal (fast; answers ArcGIS JSON only - f=geojson fails there).
            # Fallback: the older maps.six service (GeoJSON; was timing out at 40 s on 23 Sep 2026).
            "sources": [
                ("https://portal.spatial.nsw.gov.au/server/rest/services/NSW_Land_Parcel_Property_Theme/FeatureServer/8", "json"),
                ("https://maps.six.nsw.gov.au/arcgis/rest/services/public/NSW_Cadastre/MapServer/9", "geojson"),
            ],
            "url": "https://portal.spatial.nsw.gov.au/server/rest/services/NSW_Land_Parcel_Property_Theme/FeatureServer/8",
            "id_field": "lotidstring",
            "label": "Lot {lotnumber} {planlabel}",
        },
        # key, dxf layer, title, service url, layer id, label, on-site assertion
        "layers": [
            dict(key="ZONE", dxf="AI_ZONE", title="Land zoning", url=_PRIN, id=19,
                 label="{SYM_CODE}",
                 fact="zoned {SYM_CODE} {LAY_CLASS}",
                 value="{SYM_CODE}"),
            dict(key="HOB", dxf="AI_HOB", title="Height of buildings", url=_PRIN, id=14,
                 label="{MAX_B_H_M} m",
                 fact="mapped maximum building height {MAX_B_H} {UNITS} (map code {SYM_CODE}, {LEGIS_REF_CLAUSE})",
                 value="{MAX_B_H_M} m"),
            dict(key="FSR", dxf="AI_FSR", title="Floor space ratio", url=_PRIN, id=11,
                 label="{FSR}:1",
                 fact="mapped maximum floor space ratio {FSR}:1 (map code {SYM_CODE}, {LEGIS_REF_CLAUSE})",
                 value="{FSR}:1"),
            dict(key="ADDCTRL", dxf="AI_ADDCTRL", title="Additional controls (FSR / height areas)", url=_PRIN, id=[10, 13],
                 label="{LABEL}",
                 fact="inside additional controls area {LABEL} ({LEGIS_REF_CLAUSE}) - this clause can vary the mapped FSR or height; read it",
                 value="{LABEL}"),
            dict(key="HERITAGE", dxf="AI_HERITAGE", title="Heritage (LEP items and conservation areas)", url=_PRIN, id=16,
                 label="{H_ID}",
                 fact="{LAY_CLASS}: {H_NAME} ({H_ID}, {SIG} significance, {LEGIS_REF_CLAUSE})",
                 value="{H_ID}"),
            dict(key="SHR", dxf="AI_HERITAGE", title="State Heritage Register curtilage", url=_PRIN, id=221,
                 label="SHR {LISTINGNO}",
                 fact="State Heritage Register curtilage: {ITEMNAME} (listing {LISTINGNO})",
                 value="SHR {LISTINGNO}"),
            dict(key="FLOOD", dxf="AI_FLOOD", title="Flood planning", url=_HAZ, id=230,
                 label="FPA",
                 fact="mapped flood planning land ({LAY_CLASS})",
                 value="{LAY_CLASS}"),
            dict(key="BUSHFIRE", dxf="AI_BUSHFIRE", title="Bushfire prone land", url=_HAZ, id=229,
                 label="{d_Category}",
                 fact="mapped bushfire prone land ({d_Category})",
                 value="{d_Category}"),
            dict(key="ASS", dxf="AI_ASS", title="Acid sulfate soils", url=_PROT, id=234,
                 label="ASS {LAY_CLASS}",
                 fact="acid sulfate soils {LAY_CLASS}",
                 value="{LAY_CLASS}"),
        ],
    },
    # Explicit, empty.  Fill one of these in and the tool picks it up.
    "VIC": None,
    "QLD": None,
    "NZ": None,
    "SG": None,
    "HK": None,
}

# Layers the tool always writes (universal - every jurisdiction).
UNIVERSAL_DXF = ["AI_SITE", "AI_CADASTRE", "AI_SUNPATH", "AI_WIND", "AI_NORTH", "AI_RADIUS"]


# ---------------------------------------------------------------------------
# 2. Small helpers
# ---------------------------------------------------------------------------
class SafeDict(dict):
    def __missing__(self, k):
        return "?"


def fmt(template: str, props: dict) -> str:
    clean = {k: ("" if v is None else v) for k, v in (props or {}).items()}
    out = template.format_map(SafeDict(clean))
    out = re.sub(r"\(\s*,", "(", out).replace(", )", ")").replace("()", "")
    return re.sub(r"\s{2,}", " ", out).strip()


def slugify(s: str) -> str:
    s = re.sub(r"[^A-Za-z0-9]+", "_", s.strip().lower()).strip("_")
    return s[:48] or "site"


def log(msg: str):
    print(msg, flush=True)


# ---------------------------------------------------------------------------
# 3. Network with cache.  One retry, then give up on that request.
#    A host that fails twice is skipped for the rest of the run (no retry storms).
# ---------------------------------------------------------------------------
class Net:
    def __init__(self, cache_dir: Path, offline: bool = False, timeout: int = 30):
        self.cache = Path(cache_dir)
        self.cache.mkdir(parents=True, exist_ok=True)
        self.offline = offline
        self.timeout = timeout
        self.s = requests.Session()
        self.s.headers["User-Agent"] = UA
        self.bad_hosts: dict[str, int] = {}

    def _key(self, url: str) -> Path:
        return self.cache / (hashlib.sha1(url.encode()).hexdigest()[:20] + ".json")

    def get_json(self, url: str, params: dict | None = None):
        """Returns (data, full_url, retrieved_iso, from_cache)."""
        full = requests.Request("GET", url, params=params).prepare().url
        f = self._key(full)
        if f.exists():
            rec = json.loads(f.read_text(encoding="utf-8"))
            return rec["data"], full, rec["retrieved"], True
        if self.offline:
            raise RuntimeError("offline mode and this request is not in the cache")
        host = re.sub(r"^https?://([^/]+).*$", r"\1", full)
        if self.bad_hosts.get(host, 0) >= 2:
            raise RuntimeError(f"{host} failed twice earlier in this run - skipped")
        last = None
        for attempt in range(2):
            try:
                if sys.platform == "emscripten":      # running in a browser (Pyodide, inside a worker)
                    from js import XMLHttpRequest
                    xhr = XMLHttpRequest.new()
                    xhr.open("GET", full, False)
                    try:
                        xhr.timeout = self.timeout * 1000   # allowed for synchronous requests in a worker
                    except Exception:
                        pass
                    xhr.send(None)
                    if xhr.status == 0 or xhr.status >= 400:
                        raise RuntimeError(f"HTTP {xhr.status or 'no answer'}")
                    data = json.loads(xhr.responseText)
                    if isinstance(data, dict) and data.get("error"):
                        raise RuntimeError(f"service error: {str(data['error'])[:200]}")
                    stamp = dt.datetime.now().isoformat(timespec="seconds")
                    f.write_text(json.dumps({"url": full, "retrieved": stamp, "data": data}), encoding="utf-8")
                    return data, full, stamp, False
                r = self.s.get(full, timeout=self.timeout)
                if r.status_code >= 500:
                    raise RuntimeError(f"HTTP {r.status_code}")
                r.raise_for_status()
                data = r.json()
                if isinstance(data, dict) and data.get("error"):
                    raise RuntimeError(f"service error: {str(data['error'])[:200]}")
                stamp = dt.datetime.now().isoformat(timespec="seconds")
                f.write_text(json.dumps({"url": full, "retrieved": stamp, "data": data}), encoding="utf-8")
                return data, full, stamp, False
            except Exception as e:  # noqa
                last = e
                if attempt == 0:
                    time.sleep(3)
        self.bad_hosts[host] = self.bad_hosts.get(host, 0) + 1
        raise RuntimeError(f"{type(last).__name__}: {str(last)[:200]}")


def _esri_to_geojson(f):
    """ArcGIS JSON polygon (rings, outer clockwise) -> GeoJSON feature."""
    rings = (f.get("geometry") or {}).get("rings") or []
    outers, holes = [], []
    for r in rings:
        a = sum(x0 * y1 - x1 * y0 for (x0, y0), (x1, y1) in zip(r[:-1], r[1:]))
        (outers if a < 0 else holes).append(r)
    polys = [[o] for o in outers]
    for h in holes:
        hp = Polygon(h)
        for p in polys:
            if Polygon(p[0]).contains(hp.representative_point()):
                p.append(h)
                break
    geom = ({"type": "Polygon", "coordinates": polys[0]} if len(polys) == 1 else
            {"type": "MultiPolygon", "coordinates": polys})
    return {"type": "Feature", "properties": f.get("attributes") or {}, "geometry": geom}


def arcgis_query(net: Net, layer_url: str, env_lonlat, out_fields="*", depth=0, max_depth=3, fmt_="geojson"):
    """Envelope query returning GeoJSON features (GDA94 lon/lat, EPSG:4283).
    No resultOffset paging: NSW Spatial Services rejects it (HTTP 400, checked 23 Sep 2026).
    If the server caps the answer, the envelope is split into quarters and asked again."""
    params = {
        "where": "1=1",
        "geometry": ",".join(f"{v:.7f}" for v in env_lonlat),
        "geometryType": "esriGeometryEnvelope",
        "inSR": "4283", "outSR": "4283",
        "spatialRel": "esriSpatialRelIntersects",
        "outFields": out_fields, "returnGeometry": "true",
        "f": fmt_,
    }
    data, full, stamp, _ = net.get_json(layer_url.rstrip("/") + "/query", params)
    feats = data.get("features", []) or []
    if fmt_ == "json":
        feats = [_esri_to_geojson(f) for f in feats if (f.get("geometry") or {}).get("rings")]
    more = data.get("exceededTransferLimit") or (data.get("properties") or {}).get("exceededTransferLimit")
    if not more:
        return feats, [full], stamp
    if depth >= max_depth:
        log(f"    note: {layer_url.split('/')[-3]}/{layer_url.split('/')[-1]} capped at {len(feats)} features - study radius too big")
        return feats, [full], stamp
    x0, y0, x1, y1 = env_lonlat
    xm, ym = (x0 + x1) / 2, (y0 + y1) / 2
    seen, out, urls = set(), [], []
    for q in [(x0, y0, xm, ym), (xm, y0, x1, ym), (x0, ym, xm, y1), (xm, ym, x1, y1)]:
        f_, u_, stamp = arcgis_query(net, layer_url, q, out_fields, depth + 1, max_depth, fmt_)
        urls += u_
        for f in f_:
            p = f.get("properties") or {}
            key = p.get("OBJECTID") or p.get("objectid") or p.get("lotidstring") or json.dumps(f.get("geometry"))[:200]
            if key not in seen:
                seen.add(key)
                out.append(f)
    return out, [full] + urls, stamp


# ---------------------------------------------------------------------------
# 4. Local metric frame centred on the pin (transverse Mercator, GRS80, metres).
# ---------------------------------------------------------------------------
class Frame:
    def __init__(self, lat: float, lon: float):
        self.lat, self.lon = lat, lon
        self.proj = (f"+proj=tmerc +lat_0={lat} +lon_0={lon} +k=1 +x_0=0 +y_0=0 "
                     f"+ellps=GRS80 +units=m +no_defs")
        self.fwd = Transformer.from_crs("EPSG:4283", self.proj, always_xy=True)
        self.inv = Transformer.from_crs(self.proj, "EPSG:4283", always_xy=True)

    def xy(self, lon, lat):
        return self.fwd.transform(lon, lat)

    def geom(self, g):
        from shapely.ops import transform
        return transform(lambda x, y, z=None: self.fwd.transform(x, y), g)

    def envelope_lonlat(self, r):
        pts = [self.inv.transform(x, y) for x, y in [(-r, -r), (r, -r), (r, r), (-r, r)]]
        xs, ys = zip(*pts)
        return (min(xs), min(ys), max(xs), max(ys))


# ---------------------------------------------------------------------------
# 5. Register
# ---------------------------------------------------------------------------
REG_COLS = ["layer", "assertion", "source URL", "retrieved date", "confidence", "verdict", "note"]


class Register:
    def __init__(self):
        self.rows: list[dict] = []

    def add(self, layer, assertion, source="", retrieved="", confidence="MEDIUM"):
        self.rows.append({"layer": layer, "assertion": assertion, "source URL": source,
                          "retrieved date": retrieved, "confidence": confidence,
                          "verdict": "", "note": ""})

    def write(self, path: Path):
        with open(path, "w", newline="", encoding="utf-8-sig") as f:
            w = csv.DictWriter(f, fieldnames=REG_COLS)
            w.writeheader()
            w.writerows(self.rows)


# ---------------------------------------------------------------------------
# 6. Geocoding (only if the student has no pin).  The archiMap pin is preferred.
# ---------------------------------------------------------------------------
_ROAD_TYPES = r"(street|st|avenue|ave|av|road|rd|lane|ln|place|pl|parade|pde|drive|dr|crescent|cres|cr|court|ct|terrace|tce|highway|hwy|boulevard|blvd|way|close|cl|circuit|cct|square|sq|grove|gr|esplanade|esp)"


def geocode(net: Net, address: str, jurisdiction: str):
    """Returns (lat, lon, source_url, method, retrieved)."""
    if jurisdiction == "NSW":
        m = re.match(r"^\s*(\d+[A-Za-z]?)(?:-\d+)?\s+(.+?),\s*(.+?)(?:\s+NSW)?(?:\s+\d{4})?\s*$", address, re.I)
        if m:
            num, road, suburb = m.groups()
            road_name = re.sub(rf"\s+{_ROAD_TYPES}\b.*$", "", road, flags=re.I).strip()
            try:
                data, url, stamp, _ = net.get_json(
                    "https://maps.six.nsw.gov.au/services/public/Address_Location",
                    {"projection": "EPSG:4283", "houseNumber": num, "roadName": road_name, "suburb": suburb})
                adds = (data.get("addressResult") or {}).get("addresses") or []
                if adds:
                    a = adds[0]
                    lon, lat = a.get("addressPoint", {}).get("centreX"), a.get("addressPoint", {}).get("centreY")
                    if lon is None:
                        lon, lat = a.get("longitude"), a.get("latitude")
                    if lon is not None:
                        return float(lat), float(lon), url, "NSW Spatial Services address point", stamp
            except Exception as e:  # fall through to OSM
                log(f"  NSW geocoder failed ({e}); trying OpenStreetMap")
    data, url, stamp, _ = net.get_json("https://nominatim.openstreetmap.org/search",
                                       {"q": address, "format": "jsonv2", "limit": 1})
    if not data:
        raise RuntimeError("address not found - use the pin coordinates from archiMap instead")
    return float(data[0]["lat"]), float(data[0]["lon"]), url, "OpenStreetMap Nominatim (street-level at best)", stamp


# ---------------------------------------------------------------------------
# 7. archiMap DXF ingest + alignment
# ---------------------------------------------------------------------------
BASE_ROLES = [  # regex on archiMap layer names -> role used by the renderer
    ("buildings", r"build|footprint|bldg"),
    ("roads", r"road|street|highway|rail|path"),
    ("green", r"green|park|veg|tree|grass"),
    ("water", r"water|river|sea|coast"),
    ("terrain", r"terrain|contour|topo|elev"),
    ("cadastre", r"cadast|parcel|lot|property|boundar"),
    ("noise", r"noise"),
    ("wind", r"wind"),
]
INSUNITS = {0: "unitless", 1: "inches", 2: "feet", 4: "millimetres", 5: "centimetres", 6: "metres", 14: "decimetres"}
UNIT_TO_M = {4: 0.001, 5: 0.01, 6: 1.0, 14: 0.1, 1: 0.0254, 2: 0.3048}


def role_of(layer: str) -> str:
    for role, rx in BASE_ROLES:
        if re.search(rx, layer, re.I):
            return role
    return "other"


def dxf_polylines(doc, flatten=0.25, layers=None):
    """All linework in modelspace as {layer: [Nx2 arrays]}.  Blocks exploded."""
    out: dict[str, list] = {}

    def take(e, layer):
        try:
            p = ezpath.make_path(e)
        except Exception:
            return
        try:
            for sub in p.sub_paths():
                v = np.array([(q.x, q.y) for q in sub.flattening(flatten)])
                if len(v) >= 2:
                    out.setdefault(layer, []).append(v)
        except Exception:
            pass

    for e in doc.modelspace():
        lay = e.dxf.get("layer", "0")
        if layers and lay not in layers:
            continue
        if e.dxftype() == "INSERT":
            try:
                for ve in e.virtual_entities():
                    take(ve, lay if ve.dxf.get("layer", "0") == "0" else ve.dxf.layer)
            except Exception:
                pass
        elif e.dxftype() in ("TEXT", "MTEXT", "DIMENSION", "POINT"):
            continue
        else:
            take(e, lay)
    return out


def inspect_dxf(path: Path) -> dict:
    try:
        doc = ezdxf.readfile(str(path))
        recovered = False
    except Exception:
        from ezdxf import recover
        doc, auditor = recover.readfile(str(path))
        recovered = True
    msp = doc.modelspace()
    counts: dict[str, dict] = {}
    for e in msp:
        lay = e.dxf.get("layer", "0")
        counts.setdefault(lay, {})
        counts[lay][e.dxftype()] = counts[lay].get(e.dxftype(), 0) + 1
    ext = ezbbox.extents(msp, fast=True)
    zs = []
    for e in list(msp)[:5000]:
        try:
            if e.dxftype() == "LWPOLYLINE":
                zs.append(e.dxf.elevation)
            elif e.dxftype() == "POLYLINE":
                zs += [v.dxf.location.z for v in e.vertices][:5]
            elif e.dxftype() == "LINE":
                zs += [e.dxf.start.z]
        except Exception:
            pass
    ins = doc.header.get("$INSUNITS", 0)
    info = {
        "file": str(path), "dxf_version": doc.dxfversion, "recovered": recovered,
        "insunits": ins, "units": INSUNITS.get(ins, str(ins)),
        "extmin": list(ext.extmin)[:3] if ext.has_data else None,
        "extmax": list(ext.extmax)[:3] if ext.has_data else None,
        "size": [ext.size.x, ext.size.y] if ext.has_data else None,
        "z_range": [float(min(zs)), float(max(zs))] if zs else None,
        "layers": {k: {"role": role_of(k), "entities": v} for k, v in sorted(counts.items())},
    }
    info["frame_guess"] = guess_frame(info)
    return info, doc


def guess_frame(info) -> str:
    if not info["extmin"]:
        return "empty"
    cx = (info["extmin"][0] + info["extmax"][0]) / 2
    cy = (info["extmin"][1] + info["extmax"][1]) / 2
    if abs(cx) <= 180 and abs(cy) <= 90 and max(info["size"]) < 1:
        return "lonlat"
    if 1.0e5 < cx < 9.0e5 and 1.0e6 < abs(cy) < 1.0e7:
        return "utm"          # MGA / UTM eastings-northings
    if abs(cx) > 1.0e6 and abs(cy) > 1.0e5:
        return "webmercator"
    return "local"


def initial_transform(frame_kind, info, fr: Frame, lat, lon):
    """Function mapping local metres (pin at 0,0) -> DXF coordinates, before fitting."""
    k = UNIT_TO_M.get(info["insunits"], 1.0)
    if frame_kind == "utm":
        zone = int((lon + 180) // 6) + 1
        south = lat < 0
        epsg = (28300 if south else 32600) + zone   # MGA94 in AU; the fit absorbs datum offsets
        if not (-44 < lat < -9 and 112 < lon < 154):
            epsg = (32700 if south else 32600) + zone
        t = Transformer.from_crs(fr.proj, f"EPSG:{epsg}", always_xy=True)
        return lambda x, y: t.transform(x, y), f"UTM/MGA zone {zone} (EPSG:{epsg})"
    if frame_kind == "webmercator":
        t = Transformer.from_crs(fr.proj, "EPSG:3857", always_xy=True)
        return lambda x, y: t.transform(x, y), "Web Mercator (EPSG:3857)"
    if frame_kind == "lonlat":
        return lambda x, y: fr.inv.transform(x, y), "longitude/latitude degrees"
    cx = (info["extmin"][0] + info["extmax"][0]) / 2
    cy = (info["extmin"][1] + info["extmax"][1]) / 2
    return (lambda x, y: (np.asarray(x) / k + cx, np.asarray(y) / k + cy),
            f"local drawing frame, {INSUNITS.get(info['insunits'], 'unknown units')}, pin assumed at drawing centre")


def _densify(lines, step):
    pts = []
    for v in lines:
        for a, b in zip(v[:-1], v[1:]):
            d = float(np.hypot(*(b - a)))
            n = max(1, int(d / step))
            t = np.linspace(0, 1, n, endpoint=False)[:, None]
            pts.append(a + (b - a) * t)
        pts.append(v[-1:])
    return np.vstack(pts) if pts else np.zeros((0, 2))


_TREES: dict = {}


def _nearest(a, b, chunk=2000):
    """For each point in a, distance and index of nearest in b (KD-tree if scipy, else brute force)."""
    try:
        from scipy.spatial import cKDTree
        key = (id(b), b.shape)
        if key not in _TREES:
            _TREES.clear()
            _TREES[key] = cKDTree(b)
        d, idx = _TREES[key].query(a)
        return d, idx
    except ImportError:
        pass
    d = np.empty(len(a))
    idx = np.empty(len(a), dtype=int)
    for i in range(0, len(a), chunk):
        aa = a[i:i + chunk]
        dd = ((aa[:, None, :] - b[None, :, :]) ** 2).sum(-1)
        j = dd.argmin(1)
        idx[i:i + chunk] = j
        d[i:i + chunk] = np.sqrt(dd[np.arange(len(aa)), j])
    return d, idx


def _grid_subsample(p, cell):
    if len(p) == 0:
        return p
    keys = np.floor(p / cell).astype(np.int64)
    _, keep = np.unique(keys, axis=0, return_index=True)
    return p[np.sort(keep)]


def fit_rigid(src, dst, allow_rotation=True, iters=30, start_gate=25.0, end_gate=2.0, R=None, t=None):
    """Trimmed ICP: find R,t so src @ R.T + t ~ dst.  Returns (R, t, residual stats)."""
    R = np.eye(2) if R is None else R
    t = np.zeros(2) if t is None else t
    gate = start_gate
    for i in range(iters):
        cur = src @ R.T + t
        d, j = _nearest(cur, dst)
        m = d < gate
        if m.sum() < 12:
            gate *= 1.5
            continue
        A, B = cur[m], dst[j[m]]
        ca, cb = A.mean(0), B.mean(0)
        Ri = np.eye(2)
        if allow_rotation:
            H = (A - ca).T @ (B - cb)
            U, _, Vt = np.linalg.svd(H)
            Ri = Vt.T @ U.T
            if np.linalg.det(Ri) < 0:
                Vt[-1] *= -1
                Ri = Vt.T @ U.T
        ti = cb - Ri @ ca
        R, t = Ri @ R, Ri @ t + ti
        gate = max(end_gate, gate * 0.75)
    return R, t, _fit_stats(src @ R.T + t, dst, R, t)


def _fit_stats(cur, dst, R, t, inlier=5.0):
    d, _ = _nearest(cur, dst)
    inl = d[d < inlier]
    return {
        "median_m": float(np.median(inl)) if len(inl) else float("nan"),
        "median_all_m": float(np.median(d)) if len(d) else float("nan"),
        "p90_inliers_m": float(np.percentile(inl, 90)) if len(inl) else float("nan"),
        "inlier_share": float(len(inl) / max(1, len(d))),
        "rotation_deg": float(math.degrees(math.atan2(R[1, 0], R[0, 0]))),
        "shift_m": [float(t[0]), float(t[1])],
        "n_points": int(len(d)),
    }


def _line_residual(pts, lines, lo, hi, n=4000, inlier=None):
    """Point-to-linework distance (not point-to-point), so sampling density does not inflate it."""
    import shapely
    from shapely.geometry import MultiLineString
    keep = [v for v in lines if np.any((v[:, 0] > lo[0]) & (v[:, 0] < hi[0]) & (v[:, 1] > lo[1]) & (v[:, 1] < hi[1]))]
    if not keep:
        return {}
    ml = MultiLineString([v for v in keep if len(v) > 1])
    if len(pts) > n:
        pts = pts[np.random.default_rng(0).choice(len(pts), n, replace=False)]
    d = shapely.distance(shapely.points(pts), ml)
    return {"_d_line": d}


def fit_progressive(src, dst, centre, unit, R0=np.eye(2)):
    """Coarse-to-fine by distance from the pin: fit the neighbourhood first, where a small
    rotation error cannot yet throw points onto the wrong lot, then widen."""
    R = R0.copy()
    t = centre - R0 @ centre
    r_all = np.hypot(*(src - centre).T)
    for rad, g0, g1 in [(60, 20, 2), (120, 6, 1.5), (250, 4, 1), (1e9, 3, 0.75)]:
        sel = src[r_all < rad / unit]
        if len(sel) < 30:
            continue
        R, t, _ = fit_rigid(sel, dst, start_gate=g0 / unit, end_gate=g1 / unit, R=R, t=t, iters=20)
    return R, t, _fit_stats(src @ R.T + t, dst, R, t, inlier=5.0 / unit)


class Alignment:
    """Maps local metres (pin at 0,0) into the archiMap DXF's own coordinates."""

    def __init__(self, base_fn, R=np.eye(2), t=np.zeros(2), unit_m=1.0, label="", stats=None, ok=True):
        self.base_fn, self.R, self.t, self.unit_m = base_fn, R, t, unit_m
        self.label, self.stats, self.ok = label, stats or {}, ok

    def pt(self, x, y):
        X, Y = self.base_fn(np.asarray(x, float), np.asarray(y, float))
        P = np.stack([np.asarray(X, float).ravel(), np.asarray(Y, float).ravel()], 1)
        return P @ self.R.T + self.t

    def geom(self, g):
        from shapely.ops import transform

        def f(x, y, z=None):
            P = self.pt(x, y)
            return P[:, 0], P[:, 1]
        return transform(f, g)

    @property
    def rotation_deg(self):
        return math.degrees(math.atan2(self.R[1, 0], self.R[0, 0]))

    def scale_to_m(self):
        """Drawing units per metre, measured."""
        a = self.pt([0.0, 100.0], [0.0, 0.0])
        return float(np.hypot(*(a[1] - a[0]))) / 100.0


def identity_alignment():
    return Alignment(lambda x, y: (x, y), label="tool's own local frame: metres, pin at 0,0, north = +Y")


def align_to_dxf(doc, info, fr: Frame, lat, lon, cad_local: list, reg: Register, cad_layer=None):
    kind = info["frame_guess"]
    base_fn, label = initial_transform(kind, info, fr, lat, lon)
    lines = dxf_polylines(doc, flatten=0.2)
    cand = [cad_layer] if cad_layer else [k for k in lines if role_of(k) == "cadastre"]
    if not cand:
        cand = [k for k in lines if role_of(k) == "buildings"]
        note = "no cadastre layer found in the DXF - fitted against building footprints instead (weaker)"
    else:
        note = f"fitted against DXF layer(s): {', '.join(cand)}"
    dst = [v for k in cand for v in lines.get(k, [])]
    if not dst or not cad_local:
        reg.add("ALIGNMENT", "Could not align: no matching linework in the archiMap DXF or no cadastre pulled. "
                "Pulled layers written to a separate DXF in the tool's local frame instead.", "", TODAY, "NOT CHECKED")
        return None
    # pulled cadastre -> initial DXF coordinates
    src_lines = []
    for g in cad_local:
        for ring in _rings(g):
            c = np.asarray(ring.coords)
            X, Y = base_fn(c[:, 0], c[:, 1])
            src_lines.append(np.stack([np.ravel(X), np.ravel(Y)], 1))
    unit = UNIT_TO_M.get(info["insunits"], 1.0)
    step = 1.0 / unit
    src = _grid_subsample(_densify(src_lines, step), 1.5 / unit)
    dstp = _grid_subsample(_densify(dst, 0.3 * step), 0.25 / unit)
    # keep only the overlap zone to stop the fit being pulled by edges
    lo, hi = dstp.min(0), dstp.max(0)
    src = src[(src[:, 0] > lo[0]) & (src[:, 0] < hi[0]) & (src[:, 1] > lo[1]) & (src[:, 1] < hi[1])]
    if len(src) < 50:
        reg.add("ALIGNMENT", f"Could not align: pulled cadastre does not overlap the DXF ({label}). "
                "Check the pin coordinates.", "", TODAY, "NOT CHECKED")
        return None
    m = 40.0 / unit
    lo2, hi2 = src.min(0) - m, src.max(0) + m
    dstp = dstp[(dstp[:, 0] > lo2[0]) & (dstp[:, 0] < hi2[0]) & (dstp[:, 1] > lo2[1]) & (dstp[:, 1] < hi2[1])]
    # rotation search: straight first; only widen if the straight fit is poor
    centre = np.ravel(np.c_[base_fn(np.array([0.0]), np.array([0.0]))])
    best = None
    for rot0 in [0, 5, -5, 10, -10, 20, -20, 30, -30, 45, -45, 60, -60, 90, -90, 135, -135, 180]:
        th = math.radians(rot0)
        R0 = np.array([[math.cos(th), -math.sin(th)], [math.sin(th), math.cos(th)]])
        Rt, tt, st = fit_progressive(src, dstp, centre, unit, R0)
        score = (st["inlier_share"], -st["median_m"])
        if best is None or score > best[4]:
            best = (Rt, tt, st, rot0, score)
        if st["median_m"] * unit < 0.5 and st["inlier_share"] > 0.6:
            break
    Rt, tt, st, rot0, _ = best
    st.update(_line_residual(src @ Rt.T + tt, dst, lo2, hi2))
    if "_d_line" in st:
        d = st.pop("_d_line") * unit
        inl = d[d < 5.0]
        st["median_m"] = float(np.median(inl)) if len(inl) else float("nan")
        st["p90_inliers_m"] = float(np.percentile(inl, 90)) if len(inl) else float("nan")
        st["inlier_share"] = float(len(inl) / max(1, len(d)))
        st["median_all_m"] = float(np.median(d))
    else:
        st = {k: (v * unit if k in ("median_m", "median_all_m", "p90_inliers_m") else v) for k, v in st.items()}
    al = Alignment(base_fn, Rt, tt, unit, label, st)
    res = st["median_m"]
    al.ok = bool(res <= 1.0 and st["inlier_share"] > 0.5)
    verdict = "within" if al.ok else "WORSE THAN"
    reg.add("ALIGNMENT",
            f"Pulled layers fitted to the archiMap drawing ({label}; {note}). Residual on matched boundaries: median "
            f"{res:.2f} m, 90th percentile {st['p90_inliers_m']:.2f} m; {st['inlier_share']*100:.0f}% of pulled cadastre "
            f"points found a match within 5 m; rotation applied {al.rotation_deg:.2f} deg. "
            f"Fit is {verdict} the 1 m tolerance." + ("" if al.ok else
            " Pulled layers written to a separate DXF, not merged - do not trust their position against the archiMap base."),
            "", TODAY, "LOW")
    return al


def _rings(g):
    if g.is_empty:
        return []
    if isinstance(g, Polygon):
        return [g.exterior] + list(g.interiors)
    if isinstance(g, MultiPolygon):
        return [r for p in g.geoms for r in _rings(p)]
    if hasattr(g, "geoms"):
        return [r for p in g.geoms for r in _rings(p)]
    return []


def _polys(g):
    if g.is_empty:
        return []
    if isinstance(g, Polygon):
        return [g]
    if hasattr(g, "geoms"):
        return [p for q in g.geoms for p in _polys(q)]
    return []


# ---------------------------------------------------------------------------
# 8. Pull statutory layers
# ---------------------------------------------------------------------------
def pull_cadastre(net, adapter, fr: Frame, radius, reg):
    cfg = adapter["cadastre"]
    env = fr.envelope_lonlat(radius * 1.15)
    last = None
    for url, fmt_ in cfg.get("sources", [(cfg["url"], "geojson")]):
        try:
            feats, urls, stamp = arcgis_query(net, url, env, fmt_=fmt_,
                                              out_fields="lotidstring,lotnumber,planlabel,sectionnumber,planlotarea,planlotareaunits")
            break
        except Exception as e:
            last = e
            log(f"  cadastre source failed ({url.split('/services/')[1][:40]}...): {e}; trying the next one")
    else:
        raise RuntimeError(f"every cadastre source failed: {last}")
    out = []
    for f in feats:
        try:
            g = fr.geom(shape(f["geometry"]))
            if not g.is_valid:
                g = g.buffer(0)
            out.append((g, f.get("properties") or {}))
        except Exception:
            continue
    return out, urls[0] if urls else cfg["url"], stamp


def pick_site(cad, pin_xy=(0.0, 0.0), lot_ids=None):
    pin = Point(pin_xy)
    if lot_ids:
        sel = [(g, p) for g, p in cad if str(p.get("lotidstring")) in lot_ids]
        if sel:
            return unary_union([g for g, _ in sel]), [p for _, p in sel], "lots named by student"
    hit = [(g, p) for g, p in cad if g.contains(pin)]
    if hit:
        hit.sort(key=lambda gp: gp[0].area)   # smallest containing lot (strata / stratum overlaps)
        return hit[0][0], [hit[0][1]], "lot containing the pin"
    near = sorted(cad, key=lambda gp: gp[0].distance(pin))
    if near and near[0][0].distance(pin) < 20:
        return near[0][0], [near[0][1]], f"nearest lot, {near[0][0].distance(pin):.1f} m from pin (pin not inside any lot)"
    return pin.buffer(10), [], "no lot found near the pin - a 10 m circle stands in for the site"


def pull_statutory(net, adapter, fr, radius, site, reg, out_dir: Path):
    circle = Point(0, 0).buffer(radius, 128)
    env = fr.envelope_lonlat(radius)
    layers = {}
    gj_dir = out_dir
    for L in adapter["layers"]:
        ids = L["id"] if isinstance(L["id"], list) else [L["id"]]
        feats, src = [], ""
        stamp = ""
        try:
            for i in ids:
                f_, urls, stamp = arcgis_query(net, f"{L['url']}/{i}", env)
                dom, dates = layer_schema(net, f"{L['url']}/{i}")
                for f in f_:
                    decode_props(f.setdefault("properties", {}) or {}, dom, dates)
                feats += f_
                src = src or (urls[0] if urls else "")
        except Exception as e:
            reg.add(L["key"], f"NOT CHECKED - {L['title']}: the service did not answer ({e}). "
                    "Absence here is not an all-clear; check the NSW Planning Portal Spatial Viewer by hand.",
                    f"{L['url']}/{ids[0]}", TODAY, "NOT CHECKED")
            log(f"  {L['key']:9s} FAILED  {e}")
            continue
        items = []
        for f in feats:
            try:
                g = fr.geom(shape(f["geometry"]))
                if not g.is_valid:
                    g = g.buffer(0)
                gc = g.intersection(circle)
                if gc.is_empty:
                    continue
                items.append({"geom": gc, "props": f.get("properties") or {}})
            except Exception:
                continue
        layers[L["key"]] = {"cfg": L, "items": items, "source": src, "retrieved": stamp}
        # GeoJSON per layer (lon/lat, GDA94) - clipped to the analysis circle
        fc = {"type": "FeatureCollection", "name": L["dxf"],
              "features": [{"type": "Feature", "properties": it["props"],
                            "geometry": mapping(_to_lonlat(fr, it["geom"]))} for it in items]}
        (gj_dir / f"{L['key'].lower()}.geojson").write_text(json.dumps(fc), encoding="utf-8")
        _register_layer(reg, L, items, site, src, stamp)
        log(f"  {L['key']:9s} {len(items):3d} features")
    return layers


def layer_schema(net, layer_url):
    """Coded-value domains and date fields for one layer.  GeoJSON output returns raw domain codes
    (e.g. height map code 46 instead of J2) and dates as epoch milliseconds, so both are decoded here."""
    try:
        data, _, _, _ = net.get_json(layer_url, {"f": "json"})
    except Exception:
        return {}, set()
    dom, dates = {}, set()
    for f in data.get("fields") or []:
        d = f.get("domain") or {}
        if d.get("type") == "codedValue":
            dom[f["name"]] = {str(c.get("code")): c.get("name") for c in d.get("codedValues", [])}
        if f.get("type") == "esriFieldTypeDate":
            dates.add(f["name"])
    if data.get("subtypeField") and data.get("types"):          # e.g. height map code 46 -> "J2"
        dom[data["subtypeField"]] = {str(t.get("id")): t.get("name") for t in data["types"]}
    return dom, dates


def decode_props(p, dom, dates):
    for k, v in list(p.items()):
        if v is None:
            continue
        if k in dom and str(v) in dom[k]:
            p[k] = dom[k][str(v)]
        elif (k in dates or k.upper().endswith("_DATE")) and isinstance(v, (int, float)) and v > 1e11:
            p[k] = dt.datetime.fromtimestamp(v / 1000, dt.timezone.utc).date().isoformat()
    return p


def flood_coverage_note(net, layers, site, reg):
    """Most NSW councils do not publish a flood planning map on the Planning Portal (checked 23 Sep 2026:
    only about a dozen LGAs do). 'Nothing mapped' must not read as 'not flood affected'."""
    fl = layers.get("FLOOD")
    if not fl or any(it["geom"].intersects(site) for it in fl["items"]):
        return
    lga = next((it["props"].get("LGA_NAME") for it in layers.get("ZONE", {}).get("items", [])
                if it["geom"].intersects(site) and it["props"].get("LGA_NAME")), None)
    if not lga:
        return
    url = f"{fl['cfg']['url']}/{fl['cfg']['id']}/query"
    try:
        data, full, stamp, _ = net.get_json(url, {"where": f"LGA_NAME='{lga}'", "returnCountOnly": "true", "f": "json"})
        n = int(data.get("count", 0))
    except Exception as e:
        reg.add("FLOOD", f"Could not check whether {lga.title()} publishes a flood planning map ({e}).", url, TODAY, "NOT CHECKED")
        return
    if n == 0:
        reg.add("FLOOD", f"{lga.title()} council publishes NO flood planning map on the NSW Planning Portal, so the empty "
                "flood layer above says nothing about this site. Flood affectation must come from the council's flood "
                "study / flood planning level, or the section 10.7 certificate.", full, stamp, "NOT CHECKED")
    else:
        reg.add("FLOOD", f"{lga.title()} does publish a flood planning map ({n} mapped areas); none of them covers the site.",
                full, stamp, "MEDIUM")


def _to_lonlat(fr, g):
    from shapely.ops import transform
    return transform(lambda x, y, z=None: fr.inv.transform(x, y), g)


def _register_layer(reg, L, items, site, src, stamp):
    on = [it for it in items if it["geom"].intersection(site).area > max(1.0, 0.005 * site.area)]
    near = [it for it in items if it not in on]
    if on:
        vals = []
        for it in on:
            p = it["props"]
            share = it["geom"].intersection(site).area / site.area * 100
            inst = fmt("{EPI_NAME}", p)
            amend = fmt("{AMENDMENT}", p)
            cur = fmt("{CURRENCY_DATE}", p)
            if cur in ("?", "") and p.get("LastUpdate"):
                cur = str(p.get("LastUpdate"))
            tail = "; ".join(x for x in [inst if inst != "?" else "", amend if amend != "?" else "",
                                         f"currency {cur}" if cur not in ("?", "") else ""] if x)
            a = f"Site is {fmt(L['fact'], p)}" + (f" - {tail}" if tail else "")
            if share < 99:
                a += f". Covers about {share:.0f}% of the site (computed)"
            conf = "HIGH" if (cur not in ("?", "") and (inst not in ("?", "") or L["key"] == "BUSHFIRE")) else "MEDIUM"
            reg.add(L["key"], a, src, stamp, conf)
            vals.append(fmt(L["value"], p))
        if len(set(vals)) > 1:
            reg.add(L["key"], f"Site is split across {len(set(vals))} {L['title'].lower()} values: {', '.join(sorted(set(vals)))}. "
                    "Check which boundary governs which part.", src, stamp, "LOW")
    else:
        reg.add(L["key"], f"Checked: no {L['title'].lower()} mapped over the site in this dataset. "
                "Absent is not the same as unconstrained.", src, stamp, "MEDIUM")
    if near:
        if L["key"] in ("HERITAGE", "SHR"):
            for it in sorted(near, key=lambda it: it["geom"].distance(site))[:12]:
                d = it["geom"].distance(site)
                reg.add(L["key"], f"Within the study radius, {d:.0f} m from the site boundary (computed): "
                        f"{fmt(L['fact'], it['props'])}", src, stamp, "LOW")
        else:
            vals = sorted({fmt(L["value"], it["props"]) for it in near})
            reg.add(L["key"], f"Other {L['title'].lower()} values inside the study radius (drawn, not asserted one by one): "
                    f"{', '.join(vals[:15])}{' ...' if len(vals) > 15 else ''}", src, stamp, "MEDIUM")


# ---------------------------------------------------------------------------
# 9. Climate (Open-Meteo, no key) and sun (computed)
# ---------------------------------------------------------------------------
def pull_climate(net, lat, lon, reg, years=5):
    end = dt.date(dt.date.today().year - 1, 12, 31)
    start = dt.date(end.year - years + 1, 1, 1)
    params = {"latitude": f"{lat:.4f}", "longitude": f"{lon:.4f}",
              "start_date": start.isoformat(), "end_date": end.isoformat(),
              "hourly": "wind_speed_10m,wind_direction_10m",
              "daily": "temperature_2m_max,temperature_2m_min,precipitation_sum",
              "timezone": "auto", "wind_speed_unit": "ms"}
    data, url, stamp, _ = net.get_json("https://archive-api.open-meteo.com/v1/archive", params)
    d = data["daily"]
    months = np.array([int(t[5:7]) for t in d["time"]])
    tmax = np.array(d["temperature_2m_max"], float)
    tmin = np.array(d["temperature_2m_min"], float)
    rain = np.array(d["precipitation_sum"], float)
    yrs = np.array([int(t[:4]) for t in d["time"]])
    clim = {"months": list(range(1, 13)),
            "tmax": [float(np.nanmean(tmax[months == m])) for m in range(1, 13)],
            "tmin": [float(np.nanmean(tmin[months == m])) for m in range(1, 13)],
            "rain": [float(np.nansum(rain[months == m]) / len(set(yrs))) for m in range(1, 13)]}
    h = data["hourly"]
    hm = np.array([int(t[5:7]) for t in h["time"]])
    ws = np.array(h["wind_speed_10m"], dtype=float)
    wd = np.array(h["wind_direction_10m"], dtype=float)
    ok = ~np.isnan(ws) & ~np.isnan(wd)
    clim["wind"] = {"all": _rose(wd[ok], ws[ok]),
                    "Dec-Feb": _rose(wd[ok & np.isin(hm, [12, 1, 2])], ws[ok & np.isin(hm, [12, 1, 2])]),
                    "Jun-Aug": _rose(wd[ok & np.isin(hm, [6, 7, 8])], ws[ok & np.isin(hm, [6, 7, 8])])}
    glat, glon = data.get("latitude"), data.get("longitude")
    gdist = _haversine(lat, lon, glat, glon) if glat is not None else float("nan")
    src_note = (f"Open-Meteo historical weather API (reanalysis), {start.year}-{end.year}, "
                f"grid cell at {glat:.3f}, {glon:.3f}, {gdist/1000:.1f} km from the pin, elevation {data.get('elevation')} m")
    hot = int(np.argmax(clim["tmax"])) + 1
    cold = int(np.argmin(clim["tmin"])) + 1
    reg.add("CLIMATE", f"Warmest month {_mon(hot)}: mean daily max {clim['tmax'][hot-1]:.1f} C. "
            f"Coolest month {_mon(cold)}: mean daily min {clim['tmin'][cold-1]:.1f} C. ({src_note})", url, stamp, "MEDIUM")
    reg.add("CLIMATE", f"Mean annual rainfall about {sum(clim['rain']):.0f} mm; wettest month "
            f"{_mon(int(np.argmax(clim['rain']))+1)}. ({src_note})", url, stamp, "MEDIUM")
    for k in ("Dec-Feb", "Jun-Aug"):
        r = clim["wind"][k]
        reg.add("WIND", f"{k}: most frequent wind from the {r['prevailing']} ({r['prev_share']:.0f}% of hours), "
                f"mean speed {r['mean']:.1f} m/s at 10 m. Reanalysis grid, not a site measurement: sea breezes, "
                f"terrain and buildings are not resolved. ({src_note})", url, stamp, "MEDIUM")
    reg.add("CLIMATE", "Open-Meteo data is CC BY 4.0 - keep the credit line on any diagram that uses it.",
            "https://open-meteo.com/en/licence", TODAY, "HIGH")
    return clim


_DIRS16 = ["N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"]


def _rose(wd, ws, n=16, bins=(0.5, 3, 6, 9, 100)):
    if len(wd) == 0:
        return {"freq": [[0] * n for _ in range(len(bins) - 1)], "calm": 0, "prevailing": "?", "prev_share": 0, "mean": 0}
    sec = ((wd + 360 / n / 2) % 360 // (360 / n)).astype(int)
    calm = ws < bins[0]
    freq = []
    for lo, hi in zip(bins[:-1], bins[1:]):
        m = (ws >= lo) & (ws < hi)
        freq.append([float((m & (sec == i)).sum()) / len(wd) * 100 for i in range(n)])
    tot = np.sum(freq, 0)
    p = int(np.argmax(tot))
    return {"freq": freq, "calm": float(calm.mean() * 100), "prevailing": _DIRS16[p],
            "prev_share": float(tot[p]), "mean": float(np.mean(ws)), "bins": list(bins)}


def _mon(m):
    return ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"][m - 1]


def _haversine(lat1, lon1, lat2, lon2):
    R = 6371000.0
    p1, p2 = math.radians(lat1), math.radians(lat2)
    dp, dl = p2 - p1, math.radians(lon2 - lon1)
    a = math.sin(dp / 2) ** 2 + math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
    return 2 * R * math.asin(math.sqrt(a))


def sun_position(lat, lon, when_utc: np.ndarray):
    """NOAA solar position. when_utc: array of datetime64[m]. Returns altitude, azimuth (deg, azimuth from true north)."""
    t = (when_utc - np.datetime64("2000-01-01T12:00")) / np.timedelta64(1, "D")
    jc = t / 36525.0
    L0 = (280.46646 + jc * (36000.76983 + jc * 0.0003032)) % 360
    M = 357.52911 + jc * (35999.05029 - 0.0001537 * jc)
    e = 0.016708634 - jc * (0.000042037 + 0.0000001267 * jc)
    Mr = np.radians(M)
    C = np.sin(Mr) * (1.914602 - jc * (0.004817 + 0.000014 * jc)) + np.sin(2 * Mr) * (0.019993 - 0.000101 * jc) + np.sin(3 * Mr) * 0.000289
    true_long = L0 + C
    omega = 125.04 - 1934.136 * jc
    app_long = true_long - 0.00569 - 0.00478 * np.sin(np.radians(omega))
    eps0 = 23 + (26 + (21.448 - jc * (46.815 + jc * (0.00059 - jc * 0.001813))) / 60) / 60
    eps = eps0 + 0.00256 * np.cos(np.radians(omega))
    decl = np.degrees(np.arcsin(np.sin(np.radians(eps)) * np.sin(np.radians(app_long))))
    y = np.tan(np.radians(eps / 2)) ** 2
    L0r = np.radians(L0)
    eqt = 4 * np.degrees(y * np.sin(2 * L0r) - 2 * e * np.sin(Mr) + 4 * e * y * np.sin(Mr) * np.cos(2 * L0r)
                         - 0.5 * y * y * np.sin(4 * L0r) - 1.25 * e * e * np.sin(2 * Mr))
    mins = (when_utc - when_utc.astype("datetime64[D]")) / np.timedelta64(1, "m")
    tst = (mins + eqt + 4 * lon) % 1440
    ha = tst / 4 - 180
    ha = np.where(ha < -180, ha + 360, ha)
    lr, dr, har = math.radians(lat), np.radians(decl), np.radians(ha)
    cosz = np.sin(lr) * np.sin(dr) + np.cos(lr) * np.cos(dr) * np.cos(har)
    zen = np.degrees(np.arccos(np.clip(cosz, -1, 1)))
    alt = 90 - zen
    az = (np.degrees(np.arctan2(np.sin(har), np.cos(har) * np.sin(lr) - np.tan(dr) * np.cos(lr))) + 180) % 360
    return alt, az


def sun_curves(lat, lon, utc_offset_h, year):
    """Sun-path curves for the 21st of each solstice/equinox month, and hour lines."""
    dates = {"21 Dec": (12, 21), "21 Mar / 23 Sep": (3, 21), "21 Jun": (6, 21),
             "21 Jan/Nov": (1, 21), "21 Feb/Oct": (2, 21), "21 Apr/Aug": (4, 21), "21 May/Jul": (5, 21)}
    curves = {}
    for name, (m, d) in dates.items():
        base = np.datetime64(f"{year}-{m:02d}-{d:02d}T00:00") - np.timedelta64(int(utc_offset_h * 60), "m")
        tt = base + np.arange(0, 24 * 60, 5).astype("timedelta64[m]")
        alt, az = sun_position(lat, lon, tt)
        curves[name] = (alt, az, (np.arange(0, 24 * 60, 5) / 60.0))
    hours = {}
    for hr in range(4, 21):
        alts, azs = [], []
        for m in range(1, 13):
            base = np.datetime64(f"{year}-{m:02d}-21T{hr:02d}:00") - np.timedelta64(int(utc_offset_h * 60), "m")
            a, z = sun_position(lat, lon, np.array([base]))
            alts.append(a[0]); azs.append(z[0])
        hours[hr] = (np.array(alts), np.array(azs))
    return curves, hours


def register_sun(reg, lat, lon, utc_offset_h, year):
    rows = []
    for label, (m, d) in [("winter solstice", (6, 21) if lat < 0 else (12, 21)),
                          ("summer solstice", (12, 21) if lat < 0 else (6, 21)),
                          ("equinox", (3, 21))]:
        base = np.datetime64(f"{year}-{m:02d}-{d:02d}T00:00") - np.timedelta64(int(utc_offset_h * 60), "m")
        tt = base + np.arange(0, 1440).astype("timedelta64[m]")
        alt, az = sun_position(lat, lon, tt)
        i = int(np.argmax(alt))
        up = np.where(alt > -0.833)[0]
        rise = up[0] / 60 if len(up) else float("nan")
        sets = up[-1] / 60 if len(up) else float("nan")
        rows.append(f"{label} ({d} {_mon(m)}): highest sun {alt[i]:.1f} deg above the horizon at "
                    f"{int(i/60):02d}:{i%60:02d}; up from about {int(rise):02d}:{int((rise%1)*60):02d} "
                    f"to {int(sets):02d}:{int((sets%1)*60):02d} (UTC{utc_offset_h:+g}, standard time)")
    for r in rows:
        reg.add("SUN", "Computed from latitude/longitude (NOAA solar position equations): " + r,
                "https://gml.noaa.gov/grad/solcalc/calcdetails.html", TODAY, "MEDIUM")


# ---------------------------------------------------------------------------
# 10. Style card
# ---------------------------------------------------------------------------
DEFAULT_STYLE = {
    "name": "Placeholder - replace with your own",
    "type": {"family": "DejaVu Sans", "size": 7.0, "title_size": 11.0},
    "palette": {"background": "#F3F0E8", "ink": "#1D1D1B", "base": "#A7A399", "site": "#C8102E",
                "grid": "#D9D4C7"},
    "layers": {"ZONE": "#D8C9A8", "HOB": "#5B7A99", "FSR": "#7F9C6B", "ADDCTRL": "#9C6B98",
               "HERITAGE": "#8C4A2F", "SHR": "#8C4A2F", "FLOOD": "#2F6690", "BUSHFIRE": "#D1495B",
               "ASS": "#B08D57", "SUN": "#E0A100", "WIND": "#3A6EA5"},
    "zones": {"R": "#E9C9B5", "B": "#8FB3D9", "E": "#8FB3D9", "MU": "#B39DDB", "IN": "#C5B0D5",
              "SP": "#F2E394", "RE": "#A8D5A2", "C": "#7FB77E", "W": "#9EC9E2", "RU": "#E6DDB8", "SU": "#CCCCCC"},
    "weights": {"base": 0.15, "cadastre": 0.12, "statutory": 0.35, "site": 0.9, "axon_drop": 0.2},
    "hatch": {"ZONE": "solid", "HOB": "none", "FSR": "none", "ADDCTRL": "..", "HERITAGE": "////",
              "SHR": "xxxx", "FLOOD": "\\\\\\\\", "BUSHFIRE": "xx", "ASS": "..."},
    "fill_alpha": 0.55,
    "axon": {"angle": 45.0, "foreshorten": 0.6, "gap_factor": 0.8},
    "plan_layers": ["ZONE", "HERITAGE", "SHR", "FLOOD", "BUSHFIRE", "HOB"],
    "axon_layers": ["BASE", "CADASTRE", "ZONE", "HOB", "FSR", "HERITAGE", "HAZARD", "CLIMATE"],
    "references": {"images": ["ref_1.jpg", "ref_2.jpg", "ref_3.jpg"],
                   "words": "measured, archival, warm paper ground, hairline ink, restrained colour"},
    "dxf": {"hatches": False, "text_height_m": 2.0},
}


def load_style(path: Path | None) -> dict:
    st = json.loads(json.dumps(DEFAULT_STYLE))
    if path and Path(path).exists():
        if tomllib is None:
            raise RuntimeError("style card needs Python 3.11+ (tomllib)")
        with open(path, "rb") as f:
            user = tomllib.load(f)
        for tbl in list(user.values()):          # forgive top-level keys typed under a [table]
            if isinstance(tbl, dict):
                for k in ("fill_alpha", "plan_layers", "axon_layers", "name"):
                    if k in tbl and not isinstance(DEFAULT_STYLE.get(k), dict):
                        user.setdefault(k, tbl.pop(k))
        for k, v in user.items():
            if isinstance(v, dict) and isinstance(st.get(k), dict):
                st[k].update(v)
            else:
                st[k] = v
    fams = {f.name for f in matplotlib.font_manager.fontManager.ttflist}
    fam = st["type"]["family"]
    st["type"]["_family_used"] = fam if fam in fams else "DejaVu Sans"
    return st


def prompt_suffix(st) -> str:
    pal = st["palette"]
    lay = ", ".join(f"{k.lower()} {v}" for k, v in st["layers"].items())
    refs = ", ".join(st["references"]["images"])
    return (f"Style: {st['references']['words']}. Ground {pal['background']}, ink {pal['ink']}, "
            f"context linework {pal['base']}, site outline {pal['site']}. Layer colours: {lay}. "
            f"Line weights from {st['weights']['base']} mm hairline to {st['weights']['site']} mm for the site. "
            f"Hatching: " + ", ".join(f"{k.lower()} '{v}'" for k, v in st["hatch"].items() if v not in ("none",)) +
            f". Typeface: {st['type']['family']}. Match the attached reference images ({refs}). "
            f"Do not invent data, labels or north points that are not in the source drawing.")


def _mpl_hatch(h):
    return None if h in ("none", "solid", "", None) else h


def zone_colour(st, code):
    """Colour by zone family (R, B, E, IN ...), stepped darker by the zone number (R1 lighter than R4)."""
    code = str(code or "").upper()
    for pre in sorted(st["zones"], key=len, reverse=True):
        if code.startswith(pre):
            base = matplotlib.colors.to_rgb(st["zones"][pre])
            m = re.search(r"(\d)", code[len(pre):])
            k = min(3, max(-1, int(m.group(1)) - 2)) * 0.06 if m else 0.0
            return tuple(max(0.0, min(1.0, c * (1 - k))) for c in base)
    return st["layers"]["ZONE"]


def _lw(mm):  # mm at print -> points
    return mm * 72 / 25.4


# ---------------------------------------------------------------------------
# 11. DXF write-back
# ---------------------------------------------------------------------------
def _true_color(hexs):
    h = hexs.lstrip("#")
    return ezdxf.colors.rgb2int((int(h[0:2], 16), int(h[2:4], 16), int(h[4:6], 16)))


def write_dxf(out_path, base_doc, al: Alignment, site, cad, layers, st, sun, clim, radius, lat, lon, utc_offset):
    doc = base_doc if base_doc is not None else ezdxf.new("R2018", setup=True)
    if base_doc is None:
        doc.header["$INSUNITS"] = 6
    msp = doc.modelspace()
    unit = 1.0 / al.scale_to_m()    # metres per drawing unit (approx)
    th = st["dxf"]["text_height_m"] / unit
    colours = dict(st["layers"], SITE=st["palette"]["site"], CADASTRE=st["palette"]["base"],
                   NORTH=st["palette"]["ink"], RADIUS=st["palette"]["base"])

    def layer(name, key):
        if name not in doc.layers:
            L = doc.layers.add(name)
            L.true_color = _true_color(colours.get(key, "#000000"))
        return name

    def poly(g, lay, closed=True):
        for ring in _rings(g):
            c = np.asarray(ring.coords)[:, :2]
            P = al.pt(c[:, 0], c[:, 1])
            msp.add_lwpolyline(P.tolist(), close=closed, dxfattribs={"layer": lay})

    def text(s, x, y, lay, h=th, rot=0.0):
        P = al.pt([x], [y])[0]
        t = msp.add_text(s, height=h, rotation=rot, dxfattribs={"layer": lay})
        t.set_placement((P[0], P[1]), align=ezdxf.enums.TextEntityAlignment.MIDDLE_CENTER)

    rot = al.rotation_deg
    layer("AI_RADIUS", "RADIUS")
    poly(Point(0, 0).buffer(radius, 128), "AI_RADIUS")
    if base_doc is None or not al.ok:
        layer("AI_CADASTRE", "CADASTRE")
        for g, _ in cad:
            poly(g.intersection(Point(0, 0).buffer(radius, 128)), "AI_CADASTRE")
    layer("AI_SITE", "SITE")
    poly(site, "AI_SITE")
    for key, Ld in layers.items():
        name = layer(Ld["cfg"]["dxf"], key)
        for it in Ld["items"]:
            poly(it["geom"], name)
            if st["dxf"].get("hatches") and st["hatch"].get(key, "none") != "none":
                for p in _polys(it["geom"]):
                    hch = msp.add_hatch(dxfattribs={"layer": name})
                    hch.set_pattern_fill("ANSI31", scale=1.0 / unit * 0.5)
                    ext = al.pt(*np.asarray(p.exterior.coords)[:, :2].T)
                    hch.paths.add_polyline_path(ext.tolist(), is_closed=True)
            rp = it["geom"].representative_point()
            lab = fmt(Ld["cfg"]["label"], it["props"])
            if lab and lab != "?":
                text(lab, rp.x, rp.y, name, rot=rot)
    # north arrow (true north), outside the radius at 45 deg
    layer("AI_NORTH", "NORTH")
    nx, ny = radius * 0.8, radius * 0.8
    L = radius * 0.08
    arrow = Polygon([(nx, ny + L), (nx - L * 0.35, ny - L * 0.6), (nx, ny - L * 0.3), (nx + L * 0.35, ny - L * 0.6)])
    poly(arrow, "AI_NORTH")
    text("N (true)", nx, ny + L * 1.5, "AI_NORTH", rot=rot)
    # sun path (stereographic, horizontal plane, centred on the site)
    c = site.centroid
    R = max(25.0, min(radius * 0.35, math.sqrt(site.area) * 2.5))
    if sun is not None:
        lay_s = layer("AI_SUNPATH", "SUN")
        curves, hours = sun
        poly(Point(c.x, c.y).buffer(R, 128), lay_s)
        for alt in (30, 60):
            poly(Point(c.x, c.y).buffer(R * math.tan(math.radians(90 - alt) / 2), 128), lay_s)
        for name, (alt, az, _) in curves.items():
            m = alt > 0
            if m.sum() < 2:
                continue
            r = R * np.tan(np.radians(90 - alt[m]) / 2)
            xs = c.x + r * np.sin(np.radians(az[m]))
            ys = c.y + r * np.cos(np.radians(az[m]))
            P = al.pt(xs, ys)
            msp.add_lwpolyline(P.tolist(), dxfattribs={"layer": lay_s})
            k = int(np.argmax(alt[m]))
            text(name, xs[k], ys[k] - th * unit * 0.9, lay_s, h=th * 0.6, rot=rot)
        for hr, (alt, az) in hours.items():
            m = alt > 0
            if m.sum() < 2:
                continue
            r = R * np.tan(np.radians(90 - alt[m]) / 2)
            xs, ys = c.x + r * np.sin(np.radians(az[m])), c.y + r * np.cos(np.radians(az[m]))
            msp.add_lwpolyline(al.pt(xs, ys).tolist(), dxfattribs={"layer": lay_s})
            text(f"{hr}", xs[0], ys[0], lay_s, h=th * 0.6, rot=rot)
    if clim is not None:
        lay_w = layer("AI_WIND", "WIND")
        rose = clim["wind"]["all"]
        tot = np.sum(rose["freq"], 0)
        mx = max(tot.max(), 1e-6)
        n = len(tot)
        for i, f in enumerate(tot):
            a0 = math.radians(i * 360 / n - 360 / n / 2 * 0.8)
            a1 = math.radians(i * 360 / n + 360 / n / 2 * 0.8)
            rr = R * 1.25 * f / mx
            pts = [(c.x, c.y)] + [(c.x + rr * math.sin(a), c.y + rr * math.cos(a)) for a in np.linspace(a0, a1, 6)]
            poly(Polygon(pts), lay_w)
        text(f"wind from (all hours): {rose['prevailing']}", c.x, c.y - R * 1.4, lay_w, h=th * 0.7, rot=rot)
    doc.saveas(out_path)


# ---------------------------------------------------------------------------
# 12. Diagrams
# ---------------------------------------------------------------------------
def _fig(st, w=210, h=297):
    fig = plt.figure(figsize=(w / 25.4, h / 25.4), dpi=200)
    fig.patch.set_facecolor(st["palette"]["background"])
    plt.rcParams["font.family"] = st["type"]["_family_used"]
    plt.rcParams["font.size"] = st["type"]["size"]
    return fig


def _draw_geom(ax, g, T, fc, ec, lw, hatch=None, alpha=1.0, z=1):
    for p in _polys(g):
        ext = T(np.asarray(p.exterior.coords)[:, :2])
        verts = [ext] + [T(np.asarray(r.coords)[:, :2]) for r in p.interiors]
        codes = []
        allv = []
        for v in verts:
            codes += [MplPath.MOVETO] + [MplPath.LINETO] * (len(v) - 2) + [MplPath.CLOSEPOLY]
            allv.append(v)
        path = MplPath(np.vstack(allv), codes)
        ax.add_patch(PathPatch(path, facecolor=fc if fc else "none", edgecolor=ec, lw=lw,
                               hatch=hatch, alpha=alpha, zorder=z))


def base_in_local(base_lines, al: Alignment | None, radius):
    """archiMap linework mapped back into local metres, clipped (by bbox) to the radius."""
    if not base_lines or al is None:
        return {}
    Rinv = np.linalg.inv(al.R)
    # invert: dxf -> after base_fn.  Only possible cleanly for linear base_fn; do numerically.
    probe = np.array([[0, 0], [100, 0], [0, 100]], float)
    Q = al.pt(probe[:, 0], probe[:, 1])
    A = np.linalg.lstsq(np.c_[probe, np.ones(3)], Q, rcond=None)[0]   # local -> dxf affine (3x2)
    M = A[:2].T
    b = A[2]
    Minv = np.linalg.inv(M)
    out = {}
    for lay, lines in base_lines.items():
        role = role_of(lay)
        keep = []
        for v in lines:
            loc = (v - b) @ Minv.T
            if np.any(np.hypot(loc[:, 0], loc[:, 1]) < radius * 1.05):
                keep.append(loc)
        if keep:
            out.setdefault(role, []).extend(keep)
    return out


def plan_diagram(path_stem, st, site, cad, layers, base_local, radius, title, subtitle, credits):
    fig = _fig(st, 297, 297)
    ax = fig.add_axes([0.05, 0.12, 0.9, 0.8])
    ax.set_facecolor(st["palette"]["background"])
    T = lambda a: a
    clipc = MplPolygon(np.c_[radius * np.cos(np.linspace(0, 2 * np.pi, 256)),
                             radius * np.sin(np.linspace(0, 2 * np.pi, 256))], closed=True, transform=ax.transData)
    handles = []
    # statutory fills first
    for key in st["plan_layers"]:
        if key not in layers:
            continue
        Ld = layers[key]
        col = st["layers"].get(key, "#888888")
        hatch = _mpl_hatch(st["hatch"].get(key))
        for it in Ld["items"]:
            if key == "ZONE":
                fc = zone_colour(st, it["props"].get("SYM_CODE"))
                _draw_geom(ax, it["geom"], T, fc, st["palette"]["background"], _lw(0.1), None, st["fill_alpha"], 1)
            elif key == "HOB":
                _draw_geom(ax, it["geom"], T, None, col, _lw(st["weights"]["statutory"]), None, 1, 3)
            else:
                plt.rcParams["hatch.color"] = col
                _draw_geom(ax, it["geom"], T, None, col, _lw(st["weights"]["statutory"]), hatch, 0.9, 2)
        if Ld["items"]:
            handles.append((key, col, Ld["cfg"]["title"]))
    # base linework
    if base_local:
        for role, lines in base_local.items():
            w = st["weights"]["base"] * (2.0 if role == "buildings" else 1.0)
            lc = LineCollection(lines, colors=st["palette"]["base"] if role != "buildings" else st["palette"]["ink"],
                                linewidths=_lw(w), zorder=4, alpha=0.9 if role == "buildings" else 0.6)
            lc.set_clip_path(clipc)
            ax.add_collection(lc)
    else:
        lines = [np.asarray(r.coords) for g, _ in cad for r in _rings(g)]
        lc = LineCollection(lines, colors=st["palette"]["base"], linewidths=_lw(st["weights"]["cadastre"]), zorder=4)
        lc.set_clip_path(clipc)
        ax.add_collection(lc)
    # labels
    for key in ("ZONE", "HOB"):
        if key in layers:
            for it in layers[key]["items"]:
                if it["geom"].area < (radius * 0.08) ** 2:
                    continue
                rp = it["geom"].representative_point()
                lab = fmt(layers[key]["cfg"]["label"], it["props"])
                ax.text(rp.x, rp.y - (radius * 0.03 if key == "HOB" else 0), lab, ha="center", va="center",
                        fontsize=st["type"]["size"] * (1.1 if key == "ZONE" else 0.85),
                        color=st["palette"]["ink"] if key == "ZONE" else st["layers"]["HOB"], zorder=6)
    _draw_geom(ax, site, T, None, st["palette"]["site"], _lw(st["weights"]["site"]), None, 1, 8)
    ax.add_patch(MplPolygon(clipc.get_xy(), closed=True, fill=False, ec=st["palette"]["ink"], lw=_lw(0.25), zorder=9))
    _north_scale(ax, st, radius)
    ax.set_xlim(-radius * 1.08, radius * 1.08)
    ax.set_ylim(-radius * 1.08, radius * 1.08)
    ax.set_aspect("equal")
    ax.axis("off")
    _legend(fig, st, handles)
    _titles(fig, st, title, subtitle, credits)
    _save(fig, path_stem)


def _north_scale(ax, st, radius):
    x, y, L = radius * 0.92, radius * 0.92, radius * 0.07
    ax.add_patch(MplPolygon([(x, y + L), (x - L * 0.35, y - L * 0.6), (x, y - L * 0.3), (x + L * 0.35, y - L * 0.6)],
                            closed=True, fc=st["palette"]["ink"], ec="none", zorder=10))
    ax.text(x, y + L * 1.4, "N", ha="center", va="bottom", fontsize=st["type"]["size"] * 1.2, color=st["palette"]["ink"])
    nice = [10, 20, 25, 50, 100, 200, 250, 500]
    sb = min(nice, key=lambda n: abs(n - radius * 0.4))
    x0, y0 = -radius, -radius * 1.02
    for i, (a, b) in enumerate([(0, sb / 2), (sb / 2, sb)]):
        ax.add_patch(plt.Rectangle((x0 + a, y0), b - a, radius * 0.012, fc=st["palette"]["ink"] if i == 0 else "none",
                                   ec=st["palette"]["ink"], lw=_lw(0.2), zorder=10))
    for v in (0, sb / 2, sb):
        ax.text(x0 + v, y0 - radius * 0.02, f"{v:g}", ha="center", va="top", fontsize=st["type"]["size"] * 0.8,
                color=st["palette"]["ink"])
    ax.text(x0 + sb + radius * 0.02, y0, "m", va="bottom", fontsize=st["type"]["size"] * 0.8, color=st["palette"]["ink"])


def _legend(fig, st, handles):
    ax = fig.add_axes([0.05, 0.02, 0.9, 0.08])
    ax.axis("off")
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    items = [("SITE", st["palette"]["site"], "Site (cadastral lot)")] + handles
    for i, (k, col, t) in enumerate(items):
        cx, cy = (i % 4) * 0.25, 0.75 - (i // 4) * 0.35
        hatch = _mpl_hatch(st["hatch"].get(k)) if k not in ("ZONE", "SITE", "HOB") else None
        fc = col if k == "ZONE" else "none"
        plt.rcParams["hatch.color"] = col
        ax.add_patch(plt.Rectangle((cx, cy - 0.08), 0.03, 0.16, fc=fc, ec=col, hatch=hatch,
                                   lw=_lw(0.6 if k == "SITE" else 0.3), alpha=st["fill_alpha"] if k == "ZONE" else 1))
        ax.text(cx + 0.04, cy, t, va="center", fontsize=st["type"]["size"], color=st["palette"]["ink"])


def _titles(fig, st, title, subtitle, credits):
    fig.text(0.05, 0.965, title, fontsize=st["type"]["title_size"], color=st["palette"]["ink"], weight="bold")
    fig.text(0.05, 0.945, subtitle, fontsize=st["type"]["size"], color=st["palette"]["ink"])
    fig.text(0.95, 0.008, credits, fontsize=st["type"]["size"] * 0.7, color=st["palette"]["base"], ha="right")


def _save(fig, stem):
    fig.savefig(f"{stem}.svg", facecolor=fig.get_facecolor())
    fig.savefig(f"{stem}.png", dpi=200, facecolor=fig.get_facecolor())
    plt.close(fig)


def exploded_diagram(path_stem, st, site, cad, layers, base_local, radius, sun, clim, title, subtitle, credits):
    """Stacked 2.5D exploded axonometric (plan-oblique): every layer is a plane."""
    ang = math.radians(st["axon"]["angle"])
    k = st["axon"]["foreshorten"]
    stack = []
    for key in st["axon_layers"]:
        if key == "BASE":
            stack.append(("BASE", "Context (archiMap)" if base_local else "Context (cadastre)"))
        elif key == "CADASTRE":
            stack.append(("CADASTRE", "Cadastre and site"))
        elif key == "HAZARD":
            if any(h in layers for h in ("FLOOD", "BUSHFIRE", "ASS")):
                stack.append(("HAZARD", "Hazards: flood, bushfire, acid sulfate"))
        elif key == "CLIMATE":
            if sun is not None or clim is not None:
                stack.append(("CLIMATE", "Sun path and wind"))
        elif key in layers:
            stack.append((key, layers[key]["cfg"]["title"]))
    n = len(stack)
    gap = radius * 2 * k * st["axon"]["gap_factor"]
    xl = (-radius * 1.12, radius * 3.1)
    yl = (-radius * k * 1.12, gap * (n - 1) + radius * k * 1.45)
    W = 297.0
    Hmm = W * 0.94 / 0.86 * (yl[1] - yl[0]) / (xl[1] - xl[0])
    fig = _fig(st, W, max(210.0, Hmm))
    ax = fig.add_axes([0.03, 0.05, 0.94, 0.86])
    ax.set_facecolor(st["palette"]["background"])

    def T(a, z):
        a = np.asarray(a, float)
        x = a[:, 0] * math.cos(ang) - a[:, 1] * math.sin(ang)
        y = (a[:, 0] * math.sin(ang) + a[:, 1] * math.cos(ang)) * k + z
        return np.c_[x, y]

    circ = np.c_[radius * np.cos(np.linspace(0, 2 * np.pi, 256)), radius * np.sin(np.linspace(0, 2 * np.pi, 256))]
    site_pts = np.asarray(site.exterior.coords) if isinstance(site, Polygon) else np.asarray(site.convex_hull.exterior.coords)
    ink, bg = st["palette"]["ink"], st["palette"]["background"]
    zs = []
    for i, (key, label) in enumerate(stack):
        z = gap * (n - 1 - i)
        zs.append(z)
        zo = 10 + (n - i) * 10
        disc = T(circ, z)
        ax.add_patch(MplPolygon(disc, closed=True, fc=bg, ec=ink, lw=_lw(0.25), zorder=zo, alpha=0.96))
        clip = MplPolygon(disc, closed=True, transform=ax.transData)
        TT = lambda a, z=z: T(a, z)
        if key == "BASE":
            src = base_local.items() if base_local else [("cadastre", [np.asarray(r.coords) for g, _ in cad for r in _rings(g)])]
            for role, lines in src:
                lc = LineCollection([TT(v) for v in lines], colors=ink if role == "buildings" else st["palette"]["base"],
                                    linewidths=_lw(st["weights"]["base"] * (1.6 if role == "buildings" else 1)), zorder=zo + 1)
                lc.set_clip_path(clip)
                ax.add_collection(lc)
        elif key == "CADASTRE":
            lc = LineCollection([TT(np.asarray(r.coords)) for g, _ in cad for r in _rings(g)], colors=st["palette"]["base"],
                                linewidths=_lw(st["weights"]["cadastre"]), zorder=zo + 1)
            lc.set_clip_path(clip)
            ax.add_collection(lc)
        elif key == "HAZARD":
            for hk in ("FLOOD", "BUSHFIRE", "ASS"):
                if hk in layers:
                    col = st["layers"][hk]
                    plt.rcParams["hatch.color"] = col
                    for it in layers[hk]["items"]:
                        _draw_geom(ax, it["geom"], TT, None, col, _lw(st["weights"]["statutory"]),
                                   _mpl_hatch(st["hatch"].get(hk)), 0.9, zo + 1)
        elif key == "CLIMATE":
            c = site.centroid
            R = radius * 0.55
            if sun is not None:
                curves, hours = sun
                for nm, (alt, az, _) in curves.items():
                    m = alt > 0
                    r = R * np.tan(np.radians(90 - alt[m]) / 2)
                    pts = np.c_[c.x + r * np.sin(np.radians(az[m])), c.y + r * np.cos(np.radians(az[m]))]
                    if len(pts) > 1:
                        ax.plot(*TT(pts).T, color=st["layers"]["SUN"], lw=_lw(0.35), zorder=zo + 2)
            if clim is not None:
                tot = np.sum(clim["wind"]["all"]["freq"], 0)
                mx = max(tot.max(), 1e-6)
                for j, f in enumerate(tot):
                    a0, a1 = math.radians(j * 22.5 - 9), math.radians(j * 22.5 + 9)
                    rr = R * 0.9 * f / mx
                    pts = np.array([(c.x, c.y)] + [(c.x + rr * math.sin(a), c.y + rr * math.cos(a)) for a in np.linspace(a0, a1, 6)])
                    ax.add_patch(MplPolygon(TT(pts), closed=True, fc=st["layers"]["WIND"], ec="none", alpha=0.5, zorder=zo + 1))
        else:
            col = st["layers"].get(key, "#888888")
            for it in layers[key]["items"]:
                if key == "ZONE":
                    _draw_geom(ax, it["geom"], TT, zone_colour(st, it["props"].get("SYM_CODE")), bg, _lw(0.1), None,
                               st["fill_alpha"], zo + 1)
                elif key in ("HOB", "FSR"):
                    vals = [float(re.sub(r"[^0-9.]", "", str(fmt(layers[key]["cfg"]["value"], x["props"]))) or 0)
                            for x in layers[key]["items"]]
                    v = float(re.sub(r"[^0-9.]", "", str(fmt(layers[key]["cfg"]["value"], it["props"]))) or 0)
                    lo, hi = (min(vals), max(vals)) if vals else (0, 1)
                    a = 0.15 + 0.65 * ((v - lo) / (hi - lo) if hi > lo else 0.5)
                    _draw_geom(ax, it["geom"], TT, col, bg, _lw(0.1), None, a, zo + 1)
                else:
                    plt.rcParams["hatch.color"] = col
                    _draw_geom(ax, it["geom"], TT, None, col, _lw(st["weights"]["statutory"]),
                               _mpl_hatch(st["hatch"].get(key)), 0.9, zo + 1)
        # site on every plane
        ax.add_patch(MplPolygon(TT(site_pts[:, :2]), closed=True, fill=False, ec=st["palette"]["site"],
                                lw=_lw(st["weights"]["site"] * 0.7), zorder=zo + 5))
        # label, right-hand side with leader
        edge = T(np.array([[radius * math.cos(math.radians(-45 + 90 - st['axon']['angle'] + 45)),
                            radius * math.sin(math.radians(-45 + 90 - st['axon']['angle'] + 45))]]), z)[0]
        lx = radius * 1.25
        ax.plot([edge[0], lx], [edge[1], z], color=ink, lw=_lw(0.15), zorder=zo + 6)
        ax.text(lx + radius * 0.03, z, f"{i+1:02d}  {label}", va="center", fontsize=st["type"]["size"] * 1.1,
                color=ink, zorder=zo + 6)
    # drop lines through the site corners
    corners = site_pts[:-1] if len(site_pts) <= 9 else site.minimum_rotated_rectangle.exterior.coords[:-1]
    for p in np.asarray(corners)[:, :2]:
        top = T(p[None, :], zs[0])[0]
        bot = T(p[None, :], zs[-1])[0]
        ax.plot([top[0], bot[0]], [top[1], bot[1]], color=st["palette"]["site"], lw=_lw(st["weights"]["axon_drop"]),
                ls=(0, (4, 3)), zorder=1000)
    # north arrow on the top plane
    na = T(np.array([[0, radius * 1.1], [0, radius * 1.3]]), zs[0])
    ax.annotate("", xy=na[1], xytext=na[0], arrowprops=dict(arrowstyle="-|>", color=ink, lw=_lw(0.3)), zorder=1001)
    ax.text(*na[1], " N", color=ink, fontsize=st["type"]["size"], zorder=1001)
    ax.set_xlim(*xl)
    ax.set_ylim(*yl)
    ax.set_aspect("equal", adjustable="box", anchor="W")
    ax.axis("off")
    _titles(fig, st, title, subtitle, credits)
    _save(fig, path_stem)


def climate_diagram(path_stem, st, lat, lon, sun, clim, title, credits):
    fig = _fig(st, 297, 210)
    ink = st["palette"]["ink"]
    if sun is not None:
        ax = fig.add_axes([0.03, 0.12, 0.36, 0.72], projection="polar")
        ax.set_facecolor(st["palette"]["background"])
        ax.set_theta_zero_location("N")
        ax.set_theta_direction(-1)
        curves, hours = sun
        for nm, (alt, az, _) in curves.items():
            m = alt > 0
            r = np.tan(np.radians(90 - alt[m]) / 2)
            main = nm in ("21 Dec", "21 Jun", "21 Mar / 23 Sep")
            ax.plot(np.radians(az[m]), r, color=st["layers"]["SUN"], lw=_lw(0.5 if main else 0.2))
            if main and m.any():
                j = int(np.argmax(alt[m]))
                ax.text(np.radians(az[m][j]), r[j], nm, fontsize=st["type"]["size"] * 0.8, color=ink, ha="center", va="bottom")
        for hr, (alt, az) in hours.items():
            m = alt > 0
            if m.sum() > 1:
                r = np.tan(np.radians(90 - alt[m]) / 2)
                ax.plot(np.radians(az[m]), r, color=st["palette"]["base"], lw=_lw(0.15))
                ax.text(np.radians(az[m][0]), r[0], f"{hr}", fontsize=st["type"]["size"] * 0.7, color=ink)
        ax.set_rmax(1.0)
        ax.set_yticks([math.tan(math.radians(90 - a) / 2) for a in (60, 30)])
        ax.set_yticklabels(["60°", "30°"], fontsize=st["type"]["size"] * 0.7, color=ink)
        ax.set_xticks(np.radians([0, 90, 180, 270]))
        ax.set_xticklabels(["N", "E", "S", "W"], color=ink)
        ax.grid(color=st["palette"]["grid"], lw=_lw(0.15))
        ax.set_title("Sun path (stereographic, standard time)", fontsize=st["type"]["size"] * 1.1, color=ink, pad=14)
    if clim is not None:
        cols = plt.get_cmap("Blues")(np.linspace(0.35, 0.95, 4))
        for j, key in enumerate(("Dec-Feb", "Jun-Aug")):
            ax = fig.add_axes([0.44 + j * 0.25, 0.44, 0.17, 0.38], projection="polar")
            ax.set_facecolor(st["palette"]["background"])
            ax.set_theta_zero_location("N")
            ax.set_theta_direction(-1)
            r = clim["wind"][key]
            th = np.radians(np.arange(16) * 22.5)
            bottom = np.zeros(16)
            labels = ["0.5-3", "3-6", "6-9", "9+ m/s"]
            for b, f in enumerate(r["freq"]):
                ax.bar(th, f, width=np.radians(18), bottom=bottom, color=cols[b], edgecolor="none", label=labels[b])
                bottom += np.array(f)
            ax.set_xticks(np.radians([0, 90, 180, 270]))
            ax.set_xticklabels(["N", "E", "S", "W"], color=ink, fontsize=st["type"]["size"] * 0.8)
            ax.tick_params(axis="y", labelsize=st["type"]["size"] * 0.6, colors=ink)
            ax.grid(color=st["palette"]["grid"], lw=_lw(0.15))
            ax.set_title(f"Wind from, {key}\n(prevailing {r['prevailing']})", fontsize=st["type"]["size"], color=ink, pad=10)
            if j == 1:
                ax.legend(loc="lower left", bbox_to_anchor=(1.05, 0.0), fontsize=st["type"]["size"] * 0.7, frameon=False)
        ax = fig.add_axes([0.45, 0.1, 0.5, 0.22])
        ax.set_facecolor(st["palette"]["background"])
        m = np.arange(1, 13)
        ax.bar(m, clim["rain"], color=st["layers"]["FLOOD"], alpha=0.35, width=0.6)
        ax.set_ylabel("rain mm/month", color=ink, fontsize=st["type"]["size"] * 0.8)
        ax2 = ax.twinx()
        ax2.fill_between(m, clim["tmin"], clim["tmax"], color=st["layers"]["BUSHFIRE"], alpha=0.25, lw=0)
        ax2.plot(m, clim["tmax"], color=st["layers"]["BUSHFIRE"], lw=_lw(0.4))
        ax2.plot(m, clim["tmin"], color=st["layers"]["BUSHFIRE"], lw=_lw(0.4))
        ax2.set_ylabel("mean daily min-max °C", color=ink, fontsize=st["type"]["size"] * 0.8)
        ax.set_xticks(m)
        ax.set_xticklabels([_mon(i)[0] for i in m])
        for a in (ax, ax2):
            a.tick_params(colors=ink, labelsize=st["type"]["size"] * 0.7)
            for s in a.spines.values():
                s.set_color(st["palette"]["base"])
    _titles(fig, st, title, f"{lat:.5f}, {lon:.5f}", credits)
    _save(fig, path_stem)


# ---------------------------------------------------------------------------
# 13. The run
# ---------------------------------------------------------------------------
def run(args):
    t0 = time.time()
    out_root = Path(args.out)
    reg = Register()
    net = Net(Path(args.cache or out_root / "_cache"), offline=args.offline)
    jur = (args.jurisdiction or "NSW").upper()
    st = load_style(Path(args.style) if args.style else None)

    # --- pin
    if args.lat is not None and args.lon is not None:
        lat, lon = float(args.lat), float(args.lon)
        reg.add("SITE", f"Pin supplied by student at {lat:.6f}, {lon:.6f} (from archiMap or a map). "
                "Check it sits on the right lot in the diagrams.", "", TODAY, "LOW")
    else:
        lat, lon, gurl, method, stamp = geocode(net, args.address, jur)
        reg.add("SITE", f"Address '{args.address}' located at {lat:.6f}, {lon:.6f} by {method}.", gurl, stamp,
                "MEDIUM" if "NSW" in method else "LOW")
    slug = slugify(args.name or args.address or f"{lat:.4f}_{lon:.4f}")
    out = out_root / slug
    gj = out / f"layers_{slug}"
    gj.mkdir(parents=True, exist_ok=True)
    fr = Frame(lat, lon)
    radius = float(args.radius)
    log(f"Site Layers-overlay {VERSION}  |  {args.address or slug}  |  {lat:.6f}, {lon:.6f}  |  {jur}  |  radius {radius:.0f} m")

    adapter = ADAPTERS.get(jur, "missing")
    cad, site, site_how = [], Point(0, 0).buffer(10), "no cadastre"
    layers = {}
    if adapter in (None, "missing"):
        reg.add("STATUTORY", f"Statutory layers not available for jurisdiction '{jur}' in this tool. The adapter slot is "
                "empty: zoning, height, density, heritage and hazards must be found and verified by hand from the "
                "local planning authority.", "", TODAY, "NOT CHECKED")
        log(f"  no statutory adapter for {jur} - universal layers only")
    else:
        try:
            cad, cad_src, cad_stamp = pull_cadastre(net, adapter, fr, radius, reg)
            site, props, site_how = pick_site(cad, lot_ids=args.lots.split(",") if args.lots else None)
            lots = "; ".join(fmt(adapter["cadastre"]["label"], p) for p in props) or "no lot"
            plan_area = "; ".join(f"{p.get('planlotarea')} {p.get('planlotareaunits') or ''}".strip()
                                  for p in props if p.get("planlotarea"))
            reg.add("SITE", f"Site taken as {lots} ({site_how}). Area computed from the cadastre polygon about "
                    f"{site.area:.0f} m2" + (f"; area recorded on the plan: {plan_area}" if plan_area else
                    "; no plan area recorded in the cadastre") + ". Neither is title area or survey.",
                    cad_src, cad_stamp, "LOW")
            log(f"  CADASTRE  {len(cad):3d} lots; site = {lots} ({site_how})")
        except Exception as e:
            reg.add("SITE", f"NOT CHECKED - cadastre did not answer ({e}). The site outline is a 10 m circle at the pin.",
                    adapter["cadastre"]["url"], TODAY, "NOT CHECKED")
            log(f"  CADASTRE  FAILED {e}")
        layers = pull_statutory(net, adapter, fr, radius, site, reg, gj)
        flood_coverage_note(net, layers, site, reg)
        if cad:
            fc = {"type": "FeatureCollection", "name": "AI_CADASTRE", "features": [
                {"type": "Feature", "properties": p, "geometry": mapping(_to_lonlat(fr, g))} for g, p in cad]}
            (gj / "cadastre.geojson").write_text(json.dumps(fc), encoding="utf-8")
    # --- standing gaps: always in the register
    for a in ["Title area, survey, easements, covenants and restrictions (the s88B instrument) - not checked by this tool.",
              "Planning certificate (NSW: section 10.7) - not checked; only the council can certify it.",
              "Development control plan (setbacks, landscape, parking) - not read by this tool.",
              "Every statutory layer above is a point-in-time copy of a map service. Instruments change; check currency."]:
        reg.add("GAP", a, "", TODAY, "NOT CHECKED")

    # --- climate + sun
    clim = sun = None
    utc = args.utc_offset if args.utc_offset is not None else round(lon / 15)
    try:
        sun = sun_curves(lat, lon, utc, dt.date.today().year)
        register_sun(reg, lat, lon, utc, dt.date.today().year)
        log("  SUN       computed")
    except Exception as e:
        reg.add("SUN", f"NOT CHECKED - sun path computation failed ({e})", "", TODAY, "NOT CHECKED")
    try:
        clim = pull_climate(net, lat, lon, reg)
        (gj / "climate.json").write_text(json.dumps(clim), encoding="utf-8")
        log(f"  CLIMATE   ok, prevailing wind {clim['wind']['all']['prevailing']}")
    except Exception as e:
        reg.add("CLIMATE", f"NOT CHECKED - Open-Meteo did not answer ({e}). Use BoM climate statistics by hand.",
                "https://archive-api.open-meteo.com/v1/archive", TODAY, "NOT CHECKED")
        log(f"  CLIMATE   FAILED {e}")

    # --- archiMap DXF + alignment
    base_doc, al, base_local = None, None, {}
    if args.dxf:
        try:
            info, base_doc = inspect_dxf(Path(args.dxf))
            (out / "archimap_inspect.json").write_text(json.dumps(info, indent=1, default=str), encoding="utf-8")
            log(f"  DXF       {len(info['layers'])} layers, units {info['units']}, frame guess {info['frame_guess']}")
            reg.add("BASE", f"Base drawing: archiMap export '{Path(args.dxf).name}', {len(info['layers'])} layers "
                    f"({', '.join(list(info['layers'])[:10])}{'...' if len(info['layers']) > 10 else ''}), units "
                    f"{info['units']}. archiMap's own layers are passed through untouched and are not verified by this tool.",
                    "https://masslabs-archi.com/", TODAY, "LOW")
            if cad:
                al = align_to_dxf(base_doc, info, fr, lat, lon, [g for g, _ in cad], reg, args.cad_layer)
        except Exception as e:
            reg.add("BASE", f"Could not read the archiMap DXF ({e}). Pulled layers written to their own DXF.",
                    "", TODAY, "NOT CHECKED")
            log(f"  DXF       FAILED {e}")
            base_doc = None
    merged = base_doc is not None and al is not None and al.ok
    dxf_path = out / (f"site_{slug}.dxf" if merged or not args.dxf else f"site_{slug}_AI_only.dxf")
    try:
        write_dxf(dxf_path, base_doc if merged else None, al if merged else identity_alignment(),
                  site, cad, layers, st, sun, clim, radius, lat, lon, utc)
        log(f"  wrote {dxf_path.name}" + ("  (merged into archiMap drawing)" if merged else ""))
    except Exception as e:
        log(f"  DXF write FAILED {e}")
        traceback.print_exc()
    if merged:
        allb = dxf_polylines(base_doc, flatten=0.5)
        base_local = base_in_local({k: v for k, v in allb.items() if not k.upper().startswith("AI_")}, al, radius)

    # --- diagrams
    addr = args.address or f"{lat:.5f}, {lon:.5f}"
    credits = (f"Statutory: {ADAPTERS[jur]['name'] if isinstance(ADAPTERS.get(jur), dict) else 'none for ' + jur}, retrieved {TODAY}. "
               f"Climate: Open-Meteo.com (CC BY 4.0). Base: {'archiMap export' if base_local else 'NSW cadastre'}. "
               f"Site Layers-overlay {VERSION}. Not a planning certificate.")
    for nm, fn in [("diagram_plan", lambda: plan_diagram(out / "diagram_plan", st, site, cad, layers, base_local, radius,
                                                         f"Statutory overlay - {addr}", f"Study radius {radius:.0f} m. Style: {st['name']}", credits)),
                   ("diagram_exploded", lambda: exploded_diagram(out / "diagram_exploded", st, site, cad, layers, base_local, radius,
                                                                 sun, clim, f"Site layers - {addr}", "Exploded layer axonometric (stacked 2.5D)", credits)),
                   ("climate", lambda: climate_diagram(out / "climate", st, lat, lon, sun, clim, f"Sun and wind - {addr}", credits))]:
        try:
            fn()
            log(f"  wrote {nm}.svg/.png")
        except Exception as e:
            log(f"  {nm} FAILED {e}")
            traceback.print_exc()
    (out / "prompt_suffix.txt").write_text(prompt_suffix(st), encoding="utf-8")
    reg.write(out / f"register_{slug}.csv")
    log(f"  wrote register_{slug}.csv ({len(reg.rows)} rows, verdict column empty - that is your job)")
    zpath = out_root / f"site_layers_overlay_{slug}.zip"
    with zipfile.ZipFile(zpath, "w", zipfile.ZIP_DEFLATED) as z:
        for p in out.rglob("*"):
            if p.is_file():
                z.write(p, p.relative_to(out_root))
    log(f"Done in {time.time()-t0:.0f} s  ->  {zpath}")
    return out, reg


def main(argv=None):
    ap = argparse.ArgumentParser(description="Site Layers-overlay - statutory + climate layers for an archiMap DXF")
    ap.add_argument("--address", help="street address (used for labels, and to find the pin if no lat/lon)")
    ap.add_argument("--lat", type=float, help="pin latitude (copy from archiMap)")
    ap.add_argument("--lon", type=float, help="pin longitude (copy from archiMap)")
    ap.add_argument("--jurisdiction", default="NSW", help="NSW | VIC | QLD | NZ | SG | HK | other")
    ap.add_argument("--dxf", help="archiMap DXF export (optional - without it you get the pulled layers on their own)")
    ap.add_argument("--style", help="style card .toml")
    ap.add_argument("--radius", default=200, type=float, help="study radius in metres (default 200)")
    ap.add_argument("--lots", help="comma-separated lot ids if the site is several lots, e.g. 10//DP7949,11//DP7949")
    ap.add_argument("--cad-layer", help="name of the archiMap cadastre layer, if auto-detect picks the wrong one")
    ap.add_argument("--utc-offset", type=float, help="standard-time UTC offset for the sun path (default from longitude)")
    ap.add_argument("--name", help="short name for output files")
    ap.add_argument("--out", default="site_layers_overlay_out")
    ap.add_argument("--cache", help="cache folder (default <out>/_cache)")
    ap.add_argument("--offline", action="store_true", help="use cached responses only (class demo)")
    ap.add_argument("--inspect", help="only inspect an archiMap DXF and print what it contains")
    a = ap.parse_args(argv)
    if a.inspect:
        info, _ = inspect_dxf(Path(a.inspect))
        print(json.dumps(info, indent=1, default=str))
        return
    if a.lat is None and not a.address:
        ap.error("give --lat and --lon (preferred) or --address")
    run(a)


if __name__ == "__main__":
    main()