"""Build a Map 2/3 model for Nellmapius Ext. 25 from the 100_ROADS and 400_KERBING drawings.

What the supplied files give: the kerb and road-edge lines, stand lines, road names and widths, plus four
survey benchmarks (SVR detail contour plan). They give no design levels, no pavement build-up and no
stormwater or sewer data in a form that can be placed, so the model has road areas, kerbs and context only,
over an indicative plane through the four benchmarks.

Requires numpy shapely ezdxf scipy. Coordinates: local east/north metres = CAD (x=-Lo Y, y=-Lo X) minus origin.

  python3 -I extract_nellmapius.py 100_ROADS.dxf 400_KERBING.dxf OUT_DIR
"""
import base64
import hashlib
import io
import json
import math
import re
import sys
from collections import defaultdict
from pathlib import Path

import numpy as np
import shapely
from scipy.spatial import cKDTree
from shapely.geometry import LineString, MultiPoint, Point, Polygon, box
from shapely.ops import polygonize, unary_union

# Survey benchmarks printed on the SVR detail contour plan (Lo29 Y, X, level), Sep 2024
BENCHMARKS = [("ST18", 65305.974, 2847299.892, 1298.993, "12mm peg next to kerb"), ("BM3", 65203.736, 2847335.083, 1300.044, "hole in kerb"),
              ("ST17", 65152.153, 2847435.083, 1302.880, "12mm peg in paving"), ("BM2", 65182.192, 2847332.670, 1300.406, "hole in kerb")]
KERB_LAYERS = {"Kerb_300 Mountable": "300 mm mountable kerb", "Kerb_Edge Beam": "150 mm edge beam", "Kerb_Fig7_Semi-mount": "Fig. 7 semi-mountable kerb"}
NAMES = {"IDELO CRESCENT": "Idelo Crescent", "UMNGAMO STREET": "Umngamo Street", "IGWALI STREET": "Igwali Street"}   # roads built in this phase (cross-section titles)
SECTION = 20.0


def read_dxf(path, layers):
    """TEXT/MTEXT/lines of the wanted layers; repairs LibreDWG's multiline text and unsupported payloads."""
    import ezdxf
    from ezdxf.path import make_path
    lines = Path(path).read_text(errors="replace").splitlines()
    pairs, i = [], 0
    while i < len(lines):
        try:
            code = int(lines[i])
        except ValueError:
            pairs[-1] = (pairs[-1][0], pairs[-1][1] + r"\P" + lines[i]); i += 1; continue
        pairs.append((code, lines[i + 1])); i += 2
    clean, skip = [], False
    for code, value in pairs:
        if code == 0:
            skip = False
        if code == 101 and value == "Embedded Object":
            skip = True
        if not skip:
            clean.extend((str(code), value))
    doc = ezdxf.read(io.StringIO("\n".join(clean) + "\n"))
    out = defaultdict(list)
    for e in doc.modelspace():
        if e.dxf.layer not in layers:
            continue
        k = e.dxftype()
        if k in ("LINE", "ARC", "LWPOLYLINE", "POLYLINE", "CIRCLE"):
            out[e.dxf.layer].append({"type": k, "pts": [list(p)[:2] for p in make_path(e).flattening(.05)]})
        elif k in ("TEXT", "MTEXT"):
            out[e.dxf.layer].append({"type": k, "text": e.plain_text() if k == "MTEXT" else e.dxf.text, "pos": list(e.dxf.insert)[:2]})
    return out


def plane_fit(points):
    A = np.array([[1, e, n] for e, n, _ in points]); z = np.array([p[2] for p in points])
    coef, *_ = np.linalg.lstsq(A, z, rcond=None)
    return coef, z - A @ coef


def build(roads, kerbs, out_dir, source_files):
    edge = [LineString(r["pts"]) for r in kerbs["_X25 ROAD EDGE"] if len(r["pts"]) > 2]
    closed = [l for l in edge if math.dist(l.coords[0], l.coords[-1]) < 0.05]
    cx0, cy0, cx1, cy1 = unary_union(closed).bounds
    window = box(cx0 - 120, cy0 - 120, cx1 + 120, cy1 + 120)
    origin = (round((cx0 + cx1) / 2, -1), round((cy0 + cy1) / 2, -1))
    local = lambda pts: [[round(p[0] - origin[0], 3), round(p[1] - origin[1], 3)] for p in pts]
    inside = lambda pts: window.contains(LineString(pts)) if len(pts) > 1 else False

    # ---- road faces: the ring roads are the space between the closed edge lines (kerbs lie on them); the Igwali spur is closed by a hull
    kerb_lines = [(k, LineString(r["pts"])) for k in KERB_LAYERS for r in kerbs[k] if len(r["pts"]) > 1 and inside(r["pts"])]
    faces = list(polygonize(unary_union([g for _, g in kerb_lines] + closed)))
    ring = [f for f in faces if f.area > 20 and f.area / f.length < 8]
    blocks = [f for f in faces if f.area > 20 and f.area / f.length >= 8]
    road = unary_union(ring)   # the 932 m outer edge line includes the Igwali Street spur, so the face already has it
    report = {"roadFaces": [round(f.area, 1) for f in ring], "blocks": [round(b.area) for b in blocks]}

    # ---- 20 m sections: along each block's boundary (ring roads) and along the spur axis
    labels = []   # (name, x, y)
    for r in roads["_NEW ROAD NAMES"] + roads["ROAD NAMES"]:
        t = re.sub(r"[()]|\d+m", "", r["text"]).strip().upper()
        if t in NAMES:
            labels.append((NAMES[t], r["pos"][0], r["pos"][1]))
    block_rings = [b.exterior for b in blocks]
    block_minx = min(b.bounds[0] for b in blocks)
    cell = 1.0
    minx, miny, maxx, maxy = road.bounds
    cells = [(x + .5, y + .5) for x in np.arange(math.floor(minx), math.ceil(maxx), cell) for y in np.arange(math.floor(miny), math.ceil(maxy), cell) if road.contains(Point(x + .5, y + .5))]

    spur_cells = [(x, y) for x, y in cells if x < block_minx - 2]
    west = min(spur_cells, key=lambda c: c[0])
    near_j = [c[1] for c in spur_cells if block_minx - 3 < c[0] < block_minx - 2] or [west[1]]
    axis = LineString([(block_minx - 2, float(np.mean(near_j))), west])
    n_spur = max(1, int(axis.length // SECTION))

    def chunk_of(x, y):
        if x < block_minx - 2:
            return ("spur", min(n_spur - 1, int(axis.project(Point(x, y)) // SECTION)))
        p = Point(x, y)
        k = min(range(len(block_rings)), key=lambda i: block_rings[i].distance(p))
        return (f"loop{k}", int(block_rings[k].project(p) // SECTION))

    groups = defaultdict(list)
    for x, y in cells:
        groups[chunk_of(x, y)].append((x, y))
    exact = road.area / max(1, len(cells))
    sections, areas, chunk_geo = [], [], {}
    for (kind, idx), pts in sorted(groups.items()):
        cx, cy = np.mean([p[0] for p in pts]), np.mean([p[1] for p in pts])
        name = (min(labels, key=lambda l: math.hypot(l[1] - cx, l[2] - cy))[0] if labels else "Road") if kind != "spur" else "Igwali Street"
        if kind.startswith("loop") and name == "Igwali Street":
            name = "Idelo Crescent"
        tag = kind if kind != "spur" else "spur"
        f, t = (idx * SECTION, (idx + 1) * SECTION)
        key = f"{re.sub('[^a-z0-9]+', '-', name.lower())}-{tag}-ch{int(f)}-{int(t)}"
        poly = unary_union([box(x - .5, y - .5, x + .5, y + .5) for x, y in pts]).buffer(0)
        parts = list(poly.geoms) if poly.geom_type == "MultiPolygon" else [poly]
        sections.append({"key": key, "name": f"{name} ({'spur' if kind == 'spur' else 'loop ' + 'AB'[int(kind[-1])] if len(blocks) > 1 else 'loop'})", "from": f, "to": t})
        chunk_geo[key] = poly
        for part in parts:
            areas.append({"kind": "road_surface", "section": key, "area": round(part.area * exact / cell ** 2, 3), "rings": [local(list(part.simplify(.05).exterior.coords))]})
    assert abs(sum(a["area"] for a in areas) - road.area) < 1, (sum(a["area"] for a in areas), road.area)

    # ---- kerbs: each metre of kerb line goes to the section it sits in
    tree = cKDTree(np.array([(g.centroid.x, g.centroid.y) for g in chunk_geo.values()])); keys = list(chunk_geo)
    owner = lambda x, y: keys[int(tree.query([x, y])[1])] if False else min(keys, key=lambda k: chunk_geo[k].distance(Point(x, y)))
    kerb_rows, total_kerb = [], 0.0
    for lname, g in kerb_lines:
        runs, cur, cur_key = [], [], None
        n = max(1, int(round(g.length)))
        pts = [g.interpolate(i * g.length / n) for i in range(n + 1)]
        for a, b in zip(pts, pts[1:]):
            k = owner((a.x + b.x) / 2, (a.y + b.y) / 2)
            if k != cur_key and cur:
                runs.append((cur_key, cur)); cur = []
            cur_key = k; cur = cur + [a.coords[0]] if not cur else cur; cur.append(b.coords[0])
        if cur:
            runs.append((cur_key, cur))
        for k, line in runs:
            L = LineString(line).length; total_kerb += L
            kerb_rows.append({"section": k, "length": round(L, 3), "pts": local(line)})

    context = [{"kind": "tar_edge", "pts": local(list(l.coords))} for l in closed if window.contains(l)]
    context += [{"kind": "other_road", "pts": local(r["pts"])} for r in kerbs["KERB_FIG 10"] if inside(r["pts"])]
    context += [{"kind": "marking", "pts": local(r["pts"])} for r in roads["DCE-RMRS-LINE-W-SOLID"] if inside(r["pts"])]
    context += [{"kind": "other_road", "pts": local(r["pts"])} for r in roads["2025_Calc Plan"] if inside(r["pts"])]   # stand lines
    pegs = []
    for name, Y, X, z, note in BENCHMARKS:
        e, n = local([(-Y, -X)])[0]
        pegs.append({"name": name, "e": e, "n": n, "z": z, "Y": Y, "X": X})
    notes = [
        "Nellmapius Ext. 25 (SVR / Civil Concepts drawings 100_ROADS and 400_KERBING). Roads read from the closed road-edge lines and kerb lines; widths per the cross-section titles: Idelo Crescent and Umngamo Street 5.0 m in 13 m reserve, Igwali Street 8.0 m in 20 m reserve.",
        "20 m sections run along the inner edge of each ring road and along the Igwali Street spur; road names come from the nearest street-name label, so check them.",
        "No pavement build-up was supplied: each section has one 'road surface' item to tick off. No design levels or contours were supplied.",
        "The ground is an indicative plane through the four survey benchmarks (Sep 2024). It is not a survey surface: do not use it for earthworks, cut and fill or levels.",
        "No stormwater, sewer or water pipes are placed: the sewer long section supplied is sheet 3 of 3 and the services layout carries no levels.",
        f"Kerb lines measured once: {total_kerb:.0f} m (300 mm mountable, 150 mm edge beam, Fig. 7 semi-mountable).",
    ]
    road_model = {"project": "Nellmapius Ext. 25", "origin": {"cadX": origin[0], "cadY": origin[1]},
                  "builds": {"road_surface": {"name": "Road surface", "color": 0x6B7280, "layers": [{"id": "road-surface", "name": "Road surface (layers not supplied)", "thk": 0.1, "unit": "m²", "color": 0x6B7280}]}},
                  "sections": sections, "areas": areas, "kerbs": kerb_rows, "kerbQuantityFactor": 1, "context": context, "contours": [], "pegs": pegs, "water": [], "buffer": [], "notes": notes}

    # ---- indicative ground: a plane through the benchmarks
    pts3 = [(p["e"], p["n"], p["z"]) for p in pegs]
    coef, res = plane_fit(pts3)
    x0, y0 = math.floor(minx - origin[0] - 40), math.floor(miny - origin[1] - 40)
    nx, ny = math.ceil(maxx - origin[0] + 40) - x0 + 1, math.ceil(maxy - origin[1] + 40) - y0 + 1
    gx, gy = np.meshgrid(x0 + np.arange(nx), y0 + np.arange(ny))
    Z = coef[0] + coef[1] * gx + coef[2] * gy
    zmin = round(float(Z.min()) - 1, 1)
    q = np.round((Z - zmin) * 1000).astype("<i2")
    surface = {"x0": float(x0), "y0": float(y0), "cell": 1, "nx": int(nx), "ny": int(ny), "zmin": zmin, "z": base64.b64encode(q.tobytes()).decode(),
               "encoding": "int16 little-endian, mm above zmin, -1 = outside supplied contour coverage, row-major",
               "source": "Indicative plane through four survey benchmarks (SVR Land Surveyors, Sep 2024). Not a survey surface.",
               "check": {"pegs": 4, "medianAbs": round(float(np.median(np.abs(res))), 3), "p90Abs": round(float(np.percentile(np.abs(res), 90)), 3)}}
    earthworks = {"source": "None supplied", "method": "No design levels or cross-sections supplied for Nellmapius Ext. 25.", "chainageOffset": 0, "chainageEnd": 0, "sections": [], "volumes": [], "alignment": [],
                  "totals": {"cut": 0, "fill": 0}, "validation": {"samples": 0, "meanAbs": 0, "median": 0, "p90Abs": 0, "note": "No earthworks data."}}
    services = {"planRoutes": [], "stormwater": {"structures": [], "pipes": [], "notes": ["Not placed: no levels supplied."]}, "sewer": {"manholes": [], "pipes": [], "notes": ["Not placed: sewer long section sheet 3 of 3 only."]},
                "water": {"points": [], "pipeLabels": {}}, "holds": []}
    schedules = {"signs": [], "signCounts": {}, "markingLabels": [], "markingCounts": {}, "studs": [], "studCounts": {}, "markings": {}, "paintAreas": {}, "waterFittings": {}, "notes": ["Signs, studs and fittings are not quantified."]}
    out = Path(out_dir); out.mkdir(parents=True, exist_ok=True)
    files = {"bushveld": road_model, "services": services, "earthworks": earthworks, "schedules": schedules, "surface": surface}
    hashes = {}
    for name, data in files.items():
        text = json.dumps(data, separators=(",", ":")); (out / f"{name}.json").write_text(text); hashes[f"{name}.json"] = hashlib.sha256(text.encode()).hexdigest()
    report.update({"spurSections": n_spur, "sections": len(sections), "areaPieces": len(areas), "roadAreaM2": round(road.area, 1), "kerbM": round(total_kerb, 1), "kerbRows": len(kerb_rows), "contextLines": len(context),
                   "planeResidualsM": [round(float(r), 3) for r in res], "origin": origin, "sha256": hashes, "sources": source_files})
    (out / "validation.json").write_text(json.dumps(report, indent=2))
    (out / "sources.json").write_text(json.dumps({"sources": source_files, "excluded": "Cedar Gate (C3277) drawings are a different project."}, indent=2))
    print(json.dumps({k: report[k] for k in ("sections", "spurSections", "roadAreaM2", "kerbM", "planeResidualsM")}))


def main(roads_dxf, kerbs_dxf, out_dir):
    roads = read_dxf(roads_dxf, {"_NEW ROAD NAMES", "ROAD NAMES", "DCE-RMRS-LINE-W-SOLID", "2025_Calc Plan"})
    kerbs = read_dxf(kerbs_dxf, {"_X25 ROAD EDGE", "KERB_FIG 10", *KERB_LAYERS})
    src = [{"file": Path(p).name, "sha256": hashlib.sha256(Path(p).read_bytes()).hexdigest()} for p in (roads_dxf, kerbs_dxf)]
    build(roads, kerbs, out_dir, src)


if __name__ == "__main__":
    main(*sys.argv[1:4])