#!/usr/bin/env python3
"""
SITE LAYERS-OPEN  -  address in, layered DXF and site diagrams out, from open data.
ELVTR AI for Architects (APAC 2026) - Assignment 3 student tool.

    python site_layers-open.py "388 George St, Sydney NSW" --radius 300

Nothing to buy, no keys. Needs:  pip install ezdxf shapely matplotlib numpy requests
Every fact the tool pulls or infers lands in a verification register with an
empty verdict column. The tool does not know if it is right. You check.

Data: (c) OpenStreetMap contributors (ODbL) - NSW Spatial Services and NSW
Department of Planning (CC BY 4.0) - Open-Meteo (CC BY 4.0, ERA5 reanalysis)
- Terrain Tiles on AWS (Mapzen/Tilezen; SRTM and others).
"""
from __future__ import annotations

VERSION = "1.0"
TOOL = "Site Layers-open"

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

import numpy as np
import requests
import ezdxf
from shapely.geometry import (Point, Polygon, LineString, MultiPolygon, MultiLineString,
                              GeometryCollection, mapping, shape)
from shapely.ops import unary_union, polygonize, transform as shp_transform
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import PathPatch, Circle as MplCircle, Polygon as MplPolygon
from matplotlib.path import Path as MplPath

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

USER_AGENT = f"SiteLayers-open/{VERSION} (ELVTR AI for Architects student tool; open data only)"
TODAY = dt.date.today().isoformat()


# --------------------------------------------------------------------------------------
# small helpers
# --------------------------------------------------------------------------------------
def slugify(s: str, n: int = 40) -> str:
    s = re.sub(r"[^A-Za-z0-9]+", "-", s).strip("-").lower()
    return (s[:n].rstrip("-")) or "site"


def epoch_ms_to_date(v):
    try:
        if v is None:
            return ""
        return dt.datetime.utcfromtimestamp(float(v) / 1000.0).date().isoformat()
    except Exception:
        return ""


def fnum(v):
    """Parse a number out of an OSM-style tag ('12', '12 m', '12.5m', '3;4')."""
    if v is None:
        return None
    m = re.search(r"-?\d+(?:\.\d+)?", str(v).replace(",", "."))
    return float(m.group(0)) if m else None


def fmt_m(v, nd=1):
    return "" if v is None else f"{v:.{nd}f}"


class Log:
    def __init__(self, path: Path | None = None, quiet=False):
        self.lines = []
        self.path = path
        self.quiet = quiet

    def __call__(self, *a):
        s = " ".join(str(x) for x in a)
        self.lines.append(f"{dt.datetime.now().strftime('%H:%M:%S')}  {s}")
        if not self.quiet:
            print(s, flush=True)

    def save(self):
        if self.path:
            self.path.write_text("\n".join(self.lines) + "\n", encoding="utf-8")


# --------------------------------------------------------------------------------------
# network: one polite attempt, cache everything, never a retry storm
# --------------------------------------------------------------------------------------
class NotCached(Exception):
    pass


class Net:
    """Every response is cached on disk. --offline replays the cache only (the class demo).
    Policy: one attempt per URL; on 429/502/503/504 wait once (max 20 s) and try once more.
    Then give up and let the layer fail on its own - the run carries on."""

    RETRY_CODES = {502, 503, 504}

    def __init__(self, cache_dir: Path, offline=False, log=print):
        self.cache_dir = Path(cache_dir)
        self.cache_dir.mkdir(parents=True, exist_ok=True)
        self.offline = offline
        self.log = log
        self.s = requests.Session()
        self.s.headers.update({"User-Agent": USER_AGENT})
        self._last_nominatim = 0.0

    def _key(self, method, url, params, data):
        raw = json.dumps([method, url, sorted((params or {}).items()), data or ""], sort_keys=True)
        return hashlib.sha1(raw.encode("utf-8")).hexdigest()[:20]

    def fetch(self, url, params=None, data=None, timeout=60, binary=False, label=""):
        """Returns (payload, meta). payload is bytes if binary else text. meta has 'retrieved', 'url'."""
        method = "POST" if data is not None else "GET"
        key = self._key(method, url, params, data)
        body_p = self.cache_dir / f"{key}.bin"
        meta_p = self.cache_dir / f"{key}.json"
        if body_p.exists() and meta_p.exists():
            meta = json.loads(meta_p.read_text(encoding="utf-8"))
            b = body_p.read_bytes()
            return (b if binary else b.decode("utf-8", "replace")), meta
        if self.offline:
            raise NotCached(f"not in cache (offline mode): {label or url}")
        if "nominatim.openstreetmap.org" in url:  # usage policy: max 1 request per second
            wait = 1.1 - (time.time() - self._last_nominatim)
            if wait > 0:
                time.sleep(wait)
            self._last_nominatim = time.time()
        attempt = 0
        while True:
            attempt += 1
            try:
                status, final_url, content, headers = self._request(method, url, params, data, timeout)
            except Exception as e:
                raise RuntimeError(f"{label or url}: network error ({type(e).__name__})") from None
            if status in self.RETRY_CODES and attempt == 1:
                ra = min(fnum(headers.get("retry-after")) or 8, 20)
                self.log(f"   {label}: server busy ({status}); waiting {ra:.0f} s, trying once more")
                time.sleep(ra)
                continue
            if status != 200:
                raise RuntimeError(f"{label or url}: HTTP {status}")
            break
        meta = {"url": final_url, "retrieved": dt.datetime.now().astimezone().isoformat(timespec="seconds"),
                "status": status, "headers": headers}
        body_p.write_bytes(content)
        meta_p.write_text(json.dumps(meta, indent=1), encoding="utf-8")
        return (content if binary else content.decode("utf-8", "replace")), meta

    KEEP_HEADERS = ("content-type", "x-amz-meta-x-imagery-sources", "last-modified", "retry-after")

    def _request(self, method, url, params, data, timeout):
        """(status, final_url, bytes, headers). CPython uses requests; a browser (Pyodide) uses sync XHR."""
        if sys.platform == "emscripten":
            from urllib.parse import urlencode
            from pyodide.code import run_js
            full = url + (("&" if "?" in url else "?") + urlencode(params) if params else "")
            fn = run_js("""(m,u,b,ms)=>{const x=new XMLHttpRequest();x.open(m,u,false);
              try{x.timeout=ms;}catch(e){}
              x.overrideMimeType('text/plain; charset=x-user-defined');
              if(b!==null){x.setRequestHeader('Content-Type','application/x-www-form-urlencoded');}
              x.send(b);const t=x.responseText;const a=new Uint8Array(t.length);
              for(let i=0;i<t.length;i++){a[i]=t.charCodeAt(i)&255;}
              return [x.status,x.responseURL||u,a,x.getAllResponseHeaders()];}""")
            st, fu, arr, hdr = fn(method, full, urlencode(data) if data is not None else None, int(timeout * 1000)).to_py()
            heads = {}
            for line in str(hdr).splitlines():
                if ":" in line:
                    k, v = line.split(":", 1)
                    if k.strip().lower() in self.KEEP_HEADERS:
                        heads[k.strip().lower()] = v.strip()
            return int(st), str(fu), bytes(arr), heads
        if method == "POST":
            r = self.s.post(url, data=data, timeout=timeout)
        else:
            r = self.s.get(url, params=params, timeout=timeout)
        heads = {k.lower(): v for k, v in r.headers.items() if k.lower() in self.KEEP_HEADERS}
        return r.status_code, r.url, r.content, heads

    def json(self, url, params=None, data=None, timeout=60, label=""):
        txt, meta = self.fetch(url, params=params, data=data, timeout=timeout, label=label)
        try:
            obj = json.loads(txt)
        except Exception:
            raise RuntimeError(f"{label or url}: reply was not JSON ({txt[:80]!r})") from None
        if isinstance(obj, dict) and "error" in obj and not obj.get("features"):
            raise RuntimeError(f"{label or url}: service error {obj['error']}")
        return obj, meta


# --------------------------------------------------------------------------------------
# verification register
# --------------------------------------------------------------------------------------
REGISTER_COLUMNS = ["id", "layer", "assertion", "source", "source_url", "retrieved",
                    "confidence", "how_to_check", "verdict", "note"]


class Register:
    """One row per pulled or inferred fact. 'verdict' and 'note' are left blank for the student."""

    def __init__(self):
        self.rows = []

    def add(self, layer, assertion, source, url="", retrieved="", confidence="medium", how=""):
        self.rows.append({
            "id": f"R{len(self.rows) + 1:03d}", "layer": layer, "assertion": assertion,
            "source": source, "source_url": url, "retrieved": (retrieved or "")[:10],
            "confidence": confidence, "how_to_check": how, "verdict": "", "note": ""})

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


# --------------------------------------------------------------------------------------
# local metric frame centred on the pin (x east, y north, metres, true north up)
# Accurate to ~5 mm at 300 m and ~5 cm at 1 km - far inside open-data accuracy.
# --------------------------------------------------------------------------------------
class LocalFrame:
    def __init__(self, lat0, lon0):
        self.lat0, self.lon0 = lat0, lon0
        a = 6378137.0
        f = 1 / 298.257223563
        e2 = f * (2 - f)
        phi = math.radians(lat0)
        s = math.sin(phi)
        N = a / math.sqrt(1 - e2 * s * s)
        M = a * (1 - e2) / (1 - e2 * s * s) ** 1.5
        self.kx = math.radians(1) * N * math.cos(phi)
        self.ky = math.radians(1) * M

    def fwd(self, lon, lat):
        return (np.asarray(lon) - self.lon0) * self.kx, (np.asarray(lat) - self.lat0) * self.ky

    def inv(self, x, y):
        return np.asarray(x) / self.kx + self.lon0, np.asarray(y) / self.ky + self.lat0

    def to_local(self, geom):
        return shp_transform(lambda x, y, z=None: self.fwd(x, y), geom)

    def to_lonlat(self, geom):
        return shp_transform(lambda x, y, z=None: self.inv(x, y), geom)

    def bbox_lonlat(self, r):
        lon1, lat1 = self.inv(-r, -r)
        lon2, lat2 = self.inv(r, r)
        return float(lon1), float(lat1), float(lon2), float(lat2)


# --------------------------------------------------------------------------------------
# layer container
# --------------------------------------------------------------------------------------
class Layer:
    """A named set of features in the local frame. status: ok | empty | unavailable | failed | skipped"""

    def __init__(self, key, dxf, title, family="CTX"):
        self.key, self.dxf, self.title, self.family = key, dxf, title, family
        self.features = []  # list of (shapely geom in local metres, props dict)
        self.status = "skipped"
        self.note = ""
        self.source = ""
        self.url = ""
        self.retrieved = ""

    def add(self, geom, props=None):
        if geom is None or geom.is_empty:
            return
        self.features.append((geom, props or {}))

    def __len__(self):
        return len(self.features)


def clip_to(geom, circle):
    try:
        if not geom.is_valid:
            geom = geom.buffer(0) if geom.geom_type in ("Polygon", "MultiPolygon") else geom
        g = geom.intersection(circle)
        return None if g.is_empty else g
    except Exception:
        return None


def iter_parts(geom):
    """Yield simple parts (Polygon / LineString / Point) of any geometry."""
    if geom is None or geom.is_empty:
        return
    if geom.geom_type in ("Polygon", "LineString", "Point", "LinearRing"):
        yield geom
    elif hasattr(geom, "geoms"):
        for g in geom.geoms:
            yield from iter_parts(g)


# --------------------------------------------------------------------------------------
# geocoding and jurisdiction
# --------------------------------------------------------------------------------------
NOMINATIM = "https://nominatim.openstreetmap.org"
PHOTON = "https://photon.komoot.io/api/"

AU_STATES = {"new south wales": "AU-NSW", "victoria": "AU-VIC", "queensland": "AU-QLD",
             "western australia": "AU-WA", "south australia": "AU-SA", "tasmania": "AU-TAS",
             "australian capital territory": "AU-ACT", "northern territory": "AU-NT"}


def jurisdiction_from(country_code, state):
    cc = (country_code or "").lower()
    st = (state or "").lower()
    if cc == "au":
        return AU_STATES.get(st, "AU-" + (st[:3].upper() or "?"))
    return {"nz": "NZ", "sg": "SG", "hk": "HK", "my": "MY", "id": "ID", "ph": "PH", "jp": "JP",
            "kr": "KR", "cn": "CN", "in": "IN", "th": "TH", "vn": "VN"}.get(cc, cc.upper() or "UNKNOWN")


NSW_ADDRESS = "https://maps.six.nsw.gov.au/services/public/Address_Location"
ROAD_TYPES = (r"(street|st|avenue|ave|av|road|rd|lane|ln|place|pl|parade|pde|drive|dr|crescent|cres|cr|court|ct|"
              r"terrace|tce|highway|hwy|boulevard|blvd|way|close|cl|circuit|cct|square|sq|grove|gr|esplanade|esp)")


def geocode_nsw(net: Net, query: str, reg: Register):
    """NSW Spatial Services address points: lands on the property, not the street."""
    m = re.match(r"^\s*(?:unit\s*\S+\s*/\s*|\d+[A-Za-z]?\s*/\s*)?(\d+[A-Za-z]?)(?:\s*-\s*\d+[A-Za-z]?)?\s+(.+?),\s*(.+?)"
                 r"(?:,?\s+(?:NSW|New South Wales))(?:\s+\d{4})?\s*$", query, re.I)
    if not m:
        return None
    num, road, suburb = m.groups()
    road_name = re.sub(rf"\s+{ROAD_TYPES}\.?$", "", road.strip(), flags=re.I).strip()
    obj, meta = net.json(NSW_ADDRESS, params={"projection": "EPSG:4326", "houseNumber": num, "roadName": road_name,
                                              "suburb": suburb.strip()}, label="geocode (NSW address points)")
    adds = (obj.get("addressResult") or {}).get("addresses") or []
    if not adds:
        return None
    a = adds[0]
    ap = a.get("addressPoint") or {}
    lon, lat = ap.get("centreX", a.get("longitude")), ap.get("centreY", a.get("latitude"))
    if lon is None or lat is None:
        return None
    label = a.get("formattedAddress") or a.get("address") or query
    reg.add("PIN", f"Address '{query}' found as NSW address point {float(lat):.6f}, {float(lon):.6f} ({label})",
            "NSW Spatial Services address location service", meta["url"], meta["retrieved"], "medium",
            "Check the pin sits on your lot. Large sites have one address point - it may not be where you expect.")
    return dict(lat=float(lat), lon=float(lon), label=str(label), jurisdiction="AU-NSW", source="NSW address points")


def geocode(net: Net, query: str, reg: Register, log):
    """Address or 'lat, lon'. Returns dict(lat, lon, label, jurisdiction)."""
    m = re.match(r"^\s*(-?\d+(?:\.\d+)?)\s*[, ]\s*(-?\d+(?:\.\d+)?)\s*$", query)
    if m:
        lat, lon = float(m.group(1)), float(m.group(2))
        if abs(lat) > 90 and abs(lon) <= 90:
            lat, lon = lon, lat
        label, cc, state = f"{lat:.6f}, {lon:.6f}", "", ""
        try:
            obj, meta = net.json(f"{NOMINATIM}/reverse", params={"lat": lat, "lon": lon, "format": "jsonv2",
                                                                  "zoom": 18, "addressdetails": 1},
                                 label="reverse geocode")
            label = obj.get("display_name", label)
            ad = obj.get("address", {})
            cc, state = ad.get("country_code", ""), ad.get("state", "")
        except Exception as e:
            log(f"   reverse geocode failed ({e}); jurisdiction must be given with --jurisdiction")
        reg.add("PIN", f"Pin placed from coordinates you typed: {lat:.6f}, {lon:.6f} (nearest address: {label})",
                "your input", "", TODAY, "high",
                "Open the coordinates in any web map and confirm the pin is on your site.")
        return dict(lat=lat, lon=lon, label=label, jurisdiction=jurisdiction_from(cc, state), source="input")

    tried = []
    if re.search(r"\bNSW\b|New South Wales", query, re.I):
        try:
            res = geocode_nsw(net, query, reg)
            if res:
                return res
            tried.append("NSW address points: no match")
        except Exception as e:
            tried.append(f"NSW address points: {e}")
    try:
        obj, meta = net.json(f"{NOMINATIM}/search", params={"q": query, "format": "jsonv2", "limit": 1,
                                                             "addressdetails": 1}, label="geocode (Nominatim)")
        if obj:
            h = obj[0]
            ad = h.get("address", {})
            res = dict(lat=float(h["lat"]), lon=float(h["lon"]), label=h.get("display_name", query),
                       jurisdiction=jurisdiction_from(ad.get("country_code"), ad.get("state")), source="Nominatim")
            precise = bool(ad.get("house_number")) or h.get("category") in ("building",) or h.get("type") in ("house",)
            street_words = [w for w in re.findall(r"[a-z]{3,}", query.split(",")[0].lower())
                            if not re.fullmatch(ROAD_TYPES, w)]
            mismatch = street_words and not all(w in res["label"].lower() for w in street_words)
            if mismatch:
                log(f"   WARNING: the address found ({res['label']}) does not match the street you typed. "
                    "Check it, or rerun with 'lat, lon'.")
            reg.add("PIN", (f"Address '{query}' geocoded to {res['lat']:.6f}, {res['lon']:.6f} - {res['label']}"
                            + (" - THE STREET DOES NOT MATCH WHAT YOU TYPED" if mismatch else "")),
                    "OpenStreetMap Nominatim", meta["url"], meta["retrieved"],
                    "low" if (mismatch or not precise) else "medium",
                    "Check the pin sits on your lot, not the street or a neighbour. If not, rerun with 'lat, lon'.")
            return res
        tried.append("Nominatim: no match")
    except Exception as e:
        tried.append(f"Nominatim: {e}")
    try:
        obj, meta = net.json(PHOTON, params={"q": query, "limit": 1}, label="geocode (Photon)")
        feats = obj.get("features", [])
        if feats:
            f = feats[0]
            lon, lat = f["geometry"]["coordinates"]
            p = f.get("properties", {})
            label = ", ".join(str(p[k]) for k in ("housenumber", "street", "city", "state", "country") if p.get(k))
            res = dict(lat=lat, lon=lon, label=label or query,
                       jurisdiction=jurisdiction_from(p.get("countrycode"), p.get("state")), source="Photon")
            reg.add("PIN", f"Address '{query}' geocoded to {lat:.6f}, {lon:.6f} - {res['label']}",
                    "Photon (komoot, OSM-based)", meta["url"], meta["retrieved"], "low",
                    "Check the pin sits on your lot. If not, rerun with 'lat, lon'.")
            return res
        tried.append("Photon: no match")
    except Exception as e:
        tried.append(f"Photon: {e}")
    raise SystemExit("Could not find that address (" + "; ".join(tried) +
                     "). Try a simpler form, or give coordinates as 'lat, lon'.")


# --------------------------------------------------------------------------------------
# OPENSTREETMAP via Overpass  (buildings, roads, paths, rail, water, trees, land use)
# --------------------------------------------------------------------------------------
# Checked 23 Sep 2026: overpass-api.de was intermittently "too busy" (HTTP 504); its sister server z. answered in
# 3 s; maps.mail.ru answered in 18 s; kumi.systems and private.coffee timed out. Each is tried ONCE, in this order.
OVERPASS = ["https://overpass-api.de/api/interpreter",
            "https://z.overpass-api.de/api/interpreter",
            "https://maps.mail.ru/osm/tools/overpass/api/interpreter",
            "https://overpass.kumi.systems/api/interpreter"]
OVERPASS_TIMEOUT = [90, 90, 60, 45]

RAIL_KINDS = {"rail", "light_rail", "subway", "tram", "narrow_gauge", "monorail", "funicular", "preserved"}
PATH_KINDS = {"footway", "path", "pedestrian", "cycleway", "steps", "corridor", "bridleway", "track", "elevator"}
SKIP_HIGHWAY = {"proposed", "construction", "platform", "bus_stop", "street_lamp", "traffic_signals",
                "crossing", "give_way", "stop", "turning_circle", "milestone", "elevator"}
GREEN_LANDUSE = {"grass", "recreation_ground", "forest", "meadow", "village_green", "cemetery", "allotments",
                 "orchard", "vineyard", "greenfield"}
WOOD = {"wood", "scrub", "heath"}


def overpass_query(lat, lon, r):
    a = f"(around:{r},{lat:.7f},{lon:.7f})"
    return (f"[out:json][timeout:90];("
            f"way[building]{a};relation[building][type=multipolygon]{a};"
            f"way[highway]{a};way[railway]{a};"
            f"way[natural=water]{a};relation[natural=water]{a};way[waterway]{a};way[natural=coastline]{a};"
            f"node[natural=tree]{a};way[natural~\"^(wood|scrub|heath|tree_row)$\"]{a};"
            f"way[landuse]{a};relation[landuse][type=multipolygon]{a};"
            f"way[leisure~\"^(park|garden|pitch|playground|nature_reserve|common)$\"]{a};"
            f");out geom;")


def overpass_wait(net: Net, ep, log):
    """overpass-api.de rations queries per address (4 slots). When it refuses (406/429/504), ask its
    status page how long until a slot frees, wait that once (max 45 s). Never a loop."""
    try:
        st, _, body, _ = net._request("GET", ep.replace("/interpreter", "/status"), None, None, 15)
        txt = body.decode("utf-8", "replace")
        if re.search(r"(\d+) slots? available now", txt) and not txt.count("0 slots available now"):
            return 3
        secs = [int(x) for x in re.findall(r"in (\d+) seconds", txt)]
        return min(min(secs) + 2, 45) if secs else 15
    except Exception:
        return 15


def fetch_osm(net: Net, lat, lon, r, log):
    q = overpass_query(lat, lon, r)
    errors = []
    for i, ep in enumerate(OVERPASS):
        host = ep.split('/')[2]
        try:
            try:
                obj, meta = net.json(ep, data={"data": q}, timeout=OVERPASS_TIMEOUT[i], label=f"Overpass {host}")
            except RuntimeError as e:
                if i == 0 and re.search(r"HTTP (406|429|504)", str(e)):
                    w = overpass_wait(net, ep, log)
                    log(f"   Overpass {host} busy for this address; waiting {w} s once")
                    time.sleep(w)
                    obj, meta = net.json(ep, data={"data": q}, timeout=OVERPASS_TIMEOUT[i], label=f"Overpass {host}")
                else:
                    raise
            if "elements" not in obj:
                raise RuntimeError("no elements in reply")
            remark = obj.get("remark", "")
            if remark and "runtime error" in remark.lower():
                raise RuntimeError(remark[:120])
            return obj, meta, ep
        except NotCached:
            raise
        except Exception as e:
            errors.append(str(e))
            log(f"   {e} - trying next Overpass server")
    raise RuntimeError("all Overpass servers failed: " + " | ".join(errors))


def _way_coords(el):
    return [(p["lon"], p["lat"]) for p in el.get("geometry", []) if p]


def _relation_polygon(el):
    outers, inners = [], []
    for m in el.get("members", []):
        if m.get("type") != "way" or not m.get("geometry"):
            continue
        coords = [(p["lon"], p["lat"]) for p in m["geometry"] if p]
        if len(coords) < 2:
            continue
        (inners if m.get("role") == "inner" else outers).append(LineString(coords))
    if not outers:
        return None
    polys_out = list(polygonize(unary_union(outers)))
    polys_in = list(polygonize(unary_union(inners))) if inners else []
    if not polys_out:
        return None
    g = unary_union(polys_out)
    if polys_in:
        g = g.difference(unary_union(polys_in))
    return g


def _el_geom(el):
    if el["type"] == "node":
        return Point(el["lon"], el["lat"])
    if el["type"] == "way":
        c = _way_coords(el)
        if len(c) < 2:
            return None
        return LineString(c), (len(c) >= 4 and c[0] == c[-1])
    if el["type"] == "relation":
        return _relation_polygon(el)
    return None


def height_of(tags):
    """(height_m, basis). Never invents: no tag, no height."""
    h = fnum(tags.get("height"))
    if h is not None and 0 < h < 1000:
        return h, "height tag"
    lv = fnum(tags.get("building:levels"))
    if lv is not None and 0 < lv < 200:
        rl = fnum(tags.get("roof:levels")) or 0
        n = int(lv + rl)
        return round((lv + rl) * 3.2, 1), f"{n} storey{'s' if n != 1 else ''} x 3.2 m (estimate)"
    return None, "no height in OSM"


def build_osm_layers(osm, frame: LocalFrame, r, layers: dict):
    circle = Point(0, 0).buffer(r, 128)
    for el in osm.get("elements", []):
        tags = el.get("tags", {}) or {}
        g = _el_geom(el)
        if g is None:
            continue
        closed = False
        if isinstance(g, tuple):
            g, closed = g
        try:
            g = frame.to_local(g)
        except Exception:
            continue
        oid = f"{el['type'][0]}{el['id']}"
        name = tags.get("name", "")

        if "building" in tags and el["type"] in ("way", "relation"):
            poly = Polygon(g.coords) if (el["type"] == "way" and closed) else (g if el["type"] == "relation" else None)
            if poly is None or poly.is_empty:
                continue
            if not poly.is_valid:
                poly = poly.buffer(0)
            if poly.is_empty or not poly.centroid.within(circle):
                continue  # keep whole buildings whose centre is inside the radius; never cut a building
            h, basis = height_of(tags)
            layers["buildings"].add(poly, dict(osm_id=oid, name=name, building=tags.get("building", ""),
                                               height_m=h, height_basis=basis,
                                               levels=tags.get("building:levels", ""),
                                               height_tag=tags.get("height", "")))
            continue

        if "highway" in tags and el["type"] == "way":
            kind = tags["highway"]
            if kind in SKIP_HIGHWAY:
                continue
            c = clip_to(g, circle)
            if c is None:
                continue
            props = dict(osm_id=oid, name=name, kind=kind, tunnel=tags.get("tunnel", ""),
                         bridge=tags.get("bridge", ""), lanes=tags.get("lanes", ""))
            (layers["paths"] if kind in PATH_KINDS else layers["roads"]).add(c, props)
            continue

        if "railway" in tags and el["type"] == "way":
            kind = tags["railway"]
            if kind not in RAIL_KINDS:
                continue
            c = clip_to(g, circle)
            if c is not None:
                layers["rail"].add(c, dict(osm_id=oid, name=name, kind=kind, tunnel=tags.get("tunnel", ""),
                                           bridge=tags.get("bridge", "")))
            continue

        nat = tags.get("natural", "")
        if nat == "water" or "waterway" in tags or nat == "coastline":
            if nat == "water" and (closed or el["type"] == "relation"):
                gg = Polygon(g.coords) if el["type"] == "way" else g
                gg = gg if gg.is_valid else gg.buffer(0)
            else:
                gg = g
            c = clip_to(gg, circle)
            if c is not None:
                layers["water"].add(c, dict(osm_id=oid, name=name,
                                            kind=tags.get("waterway") or tags.get("water") or nat))
            continue

        if nat == "tree" and el["type"] == "node":
            if g.within(circle):
                crown = fnum(tags.get("diameter_crown"))
                layers["trees"].add(g, dict(osm_id=oid, kind="tree", species=tags.get("species", "") or tags.get("genus", ""),
                                            crown_m=crown or ""))
            continue
        if nat == "tree_row" and el["type"] == "way":
            c = clip_to(g, circle)
            if c is not None:
                layers["trees"].add(c, dict(osm_id=oid, kind="tree_row"))
            continue
        if nat in WOOD and (closed or el["type"] == "relation"):
            gg = Polygon(g.coords) if el["type"] == "way" else g
            c = clip_to(gg if gg.is_valid else gg.buffer(0), circle)
            if c is not None:
                layers["trees"].add(c, dict(osm_id=oid, kind=nat))
            continue

        lu = tags.get("landuse") or tags.get("leisure")
        if lu and (closed or el["type"] == "relation"):
            gg = Polygon(g.coords) if el["type"] == "way" else g
            if gg is None:
                continue
            c = clip_to(gg if gg.is_valid else gg.buffer(0), circle)
            if c is not None:
                green = lu in GREEN_LANDUSE or "leisure" in tags
                layers["landuse"].add(c, dict(osm_id=oid, name=name, kind=lu, green="yes" if green else ""))


# --------------------------------------------------------------------------------------
# ARCGIS REST helpers (NSW services)
# --------------------------------------------------------------------------------------
def esri_to_shapely(geom):
    """Esri JSON geometry (rings / paths / x,y) in lon/lat -> shapely."""
    if not geom:
        return None
    if "rings" in geom:
        polys = []
        for ring in geom["rings"]:
            if len(ring) < 4:
                continue
            p = Polygon(ring)
            if not p.is_valid:
                p = p.buffer(0)
            polys.append(p)
        if not polys:
            return None
        # outer rings contain holes: xor-assemble
        out = polys[0]
        for p in polys[1:]:
            out = out.symmetric_difference(p)
        return out
    if "paths" in geom:
        ls = [LineString(p) for p in geom["paths"] if len(p) >= 2]
        return ls[0] if len(ls) == 1 else MultiLineString(ls) if ls else None
    if "x" in geom:
        return Point(geom["x"], geom["y"])
    return None


def arcgis_query(net: Net, layer_url, frame: LocalFrame, r=None, point=None, out_fields="*",
                 fmt="geojson", label="", max_pages=10, extra=None):
    """Envelope (radius r around the pin) or point query. Returns (list of (geom lonlat, props), meta)."""
    params = {"inSR": 4326, "outSR": 4326, "spatialRel": "esriSpatialRelIntersects",
              "outFields": out_fields, "returnGeometry": "true", "f": fmt}
    if point is not None:
        params.update(geometry=f"{point[0]},{point[1]}", geometryType="esriGeometryPoint")
    else:
        x1, y1, x2, y2 = frame.bbox_lonlat(r)
        params.update(geometry=f"{x1:.7f},{y1:.7f},{x2:.7f},{y2:.7f}", geometryType="esriGeometryEnvelope")
    if extra:
        params.update(extra)
    out, meta0 = [], None
    offset = 0
    for page in range(max_pages):
        p = dict(params)
        if offset:
            p["resultOffset"] = offset
        obj, meta = net.json(layer_url.rstrip("/") + "/query", params=p, timeout=90, label=label)
        meta0 = meta0 or meta
        feats = obj.get("features", [])
        for f in feats:
            if fmt == "geojson":
                g = shape(f["geometry"]) if f.get("geometry") else None
                props = f.get("properties", {}) or {}
            else:
                g = esri_to_shapely(f.get("geometry"))
                props = f.get("attributes", {}) or {}
            if g is not None and not g.is_empty:
                out.append((g, props))
        more = obj.get("exceededTransferLimit") or (obj.get("properties", {}) or {}).get("exceededTransferLimit")
        if not more or not feats:
            break
        offset += len(feats)
    return out, meta0


def layer_schema(net: Net, layer_url):
    """Coded-value domains and subtype names. GeoJSON returns raw codes (a height map code 46 instead
    of 'J2'), so they are decoded here (found by the overlay build, 23 Sep 2026)."""
    try:
        obj, _ = net.json(layer_url, params={"f": "json"}, timeout=60, label="layer schema")
    except NotCached:
        raise
    except Exception:
        return {}
    dom = {}
    for f in obj.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 obj.get("subtypeField") and obj.get("types"):
        dom[obj["subtypeField"]] = {str(t.get("id")): t.get("name") for t in obj["types"]}
    return dom


def decode(props, dom):
    for k, v in list(props.items()):
        if v is not None and k in dom and str(v) in dom[k]:
            props[k] = dom[k][str(v)]
    return props


# --------------------------------------------------------------------------------------
# JURISDICTION ADAPTERS
# One row per jurisdiction. NSW is filled. The others are explicit empty slots: they
# produce a register row saying "not available for this jurisdiction" and the run goes on.
# To add one: copy the NSW block, point each layer at that state's open map service.
# --------------------------------------------------------------------------------------
NSW_PORTAL = "https://portal.spatial.nsw.gov.au/server/rest/services"
NSW_EP = "https://mapprod3.environment.nsw.gov.au/arcgis/rest/services/ePlanning"


def _hob_label(p):
    v = p.get("MAX_B_H") if p.get("MAX_B_H") is not None else p.get("MAX_B_H_M")
    if v is None and p.get("MAX_B_H_RL") is not None:
        return f"RL {p['MAX_B_H_RL']}"
    return "" if v is None else f"{v:g} {p.get('UNITS') or 'm'}"


def _fsr_label(p):
    v = p.get("FSR")
    return (f"{v:g}:1" if isinstance(v, (int, float)) else str(p.get("LABEL") or "")) if v is not None else str(p.get("LABEL") or "")


STAT_NSW = [
    # key, dxf layer, title, url, label fn, description fn
    dict(key="zone", dxf="AI_ZONE", title="Land zoning",
         url=f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/19",
         label=lambda p: str(p.get("SYM_CODE") or p.get("LABEL") or ""),
         desc=lambda p: f"{p.get('SYM_CODE') or p.get('LABEL')} {p.get('LAY_CLASS') or ''}".strip()),
    dict(key="hob", dxf="AI_HOB", title="Height of buildings (LEP)",
         url=f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/14",
         label=_hob_label, desc=lambda p: f"maximum height {_hob_label(p)}"),
    dict(key="fsr", dxf="AI_FSR", title="Floor space ratio (LEP)",
         url=f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/11",
         label=_fsr_label, desc=lambda p: f"maximum FSR {_fsr_label(p)}"),
    dict(key="addctrl", dxf="AI_ADDCTRL", title="Additional controls (height / FSR areas)",
         url=[f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/10", f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/13"],
         label=lambda p: str(p.get("LABEL") or p.get("SYM_CODE") or ""),
         desc=lambda p: f"additional controls area {p.get('LABEL') or p.get('SYM_CODE') or ''} ({p.get('LEGIS_REF_CLAUSE') or 'see LEP'})"),
    dict(key="heritage", dxf="AI_HERITAGE", title="Heritage (LEP items and conservation areas)",
         url=f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/16",
         label=lambda p: str(p.get("H_ID") or ""),
         desc=lambda p: f"{p.get('LAY_CLASS') or 'heritage'} {p.get('H_ID') or ''} '{p.get('H_NAME') or ''}' ({p.get('SIG') or '?'} significance)"),
    dict(key="shr", dxf="AI_HERITAGE_SHR", title="State Heritage Register curtilage",
         url=f"{NSW_EP}/Planning_Portal_Principal_Planning/MapServer/221",
         label=lambda p: str(p.get("LISTINGNO") or ""),
         desc=lambda p: f"SHR {p.get('LISTINGNO') or ''} '{p.get('ITEMNAME') or ''}'"),
    dict(key="flood", dxf="AI_FLOOD", title="Flood planning (LEP/SEPP map)",
         url=f"{NSW_EP}/Planning_Portal_Hazard/MapServer/230",
         label=lambda p: "FPA", desc=lambda p: f"{p.get('LAY_CLASS') or 'flood planning area'}"),
    dict(key="bushfire", dxf="AI_BUSHFIRE", title="Bush fire prone land",
         url=f"{NSW_EP}/Planning_Portal_Hazard/MapServer/229",
         label=lambda p: "BFPL", desc=lambda p: f"{p.get('d_Category') or 'bush fire prone land'}"),
    dict(key="ass", dxf="AI_ASS", title="Acid sulfate soils",
         url=f"{NSW_EP}/Planning_Portal_Protection/MapServer/234",
         label=lambda p: str(p.get("LABEL") or p.get("LAY_CLASS") or ""),
         desc=lambda p: f"acid sulfate soils {p.get('LAY_CLASS') or p.get('LABEL') or ''}"),
]

ADAPTERS = {
    "AU-NSW": dict(name="New South Wales",
                   cadastre=f"{NSW_PORTAL}/NSW_Land_Parcel_Property_Theme/FeatureServer/8",
                   contours=f"{NSW_PORTAL}/NSW_Elevation_and_Depth_Theme/FeatureServer/2",
                   statutory=STAT_NSW,
                   planning_viewer="https://www.planningportal.nsw.gov.au/spatialviewer/",
                   note="NSW Planning Portal open map services (no key)."),
    # ---- empty slots: fill when a keyless service is verified live ----
    "AU-VIC": dict(name="Victoria", cadastre=None, contours=None, statutory=[], planning_viewer="https://mapshare.vic.gov.au/vicplan/",
                   note="Slot empty. Candidate: Vicmap / VicPlan open WFS. Not verified."),
    "AU-QLD": dict(name="Queensland", cadastre=None, contours=None, statutory=[], planning_viewer="https://planning.statedevelopment.qld.gov.au/maps",
                   note="Slot empty. Candidate: QSpatial / council planning scheme maps. Not verified."),
    "NZ": dict(name="New Zealand", cadastre=None, contours=None, statutory=[], planning_viewer="",
               note="Slot empty. LINZ cadastre needs a free key; district plan zones are per council. Not verified."),
    "SG": dict(name="Singapore", cadastre=None, contours=None, statutory=[], planning_viewer="https://www.ura.gov.sg/maps/",
               note="Slot empty. URA Master Plan via data.gov.sg. Not verified."),
    "HK": dict(name="Hong Kong", cadastre=None, contours=None, statutory=[], planning_viewer="https://www.ozp.tpb.gov.hk/",
               note="Slot empty. Town Planning Board OZP / CSDI portal. Not verified."),
}


def adapter_for(j):
    return ADAPTERS.get(j) or dict(name=j, cadastre=None, contours=None, statutory=[], planning_viewer="",
                                    note="No adapter for this jurisdiction.")


# --------------------------------------------------------------------------------------
# TERRAIN: NSW contours (Spatial Services, 2 m in cities) or SRTM terrain tiles elsewhere
# --------------------------------------------------------------------------------------
TERRARIUM = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium/{z}/{x}/{y}.png"


def _tile_xy(lon, lat, z):
    n = 2 ** z
    x = (lon + 180.0) / 360.0 * n
    y = (1.0 - math.asinh(math.tan(math.radians(lat))) / math.pi) / 2.0 * n
    return x, y


def _tile_lonlat(x, y, z):
    n = 2 ** z
    lon = x / n * 360.0 - 180.0
    lat = math.degrees(math.atan(math.sinh(math.pi * (1 - 2 * y / n))))
    return lon, lat


def terrain_grid(net: Net, frame: LocalFrame, r, log):
    """Decode Terrarium PNG tiles to a height grid in local metres. Returns (X, Y, Z, sources, meta)."""
    z = 15 if r <= 600 else 14
    lon1, lat1, lon2, lat2 = frame.bbox_lonlat(r * 1.05)
    x1, y1 = _tile_xy(lon1, lat2, z)
    x2, y2 = _tile_xy(lon2, lat1, z)
    tx = range(int(x1), int(x2) + 1)
    ty = range(int(y1), int(y2) + 1)
    if len(tx) * len(ty) > 16:
        raise RuntimeError("terrain area too large")
    rows, sources, meta0 = [], set(), None
    for yy in ty:
        row = []
        for xx in tx:
            b, meta = net.fetch(TERRARIUM.format(z=z, x=xx, y=yy), binary=True, timeout=60, label=f"terrain tile {z}/{xx}/{yy}")
            meta0 = meta0 or meta
            src = meta.get("headers", {}).get("x-amz-meta-x-imagery-sources") or meta.get("headers", {}).get("X-Amz-Meta-X-Imagery-Sources")
            if src:
                sources.update(s.strip().split("/")[0] for s in src.split(","))
            img = plt.imread(io.BytesIO(b), format="png")[:, :, :3] * 255.0
            row.append(img[:, :, 0] * 256 + img[:, :, 1] + img[:, :, 2] / 256.0 - 32768)
        rows.append(np.hstack(row))
    Zg = np.vstack(rows)
    ny, nx = Zg.shape
    px = np.arange(nx) + 0.5
    py = np.arange(ny) + 0.5
    lon = np.array([_tile_lonlat(tx[0] + p / 256.0, ty[0], z)[0] for p in px])
    lat = np.array([_tile_lonlat(tx[0], ty[0] + p / 256.0, z)[1] for p in py])
    LON, LAT = np.meshgrid(lon, lat)
    X, Y = frame.fwd(LON, LAT)
    return X, Y, Zg, sources, meta0


def contour_lines(X, Y, Z, interval):
    zmin, zmax = np.nanmin(Z), np.nanmax(Z)
    levels = np.arange(math.floor(zmin / interval) * interval, zmax + interval, interval)
    fig = plt.figure()
    cs = plt.contour(X, Y, Z, levels=levels)
    out = []
    for lev, segs in zip(cs.levels, cs.allsegs):
        for seg in segs:
            if len(seg) >= 2:
                out.append((float(lev), LineString(seg)))
    plt.close(fig)
    return out


def pick_interval(relief):
    for iv in (1, 2, 5, 10, 20, 50):
        if relief / iv <= 30:
            return iv
    return 100


# --------------------------------------------------------------------------------------
# CLIMATE (Open-Meteo historical reanalysis, no key) and SUN (computed from latitude)
# --------------------------------------------------------------------------------------
OPEN_METEO_ARCHIVE = "https://archive-api.open-meteo.com/v1/archive"
SECTORS16 = ["N", "NNE", "NE", "ENE", "E", "ESE", "SE", "SSE", "S", "SSW", "SW", "WSW", "W", "WNW", "NW", "NNW"]
SPEED_BINS = [0.5, 2, 4, 6, 8, 99]  # m/s; below 0.5 = calm
SPEED_LABELS = ["0.5-2", "2-4", "4-6", "6-8", "8+"]
MONTHS = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]


def fetch_climate(net: Net, lat, lon, years=10):
    y2 = dt.date.today().year - 1
    y1 = y2 - years + 1
    params = {"latitude": round(lat, 4), "longitude": round(lon, 4),
              "start_date": f"{y1}-01-01", "end_date": f"{y2}-12-31",
              "hourly": "wind_speed_10m,wind_direction_10m",
              "daily": "temperature_2m_max,temperature_2m_min,precipitation_sum",
              "wind_speed_unit": "ms", "timezone": "auto"}
    obj, meta = net.json(OPEN_METEO_ARCHIVE, params=params, timeout=120, label="climate (Open-Meteo)")
    return obj, meta, (y1, y2)


def seasons_for(lat):
    # meteorological seasons, hemisphere-aware
    if lat < 0:
        return {"Summer (Dec-Feb)": (12, 1, 2), "Autumn (Mar-May)": (3, 4, 5),
                "Winter (Jun-Aug)": (6, 7, 8), "Spring (Sep-Nov)": (9, 10, 11)}
    return {"Winter (Dec-Feb)": (12, 1, 2), "Spring (Mar-May)": (3, 4, 5),
            "Summer (Jun-Aug)": (6, 7, 8), "Autumn (Sep-Nov)": (9, 10, 11)}


def wind_table(spd, dirn):
    """16 x 5 frequency table (% of all hours), and calm %."""
    spd = np.asarray(spd, float)
    dirn = np.asarray(dirn, float)
    ok = ~(np.isnan(spd) | np.isnan(dirn))
    spd, dirn = spd[ok], dirn[ok]
    n = len(spd)
    if n == 0:
        return np.zeros((16, 5)), 0.0, 0
    calm = float(np.sum(spd < SPEED_BINS[0])) / n * 100
    sec = (np.floor((dirn + 11.25) / 22.5).astype(int)) % 16
    tab = np.zeros((16, 5))
    for i in range(5):
        m = (spd >= SPEED_BINS[i]) & (spd < SPEED_BINS[i + 1])
        tab[:, i] = np.bincount(sec[m], minlength=16)[:16] / n * 100
    return tab, calm, n


def analyse_climate(obj, lat):
    h = obj["hourly"]
    times = np.array(h["time"], dtype="datetime64[m]")
    months = (times.astype("datetime64[M]").astype(int) % 12) + 1
    spd = np.array([np.nan if v is None else v for v in h["wind_speed_10m"]], float)
    dr = np.array([np.nan if v is None else v for v in h["wind_direction_10m"]], float)
    res = {"annual": wind_table(spd, dr), "seasons": {}}
    for name, ms in seasons_for(lat).items():
        m = np.isin(months, ms)
        res["seasons"][name] = wind_table(spd[m], dr[m])
    res["mean_speed"] = float(np.nanmean(spd))
    d = obj["daily"]
    dt_ = np.array(d["time"], dtype="datetime64[D]")
    dm = (dt_.astype("datetime64[M]").astype(int) % 12) + 1
    dy = dt_.astype("datetime64[Y]").astype(int) + 1970
    tmax = np.array([np.nan if v is None else v for v in d["temperature_2m_max"]], float)
    tmin = np.array([np.nan if v is None else v for v in d["temperature_2m_min"]], float)
    pr = np.array([np.nan if v is None else v for v in d["precipitation_sum"]], float)
    res["tmax"] = [float(np.nanmean(tmax[dm == k])) for k in range(1, 13)]
    res["tmin"] = [float(np.nanmean(tmin[dm == k])) for k in range(1, 13)]
    years = sorted(set(dy.tolist()))
    res["rain"] = [float(np.nanmean([np.nansum(pr[(dm == k) & (dy == y)]) for y in years])) for k in range(1, 13)]
    res["rain_year"] = float(sum(res["rain"]))
    res["grid"] = (obj.get("latitude"), obj.get("longitude"), obj.get("elevation"))
    return res


def prevailing(tab):
    tot = tab.sum(axis=1)
    i = int(np.argmax(tot))
    return SECTORS16[i], float(tot[i])


# ---- sun position: NOAA approximation, apparent solar time, accurate to a fraction of a degree ----
def declination(doy):
    g = 2 * math.pi / 365.0 * (doy - 1)
    return (0.006918 - 0.399912 * math.cos(g) + 0.070257 * math.sin(g) - 0.006758 * math.cos(2 * g)
            + 0.000907 * math.sin(2 * g) - 0.002697 * math.cos(3 * g) + 0.00148 * math.sin(3 * g))


def sun_pos(lat, doy, solar_hour):
    """(altitude deg, azimuth deg clockwise from true north)."""
    phi = math.radians(lat)
    dec = declination(doy)
    w = math.radians(15.0 * (solar_hour - 12.0))
    sa = math.sin(phi) * math.sin(dec) + math.cos(phi) * math.cos(dec) * math.cos(w)
    alt = math.asin(max(-1, min(1, sa)))
    az = math.atan2(math.sin(w), math.cos(w) * math.sin(phi) - math.tan(dec) * math.cos(phi)) + math.pi
    return math.degrees(alt), math.degrees(az) % 360


def day_length(lat, doy):
    phi = math.radians(lat)
    dec = declination(doy)
    c = -math.tan(phi) * math.tan(dec)
    if c <= -1:
        return 24.0
    if c >= 1:
        return 0.0
    return 2 * math.degrees(math.acos(c)) / 15.0


def key_days(lat):
    """doy for the dates that matter, labelled for the hemisphere."""
    if lat < 0:
        return [("Summer solstice 21 Dec", 355), ("Equinox 21 Mar / 23 Sep", 80), ("Winter solstice 21 Jun", 172)]
    return [("Summer solstice 21 Jun", 172), ("Equinox 21 Mar / 23 Sep", 80), ("Winter solstice 21 Dec", 355)]


def sun_curve(lat, doy, step=0.1):
    pts = []
    t = 0.0
    while t <= 24.0001:
        alt, az = sun_pos(lat, doy, t)
        if alt >= 0:
            pts.append((t, alt, az))
        t += step
    return pts


def sun_facts(lat):
    out = []
    for name, doy in key_days(lat):
        noon_alt, _ = sun_pos(lat, doy, 12.0)
        pts = sun_curve(lat, doy, 0.02)
        rise_az = pts[0][2] if pts else None
        set_az = pts[-1][2] if pts else None
        out.append(dict(name=name, doy=doy, noon_alt=noon_alt, day_len=day_length(lat, doy),
                        rise_az=rise_az, set_az=set_az))
    return out


def polar_xy(alt, az, R):
    """Equidistant sun-path projection in plan: centre = zenith, rim = horizon. North up."""
    rr = (90.0 - alt) / 90.0 * R
    a = math.radians(az)
    return rr * math.sin(a), rr * math.cos(a)


# --------------------------------------------------------------------------------------
# PIPELINE - each step fails on its own; the run always finishes
# --------------------------------------------------------------------------------------
STAT_GENERIC = [("zone", "AI_ZONE", "Land zoning"), ("hob", "AI_HOB", "Height of buildings"),
                ("fsr", "AI_FSR", "Floor space ratio"), ("addctrl", "AI_ADDCTRL", "Additional controls"), ("heritage", "AI_HERITAGE", "Heritage"),
                ("flood", "AI_FLOOD", "Flood planning"), ("bushfire", "AI_BUSHFIRE", "Bush fire prone land"),
                ("ass", "AI_ASS", "Acid sulfate soils")]


def new_layers():
    L = {}
    for key, dxf, title in [("site", "CTX_SITE", "Site lot"), ("cadastre", "CTX_CADASTRE", "Cadastre (lots)"),
                            ("buildings", "CTX_BUILDINGS", "Buildings"), ("roads", "CTX_ROADS", "Roads"),
                            ("paths", "CTX_PATHS", "Paths and cycleways"), ("rail", "CTX_RAIL", "Rail"),
                            ("water", "CTX_WATER", "Water"), ("trees", "CTX_TREES", "Trees and vegetation"),
                            ("landuse", "CTX_LANDUSE", "Land use and open space"),
                            ("contours", "CTX_CONTOURS", "Contours")]:
        L[key] = Layer(key, dxf, title, "CTX")
    for key, dxf, title in STAT_GENERIC + [("shr", "AI_HERITAGE_SHR", "State Heritage Register")]:
        L[key] = Layer(key, dxf, title, "AI")
    return L


class Ctx:
    pass


def _share(g, site):
    try:
        pct = g.intersection(site).area / max(site.area, 1e-6) * 100
    except Exception:
        return ""
    return "" if pct >= 99 else f" (about {pct:.0f}% of the lot, computed)"


def step_osm(c: Ctx):
    osm, meta, ep = fetch_osm(c.net, c.pin["lat"], c.pin["lon"], c.r, c.log)
    build_osm_layers(osm, c.frame, c.r, c.L)
    stamp = (osm.get("osm3s") or {}).get("timestamp_osm_base", "")
    src = f"OpenStreetMap via Overpass ({ep.split('/')[2]}), data as of {stamp[:10]}"
    for k in ("buildings", "roads", "paths", "rail", "water", "trees", "landuse"):
        L = c.L[k]
        L.status = "ok" if len(L) else "empty"
        L.source, L.url, L.retrieved = src, meta["url"], meta["retrieved"]
    b = c.L["buildings"]
    hs = [p["height_m"] for _, p in b.features if p["height_m"] is not None]
    tagged = sum(1 for _, p in b.features if p["height_basis"] == "height tag")
    lv = sum(1 for _, p in b.features if "storey" in p["height_basis"])
    how_osm = "Walk it, or compare with current aerial imagery (Google/Nearmap/SIX Maps). OSM is volunteer-mapped."
    c.reg.add("CTX_BUILDINGS", f"{len(b)} building footprints mapped within {c.r} m", src, meta["url"],
              meta["retrieved"], "medium", how_osm)
    if len(b):
        cov = (tagged + lv) / len(b) * 100
        c.reg.add("CTX_BUILDINGS", f"Heights known for {tagged + lv} of {len(b)} buildings ({cov:.0f}%): "
                  f"{tagged} from a height tag, {lv} estimated from storeys x 3.2 m; the rest have no height",
                  src, meta["url"], meta["retrieved"], "low" if cov < 50 else "medium",
                  "Spot-check the tallest three against a photo or the council DA record.")
        if cov < 30:
            c.reg.add("CTX_BUILDINGS", "OSM height coverage here is THIN - the massing and exploded axon under-represent the real city",
                      "tool inference", "", TODAY, "high", "Say so in your diagram caption.")
        if hs:
            i = int(np.argmax([p["height_m"] or -1 for _, p in b.features]))
            p = b.features[i][1]
            c.reg.add("CTX_BUILDINGS", f"Tallest mapped building within radius: {p['name'] or p['osm_id']} at about {p['height_m']:.0f} m ({p['height_basis']})",
                      src, f"https://www.openstreetmap.org/{ {'w': 'way', 'r': 'relation'}.get(p['osm_id'][0], 'way')}/{p['osm_id'][1:]}",
                      meta["retrieved"], "low", "Check against a published height for that building.")
    rail = c.L["rail"]
    if len(rail):
        kinds = sorted({p["kind"] for _, p in rail.features})
        tun = sum(1 for _, p in rail.features if p.get("tunnel") in ("yes", "building_passage"))
        c.reg.add("CTX_RAIL", f"Rail within radius: {', '.join(kinds)} ({tun} segment(s) in tunnel)", src, meta["url"],
                  meta["retrieved"], "medium", "Check the transport operator's network map; tunnels matter for structure and vibration.")
    else:
        c.reg.add("CTX_RAIL", "No rail mapped within radius", src, meta["url"], meta["retrieved"], "medium", how_osm)
    names = sorted({p["name"] for _, p in c.L["roads"].features if p.get("name")})
    c.reg.add("CTX_ROADS", f"{len(c.L['roads'])} road segments within radius" + (f", including {', '.join(names[:6])}" if names else ""),
              src, meta["url"], meta["retrieved"], "medium", "Road hierarchy and one-way status: check the council / state road register.")
    c.reg.add("CTX_WATER", ("Water mapped within radius: " + ", ".join(sorted({p['kind'] for _, p in c.L['water'].features})))
              if len(c.L["water"]) else "No water body or watercourse mapped within radius",
              src, meta["url"], meta["retrieved"], "medium", how_osm)
    nt = sum(1 for _, p in c.L["trees"].features if p["kind"] == "tree")
    c.reg.add("CTX_TREES", f"{nt} individual trees and {len(c.L['trees']) - nt} vegetated areas / tree rows mapped",
              src, meta["url"], meta["retrieved"], "low",
              "OSM rarely maps every street tree. Check council tree canopy data or a site visit.")


def step_cadastre(c: Ctx):
    ad = c.adapter
    site = c.L["site"]
    cad = c.L["cadastre"]
    if not ad.get("cadastre"):
        for L in (site, cad):
            L.status, L.note = "unavailable", f"no open cadastre adapter for {c.jur}"
        c.reg.add("CTX_CADASTRE", f"Lot boundaries not available for this jurisdiction ({c.jur}); OSM does not map land parcels",
                  "adapter table (empty slot)", "", TODAY, "n/a",
                  "Get the lot boundary from the title plan / survey or the state land registry viewer.")
        c.site_geom = Point(0, 0).buffer(6, 32)
        c.reg.add("CTX_SITE", "Site shown as a 6 m circle at the pin - NOT a lot boundary", "tool placeholder", "", TODAY, "high",
                  "Draw your real boundary in Rhino/Revit before using any site figure.")
        return
    url = ad["cadastre"]
    feats, meta = arcgis_query(c.net, url, c.frame, r=c.r, fmt="json",
                               out_fields="lotidstring,lotnumber,sectionnumber,planlabel", label="cadastre")
    circle = Point(0, 0).buffer(c.r, 128)
    for g, p in feats:
        gl = clip_to(c.frame.to_local(g), circle)
        if gl is not None:
            cad.add(gl, dict(lot=p.get("lotidstring", ""), plan=p.get("planlabel", "")))
    cad.status = "ok" if len(cad) else "empty"
    cad.source, cad.url, cad.retrieved = "NSW Spatial Services digital cadastre", meta["url"], meta["retrieved"]
    # the lot under the pin; if the pin is on a road, the nearest lot within 30 m
    hit, meta2 = arcgis_query(c.net, url, c.frame, point=(c.pin["lon"], c.pin["lat"]), fmt="json",
                              out_fields="lotidstring,planlabel", label="cadastre at pin")
    how = "Check lot / DP on the title search or SIX Maps; confirm the boundary shape against the DP plan."
    lots = [(c.frame.to_local(g), p) for g, p in hit]
    note = ""
    if not lots:
        cands = []
        for g, p in feats:
            gl = c.frame.to_local(g)
            d = gl.distance(Point(0, 0))
            if d <= 30:
                cands.append((d, gl, p))
        if cands:
            cands.sort(key=lambda t: t[0])
            lots = [(cands[0][1], cands[0][2])]
            note = f" (pin fell {cands[0][0]:.0f} m outside any lot - nearest lot taken)"
    if lots:
        g, p = lots[0]
        g = g if g.is_valid else g.buffer(0)
        site.add(g, dict(lot=p.get("lotidstring", ""), plan=p.get("planlabel", ""), area_m2=round(g.area, 1)))
        site.status, site.source, site.url, site.retrieved = "ok", cad.source, meta2["url"], meta2["retrieved"]
        c.site_geom = g
        c.site_label = p.get("lotidstring", "")
        c.reg.add("CTX_SITE", f"Site lot is {p.get('lotidstring', '?')}{note}", cad.source, meta2["url"], meta2["retrieved"],
                  "low" if note else "medium", how)
        c.reg.add("CTX_SITE", f"Lot area about {g.area:,.0f} m2 (computed from the digital cadastre, not survey)", "tool calculation",
                  meta2["url"], meta2["retrieved"], "medium",
                  "Compare with the area on the DP / title. Your site may be several lots - add them in CAD.")
        if len(lots) > 1 or any(p.get("lotidstring", "").startswith("CP/") for _, p in lots):
            c.reg.add("CTX_SITE", "Strata / common property at the pin - the building lot may be a strata plan",
                      cad.source, meta2["url"], meta2["retrieved"], "low", how)
    else:
        site.status, site.note = "empty", "no lot at or near the pin"
        c.site_geom = Point(0, 0).buffer(6, 32)
        c.reg.add("CTX_SITE", "No lot found at or within 30 m of the pin; site shown as a 6 m circle", cad.source,
                  meta2["url"], meta2["retrieved"], "high", "Move the pin (rerun with 'lat, lon') or draw the boundary yourself.")
    c.reg.add("CTX_CADASTRE", f"{len(cad)} lots within {c.r} m", cad.source, meta["url"], meta["retrieved"], "high",
              "Digital cadastre is indicative (can be metres out in older areas). Survey governs.")
    nb = len(c.L["buildings"])
    if nb and len(cad) and nb < 0.5 * len(cad):
        c.reg.add("CTX_BUILDINGS", f"OSM shows {nb} buildings against {len(cad)} lots - likely UNDER-MAPPED here "
                  "(missing buildings, or whole terrace rows drawn as one building)", "tool comparison of two sources",
                  "", TODAY, "medium", "Compare the footprints with an aerial photo before you use them for grain or density.")


def step_terrain(c: Ctx):
    L = c.L["contours"]
    ad = c.adapter
    circle = Point(0, 0).buffer(c.r, 128)
    how = "Compare with a survey or the state LiDAR (ELVIS: elevation.fsdf.io). Never design levels from this."
    if ad.get("contours"):
        try:
            feats, meta = arcgis_query(c.net, ad["contours"], c.frame, r=c.r, fmt="json",
                                       out_fields="elevation,verticalaccuracy", label="contours (NSW)", max_pages=15)
            for g, p in feats:
                gl = clip_to(c.frame.to_local(g), circle)
                if gl is not None and p.get("elevation") is not None:
                    L.add(gl, dict(elev=float(p["elevation"])))
            if len(L):
                els = sorted({p["elev"] for _, p in L.features})
                steps = np.diff(els)
                iv = float(np.min(steps)) if len(steps) else 0
                L.status, L.source, L.url, L.retrieved = "ok", "NSW Spatial Services topographic contours", meta["url"], meta["retrieved"]
                L.note = f"interval {iv:g} m"
                c.reg.add("CTX_CONTOURS", f"Ground within radius ranges {els[0]:g}-{els[-1]:g} m (AHD), contours at {iv:g} m interval",
                          L.source, meta["url"], meta["retrieved"], "medium", how)
                near = min(L.features, key=lambda t: t[0].distance(c.site_geom))
                c.reg.add("CTX_CONTOURS", f"Nearest contour to the site is {near[1]['elev']:g} m AHD", L.source, meta["url"],
                          meta["retrieved"], "low", how)
                c.terrain_iv = iv
                return
            c.log("   NSW contours empty here - falling back to SRTM terrain tiles")
        except NotCached:
            raise
        except Exception as e:
            c.log(f"   NSW contours failed ({e}) - falling back to SRTM terrain tiles")
    X, Y, Z, sources, meta = terrain_grid(c.net, c.frame, c.r, c.log)
    m = (X ** 2 + Y ** 2) <= (c.r * 1.05) ** 2
    relief = float(np.nanmax(Z[m]) - np.nanmin(Z[m]))
    iv = max(2, pick_interval(relief))
    for lev, ls in contour_lines(X, Y, Z, iv):
        gl = clip_to(ls, circle)
        if gl is not None:
            L.add(gl, dict(elev=lev))
    src = ", ".join(sorted(sources)) or "unknown"
    L.status = "ok" if len(L) else "empty"
    L.source = f"Terrain Tiles on AWS (Terrarium) - source DEM: {src}"
    L.url, L.retrieved = meta["url"], meta["retrieved"]
    L.note = f"interval {iv} m from a {'30 m SRTM' if 'srtm' in src else src} surface"
    zc = float(Z[np.unravel_index(np.argmin(X ** 2 + Y ** 2), X.shape)])
    c.reg.add("CTX_CONTOURS", f"Ground at the pin about {zc:.0f} m above sea level; relief within radius about {relief:.0f} m; contours drawn at {iv} m",
              L.source, meta["url"], meta["retrieved"], "low", how)
    if "srtm" in src:
        c.reg.add("CTX_CONTOURS", "Terrain here is from 30 m SRTM radar data: it smooths out anything smaller than a street block and can read roof tops as ground",
                  "Tilezen data sources list", "https://github.com/tilezen/joerd/blob/master/docs/data-sources.md", TODAY, "high", how)
    c.terrain_iv = iv


def step_statutory(c: Ctx):
    ad = c.adapter
    if not ad.get("statutory"):
        for key, dxf, title in STAT_GENERIC:
            L = c.L[key]
            L.status, L.note = "unavailable", f"not available for this jurisdiction ({c.jur})"
            c.reg.add(dxf, f"{title}: not available for this jurisdiction ({c.jur}) - adapter slot empty",
                      "adapter table", "", TODAY, "n/a",
                      f"Look it up in the local planning viewer{': ' + ad['planning_viewer'] if ad.get('planning_viewer') else ''} and add it by hand.")
        c.L["shr"].status = "unavailable"
        return
    circle = Point(0, 0).buffer(c.r, 128)
    site = c.site_geom
    near50 = site.buffer(50)
    viewer = ad.get("planning_viewer", "")
    for spec in ad["statutory"]:
        L = c.L[spec["key"]]
        L.dxf, L.title = spec["dxf"], spec["title"]
        try:
            feats, meta = [], None
            for u in (spec["url"] if isinstance(spec["url"], list) else [spec["url"]]):
                f_, m_ = arcgis_query(c.net, u, c.frame, r=c.r, fmt="geojson", label=spec["title"])
                dom = layer_schema(c.net, u)
                feats += [(g, decode(p, dom)) for g, p in f_]
                meta = meta or m_
        except NotCached:
            raise
        except Exception as e:
            L.status, L.note = "failed", str(e)
            c.reg.add(spec["dxf"], f"{spec['title']}: service did not answer ({e}) - layer missing from this run",
                      "NSW Planning Portal map service", spec["url"], TODAY, "n/a",
                      f"Rerun later, or read it off the Spatial Viewer: {viewer}")
            continue
        L.source, L.url, L.retrieved = "NSW Planning Portal (Dept of Planning) map service", meta["url"], meta["retrieved"]
        for g, p in feats:
            gl = clip_to(c.frame.to_local(g), circle)
            if gl is not None:
                pp = dict(p)
                pp["_label"] = spec["label"](p)
                pp["_desc"] = spec["desc"](p)
                pp["_epi"] = p.get("EPI_NAME") or ""
                pp["_currency"] = epoch_ms_to_date(p.get("CURRENCY_DATE") or p.get("LastUpdate") or p.get("LAST_UPDATE"))
                pp["_clause"] = p.get("LEGIS_REF_CLAUSE") or ""
                L.add(gl, pp)
        L.status = "ok" if len(L) else "empty"
        how = (f"Spatial Viewer ({viewer}) - search the address, switch on this layer; then read the clause on "
               f"legislation.nsw.gov.au. A planning certificate (s10.7) is the authoritative answer.")
        tol = max(1.0, 0.01 * site.area)  # ignore slivers: boundaries from different sources never meet exactly
        on_site = [(g, p) for g, p in L.features if g.intersects(site) and g.intersection(site).area > tol]
        if spec["key"] in ("zone", "hob", "fsr", "ass", "addctrl"):
            if on_site:
                parts = []
                for g, p in sorted(on_site, key=lambda t: -t[0].intersection(site).area):
                    share = g.intersection(site).area / max(site.area, 1e-6) * 100
                    parts.append(f"{p['_desc']}" + (f" ({share:.0f}% of site)" if len(on_site) > 1 else ""))
                p0 = on_site[0][1]
                src = f"{p0['_epi']}{', ' + p0['_clause'] if p0['_clause'] else ''} - map currency {p0['_currency'] or '?'}"
                c.reg.add(spec["dxf"], f"{spec['title']} on the site: " + "; ".join(parts), src, meta["url"], meta["retrieved"],
                          "medium", how)
            else:
                if spec["key"] == "addctrl":
                    continue
                c.reg.add(spec["dxf"], f"{spec['title']}: nothing mapped on the site" +
                          (" (this LEP may not map it, or a site-specific clause applies)" if spec["key"] in ("hob", "fsr") else ""),
                          "NSW Planning Portal map service", meta["url"], meta["retrieved"], "low", how)
        elif spec["key"] in ("heritage", "shr"):
            if on_site:
                for g, p in on_site[:5]:
                    c.reg.add(spec["dxf"], f"Site is affected by {p['_desc']}{_share(g, site)}",
                              f"{p['_epi'] or 'State Heritage Register'} - currency {p['_currency'] or '?'}",
                              meta["url"], meta["retrieved"], "medium", how + " Heritage inventory sheet: hms.heritage.nsw.gov.au")
            else:
                c.reg.add(spec["dxf"], f"No {spec['title']} mapped on the site", "NSW Planning Portal map service",
                          meta["url"], meta["retrieved"], "medium", how)
            nb = [(g, p) for g, p in L.features if g.intersects(near50) and (g, p) not in on_site]
            if nb:
                lab = "; ".join(f"{p['_label']} {p.get('H_NAME') or p.get('ITEMNAME') or ''}".strip() for _, p in nb[:6])
                c.reg.add(spec["dxf"], f"{len(nb)} {spec['title']} feature(s) within 50 m of the site: {lab}"
                          + (" ..." if len(nb) > 6 else ""), "NSW Planning Portal map service", meta["url"],
                          meta["retrieved"], "medium", "Adjoining heritage changes what you can do next to it. " + how)
            if len(L):
                c.reg.add(spec["dxf"], f"{len(L)} {spec['title']} feature(s) within {c.r} m", "NSW Planning Portal map service",
                          meta["url"], meta["retrieved"], "high", how)
        elif spec["key"] == "flood":
            if on_site:
                fu = unary_union([g for g, _ in on_site])
                c.reg.add("AI_FLOOD", "Site is inside a mapped flood planning area" + _share(fu, site) + ": "
                          + "; ".join(sorted({p["_desc"] for _, p in on_site})),
                          f"{on_site[0][1]['_epi']}", meta["url"], meta["retrieved"], "medium",
                          how + " Get the flood level from council (flood study / flood certificate).")
            else:
                lga = next((p.get("LGA_NAME") for g, p in c.L["zone"].features if g.intersects(site) and p.get("LGA_NAME")), None)
                n = None
                if lga:
                    try:
                        cnt, cm = c.net.json(spec["url"].rstrip("/") + "/query", params={
                            "where": f"LGA_NAME='{lga}'", "returnCountOnly": "true", "f": "json"}, label="flood map coverage")
                        n = int(cnt.get("count", 0))
                    except NotCached:
                        raise
                    except Exception:
                        n = None
                if n == 0:
                    L.note = f"NOT CHECKED - {lga.title()} council publishes no flood map here"
                    c.reg.add("AI_FLOOD", f"NOT CHECKED: {lga.title()} council publishes NO flood planning map on the NSW Planning Portal, "
                              "so an empty flood layer says nothing about this site", "NSW Planning Portal map service (count by council)",
                              meta["url"], meta["retrieved"], "n/a",
                              "Council flood maps / flood information request, and the s10.7 planning certificate.")
                elif n:
                    c.reg.add("AI_FLOOD", f"{lga.title()} does publish a flood planning map ({n} mapped areas); none covers the site",
                              "NSW Planning Portal map service", meta["url"], meta["retrieved"], "medium",
                              "Council flood maps; s10.7 certificate. Flood planning maps are not the only flood risk (overland flow).")
                else:
                    c.reg.add("AI_FLOOD", "No flood planning area mapped on the site IN THIS LAYER - not the same as 'not flood affected': "
                              "most NSW councils publish no flood map to the Planning Portal (checked 23 Sep 2026)",
                              "NSW Planning Portal map service", meta["url"], meta["retrieved"], "low",
                              "Council flood maps / flood information request, and the s10.7 planning certificate.")
        elif spec["key"] == "bushfire":
            if on_site:
                bu = unary_union([g for g, _ in on_site])
                c.reg.add("AI_BUSHFIRE", "Site is mapped bush fire prone land" + _share(bu, site) + ": "
                          + "; ".join(sorted({p['_desc'] for _, p in on_site})),
                          f"NSW Bush Fire Prone Land map - updated {on_site[0][1]['_currency'] or '?'}", meta["url"], meta["retrieved"],
                          "medium", "NSW RFS 'Check if you're in bush fire prone land' tool; Planning for Bush Fire Protection applies.")
            else:
                c.reg.add("AI_BUSHFIRE", "Site not mapped as bush fire prone land", "NSW Bush Fire Prone Land map (Planning Portal)",
                          meta["url"], meta["retrieved"], "medium", "NSW RFS online check; planning certificate.")
            if len(L) and not on_site:
                c.reg.add("AI_BUSHFIRE", f"Bush fire prone land mapped within {c.r} m of the site", "Planning Portal", meta["url"],
                          meta["retrieved"], "medium", "Buffer categories reach into neighbouring lots - check distance to the site.")


def step_climate(c: Ctx):
    obj, meta, (y1, y2) = fetch_climate(c.net, c.pin["lat"], c.pin["lon"])
    cl = analyse_climate(obj, c.pin["lat"])
    c.climate = cl
    c.climate_meta = meta
    glat, glon, gel = cl["grid"]
    src = f"Open-Meteo historical weather API (ERA5 reanalysis), {y1}-{y2}, grid point {glat:.3f}, {glon:.3f}"
    how = "Compare with the nearest Bureau of Meteorology station (bom.gov.au/climate/averages). Reanalysis is a 10-25 km model grid, not your street."
    d, pct = prevailing(cl["annual"][0])
    c.reg.add("AI_WIND", f"Prevailing wind (all hours, 10 m) from the {d} ({pct:.0f}% of hours); mean speed {cl['mean_speed']:.1f} m/s; calm {cl['annual'][1]:.0f}%",
              src, meta["url"], meta["retrieved"], "low", how)
    for name, (tab, calm, n) in cl["seasons"].items():
        d, pct = prevailing(tab)
        c.reg.add("AI_WIND", f"{name}: wind mostly from the {d} ({pct:.0f}% of hours)", src, meta["url"], meta["retrieved"], "low", how)
    tmax, tmin, rain = cl["tmax"], cl["tmin"], cl["rain"]
    hi, lo, wet = int(np.argmax(tmax)), int(np.argmin(tmin)), int(np.argmax(rain))
    c.reg.add("AI_CLIMATE", f"Warmest month {MONTHS[hi]} (mean daily max {tmax[hi]:.1f} C); coolest {MONTHS[lo]} (mean daily min {tmin[lo]:.1f} C)",
              src, meta["url"], meta["retrieved"], "medium", how)
    c.reg.add("AI_CLIMATE", f"Mean annual rainfall about {cl['rain_year']:.0f} mm; wettest month {MONTHS[wet]} ({rain[wet]:.0f} mm)",
              src, meta["url"], meta["retrieved"], "low", how + " Reanalysis rainfall is often out by 10-30%.")


def step_sun(c: Ctx):
    lat = c.pin["lat"]
    c.sun = sun_facts(lat)
    for f in c.sun:
        c.reg.add("AI_SUNPATH", f"{f['name']}: sun at noon {f['noon_alt']:.0f} deg above the horizon; daylight {f['day_len']:.1f} h; "
                  f"rises at azimuth {f['rise_az']:.0f} deg, sets at {f['set_az']:.0f} deg (true north)",
                  "computed from latitude (NOAA solar position approximation)", "", TODAY, "high",
                  "Check with any sun-path tool (e.g. andrewmarsh.com/apps sun path, or Ladybug). Times are SOLAR time, not clock time.")
    c.reg.add("AI_SUNPATH", "North in all drawings is TRUE north. Magnetic north differs by several degrees; grid north (MGA) by up to ~1.5 deg",
              "tool convention", "", TODAY, "high", "Check the north point on the survey you use.")


def standing_gaps(c: Ctx):
    how = "Outside this tool. Name where you checked it, or say it was not checked."
    for a in ["Title, survey, easements, covenants and restrictions on the title - not checked by this tool",
              "Planning certificate (NSW: section 10.7) - not checked; only the council issues it",
              "Development control plan (setbacks, landscape, parking, solar access) - not read by this tool",
              "Every statutory layer is a point-in-time copy of a map service; instruments change - check currency on the day you submit",
              "Layers from different sources (OSM, cadastre, planning maps) can disagree by 1-2 m or more: different datums and digitising"]:
        c.reg.add("GAP", a, "tool", "", TODAY, "n/a", how)


# --------------------------------------------------------------------------------------
# DXF - one layer scheme, documented in LAYER_SCHEME.md. Units metres, origin = pin,
# true north = +Y. CTX_ = physical context. AI_ = pulled or computed by the tool: verify.
# --------------------------------------------------------------------------------------
# name: (ACI colour, lineweight in 1/100 mm, linetype, description)
LAYER_SPEC = {
    "CTX_ORIGIN":          (7, 25, "CONTINUOUS", "Pin at 0,0; north point; data credits"),
    "CTX_RADIUS":          (8, 13, "DASHED", "Study radius"),
    "CTX_SITE":            (1, 50, "CONTINUOUS", "Site lot boundary (cadastre at the pin)"),
    "CTX_CADASTRE":        (8, 13, "CONTINUOUS", "Lot boundaries"),
    "CTX_CADASTRE_TX":     (8, 9, "CONTINUOUS", "Lot labels"),
    "CTX_BUILDINGS":       (7, 25, "CONTINUOUS", "Building footprints (OSM)"),
    "CTX_MASSING":         (7, 13, "CONTINUOUS", "3D faces: buildings with a known or estimated height only"),
    "CTX_ROADS":           (8, 18, "CONTINUOUS", "Road centrelines (tunnels dashed)"),
    "CTX_PATHS":           (9, 13, "DASHED2", "Footpaths, cycleways, steps"),
    "CTX_RAIL":            (8, 35, "CONTINUOUS", "Rail, light rail, metro (tunnels dashed)"),
    "CTX_WATER":           (5, 18, "CONTINUOUS", "Water bodies, watercourses, coastline"),
    "CTX_TREES":           (3, 13, "CONTINUOUS", "Mapped trees (circles), tree rows, woodland"),
    "CTX_LANDUSE":         (62, 9, "CONTINUOUS", "Land use and open space areas"),
    "CTX_CONTOURS":        (33, 9, "CONTINUOUS", "Contours at elevation (Z = level)"),
    "CTX_CONTOURS_INDEX":  (33, 25, "CONTINUOUS", "Index contours at elevation"),
    "CTX_CONTOURS_TX":     (33, 9, "CONTINUOUS", "Contour labels"),
    "AI_ZONE":             (30, 25, "CONTINUOUS", "Land zoning boundaries"),
    "AI_ZONE_TX":          (30, 9, "CONTINUOUS", "Zone codes"),
    "AI_HOB":              (4, 18, "CONTINUOUS", "Height of buildings control boundaries"),
    "AI_HOB_TX":           (4, 9, "CONTINUOUS", "Height control labels"),
    "AI_FSR":              (6, 18, "CONTINUOUS", "Floor space ratio boundaries"),
    "AI_FSR_TX":           (6, 9, "CONTINUOUS", "FSR labels"),
    "AI_ADDCTRL":          (6, 13, "DASHED", "Additional controls areas (LEP clauses vary height / FSR here)"),
    "AI_ADDCTRL_TX":       (6, 9, "CONTINUOUS", "Additional controls labels"),
    "AI_HERITAGE":         (210, 25, "CONTINUOUS", "LEP heritage items and conservation areas"),
    "AI_HERITAGE_TX":      (210, 9, "CONTINUOUS", "Heritage IDs"),
    "AI_HERITAGE_SHR":     (200, 35, "DASHED", "State Heritage Register curtilage"),
    "AI_FLOOD":            (150, 25, "CONTINUOUS", "Flood planning area (where mapped)"),
    "AI_BUSHFIRE":         (10, 25, "CONTINUOUS", "Bush fire prone land"),
    "AI_ASS":              (40, 13, "DASHDOT", "Acid sulfate soils classes"),
    "AI_ASS_TX":           (40, 9, "CONTINUOUS", "Acid sulfate soils class labels"),
    "AI_SUNPATH":          (2, 18, "CONTINUOUS", "Sun path diagram at the pin (diagrammatic scale)"),
    "AI_WIND":             (140, 18, "CONTINUOUS", "Wind rose at the pin (diagrammatic scale)"),
}

TX_LAYER = {"cadastre": "CTX_CADASTRE_TX", "zone": "AI_ZONE_TX", "hob": "AI_HOB_TX", "fsr": "AI_FSR_TX", "addctrl": "AI_ADDCTRL_TX",
            "heritage": "AI_HERITAGE_TX", "ass": "AI_ASS_TX"}


def _add_poly(msp, coords, layer, closed=True, elevation=0.0, linetype=None):
    pts = [(float(x), float(y)) for x, y in list(coords)[: (-1 if closed else None)]]
    if len(pts) < 2:
        return
    att = {"layer": layer}
    if elevation:
        att["elevation"] = float(elevation)
    if linetype:
        att["linetype"] = linetype
    msp.add_lwpolyline(pts, close=closed, dxfattribs=att)


def add_geom(msp, geom, layer, elevation=0.0, linetype=None, point_r=1.5):
    for part in iter_parts(geom):
        if part.geom_type == "Polygon":
            _add_poly(msp, part.exterior.coords, layer, True, elevation, linetype)
            for ring in part.interiors:
                _add_poly(msp, ring.coords, layer, True, elevation, linetype)
        elif part.geom_type in ("LineString", "LinearRing"):
            _add_poly(msp, part.coords, layer, False, elevation, linetype)
        elif part.geom_type == "Point":
            msp.add_circle((part.x, part.y), point_r, dxfattribs={"layer": layer})


def add_text(msp, text, xy, h, layer, align="MIDDLE_CENTER", rot=0.0):
    from ezdxf.enums import TextEntityAlignment
    if not text:
        return
    t = msp.add_text(str(text), height=h, rotation=rot, dxfattribs={"layer": layer})
    t.set_placement((float(xy[0]), float(xy[1])), align=getattr(TextEntityAlignment, align))


def label_point(geom):
    try:
        g = max(iter_parts(geom), key=lambda p: p.area if p.geom_type == "Polygon" else p.length)
        return g.representative_point() if g.geom_type == "Polygon" else g.interpolate(0.5, normalized=True)
    except Exception:
        return geom.representative_point()


def triangles(poly):
    """Triangulate a polygon (with holes) for 3D caps."""
    try:
        import shapely
        tris = shapely.constrained_delaunay_triangles(poly)
        return [t for t in tris.geoms]
    except Exception:
        from shapely.ops import triangulate
        return [t for t in triangulate(poly) if t.representative_point().within(poly)]


def add_massing(msp, poly, h, layer="CTX_MASSING"):
    for part in iter_parts(poly):
        if part.geom_type != "Polygon":
            continue
        part = part.simplify(0.05)
        for ring in [part.exterior] + list(part.interiors):
            c = list(ring.coords)
            for (x1, y1), (x2, y2) in zip(c[:-1], c[1:]):
                msp.add_3dface([(x1, y1, 0), (x2, y2, 0), (x2, y2, h), (x1, y1, h)], dxfattribs={"layer": layer})
        for t in triangles(part):
            tc = list(t.exterior.coords)[:3]
            msp.add_3dface([(x, y, h) for x, y in tc] + [(tc[-1][0], tc[-1][1], h)], dxfattribs={"layer": layer})
            msp.add_3dface([(x, y, 0) for x, y in tc] + [(tc[-1][0], tc[-1][1], 0)], dxfattribs={"layer": layer})


def write_dxf(c: Ctx, path: Path):
    doc = ezdxf.new("R2010", setup=True)
    doc.header["$INSUNITS"] = 6      # metres
    doc.header["$MEASUREMENT"] = 1   # metric
    doc.header["$LUNITS"] = 2
    for name, (col, lw, lt, desc) in LAYER_SPEC.items():
        lay = doc.layers.add(name, color=col, linetype=lt if lt in doc.linetypes else "CONTINUOUS")
        lay.dxf.lineweight = lw
        try:
            lay.description = desc
        except Exception:
            pass
    msp = doc.modelspace()
    r = c.r
    th = max(1.2, r / 150.0)  # text height scales with the radius

    # origin, radius, north
    msp.add_circle((0, 0), r, dxfattribs={"layer": "CTX_RADIUS"})
    msp.add_line((-3, 0), (3, 0), dxfattribs={"layer": "CTX_ORIGIN"})
    msp.add_line((0, -3), (0, 3), dxfattribs={"layer": "CTX_ORIGIN"})
    nx, ny, s = r * 0.9, r * 0.9, r * 0.05
    msp.add_lwpolyline([(nx, ny + s), (nx - s * 0.5, ny - s), (nx, ny - s * 0.5), (nx + s * 0.5, ny - s)], close=True,
                       dxfattribs={"layer": "CTX_ORIGIN"})
    add_text(msp, "N (true)", (nx, ny + s * 1.6), th, "CTX_ORIGIN")
    lines = [f"{TOOL} v{VERSION} - {c.address}",
             f"Origin 0,0 = pin at lat {c.pin['lat']:.7f}, lon {c.pin['lon']:.7f} - units metres - +Y = TRUE north",
             f"Generated {TODAY}. Jurisdiction {c.jur}. Study radius {r} m. CTX_ = context, AI_ = pulled/computed: verify against the register.",
             "Data: (c) OpenStreetMap contributors ODbL; NSW Spatial Services & NSW Dept of Planning CC BY 4.0; Open-Meteo CC BY 4.0; Terrain Tiles (SRTM et al.)"]
    for i, t in enumerate(lines):
        add_text(msp, t, (-r, -r - th * 2.2 * (i + 1)), th * 0.8, "CTX_ORIGIN", align="LEFT")

    L = c.L
    # physical
    for g, p in L["landuse"].features:
        add_geom(msp, g, "CTX_LANDUSE")
    for g, p in L["water"].features:
        add_geom(msp, g, "CTX_WATER")
    for g, p in L["trees"].features:
        cr = fnum(p.get("crown_m")) or 5.0
        add_geom(msp, g, "CTX_TREES", point_r=cr / 2)
    for key, layer in (("roads", "CTX_ROADS"), ("paths", "CTX_PATHS"), ("rail", "CTX_RAIL")):
        for g, p in L[key].features:
            lt = "DASHED" if p.get("tunnel") in ("yes", "building_passage", "culvert") else None
            add_geom(msp, g, layer, linetype=lt)
    iv = getattr(c, "terrain_iv", 2) or 2
    for g, p in L["contours"].features:
        z = p["elev"]
        idx = abs((z / (iv * 5)) - round(z / (iv * 5))) < 1e-6
        add_geom(msp, g, "CTX_CONTOURS_INDEX" if idx else "CTX_CONTOURS", elevation=z)
        if idx:
            lp = g.interpolate(0.5, normalized=True) if g.geom_type == "LineString" else label_point(g)
            add_text(msp, f"{z:g}", (lp.x, lp.y), th * 0.7, "CTX_CONTOURS_TX")
    for g, p in L["cadastre"].features:
        add_geom(msp, g, "CTX_CADASTRE")
    if r <= 400:
        for g, p in L["cadastre"].features:
            lp = label_point(g)
            add_text(msp, p.get("lot", "").split("//")[0], (lp.x, lp.y), th * 0.45, "CTX_CADASTRE_TX")
    for g, p in L["buildings"].features:
        add_geom(msp, g, "CTX_BUILDINGS")
        if p.get("height_m"):
            try:
                add_massing(msp, g, float(p["height_m"]))
            except Exception:
                pass
    for g, p in L["site"].features:
        add_geom(msp, g, "CTX_SITE")
    if not len(L["site"]) and getattr(c, "site_geom", None) is not None:
        add_geom(msp, c.site_geom, "CTX_SITE")

    # statutory
    for key in ("zone", "hob", "fsr", "addctrl", "heritage", "shr", "flood", "bushfire", "ass"):
        lay = L[key]
        for g, p in lay.features:
            add_geom(msp, g, lay.dxf)
            if key in TX_LAYER:
                lp = label_point(g)
                add_text(msp, p.get("_label", ""), (lp.x, lp.y), th * 0.9, TX_LAYER[key])

    # sun path and wind rose, drawn at the pin to a diagrammatic scale
    if getattr(c, "sun", None):
        R = r * 0.35
        for alt in (0, 30, 60):
            msp.add_circle((0, 0), (90 - alt) / 90 * R, dxfattribs={"layer": "AI_SUNPATH"})
        for az, lab in ((0, "N"), (90, "E"), (180, "S"), (270, "W")):
            x, y = polar_xy(0, az, R * 1.08)
            add_text(msp, lab, (x, y), th, "AI_SUNPATH")
        for name, doy in key_days(c.pin["lat"]):
            pts = [polar_xy(alt, az, R) for t, alt, az in sun_curve(c.pin["lat"], doy, 0.1)]
            if len(pts) > 1:
                msp.add_lwpolyline(pts, dxfattribs={"layer": "AI_SUNPATH"})
        for hr in range(5, 20):
            pts = []
            for doy in range(1, 366, 5):
                alt, az = sun_pos(c.pin["lat"], doy, hr)
                if alt > 0:
                    pts.append((doy, polar_xy(alt, az, R)))
            if len(pts) > 2:
                pts.sort(key=lambda t: declination(t[0]))
                msp.add_lwpolyline([p for _, p in pts], dxfattribs={"layer": "AI_SUNPATH", "linetype": "DASHED2"})
                x, y = pts[0][1]
                add_text(msp, f"{hr}h", (x, y), th * 0.6, "AI_SUNPATH")
    if getattr(c, "climate", None):
        tab = c.climate["annual"][0]
        tot = tab.sum(axis=1)
        Rw = r * 0.3 / max(tot.max(), 1e-6)
        for i in range(16):
            a0 = math.radians(i * 22.5 - 9)
            a1 = math.radians(i * 22.5 + 9)
            rr = tot[i] * Rw
            msp.add_lwpolyline([(0, 0), (rr * math.sin(a0), rr * math.cos(a0)), (rr * math.sin(a1), rr * math.cos(a1))],
                               close=True, dxfattribs={"layer": "AI_WIND"})
        add_text(msp, f"wind rose: petal = % of hours wind blows FROM that direction (max {tot.max():.1f}%)",
                 (0, -r * 0.33), th * 0.6, "AI_WIND")
    doc.saveas(path)
    return path


# --------------------------------------------------------------------------------------
# STYLE CARD - written once by the student. Drives every diagram and the prompt suffix.
# The tool never chooses the look; missing keys fall back to the neutral defaults below.
# --------------------------------------------------------------------------------------
DEFAULT_STYLE = {
    "meta": {"author": "", "studio": "", "title_prefix": "SITE LAYERS-OPEN"},
    "palette": {
        "background": "#F4F1EA", "ink": "#1D1D1B", "plate": "#FBFAF6",
        "buildings": "#C9C2B6", "buildings_edge": "#5E584F", "massing_top": "#E4DED3", "massing_side": "#A69E91",
        "site": "#D1462F", "roads": "#8F8A82", "paths": "#A7A198", "rail": "#2F2F2D",
        "water": "#A7BFCC", "green": "#C3CDAE", "trees": "#7E9A6A", "landuse": "#E7E2D8",
        "contours": "#A08F79", "cadastre": "#A39D93",
        "categorical": ["#E3A587", "#8FAFC4", "#C9B458", "#9DB88E", "#B394B8", "#D6897E", "#7FA7A0", "#C7A27C",
                        "#8C9BC9", "#B5B5B5"],
        "sequential_low": "#F1E7CF", "sequential_high": "#7C3A2D",
        "heritage": "#7A5230", "flood": "#2F6F9F", "bushfire": "#C0502C", "acid_sulfate": "#8A7A2A",
        "sun": "#D9A21B", "wind": "#3E6C8A"},
    "lines": {"hairline": 0.25, "fine": 0.45, "medium": 0.9, "heavy": 1.8},
    "hatch": {"heritage": "////", "shr": "xx", "flood": "....", "bushfire": "\\\\\\\\", "acid_sulfate": "--"},
    "type": {"family": "DejaVu Sans", "size": 7.0, "title_size": 12.0, "case": "upper"},
    "references": {"image_1": "", "image_2": "", "image_3": ""},
    "prompt": {"medium": "flat architectural analysis diagram, orthographic", "suffix": ""},
}


def _merge(a, b):
    out = dict(a)
    for k, v in (b or {}).items():
        out[k] = _merge(a[k], v) if isinstance(v, dict) and isinstance(a.get(k), dict) else v
    return out


def load_style(path, log=print):
    if not path:
        return json.loads(json.dumps(DEFAULT_STYLE)), "defaults (no style card given)"
    p = Path(path)
    if not p.exists():
        log(f"   style card {p} not found - using defaults")
        return json.loads(json.dumps(DEFAULT_STYLE)), "defaults (style card not found)"
    txt = p.read_text(encoding="utf-8")
    try:
        if tomllib:
            data = tomllib.loads(txt)
        else:
            data = _mini_toml(txt)
    except Exception as e:
        log(f"   style card could not be read ({e}) - using defaults. Check quotes and brackets.")
        return json.loads(json.dumps(DEFAULT_STYLE)), f"defaults (style card error: {e})"
    return _merge(DEFAULT_STYLE, data), str(p.name)


def _mini_toml(txt):
    """Enough TOML for a style card on Python < 3.11: [sections], key = "str" | number | [list]."""
    data = {}
    cur = data
    for raw in txt.splitlines():
        line = re.sub(r'\s+#[^"\]]*$', "", raw).strip()
        if not line or line.startswith("#"):
            continue
        m = re.match(r"^\[([^\]]+)\]$", line)
        if m:
            cur = data.setdefault(m.group(1).strip(), {})
            continue
        if "=" in line:
            k, v = [s.strip() for s in line.split("=", 1)]
            cur[k] = json.loads(v.replace("'", '"')) if v[:1] in '["' else float(v)
    return data


def apply_font(style, log=print):
    from matplotlib import font_manager
    fam = style["type"].get("family") or "DejaVu Sans"
    try:
        font_manager.findfont(fam, fallback_to_default=False)
        used = fam
    except Exception:
        log(f"   typeface '{fam}' is not installed here - using DejaVu Sans (install it, or change the style card)")
        used = "DejaVu Sans"
    plt.rcParams.update({"font.family": used, "font.size": float(style["type"]["size"]),
                         "svg.fonttype": "none", "hatch.linewidth": float(style["lines"]["hairline"]),
                         "axes.linewidth": 0, "savefig.facecolor": style["palette"]["background"]})
    return used


def TXT(style, s):
    case = style["type"].get("case", "upper")
    return s.upper() if case == "upper" else s.lower() if case == "lower" else s


def prompt_suffix(style):
    p, ln, h, t = style["palette"], style["lines"], style["hatch"], style["type"]
    refs = [v for v in style["references"].values() if v]
    parts = [f"Style: {style['prompt'].get('medium', '')}.",
             f"Palette only: background {p['background']}, ink {p['ink']}, buildings {p['buildings']}, "
             f"site accent {p['site']}, water {p['water']}, green {p['green']}.",
             f"Line language: hairline {ln['hairline']} pt, fine {ln['fine']} pt, heavy {ln['heavy']} pt; no gradients, no shadows unless the reference shows them.",
             f"Hatches: heritage '{h['heritage']}', flood '{h['flood']}', bushfire '{h['bushfire']}'.",
             f"Typeface {t['family']}, {t.get('case', 'upper')} case, labels only.",
             "No photorealism. Do not invent buildings, roads or labels that are not in the source drawing."]
    if refs:
        parts.append("Match the graphic register of these references: " + "; ".join(refs) + ".")
    if style["prompt"].get("suffix"):
        parts.append(style["prompt"]["suffix"].strip())
    return " ".join(parts)


# --------------------------------------------------------------------------------------
# drawing primitives (work in plan or through an axonometric transform)
# --------------------------------------------------------------------------------------
def _ring_path(coords, T):
    xs, ys = T(np.array([c[0] for c in coords]), np.array([c[1] for c in coords]))
    return list(zip(xs, ys))


def draw_poly(ax, geom, T, fc="none", ec="none", lw=0.5, hatch=None, alpha=1.0, z=1, ls="-"):
    for part in iter_parts(geom):
        if part.geom_type != "Polygon":
            continue
        verts, codes = [], []
        for ring in [part.exterior] + list(part.interiors):
            pts = _ring_path(ring.coords, T)
            if len(pts) < 3:
                continue
            verts += pts
            codes += [MplPath.MOVETO] + [MplPath.LINETO] * (len(pts) - 2) + [MplPath.CLOSEPOLY]
        if not verts:
            continue
        patch = PathPatch(MplPath(verts, codes), facecolor=fc, edgecolor=ec, lw=lw, hatch=hatch, alpha=alpha,
                          zorder=z, linestyle=ls, joinstyle="round")
        ax.add_patch(patch)


def draw_line(ax, geom, T, color, lw=0.5, z=1, ls="-", alpha=1.0):
    for part in iter_parts(geom):
        if part.geom_type == "Polygon":
            for ring in [part.exterior] + list(part.interiors):
                x, y = T(*np.array(ring.coords).T[:2])
                ax.plot(x, y, color=color, lw=lw, zorder=z, ls=ls, alpha=alpha, solid_capstyle="round")
        elif part.geom_type in ("LineString", "LinearRing"):
            x, y = T(*np.array(part.coords).T[:2])
            ax.plot(x, y, color=color, lw=lw, zorder=z, ls=ls, alpha=alpha, solid_capstyle="round")


def draw_points(ax, geoms, T, color, r=2.5, z=1, alpha=1.0):
    for g in geoms:
        circ = Point(g.x, g.y).buffer(r, 8)
        draw_poly(ax, circ, T, fc=color, ec="none", z=z, alpha=alpha)


def ident(x, y):
    return np.asarray(x), np.asarray(y)


def cat_colors(values, style):
    pal = style["palette"]["categorical"]
    vals = sorted({v for v in values if v not in (None, "")}, key=str)
    return {v: pal[i % len(pal)] for i, v in enumerate(vals)}


def seq_color(style, t):
    import matplotlib.colors as mc
    a = np.array(mc.to_rgb(style["palette"]["sequential_low"]))
    b = np.array(mc.to_rgb(style["palette"]["sequential_high"]))
    return mc.to_hex(a + (b - a) * max(0.0, min(1.0, t)))


def _num(v):
    try:
        return float(v)
    except Exception:
        return None


def stat_fill(c, key, style):
    """{feature index: (fill colour, legend label)} for a statutory layer."""
    L = c.L[key]
    out, legend = {}, {}
    if key in ("zone", "ass"):
        cmap = cat_colors([p.get("_label") for _, p in L.features], style)
        for i, (g, p) in enumerate(L.features):
            col = cmap.get(p.get("_label"), style["palette"]["categorical"][-1])
            out[i] = col
            legend[p.get("_label") or "?"] = (col, p.get("_desc", ""))
    elif key in ("hob", "fsr"):
        f = "MAX_B_H" if key == "hob" else "FSR"
        vals = [_num(p.get(f)) for _, p in L.features]
        good = [v for v in vals if v is not None]
        lo, hi = (min(good), max(good)) if good else (0, 1)
        for i, v in enumerate(vals):
            t = 0.5 if v is None or hi == lo else (v - lo) / (hi - lo)
            out[i] = seq_color(style, t)
            lab = L.features[i][1].get("_label") or "?"
            legend[lab] = (out[i], "")
        legend = dict(sorted(legend.items(), key=lambda kv: _num(re.sub(r"[^0-9.]", "", kv[0]) or 0) or 0))
    return out, legend


# --------------------------------------------------------------------------------------
# PLAN (base + statutory overlay)
# --------------------------------------------------------------------------------------
def draw_base(ax, c, style, T=ident, z0=0, faint=False):
    p, ln = style["palette"], style["lines"]
    a = 0.45 if faint else 1.0
    for g, pr in c.L["landuse"].features:
        draw_poly(ax, g, T, fc=p["green"] if pr.get("green") else p["landuse"], z=z0 + 1, alpha=a)
    for g, pr in c.L["water"].features:
        if g.geom_type in ("Polygon", "MultiPolygon"):
            draw_poly(ax, g, T, fc=p["water"], z=z0 + 2, alpha=a)
        else:
            draw_line(ax, g, T, p["water"], lw=ln["medium"] * 2, z=z0 + 2, alpha=a)
    for g, pr in c.L["trees"].features:
        if g.geom_type == "Point":
            draw_points(ax, [g], T, p["trees"], r=(fnum(pr.get("crown_m")) or 5) / 2, z=z0 + 3, alpha=0.8 * a)
        elif g.geom_type in ("Polygon", "MultiPolygon"):
            draw_poly(ax, g, T, fc=p["trees"], z=z0 + 3, alpha=0.5 * a)
        else:
            draw_line(ax, g, T, p["trees"], lw=ln["medium"] * 2, z=z0 + 3, alpha=a)
    if not faint:
        for g, pr in c.L["contours"].features:
            draw_line(ax, g, T, p["contours"], lw=ln["hairline"], z=z0 + 4)
        for g, pr in c.L["cadastre"].features:
            draw_line(ax, g, T, p["cadastre"], lw=ln["hairline"] * 0.8, z=z0 + 5)
    for g, pr in c.L["roads"].features:
        w = ln["medium"] if pr.get("kind") in ("primary", "secondary", "trunk", "motorway", "tertiary") else ln["fine"]
        draw_line(ax, g, T, p["roads"], lw=w, z=z0 + 6, ls="--" if pr.get("tunnel") == "yes" else "-", alpha=a)
    for g, pr in c.L["paths"].features:
        draw_line(ax, g, T, p["paths"], lw=ln["hairline"], z=z0 + 6, ls=(0, (1, 1.5)), alpha=a)
    for g, pr in c.L["rail"].features:
        draw_line(ax, g, T, p["rail"], lw=ln["medium"], z=z0 + 7, ls=(0, (4, 2)) if pr.get("tunnel") == "yes" else "-", alpha=a)
    for g, pr in c.L["buildings"].features:
        draw_poly(ax, g, T, fc=p["buildings"], ec=p["buildings_edge"], lw=ln["hairline"], z=z0 + 8, alpha=a)


def draw_site(ax, c, style, T=ident, z=50, lw=None):
    draw_poly(ax, c.site_geom, T, fc="none", ec=style["palette"]["site"], lw=lw or style["lines"]["heavy"], z=z)


def frame_axes(ax, c, style, pad=1.06):
    r = c.r
    ax.set_xlim(-r * pad, r * pad)
    ax.set_ylim(-r * pad, r * pad)
    ax.set_aspect("equal")
    ax.axis("off")
    ax.set_facecolor(style["palette"]["background"])


def north_and_scale(ax, c, style, size=1.0):
    r, ink = c.r, style["palette"]["ink"]
    x, y, s = r * 0.93, r * 0.86, r * 0.055 * size
    ax.add_patch(MplPolygon([(x, y + s), (x - s * 0.45, y - s * 0.8), (x, y - s * 0.4)], closed=True, fc=ink, ec=ink, lw=0.3, zorder=90))
    ax.add_patch(MplPolygon([(x, y + s), (x + s * 0.45, y - s * 0.8), (x, y - s * 0.4)], closed=True, fc="none", ec=ink, lw=0.3, zorder=90))
    ax.text(x, y + s * 1.25, "N", ha="center", va="bottom", color=ink, fontsize=style["type"]["size"], zorder=90)
    nice = [10, 20, 25, 50, 100, 200, 250, 500]
    L = min(nice, key=lambda v: abs(v - r / 3))
    x0, y0 = -r * 0.98, -r * 0.98
    for i in range(4):
        ax.add_patch(plt.Rectangle((x0 + i * L / 4, y0), L / 4, r * 0.012, fc=ink if i % 2 == 0 else "none", ec=ink, lw=0.3, zorder=90))
    ax.text(x0, y0 + r * 0.025, "0", fontsize=style["type"]["size"] * 0.8, color=ink, zorder=90)
    ax.text(x0 + L, y0 + r * 0.025, f"{L} m", fontsize=style["type"]["size"] * 0.8, color=ink, ha="center", zorder=90)


def title_block(fig, c, style, title, sub=""):
    ink = style["palette"]["ink"]
    fig.text(0.04, 0.975, TXT(style, f"{style['meta'].get('title_prefix', 'SITE LAYERS-OPEN')} - {title}"), fontsize=style["type"]["title_size"],
             color=ink, va="top", fontweight="bold")
    s2 = f"{c.address}  |  radius {c.r} m  |  {TODAY}" + (f"  |  {sub}" if sub else "")
    fig.text(0.04, 0.948, s2, fontsize=style["type"]["size"], color=ink, va="top")
    who = " / ".join(v for v in (style["meta"].get("author"), style["meta"].get("studio")) if v)
    fig.text(0.04, 0.012, ("Data: (c) OpenStreetMap contributors (ODbL); NSW Spatial Services & NSW Planning (CC BY 4.0); "
                           "Open-Meteo (CC BY 4.0). Unverified until the register says otherwise." + (f"   {who}" if who else "")),
             fontsize=style["type"]["size"] * 0.75, color=ink, alpha=0.7)


def save_fig(fig, stem: Path, dpi=200):
    fig.savefig(stem.with_suffix(".svg"))
    fig.savefig(stem.with_suffix(".png"), dpi=dpi)
    plt.close(fig)
    return [stem.with_suffix(".svg"), stem.with_suffix(".png")]


def render_plan(c, style, out: Path):
    p, ln = style["palette"], style["lines"]
    fig = plt.figure(figsize=(8.3, 9.2), facecolor=p["background"])
    ax = fig.add_axes([0.03, 0.12, 0.94, 0.8])
    frame_axes(ax, c, style)
    draw_base(ax, c, style)
    zf, zleg = stat_fill(c, "zone", style)
    for i, (g, pr) in enumerate(c.L["zone"].features):
        draw_poly(ax, g, ident, fc=zf.get(i, "none"), alpha=0.35, z=9)
        draw_line(ax, g, ident, p["ink"], lw=ln["fine"], z=20, alpha=0.6)
        lp = label_point(g)
        ax.text(lp.x, lp.y, pr.get("_label", ""), fontsize=style["type"]["size"], color=p["ink"], ha="center", va="center", zorder=30,
                fontweight="bold")
    hat = style["hatch"]
    for key, col, h in (("heritage", p["heritage"], hat["heritage"]), ("shr", p["heritage"], hat["shr"]),
                        ("flood", p["flood"], hat["flood"]), ("bushfire", p["bushfire"], hat["bushfire"])):
        for g, pr in c.L[key].features:
            draw_poly(ax, g, ident, fc="none", ec=col, lw=ln["hairline"], hatch=h, z=21)
            draw_line(ax, g, ident, col, lw=ln["fine"], z=22)
    draw_site(ax, c, style)
    ax.add_patch(MplCircle((0, 0), c.r, fc="none", ec=p["ink"], lw=ln["fine"], ls=(0, (6, 4)), zorder=40))
    north_and_scale(ax, c, style)
    # legend
    lax = fig.add_axes([0.04, 0.035, 0.92, 0.075])
    lax.axis("off")
    items = [(p["site"], "site", None, True), (p["buildings"], "buildings (OSM)", None, False)]
    items += [(col, f"{k}", None, False) for k, (col, d) in list(zleg.items())[:14]]
    for key, col, lab, h in (("heritage", p["heritage"], "heritage", hat["heritage"]), ("shr", p["heritage"], "state heritage register", hat["shr"]),
                             ("flood", p["flood"], "flood planning", hat["flood"]), ("bushfire", p["bushfire"], "bush fire prone", hat["bushfire"])):
        if len(c.L[key]):
            items.append((col, lab, h, False))
    rows = max(1, (len(items) + 6) // 7)
    for i, (col, lab, h, outline) in enumerate(items):
        x = (i % 7) / 7.0
        y = 0.85 - (i // 7) * (0.8 / max(rows - 1, 1) if rows > 1 else 0)
        lax.add_patch(plt.Rectangle((x, y - 0.1), 0.02, 0.2, fc="none" if (h or outline) else col, ec=col,
                                    hatch=h, lw=1.2 if outline else 0.4, transform=lax.transAxes))
        lax.text(x + 0.032, y, TXT(style, lab), fontsize=style["type"]["size"] * 0.9, va="center", transform=lax.transAxes, color=p["ink"])
    status = [k for k in ("zone", "heritage", "flood", "bushfire") if c.L[k].status in ("unavailable", "failed")]
    title_block(fig, c, style, "SITE AND STATUTORY PLAN",
                "statutory layers not available here: " + ", ".join(status) if status else "")
    return save_fig(fig, out / "diagram_plan")


def render_statutory(c, style, out: Path):
    p, ln = style["palette"], style["lines"]
    panels = [("zone", "Zoning"), ("hob", "Height of buildings"), ("fsr", "Floor space ratio"), ("heritage", "Heritage"),
              ("flood", "Flood planning"), ("bushfire", "Bush fire prone"), ("ass", "Acid sulfate soils"), ("_ctx", "Context")]
    fig = plt.figure(figsize=(14, 8.2), facecolor=p["background"])
    hat = style["hatch"]
    for i, (key, name) in enumerate(panels):
        ax = fig.add_axes([0.02 + (i % 4) * 0.245, 0.49 - (i // 4) * 0.44, 0.23, 0.39])
        frame_axes(ax, c, style, pad=1.02)
        ax.add_patch(MplCircle((0, 0), c.r, fc=p["plate"], ec="none", zorder=0))
        if key == "_ctx":
            draw_base(ax, c, style)
            status = f"{len(c.L['buildings'])} buildings"
        else:
            draw_base(ax, c, style, faint=True)
            L = c.L[key]
            fills, leg = stat_fill(c, key, style)
            for j, (g, pr) in enumerate(L.features):
                if key in ("zone", "ass", "hob", "fsr"):
                    draw_poly(ax, g, ident, fc=fills.get(j, "none"), alpha=0.75, z=10)
                    draw_line(ax, g, ident, p["ink"], lw=ln["hairline"], z=11, alpha=0.7)
                    lp = label_point(g)
                    if g.area > (c.r ** 2) * 0.004:
                        ax.text(lp.x, lp.y, pr.get("_label", ""), fontsize=style["type"]["size"] * 0.8, ha="center", va="center",
                                color=p["ink"], zorder=12)
                else:
                    col = {"heritage": p["heritage"], "flood": p["flood"], "bushfire": p["bushfire"]}[key]
                    draw_poly(ax, g, ident, fc="none", ec=col, hatch=hat.get(key, "//"), lw=ln["hairline"], z=10)
                    draw_line(ax, g, ident, col, lw=ln["fine"], z=11)
            if key in ("hob", "fsr"):
                for g, pr in c.L["addctrl"].features:
                    draw_line(ax, g, ident, p["ink"], lw=ln["fine"], z=13, ls=(0, (3, 2)))
            if key == "heritage":
                for g, pr in c.L["shr"].features:
                    draw_poly(ax, g, ident, fc="none", ec=p["heritage"], hatch=hat["shr"], lw=ln["hairline"], z=10)
            status = {"ok": f"{len(L)} features", "empty": "nothing mapped within radius",
                      "unavailable": "not available for this jurisdiction", "failed": "service failed - see register",
                      "skipped": "not run"}[L.status]
            if L.note and L.status == "empty":
                status = L.note
        draw_site(ax, c, style, lw=ln["medium"])
        ax.add_patch(MplCircle((0, 0), c.r, fc="none", ec=p["ink"], lw=ln["hairline"], zorder=40))
        ax.set_title(TXT(style, name), fontsize=style["type"]["size"] * 1.2, color=p["ink"], loc="left", pad=4, fontweight="bold")
        ax.text(-c.r, -c.r * 1.08, status, fontsize=style["type"]["size"] * 0.85, color=p["ink"], alpha=0.8)
    title_block(fig, c, style, "STATUTORY LAYERS")
    return save_fig(fig, out / "diagram_statutory")


# --------------------------------------------------------------------------------------
# EXPLODED LAYER AXONOMETRIC - the hero image. 2.5D: every layer on its own plate,
# plan rotated 45 deg and foreshortened, plates stacked; buildings extruded where
# a height is known (true scale); the site carried through every plate.
# --------------------------------------------------------------------------------------
AXO_ROT = math.radians(40)
AXO_K = 0.5


def axo_T(dz):
    ca, sa = math.cos(AXO_ROT), math.sin(AXO_ROT)

    def T(x, y):
        x = np.asarray(x, float)
        y = np.asarray(y, float)
        xr = x * ca - y * sa
        yr = x * sa + y * ca
        return xr, yr * AXO_K + dz
    return T


def _depth(x, y):
    return x * math.sin(AXO_ROT) + y * math.cos(AXO_ROT)


def draw_extruded(ax, c, style, dz, z0, hscale=1.0, cap=None):
    p, ln = style["palette"], style["lines"]
    T = axo_T(dz)
    import matplotlib.colors as mc
    blds = []
    site = c.site_geom
    for g, pr in c.L["buildings"].features:
        h = pr.get("height_m")
        cen = g.centroid
        on_site = g.intersects(site) and g.intersection(site).area > 0.3 * g.area
        blds.append((_depth(cen.x, cen.y), g, (min(h, cap) if (h and cap) else h), on_site))
    blds.sort(key=lambda t: -t[0])  # far first
    tint = lambda col: mc.to_hex(0.45 * np.array(mc.to_rgb(col)) + 0.55 * np.array(mc.to_rgb(p["site"])))
    z = z0
    for _, g, h, on_site in blds:
        z += 0.001
        side = tint(p["massing_side"]) if on_site else p["massing_side"]
        top = tint(p["massing_top"]) if on_site else p["massing_top"]
        if not h:
            draw_poly(ax, g, T, fc=tint(p["buildings"]) if on_site else p["buildings"], ec=p["buildings_edge"],
                      lw=ln["hairline"], z=z)
            continue
        H = h * hscale
        Tt = axo_T(dz + H)
        for part in iter_parts(g):
            if part.geom_type != "Polygon":
                continue
            ext = list(part.exterior.coords)
            walls = []
            for (x1, y1), (x2, y2) in zip(ext[:-1], ext[1:]):
                walls.append((_depth((x1 + x2) / 2, (y1 + y2) / 2), (x1, y1), (x2, y2)))
            walls.sort(key=lambda t: -t[0])
            for _, (x1, y1), (x2, y2) in walls:
                bx, by = T([x1, x2], [y1, y2])
                tx, ty = Tt([x1, x2], [y1, y2])
                ax.add_patch(MplPolygon([(bx[0], by[0]), (bx[1], by[1]), (tx[1], ty[1]), (tx[0], ty[0])], closed=True,
                                        fc=side, ec=p["buildings_edge"], lw=ln["hairline"] * 0.6, zorder=z))
            draw_poly(ax, part, Tt, fc=top, ec=p["buildings_edge"], lw=ln["hairline"], z=z + 0.0005)


def render_exploded(c, style, out: Path):
    p, ln = style["palette"], style["lines"]
    r = c.r
    gap = r * 0.55
    hat = style["hatch"]
    hmax = max([pr.get("height_m") or 0 for _, pr in c.L["buildings"].features] + [0])
    cap = r * 0.9
    bgap = min(hmax, cap) + 2 * r * AXO_K * 0.95

    def plate_status(keys):
        st = [c.L[k].status for k in keys]
        if all(s == "unavailable" for s in st):
            return "not available for this jurisdiction"
        if all(s in ("failed",) for s in st):
            return "service failed - see register"
        if all(s in ("empty", "unavailable", "failed", "skipped") for s in st):
            return "nothing mapped within radius"
        return ""

    def P_terrain(ax, T, z):
        for g, pr in c.L["contours"].features:
            draw_line(ax, g, T, p["contours"], lw=ln["fine"], z=z)
        return "TERRAIN", (c.L["contours"].note or ""), plate_status(["contours"])

    def P_land(ax, T, z):
        for g, pr in c.L["cadastre"].features:
            draw_line(ax, g, T, p["cadastre"], lw=ln["hairline"], z=z)
        return "LAND", "lots and site", plate_status(["cadastre"])

    def P_ground(ax, T, z):
        for g, pr in c.L["landuse"].features:
            draw_poly(ax, g, T, fc=p["green"] if pr.get("green") else p["landuse"], z=z)
        for g, pr in c.L["water"].features:
            if g.geom_type in ("Polygon", "MultiPolygon"):
                draw_poly(ax, g, T, fc=p["water"], z=z + 0.1)
            else:
                draw_line(ax, g, T, p["water"], lw=ln["medium"] * 2, z=z + 0.1)
        for g, pr in c.L["trees"].features:
            if g.geom_type == "Point":
                draw_points(ax, [g], T, p["trees"], r=(fnum(pr.get("crown_m")) or 5) / 2, z=z + 0.2)
            elif g.geom_type in ("Polygon", "MultiPolygon"):
                draw_poly(ax, g, T, fc=p["trees"], alpha=0.6, z=z + 0.2)
            else:
                draw_line(ax, g, T, p["trees"], lw=ln["medium"] * 2, z=z + 0.2)
        return "GROUND", "open space, trees, water", plate_status(["landuse", "water", "trees"])

    def P_move(ax, T, z):
        for g, pr in c.L["roads"].features:
            w = ln["medium"] * 1.2 if pr.get("kind") in ("primary", "secondary", "trunk", "motorway", "tertiary") else ln["fine"]
            draw_line(ax, g, T, p["roads"], lw=w, z=z, ls="--" if pr.get("tunnel") == "yes" else "-")
        for g, pr in c.L["paths"].features:
            draw_line(ax, g, T, p["paths"], lw=ln["hairline"], z=z, ls=(0, (1, 1.5)))
        for g, pr in c.L["rail"].features:
            draw_line(ax, g, T, p["rail"], lw=ln["medium"], z=z + 0.1, ls=(0, (4, 2)) if pr.get("tunnel") == "yes" else "-")
        return "MOVEMENT", "roads, paths, rail", plate_status(["roads", "paths", "rail"])

    def P_stat(key, title, sub):
        def f(ax, T, z):
            fills, leg = stat_fill(c, key, style)
            for j, (g, pr) in enumerate(c.L[key].features):
                if key in ("zone", "hob", "fsr", "ass"):
                    draw_poly(ax, g, T, fc=fills.get(j, "none"), ec=p["ink"], lw=ln["hairline"], alpha=0.85, z=z)
                    if g.area > r * r * 0.01:
                        lp = label_point(g)
                        X, Y = T([lp.x], [lp.y])
                        ax.text(X[0], Y[0], pr.get("_label", ""), fontsize=style["type"]["size"] * 0.75, ha="center",
                                va="center", color=p["ink"], zorder=z + 0.3)
            st = plate_status([key])
            return title, sub, st
        return f

    def P_heritage(ax, T, z):
        for key, h in (("heritage", hat["heritage"]), ("shr", hat["shr"])):
            for g, pr in c.L[key].features:
                draw_poly(ax, g, T, fc="none", ec=p["heritage"], hatch=h, lw=ln["hairline"], z=z)
                draw_line(ax, g, T, p["heritage"], lw=ln["fine"], z=z + 0.1)
        return "HERITAGE", "items, conservation areas, SHR", plate_status(["heritage", "shr"])

    def P_hazard(ax, T, z):
        for key, col, h in (("ass", p["acid_sulfate"], hat["acid_sulfate"]), ("flood", p["flood"], hat["flood"]),
                            ("bushfire", p["bushfire"], hat["bushfire"])):
            for g, pr in c.L[key].features:
                draw_poly(ax, g, T, fc="none", ec=col, hatch=h, lw=ln["hairline"], z=z)
                draw_line(ax, g, T, col, lw=ln["fine"], z=z + 0.1)
        st = plate_status(["flood", "bushfire", "ass"])
        if c.L["flood"].note and c.L["flood"].status == "empty":
            st = (st + " - " if st else "") + "flood: " + c.L["flood"].note
        return "HAZARD", "flood, bush fire, acid sulfate", st

    def P_climate(ax, T, z):
        if getattr(c, "sun", None):
            R = r * 0.8
            for alt in (0, 30, 60):
                draw_line(ax, Point(0, 0).buffer((90 - alt) / 90 * R, 64).exterior, T, p["sun"], lw=ln["hairline"], z=z)
            for name, doy in key_days(c.pin["lat"]):
                pts = [polar_xy(alt, az, R) for t, alt, az in sun_curve(c.pin["lat"], doy, 0.1)]
                if len(pts) > 1:
                    draw_line(ax, LineString(pts), T, p["sun"], lw=ln["medium"], z=z + 0.1)
        if getattr(c, "climate", None):
            tot = c.climate["annual"][0].sum(axis=1)
            Rw = r * 0.55 / max(tot.max(), 1e-6)
            for i in range(16):
                a0, a1 = math.radians(i * 22.5 - 9), math.radians(i * 22.5 + 9)
                rr = tot[i] * Rw
                draw_poly(ax, Polygon([(0, 0), (rr * math.sin(a0), rr * math.cos(a0)), (rr * math.sin(a1), rr * math.cos(a1))]),
                          T, fc=p["wind"], alpha=0.55, ec=p["wind"], lw=ln["hairline"], z=z + 0.2)
        st = "" if (getattr(c, "sun", None) or getattr(c, "climate", None)) else "climate step failed"
        return "CLIMATE", "sun path, annual wind rose", st

    plates = [P_terrain, P_land, P_ground, P_move, "BUILT",
              P_stat("zone", "ZONING", "land use zones"), P_stat("hob", "HEIGHT", "maximum building height"),
              P_stat("fsr", "DENSITY", "floor space ratio"), P_heritage, P_hazard, P_climate]
    # vertical positions
    dzs, dz = [], 0.0
    for i, pl in enumerate(plates):
        dzs.append(dz)
        dz += bgap if pl == "BUILT" else gap
    top = dzs[-1] + r * AXO_K
    fig_h = 2.0 + 13.0 * (top + r * AXO_K) / (r * 5.5)
    fig = plt.figure(figsize=(10, max(12, fig_h)), facecolor=p["background"])
    ax = fig.add_axes([0.0, 0.03, 1.0, 0.9])
    ax.set_aspect("equal")
    ax.axis("off")
    labx = -r * 2.25
    scen = c.site_geom.centroid
    for i, pl in enumerate(plates):
        z = 10 + i * 10
        T = axo_T(dzs[i])
        disc = Point(0, 0).buffer(r, 96)
        draw_poly(ax, disc, T, fc=p["plate"], ec=p["ink"], lw=ln["fine"], alpha=0.93, z=z)
        if pl == "BUILT":
            draw_extruded(ax, c, style, dzs[i], z + 1, cap=cap)
            title, sub = "BUILT FORM", "buildings; extruded where OSM has a height"
            n = len(c.L["buildings"])
            nh = sum(1 for _, pr in c.L["buildings"].features if pr.get("height_m"))
            sub += f" ({nh}/{n})"
            st = plate_status(["buildings"])
        else:
            title, sub, st = pl(ax, T, z + 1)
        draw_poly(ax, c.site_geom, T, fc=p["site"], alpha=0.25, ec=p["site"], lw=ln["medium"], z=(z + 0.5) if pl == "BUILT" else (z + 5))
        if st:
            X, Y = T([0], [0])
            ax.text(X[0], Y[0], TXT(style, st), ha="center", va="center", fontsize=style["type"]["size"],
                    color=p["ink"], alpha=0.7, zorder=z + 6)
        # label and leader
        ly = [dzs[i]]
        lx = [-r * 1.02]
        ax.plot([labx + r * 0.9, lx[0]], [ly[0], ly[0]], color=p["ink"], lw=ln["hairline"], zorder=z + 6)
        ax.text(labx, ly[0] + r * 0.03, TXT(style, title), fontsize=style["type"]["size"] * 1.35, fontweight="bold",
                color=p["ink"], va="bottom", zorder=z + 6)
        ax.text(labx, ly[0] - r * 0.02, sub, fontsize=style["type"]["size"] * 0.95, color=p["ink"], va="top", zorder=z + 6)
        if i < len(plates) - 1:
            T2 = axo_T(dzs[i + 1])
            X1, Y1 = T([scen.x], [scen.y])
            X2, Y2 = T2([scen.x], [scen.y])
            ax.plot([X1[0], X2[0]], [Y1[0], Y2[0]], color=p["site"], lw=ln["fine"], ls=(0, (3, 2)), zorder=z + 7)
    # north arrow for the rotated plan
    T0 = axo_T(0)
    dx, dy = T0([0, 0], [0, 1])
    vx, vy = dx[1] - dx[0], dy[1] - dy[0]
    nrm = math.hypot(vx, vy)
    vx, vy = vx / nrm * r * 0.18, vy / nrm * r * 0.18
    bx, by = r * 0.95, -r * AXO_K - r * 0.02
    ax.annotate("", xy=(bx + vx, by + vy), xytext=(bx, by), arrowprops=dict(arrowstyle="-|>", color=p["ink"], lw=0.8), zorder=500)
    ax.text(bx + vx * 1.15, by + vy * 1.15, "N", color=p["ink"], fontsize=style["type"]["size"], ha="center", va="center", zorder=500)
    ax.set_xlim(labx - r * 0.05, r * 1.2)
    ax.set_ylim(-r * AXO_K - r * 0.15, top + r * 0.1)
    title_block(fig, c, style, "EXPLODED SITE LAYERS", "heights true scale" + (f", capped at {cap:.0f} m" if hmax > cap else ""))
    return save_fig(fig, out / "diagram_exploded", dpi=180)


# --------------------------------------------------------------------------------------
# CLIMATE SHEET - sun path, annual and seasonal wind roses, monthly temperature and rain
# --------------------------------------------------------------------------------------
def _rose(ax, tab, style, title, big=False):
    p = style["palette"]
    import matplotlib.colors as mc
    theta = np.radians(np.arange(16) * 22.5)
    width = np.radians(22.5) * 0.9
    base = np.zeros(16)
    a = np.array(mc.to_rgb(p["plate"]))
    b = np.array(mc.to_rgb(p["wind"]))
    for i in range(tab.shape[1]):
        col = mc.to_hex(a + (b - a) * (0.3 + 0.7 * i / (tab.shape[1] - 1)))
        ax.bar(theta, tab[:, i], width=width, bottom=base, color=col, edgecolor=p["ink"], lw=0.15,
               label=f"{SPEED_LABELS[i]} m/s" if big else None)
        base += tab[:, i]
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    ax.set_xticks(np.radians([0, 90, 180, 270]))
    ax.set_xticklabels(["N", "E", "S", "W"], fontsize=style["type"]["size"] * 0.9)
    ax.tick_params(axis="y", labelsize=style["type"]["size"] * 0.6)
    ax.set_title(TXT(style, title), fontsize=style["type"]["size"] * (1.2 if big else 0.95), color=p["ink"], pad=8)
    ax.set_facecolor(p["plate"])
    ax.grid(color=p["ink"], alpha=0.2, lw=0.3)
    if big:
        ax.legend(loc="lower left", bbox_to_anchor=(-0.25, -0.2), fontsize=style["type"]["size"] * 0.75, frameon=False)


def render_climate(c, style, out: Path):
    p, ln = style["palette"], style["lines"]
    lat = c.pin["lat"]
    fig = plt.figure(figsize=(14, 9.6), facecolor=p["background"])
    # sun path
    ax = fig.add_axes([0.04, 0.42, 0.30, 0.44], projection="polar")
    ax.set_theta_zero_location("N")
    ax.set_theta_direction(-1)
    ax.set_ylim(0, 90)
    ax.set_yticks([0, 30, 60, 90])
    ax.set_yticklabels(["90", "60", "30", "0"], fontsize=style["type"]["size"] * 0.7)
    ax.set_xticks(np.radians(np.arange(0, 360, 45)))
    ax.set_xticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW"], fontsize=style["type"]["size"] * 0.9)
    ax.set_facecolor(p["plate"])
    ax.grid(color=p["ink"], alpha=0.2, lw=0.3)
    for (name, doy), ls in zip(key_days(lat), ["-", "--", "-"]):
        pts = sun_curve(lat, doy, 0.05)
        if pts:
            ax.plot([math.radians(a) for _, _, a in pts], [90 - al for _, al, _ in pts], color=p["sun"],
                    lw=ln["heavy"] if "Summer" in name else ln["medium"], ls=ls, label=name)
    for hr in range(5, 20):
        pts = [(sun_pos(lat, d, hr)) for d in range(1, 366, 4)]
        pts = [(math.radians(az), 90 - al) for al, az in pts if al > 0]
        if len(pts) > 2:
            ax.plot(*zip(*sorted(pts, key=lambda t: t[1])), ".", color=p["ink"], ms=0.8, alpha=0.5)
            th, rr = min(pts, key=lambda t: t[1])
            ax.text(th, rr, f"{hr}", fontsize=style["type"]["size"] * 0.7, color=p["ink"], ha="center", va="bottom")
    ax.legend(loc="upper left", bbox_to_anchor=(-0.08, 1.02), fontsize=style["type"]["size"] * 0.75, frameon=False)
    ax.set_title(TXT(style, "Sun path (solar time, true north)"), fontsize=style["type"]["size"] * 1.2, pad=10, color=p["ink"])

    if getattr(c, "climate", None):
        cl = c.climate
        ax2 = fig.add_axes([0.40, 0.42, 0.26, 0.44], projection="polar")
        y1, y2 = None, None
        _rose(ax2, cl["annual"][0], style, f"Wind, all hours (calm {cl['annual'][1]:.0f}%)", big=True)
        ax3 = fig.add_axes([0.72, 0.52, 0.26, 0.33])
        m = np.arange(12)
        ax3.bar(m, cl["rain"], color=p["water"], edgecolor="none", label="rain (mm)")
        ax3.set_ylabel("rain mm / month", fontsize=style["type"]["size"] * 0.8, color=p["ink"])
        ax3b = ax3.twinx()
        ax3b.plot(m, cl["tmax"], color=p["site"], lw=ln["medium"], marker="o", ms=2.5, label="mean daily max")
        ax3b.plot(m, cl["tmin"], color=p["wind"], lw=ln["medium"], marker="o", ms=2.5, label="mean daily min")
        ax3b.set_ylabel("deg C", fontsize=style["type"]["size"] * 0.8, color=p["ink"])
        ax3.set_xticks(m)
        ax3.set_xticklabels([s[0] for s in MONTHS], fontsize=style["type"]["size"] * 0.8)
        for a_ in (ax3, ax3b):
            a_.tick_params(labelsize=style["type"]["size"] * 0.75)
            for s in a_.spines.values():
                s.set_visible(False)
        ax3.set_facecolor(p["plate"])
        ax3b.legend(loc="upper center", bbox_to_anchor=(0.5, -0.1), ncol=2, fontsize=style["type"]["size"] * 0.75, frameon=False)
        ax3.set_title(TXT(style, "Monthly temperature and rain"), fontsize=style["type"]["size"] * 1.2, color=p["ink"])
        for i, (name, (tab, calm, n)) in enumerate(cl["seasons"].items()):
            axs = fig.add_axes([0.05 + i * 0.235, 0.07, 0.17, 0.24], projection="polar")
            _rose(axs, tab, style, name)
        gl = cl["grid"]
        sub = f"wind and weather: Open-Meteo ERA5 reanalysis grid point {gl[0]:.3f}, {gl[1]:.3f}"
    else:
        fig.text(0.55, 0.6, "Climate data did not download - see the register.", fontsize=12, color=p["ink"])
        sub = "climate step failed"
    title_block(fig, c, style, "SUN AND WIND", sub)
    return save_fig(fig, out / "diagram_climate")


# --------------------------------------------------------------------------------------
# DATA OUTPUTS
# --------------------------------------------------------------------------------------
def write_geojson(c, folder: Path):
    folder.mkdir(parents=True, exist_ok=True)
    written = []
    for key, L in c.L.items():
        if not len(L):
            continue
        feats = []
        for g, p in L.features:
            props = {k: (v if isinstance(v, (str, int, float, bool)) or v is None else str(v))
                     for k, v in p.items() if not k.startswith("_") or k in ("_label", "_desc", "_epi", "_currency")}
            props["layer"] = L.dxf
            try:
                gl = c.frame.to_lonlat(g)
                feats.append({"type": "Feature", "geometry": mapping(gl), "properties": props})
            except Exception:
                continue
        fc = {"type": "FeatureCollection", "name": L.dxf,
              "metadata": {"source": L.source, "source_url": L.url, "retrieved": L.retrieved, "status": L.status,
                           "note": L.note, "tool": f"{TOOL} v{VERSION}"},
              "features": feats}
        p = folder / f"{L.dxf}.geojson"
        p.write_text(json.dumps(fc, separators=(",", ":")), encoding="utf-8")
        written.append(p)
    return written


def write_heights(c, path: Path):
    with open(path, "w", newline="", encoding="utf-8-sig") as f:
        w = csv.writer(f)
        w.writerow(["osm_id", "name", "building", "height_m", "height_basis", "levels_tag", "height_tag",
                    "footprint_m2", "centroid_x", "centroid_y", "check"])
        for g, p in sorted(c.L["buildings"].features, key=lambda t: -(t[1].get("height_m") or 0)):
            cen = g.centroid
            w.writerow([p["osm_id"], p["name"], p["building"], fmt_m(p["height_m"]), p["height_basis"], p["levels"],
                        p["height_tag"], f"{g.area:.0f}", f"{cen.x:.1f}", f"{cen.y:.1f}", ""])


def write_about(c, path: Path, files):
    L = c.L
    lines = [f"# {TOOL} - {c.address}", "",
             f"Run {dt.datetime.now().strftime('%Y-%m-%d %H:%M')} - v{VERSION} - radius {c.r} m - jurisdiction {c.jur}"
             + (" - OFFLINE (cached data)" if c.net.offline else ""), "",
             f"Pin: {c.pin['lat']:.7f}, {c.pin['lon']:.7f} ({c.pin.get('source', '')}). DXF origin 0,0 = pin; metres; +Y = true north.", "",
             "| layer | DXF | status | features | source |", "|---|---|---|---|---|"]
    for key, lay in L.items():
        lines.append(f"| {lay.title} | {lay.dxf} | {lay.status}{(' - ' + lay.note) if lay.note else ''} | {len(lay)} | {lay.source} |")
    lines += ["", f"Register: {len(c.reg.rows)} rows. **Every row has an empty verdict. Fill it.**", "",
              "Files:", ""] + [f"- {Path(f).name}" for f in files]
    path.write_text("\n".join(lines) + "\n", encoding="utf-8")


# --------------------------------------------------------------------------------------
# ARCHIMAP COMPARISON - "what five dollars buys"
# --------------------------------------------------------------------------------------
ARCHIMAP_PUBLISHED = {  # checked on masslabs-archi.com, 23 Sep 2026
    "price": "USD 4.99 / 9.99 / 14.90 per month (1 / 2 / 3 km diameter; 10 / 15 / 20 runs a month); cheapest tier free on promotion when checked 23 Sep 2026",
    "layers": "cadastre, building footprints, building use, roads, terrain, green space, traffic noise, wind",
    "exports": "PNG, SVG, DXF, PDF; 3DM and SKP in 3D",
}
THEMES = [
    ("Building footprints", ["build", "bldg", "footprint", "house"], ["CTX_BUILDINGS"]),
    ("Building heights / 3D massing", ["height", "3d", "mass", "volume"], ["CTX_MASSING"]),
    ("Building use", ["use", "function", "usage"], []),
    ("Roads", ["road", "street", "highway", "traffic"], ["CTX_ROADS"]),
    ("Paths / cycle", ["path", "foot", "cycle", "pedestr"], ["CTX_PATHS"]),
    ("Rail", ["rail", "train", "tram", "metro"], ["CTX_RAIL"]),
    ("Water", ["water", "river", "sea", "coast"], ["CTX_WATER"]),
    ("Green space / trees", ["green", "park", "tree", "veget", "grass"], ["CTX_TREES", "CTX_LANDUSE"]),
    ("Terrain / contours", ["contour", "terrain", "topo", "elev"], ["CTX_CONTOURS", "CTX_CONTOURS_INDEX"]),
    ("Cadastre / lots", ["cadast", "parcel", "lot", "plot", "property"], ["CTX_CADASTRE", "CTX_SITE"]),
    ("Traffic noise", ["noise"], []),
    ("Wind", ["wind"], ["AI_WIND"]),
    ("Sun path", ["sun", "solar"], ["AI_SUNPATH"]),
    ("Zoning", ["zone", "zoning"], ["AI_ZONE"]),
    ("Height control", ["hob"], ["AI_HOB"]),
    ("Floor space ratio", ["fsr", "density"], ["AI_FSR"]),
    ("Heritage", ["herit"], ["AI_HERITAGE", "AI_HERITAGE_SHR"]),
    ("Flood", ["flood"], ["AI_FLOOD"]),
    ("Bush fire", ["bush", "fire"], ["AI_BUSHFIRE"]),
    ("Acid sulfate soils", ["acid", "sulf"], ["AI_ASS"]),
]


def dxf_inventory(path):
    doc = ezdxf.readfile(str(path))
    msp = doc.modelspace()
    inv = {}
    xs, ys = [], []
    for e in msp:
        lay = e.dxf.get("layer", "0")
        d = inv.setdefault(lay, {})
        d[e.dxftype()] = d.get(e.dxftype(), 0) + 1
        try:
            if e.dxftype() == "LWPOLYLINE":
                for x, y, *_ in e.get_points():
                    xs.append(x); ys.append(y)
            elif e.dxftype() in ("LINE",):
                xs += [e.dxf.start.x, e.dxf.end.x]; ys += [e.dxf.start.y, e.dxf.end.y]
        except Exception:
            pass
    units = doc.header.get("$INSUNITS", 0)
    ext = (float(min(xs)), float(min(ys)), float(max(xs)), float(max(ys))) if xs else None
    return inv, units, ext, doc


def guess_frame(ext):
    if not ext:
        return "unknown (no 2D linework found)"
    x1, y1, x2, y2 = ext
    if -180 <= x1 <= 180 and -90 <= y1 <= 90 and (x2 - x1) < 1:
        return "looks like longitude/latitude degrees"
    if 1e5 < abs(x1) < 1e6 and 1e6 < abs(y1) < 1e7:
        return "looks like projected UTM/MGA metres"
    if max(abs(x1), abs(x2), abs(y1), abs(y2)) < 1e4:
        return "looks like a local frame near 0,0"
    return "unrecognised coordinate frame"


def compare_archimap(ours: Path, theirs: Path, out: Path, style=None, log=print):
    style = style or json.loads(json.dumps(DEFAULT_STYLE))
    inv_o, u_o, ext_o, doc_o = dxf_inventory(ours)
    inv_t, u_t, ext_t, doc_t = dxf_inventory(theirs)
    UN = {0: "unitless", 1: "inches", 4: "mm", 5: "cm", 6: "metres"}

    def has(inv, words, names):
        hits = [k for k in inv if any(w in k.lower() for w in words)] if words else []
        hits += [k for k in names if k in inv]
        n = sum(sum(inv[k].values()) for k in set(hits))
        return (", ".join(sorted(set(hits))) + f" ({n} entities)") if hits else "-"

    rows = []
    for theme, words, names in THEMES:
        rows.append((theme, has(inv_t, words, []), has(inv_o, [], names)))
    md = [f"# archiMap vs Site Layers-open - {TODAY}", "",
          f"- archiMap file: `{theirs.name}` - units {UN.get(u_t, u_t)} - extents {tuple(round(v, 1) for v in ext_t) if ext_t else '-'} - {guess_frame(ext_t)}",
          f"- Site Layers-open file: `{ours.name}` - units {UN.get(u_o, u_o)} - {guess_frame(ext_o)}", "",
          f"archiMap as published (checked 23 Sep 2026): {ARCHIMAP_PUBLISHED['price']}. Layers: {ARCHIMAP_PUBLISHED['layers']}. Exports: {ARCHIMAP_PUBLISHED['exports']}.", "",
          "| theme | archiMap DXF | Site Layers-open DXF |", "|---|---|---|"]
    md += [f"| {a} | {b} | {c_} |" for a, b, c_ in rows]
    md += ["", "## Every layer in the archiMap file", "", "| layer | entities |", "|---|---|"]
    md += [f"| {k} | {', '.join(f'{t} {n}' for t, n in v.items())} |" for k, v in sorted(inv_t.items())]
    only_t = [a for a, b, c_ in rows if b != "-" and c_ == "-"]
    only_o = [a for a, b, c_ in rows if b == "-" and c_ != "-"]
    both = [a for a, b, c_ in rows if b != "-" and c_ != "-"]
    md += ["", "## Summary", "",
           f"- In both files: {', '.join(both) or 'nothing matched'}.",
           f"- Only in the archiMap file: {', '.join(only_t) or 'nothing'}.",
           f"- Only in the Site Layers-open file: {', '.join(only_o) or 'nothing'}.",
           "- Only Site Layers-open writes a verification register naming the source of each fact.", "",
           "Layer matching is by name and is a guess: check the archiMap layer list above by eye before saying any of this in class. "
           "Compare the two drawings for quality too (see compare_archimap.png) - counts do not show accuracy."]
    (out / "compare_archimap.md").write_text("\n".join(md) + "\n", encoding="utf-8")
    # side-by-side picture, each file in its own frame
    fig, axs = plt.subplots(1, 2, figsize=(14, 7.4), facecolor=style["palette"]["background"])
    for ax, doc, name in ((axs[0], doc_t, "archiMap export"), (axs[1], doc_o, "Site Layers-open (free)")):
        for e in doc.modelspace():
            try:
                if e.dxftype() == "LWPOLYLINE":
                    pts = np.array([(x, y) for x, y, *_ in e.get_points()])
                    if e.closed:
                        pts = np.vstack([pts, pts[:1]])
                    ax.plot(pts[:, 0], pts[:, 1], lw=0.2, color=style["palette"]["ink"])
                elif e.dxftype() == "LINE":
                    ax.plot([e.dxf.start.x, e.dxf.end.x], [e.dxf.start.y, e.dxf.end.y], lw=0.2, color=style["palette"]["ink"])
                elif e.dxftype() == "POLYLINE":
                    pts = np.array([(v.dxf.location.x, v.dxf.location.y) for v in e.vertices])
                    if len(pts):
                        ax.plot(pts[:, 0], pts[:, 1], lw=0.2, color=style["palette"]["ink"])
            except Exception:
                continue
        ax.set_aspect("equal"); ax.axis("off"); ax.set_title(name)
    fig.suptitle(f"Same site, two drawings - {TODAY}")
    save_fig(fig, out / "compare_archimap")
    log(f"   comparison written: compare_archimap.md / .png")
    return out / "compare_archimap.md"


# --------------------------------------------------------------------------------------
# RUN
# --------------------------------------------------------------------------------------
def run(address, radius=300, style_path=None, out_root="site_layers-open_out", offline=False, jurisdiction=None,
        compare=None, cache_dir=None, quiet=False):
    radius = int(max(100, min(1500, radius)))
    out_root = Path(out_root)
    out_root.mkdir(parents=True, exist_ok=True)
    cache = Path(cache_dir) if cache_dir else out_root / "_cache"
    if cache.suffix.lower() == ".zip" and cache.is_file():   # a zipped cache pack (e.g. the class demo): unpack once
        target = cache.with_suffix("")
        if not target.exists():
            with zipfile.ZipFile(cache) as z:
                for m in z.namelist():
                    if m.endswith("/"):
                        continue
                    (target / Path(m).name).parent.mkdir(parents=True, exist_ok=True)
                    (target / Path(m).name).write_bytes(z.read(m))
        cache = target
    c = Ctx()
    c.address, c.r = address, radius
    tmp_log = Log(quiet=quiet)
    c.log = tmp_log
    c.reg = Register()
    c.net = Net(cache, offline=offline, log=c.log)
    c.log(f"{TOOL} v{VERSION}  -  {address}  -  radius {radius} m" + ("  -  OFFLINE, from cache" if offline else ""))
    c.pin = geocode(c.net, address, c.reg, c.log)
    c.jur = jurisdiction or c.pin["jurisdiction"]
    c.adapter = adapter_for(c.jur)
    c.log(f"   pin {c.pin['lat']:.6f}, {c.pin['lon']:.6f}  -  jurisdiction {c.jur} ({c.adapter['name']})")
    slug = slugify(address)
    out = out_root / f"site_{slug}"
    out.mkdir(parents=True, exist_ok=True)
    c.log.path = out / "run_log.txt"
    c.frame = LocalFrame(c.pin["lat"], c.pin["lon"])
    c.L = new_layers()
    c.site_geom = Point(0, 0).buffer(6, 32)
    c.site_label = ""
    c.climate = None
    c.sun = None
    style, style_src = load_style(style_path, c.log)
    font = apply_font(style, c.log)

    steps = [("OpenStreetMap physical layers", step_osm, ["buildings", "roads", "paths", "rail", "water", "trees", "landuse"]),
             ("Cadastre", step_cadastre, ["cadastre", "site"]),
             ("Terrain", step_terrain, ["contours"]),
             ("Statutory layers", step_statutory, [k for k, _, _ in STAT_GENERIC] + ["shr"]),
             ("Climate", step_climate, []),
             ("Sun", step_sun, []),
             ("Standing gaps", standing_gaps, [])]
    for name, fn, keys in steps:
        t0 = time.time()
        c.log(f"-> {name}")
        try:
            fn(c)
            c.log(f"   done in {time.time() - t0:.1f} s")
        except (Exception, NotCached) as e:
            msg = str(e) or type(e).__name__
            c.log(f"   FAILED: {msg} - carrying on without it")
            for k in keys:
                if c.L[k].status in ("skipped", "ok", "empty") and not len(c.L[k]):
                    c.L[k].status, c.L[k].note = "failed", msg
            lay = {"Climate": "AI_WIND", "Sun": "AI_SUNPATH"}.get(name, c.L[keys[0]].dxf if keys else name)
            c.reg.add(lay, f"{name}: layer missing from this run ({msg})", "tool", "", TODAY, "n/a",
                      "Rerun later (servers are sometimes busy), or add this layer by hand from the source named in README.")

    files = []
    c.log("-> writing outputs")
    for label, fn in (("DXF", lambda: [write_dxf(c, out / f"site_{slug}.dxf")]),
                      ("GeoJSON", lambda: write_geojson(c, out / f"layers_{slug}")),
                      ("heights CSV", lambda: (write_heights(c, out / f"building_heights_{slug}.csv"), [out / f"building_heights_{slug}.csv"])[1]),
                      ("plan diagram", lambda: render_plan(c, style, out)),
                      ("statutory diagram", lambda: render_statutory(c, style, out)),
                      ("exploded axon", lambda: render_exploded(c, style, out)),
                      ("climate diagram", lambda: render_climate(c, style, out))):
        try:
            got = fn()
            files += [f for f in (got or [])]
            c.log(f"   {label}: ok")
        except Exception as e:
            c.log(f"   {label}: FAILED ({e})")
            c.log("      " + traceback.format_exc().strip().splitlines()[-1])
    # style card, prompt suffix, register, about
    (out / "prompt_suffix.txt").write_text(prompt_suffix(style) + "\n", encoding="utf-8")
    (out / "image_prompt_exploded_axon.txt").write_text(
        "Redraw the attached exploded site-layer axonometric as a finished presentation diagram. Keep every plate, its order, "
        "its geometry and its label exactly; change only the graphic treatment. " + prompt_suffix(style) + "\n", encoding="utf-8")
    c.reg.add("STYLE", f"Diagrams drawn with style card: {style_src}; typeface used: {font}", "your style card", "", TODAY, "high",
              "If this says 'defaults', your style card was not read.")
    reg_p = out / f"register_{slug}.csv"
    c.reg.save(reg_p)
    files += [reg_p, out / "prompt_suffix.txt", out / "image_prompt_exploded_axon.txt"]
    if compare:
        try:
            files.append(compare_archimap(out / f"site_{slug}.dxf", Path(compare), out, style, c.log))
            files.append(out / "compare_archimap.png")
        except Exception as e:
            c.log(f"   archiMap comparison FAILED ({e})")
    about = out / "ABOUT_THIS_RUN.md"
    write_about(c, about, files)
    files.append(about)
    c.log.save()
    files.append(out / "run_log.txt")
    zp = out_root / f"site_{slug}.zip"
    with zipfile.ZipFile(zp, "w", zipfile.ZIP_DEFLATED) as z:
        for f in out.rglob("*"):
            if f.is_file():
                z.write(f, f.relative_to(out_root))
    c.log(f"\nDone. {len(c.reg.rows)} facts in the register, every verdict blank. Outputs: {out}  (zip: {zp.name})")
    c.log.save()
    c.out, c.zip = out, zp
    return c


def main(argv=None):
    ap = argparse.ArgumentParser(description=f"{TOOL} v{VERSION}: address in, layered DXF and diagrams out (open data).")
    ap.add_argument("address", nargs="?", help="street address, or 'lat, lon'")
    ap.add_argument("--radius", type=int, default=300, help="study radius in metres (100-1500, default 300)")
    ap.add_argument("--style", default="style_card-open.toml", help="your style card (default style_card-open.toml)")
    ap.add_argument("--out", default="site_layers-open_out", help="output folder")
    ap.add_argument("--offline", action="store_true", help="use cached data only (for demos without network)")
    ap.add_argument("--cache", default=None, help="cache folder (default <out>/_cache)")
    ap.add_argument("--jurisdiction", default=None, help="override, e.g. AU-NSW, AU-VIC, NZ, SG")
    ap.add_argument("--compare", default=None, help="an archiMap DXF of the same site, to compare")
    ap.add_argument("--compare-only", nargs=2, metavar=("OURS_DXF", "ARCHIMAP_DXF"), help="compare two existing DXFs and stop")
    a = ap.parse_args(argv)
    if a.compare_only:
        o = Path(a.compare_only[0]).parent
        compare_archimap(Path(a.compare_only[0]), Path(a.compare_only[1]), o)
        return
    if not a.address:
        ap.error("give an address, e.g.  python site_layers-open.py \"388 George St, Sydney NSW\"")
    style = a.style if Path(a.style).exists() else None
    run(a.address, a.radius, style, a.out, a.offline, a.jurisdiction, a.compare, a.cache)


if __name__ == "__main__":
    main()