"""Extract Cedar Gate Road's 40 drawn cross-sections from C3277-912-001 Rev A.

Uses the drawing's H 1:250 / V 1:100 scales and snaps axes to the plotted grid.
Connected road design paths are separated from the independent service symbols.
Quantities are drawing-derived end-area estimates, not measured backfill or as-built volumes.
"""
import argparse
import base64
import hashlib
import json
import math
from pathlib import Path
import re
import subprocess
import tempfile

import numpy as np
from shapely.geometry import LineString, Point

from extract_cedar_gate import alignment
from extract_xs import load_svg_paths, load_words

H_SCALE = 1000/250/25.4*72
V_SCALE = -1000/100/25.4*72
STEP = .05


def read_sections(svg_path, words_path):
    paths = load_svg_paths(Path(svg_path).read_text())
    words = load_words(Path(words_path).read_text())
    black = lambda p: p['stroke'] == 'rgb(0%, 0%, 0%)'
    grey = lambda p: p['stroke'].startswith('rgb(72.')
    inside = lambda r, x, y, pad=0: r[0]-pad <= x <= r[2]+pad and r[1]-pad <= y <= r[3]+pad
    rects = []
    for p in paths:
        if black(p) and p['w'] == 17 and len(p['pts']) >= 4:
            xx, yy = zip(*p['pts'])
            if max(xx)-min(xx) > 300 and max(yy)-min(yy) > 30:
                rects.append((min(xx), min(yy), max(xx), max(yy)))
    result = []
    for i, word in enumerate(words):
        if word[0] != 'Peg':
            continue
        cx, cy = (word[1]+word[3])/2, (word[2]+word[4])/2
        strip = next((r for r in rects if inside(r, cx, cy) and r[3]-r[1] < 60), None)
        if strip is None:
            raise ValueError('Cross-section caption lacks a frame')
        plot = next((r for r in rects if abs(r[3]-strip[1]) < 2 and abs(r[0]-strip[0]) < 2 and r[3]-r[1] > 60), None)
        if plot is None:
            raise ValueError('Cross-section plot lacks a frame')
        following = [w for w in words[i+1:i+7] if abs((w[2]+w[4])/2-cy) < 3]
        chainage = next((float(w[0]) for w in following if re.fullmatch(r'\d+(\.\d+)?', w[0])), None)
        if chainage is None:
            raise ValueError('Cross-section chainage unreadable')
        gx = sorted({p['pts'][0][0] for p in paths if grey(p) and len(p['pts']) == 2 and abs(p['pts'][0][0]-p['pts'][1][0]) < .01 and inside(plot, *p['pts'][0])})
        gy = sorted({p['pts'][0][1] for p in paths if grey(p) and len(p['pts']) == 2 and abs(p['pts'][0][1]-p['pts'][1][1]) < .01 and inside(plot, *p['pts'][0])})
        offsets = [(float(w[0]), (w[1]+w[3])/2) for w in words if re.fullmatch(r'-?\d{1,2}', w[0]) and inside(strip, (w[1]+w[3])/2, (w[2]+w[4])/2)]
        snapped = [(o, min(gx, key=lambda g: abs(g-x))) for o, x in offsets if o % 5 == 0 and min(abs(g-x) for g in gx) < 2]
        if len(snapped) < 5:
            raise ValueError(f'{chainage}: insufficient offset labels')
        x0 = float(np.median([x-H_SCALE*o for o, x in snapped]))
        levels = [(float(w[0]), (w[2]+w[4])/2) for w in words if re.fullmatch(r'13\d{2}', w[0]) and inside(plot, (w[1]+w[3])/2, (w[2]+w[4])/2, 1)]
        level_pairs = [(z, min([g for g in gy if 0 < g-y < 10])) for z, y in levels if any(0 < g-y < 10 for g in gy)]
        if len(level_pairs) < 2:
            raise ValueError(f'{chainage}: insufficient level labels')
        y0 = float(np.median([y-V_SCALE*z for z, y in level_pairs]))
        to_m = lambda p: ((p[0]-x0)/H_SCALE, (p[1]-y0)/V_SCALE)
        local = []
        for p in paths:
            if black(p) and p['w'] in (6, 9, 12) and all(inside(plot, x, y, -.3) for x, y in p['pts']):
                local.append({'w': p['w'], 'pts': [to_m(p) for p in p['pts']]})
        # The 0.72 pt finished road profile is the seed. Select road layerwork,
        # kerbs, sidewalks and batters connected to it; exclude separate pipe markers.
        designs = [p for p in local if p['w'] in (6, 12)]
        connected = [p for p in designs if p['w'] == 12]
        if len(connected) != 1:
            raise ValueError(f'{chainage}: expected one finished road profile')
        remaining = [p for p in designs if p['w'] != 12]
        for _ in range(len(remaining)+1):
            geometry = [LineString(p['pts']) for p in connected]
            joined = [p for p in remaining if any(LineString(p['pts']).distance(g) < .012 for g in geometry)]
            if not joined:
                break
            connected.extend(joined)
            remaining = [p for p in remaining if p not in joined]
        ground = []
        for p in local:
            if p['w'] != 9:
                continue
            for a, b in zip(p['pts'], p['pts'][1:]):
                if abs(a[0]-b[0]) > .02:
                    ground.append(tuple(sorted((a, b))))
        if not ground:
            raise ValueError(f'{chainage}: no natural ground line')
        xs = np.arange(-15, 15+STEP/2, STEP)
        def sample(segments, aggregate):
            out = np.full(len(xs), np.nan)
            for (ax, ay), (bx, by) in segments:
                mask = (xs >= ax-1e-8) & (xs <= bx+1e-8)
                value = ay+(xs[mask]-ax)*(by-ay)/(bx-ax)
                out[mask] = np.where(np.isnan(out[mask]), value, aggregate(out[mask], value))
            return out
        ngl = sample(ground, np.maximum)
        good = np.isfinite(ngl)
        if good.sum() < 20:
            raise ValueError(f'{chainage}: insufficient natural ground')
        ngl = np.interp(xs, xs[good], ngl[good])  # bridge plotted dashes
        design_segments = [tuple(sorted((a, b))) for p in connected for a, b in zip(p['pts'], p['pts'][1:]) if abs(a[0]-b[0]) > .001]
        top = sample(design_segments, np.maximum)
        bottom = sample(design_segments, np.minimum)
        cut = np.maximum(ngl-bottom, 0)
        fill = np.maximum(bottom-ngl, 0)
        # Integrate only contiguous runs of drawn formation; do not bridge blank spans.
        def integrate(values):
            mask = np.isfinite(values[:-1]) & np.isfinite(values[1:])
            return float(np.sum((values[:-1][mask]+values[1:][mask])*.5*STEP))
        result.append({'ch': chainage, 'x': xs, 'ngl': ngl, 'top': top, 'bot': bottom,
                       'cutArea': integrate(cut), 'fillArea': integrate(fill),
                       'roadWidth': float(np.ptp(xs[np.isfinite(bottom)])), 'excludedIndependentPaths': len(remaining)})
    result.sort(key=lambda s: s['ch'])
    if len(result) != 40 or [s['ch'] for s in result] != [*range(0, 381, 10), 385.411]:
        raise ValueError('Expected all 40 drawn sections from 0 to 385.411 m')
    return result


def dem(surface):
    q = np.frombuffer(base64.b64decode(surface['z']), dtype='<i2')
    grid = np.where(q < 0, np.nan, surface['zmin']+q/1000).reshape(surface['ny'], surface['nx'])
    def at(e, n):
        x, y = (e-surface['x0'])/surface['cell'], (n-surface['y0'])/surface['cell']
        if x < 0 or y < 0 or x > surface['nx']-1 or y > surface['ny']-1:
            return None
        i, j = min(int(x), surface['nx']-2), min(int(y), surface['ny']-2)
        u, v = x-i, y-j
        values = grid[j:j+2, i:i+2].ravel(); weights = np.array([(1-u)*(1-v), u*(1-v), (1-u)*v, u*v])
        valid = np.isfinite(values)
        if not valid.any() or weights[valid].sum() == 0:
            return None
        return float(np.sum(values[valid]*weights[valid])/np.sum(weights[valid]))
    return at


def build_sections(sections, raw, folder, source):
    road = json.loads((folder/'bushveld.json').read_text())
    surface = json.loads((folder/'surface.json').read_text())
    previous = json.loads((folder/'earthworks.json').read_text())
    rows = raw['road']['RDS_LO$0$CC_RDS_NEW_INT_CL']
    line = alignment(rows, rows[0]['pts'][0], 385.411)
    ox, oy = road['origin']['cadX'], road['origin']['cadY']
    def at(ch):
        distance = min(max(ch, 0), line.length)
        point = line.interpolate(distance)
        a, b = line.interpolate(max(0, distance-.05)), line.interpolate(min(line.length, distance+.05))
        direction = np.array([b.x-a.x, b.y-a.y]); direction /= np.linalg.norm(direction)
        return np.array([point.x-ox, point.y-oy]), direction
    level = dem(surface)
    def differences(side):
        errors = []
        for s in sections:
            point, direction = at(s['ch']); normal = np.array([direction[1], -direction[0]])*side
            for offset in np.arange(-2, 2.01, .5):
                index = int(round((offset+15)/STEP)); z = s['top'][index]
                if not np.isfinite(z):
                    continue
                e, n = point+normal*offset; value = level(e, n)
                if value is not None:
                    errors.append(float(z-value))
        return np.array(errors)
    errors = {side: differences(side) for side in (1, -1)}
    side = min(errors, key=lambda side: np.median(np.abs(errors[side])))
    delta = errors[side]
    if len(delta) < 200 or np.median(np.abs(delta)) > .15 or np.percentile(np.abs(delta), 90) > .35:
        raise ValueError(f'Cross-sections do not agree with finished-level contours: {len(delta)} samples, median {np.median(np.abs(delta))}, p90 {np.percentile(np.abs(delta), 90)}')
    out = []
    for s in sections:
        point, direction = at(s['ch'])
        keep = np.arange(0, len(s['x']), 5)  # 0.25 m browser samples
        out.append({'ch': s['ch'], 'e': round(float(point[0]), 4), 'n': round(float(point[1]), 4), 'dir': direction.round(7).tolist(), 'side': side,
                    'x': s['x'][keep].round(3).tolist(), 'ngl': s['ngl'][keep].round(4).tolist(),
                    'top': [round(float(z), 4) if np.isfinite(z) else None for z in s['top'][keep]],
                    'bot': [round(float(z), 4) if np.isfinite(z) else None for z in s['bot'][keep]],
                    'cutArea': round(s['cutArea'], 4), 'fillArea': round(s['fillArea'], 4)})
    volumes = [{'from': a['ch'], 'to': b['ch'], 'cut': round((a['cutArea']+b['cutArea'])*.5*(b['ch']-a['ch']), 3),
                'fill': round((a['fillArea']+b['fillArea'])*.5*(b['ch']-a['ch']), 3)} for a, b in zip(sections, sections[1:])]
    section_keys = [s['key'] for s in road['sections'] if s.get('name') == 'Cedar Gate Road']
    for s in road['sections']:
        if s['key'] == section_keys[-1]:
            s['to'] = 385.411
    result = {'source': source, 'method': 'End-area between 10 m drawn sections (last strip 5.411 m); natural ground versus connected road formation. Service trench symbols excluded.',
              'sectionKeys': section_keys, 'chainageOffset': 0, 'chainageEnd': 385.411, 'sections': out, 'volumes': volumes,
              'alignment': [[round(float(ch), 3), *at(ch)[0].round(4).tolist()] for ch in np.arange(0, line.length, 2)],
              'totals': {key: round(sum(v[key] for v in volumes), 3) for key in ['cut', 'fill']},
              'validation': {'samples': len(delta), 'meanAbs': round(float(np.mean(np.abs(delta))), 4), 'median': round(float(np.median(delta)), 4),
                             'medianAbs': round(float(np.median(np.abs(delta))), 4), 'p90Abs': round(float(np.percentile(np.abs(delta), 90)), 4),
                             'note': 'Drawn finished road profile versus independently extracted FRL_CONTOURS; tested both offset directions.'}}
    if previous.get('groundSurface'):
        result['groundSurface'] = previous['groundSurface']
    old = 'Bulk cut/fill quantities are unavailable until the road cross-sections and formation levels are reconciled; no zero-volume claim is made.'
    replacement = 'Cut/fill covers Cedar Gate Road only, from C3277-912-001 Rev A sections 0–385.411 m. Road formation earthworks exclude service trenches and external roads; fill is design compacted fill, not pipe-trench backfill.'
    road['notes'] = [replacement if note == old else note for note in road.get('notes', [])]
    (folder/'earthworks.json').write_text(json.dumps(result, separators=(',', ':'), allow_nan=False))
    (folder/'bushveld.json').write_text(json.dumps(road, separators=(',', ':'), allow_nan=False))
    report = {'sections': len(out), 'chainage': [0, 385.411], 'totals': result['totals'], 'validation': result['validation'],
              'offset_side': side, 'widths': [[s['ch'], round(s['roadWidth'], 3)] for s in sections],
              'excluded_independent_service_paths': sum(s['excludedIndependentPaths'] for s in sections), 'scope': replacement,
              'earthworks_sha256': hashlib.sha256((folder/'earthworks.json').read_bytes()).hexdigest()}
    (folder/'earthworks-validation.json').write_text(json.dumps(report, indent=2))
    validation_path = folder/'validation.json'
    if validation_path.exists():
        validation = json.loads(validation_path.read_text()); validation['earthworks'] = report; validation['notes'] = road['notes']
        validation['sha256'] = {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in [folder/f'{f}.json' for f in ['bushveld','services','earthworks','schedules','surface']]}
        validation_path.write_text(json.dumps(validation, indent=2))
    print(json.dumps(report, indent=2))


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('pdf', type=Path)
    parser.add_argument('raw_geometry', type=Path)
    parser.add_argument('model_folder', type=Path)
    parser.add_argument('--poppler-bin', type=Path, default=Path('/usr/local/bin'))
    args = parser.parse_args()
    with tempfile.TemporaryDirectory() as tmp:
        svg, words = Path(tmp)/'sections.svg', Path(tmp)/'words.html'
        subprocess.run([str(args.poppler_bin/'pdftocairo'), '-svg', str(args.pdf), str(svg)], check=True)
        subprocess.run([str(args.poppler_bin/'pdftotext'), '-bbox', str(args.pdf), str(words)], check=True)
        build_sections(read_sections(svg, words), json.loads(args.raw_geometry.read_text()), args.model_folder, args.pdf.name)


if __name__ == '__main__':
    main()