cb4a1ce461
- Sites → Maps in the sidebar replaces the per-site Site map / Map 2 links. The Maps page shows each site's model as a plan preview (built sections solid) with % of road layerworks built and links to Site map and Map 2. - Double-clicking a road section opens its progress tick-boxes over the map (Site map saves to the server; Map 2 stays a local scenario). - Cedar Gate: 40 drawn cross-sections extracted from C3277-912-001 Rev A (extract_cedar_gate_xs.py), cut/fill limited to Cedar Gate Road, a Cut & fill workspace in Map 2 with a chainage chart, and pavement layers shown in their material colours. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
231 lines
14 KiB
Python
231 lines
14 KiB
Python
"""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()
|