223 lines
11 KiB
Python
223 lines
11 KiB
Python
"""Read the road cross-section PDF (vector) and compute cut/fill per section and between sections.
|
|
|
|
Usage:
|
|
python3 -I extract_xs.py <cross-section.pdf> <road-layout.dxf> <site-progress-test/data>
|
|
|
|
How it reads the drawing (AutoCAD plot, pdftocairo SVG + pdftotext word boxes):
|
|
* each section = a thick frame (plot area) above a strip holding "Peg Distance (m) NNN m" and the offset labels
|
|
* offsets from the -35 … 35 labels in the strip, levels from the 14xx labels snapped to the grey grid lines
|
|
* natural ground = dashed line (stroke width 3); design = solid black lines (widths 2, 6, 8: surface, kerbs,
|
|
layerworks, 1:3 side slopes). Design top = finished surface, design bottom = underside of layerworks / slope.
|
|
Earthworks are measured between natural ground and the design bottom (formation). Layerworks themselves are
|
|
in the pavement BoQ already.
|
|
"""
|
|
import base64
|
|
import json
|
|
import math
|
|
import os
|
|
import re
|
|
import subprocess
|
|
import sys
|
|
import tempfile
|
|
from html import unescape
|
|
|
|
import numpy as np
|
|
|
|
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
|
from alignment import load_alignment # noqa: E402
|
|
|
|
H_SCALE = 1000 / 500 / 25.4 * 72 # 1:500 → pt per metre (5.669)
|
|
V_SCALE = 1000 / 250 / 25.4 * 72 # 1:250 → pt per metre (11.339)
|
|
STEP = 0.1 # sampling step across the section (m)
|
|
|
|
|
|
def load_svg_paths(svg):
|
|
paths = []
|
|
for m in re.finditer(r'<path ([^>]*?)/>', svg):
|
|
attrs = dict(re.findall(r'([a-z\-]+)="([^"]*)"', m.group(1)))
|
|
if 'd' not in attrs or attrs.get('fill', 'none') != 'none':
|
|
continue
|
|
a, b, c, d, e, f = (float(v) for v in re.findall(r'[-\d.e]+', attrs.get('transform', 'matrix(1,0,0,1,0,0)')))
|
|
pts, cur = [], []
|
|
for cmd, nums in re.findall(r'([MLCZ])([^MLCZ]*)', attrs['d']):
|
|
v = [float(x) for x in nums.split()]
|
|
if cmd == 'M':
|
|
if len(cur) > 1: pts.append(cur)
|
|
cur = [(v[0], v[1])]
|
|
elif cmd == 'L':
|
|
cur.append((v[0], v[1]))
|
|
elif cmd == 'C':
|
|
cur.append((v[4], v[5]))
|
|
elif cmd == 'Z' and cur:
|
|
cur.append(cur[0])
|
|
if len(cur) > 1: pts.append(cur)
|
|
for poly in pts:
|
|
paths.append({'stroke': attrs.get('stroke', ''), 'w': float(attrs.get('stroke-width', 1)),
|
|
'pts': [(a * x + c * y + e, b * x + d * y + f) for x, y in poly]})
|
|
return paths
|
|
|
|
|
|
def load_words(html):
|
|
return [(unescape(w), float(x0), float(y0), float(x1), float(y1)) for x0, y0, x1, y1, w in
|
|
re.findall(r'<word xMin="([\d.]+)" yMin="([\d.]+)" xMax="([\d.]+)" yMax="([\d.]+)">([^<]*)</word>', html)]
|
|
|
|
|
|
def main(pdf_path, dxf_path, data_dir):
|
|
base = json.load(open(os.path.join(data_dir, 'bushveld.json')))
|
|
ox, oy = base['origin']['cadX'], base['origin']['cadY']
|
|
tmp = tempfile.mkdtemp()
|
|
subprocess.run(['pdftocairo', '-svg', pdf_path, os.path.join(tmp, 'xs.svg')], capture_output=True)
|
|
subprocess.run(['pdftotext', '-bbox', pdf_path, os.path.join(tmp, 'w.html')], capture_output=True)
|
|
paths = load_svg_paths(open(os.path.join(tmp, 'xs.svg')).read())
|
|
words = load_words(open(os.path.join(tmp, 'w.html')).read())
|
|
|
|
black = lambda p: p['stroke'] in ('rgb(0%, 0%, 0%)',)
|
|
grey = lambda p: p['stroke'].startswith('rgb(50.')
|
|
rects = []
|
|
for p in paths:
|
|
if black(p) and p['w'] >= 6 and len(p['pts']) >= 4:
|
|
xs = [q[0] for q in p['pts']]; ys = [q[1] for q in p['pts']]
|
|
if max(xs) - min(xs) > 150 and max(ys) - min(ys) > 8:
|
|
rects.append((min(xs), min(ys), max(xs), max(ys)))
|
|
inside = lambda r, x, y, pad=0: r[0] - pad <= x <= r[2] + pad and r[1] - pad <= y <= r[3] + pad
|
|
|
|
sections = []
|
|
for i, w in enumerate(words):
|
|
if w[0] != 'Peg':
|
|
continue
|
|
cx, cy = (w[1] + w[3]) / 2, (w[2] + w[4]) / 2
|
|
strip = next((r for r in rects if inside(r, cx, cy) and r[3] - r[1] < 40), None)
|
|
if not strip:
|
|
continue
|
|
plot = min((r for r in rects if abs(r[3] - strip[1]) < 2 and abs(r[0] - strip[0]) < 2 and r[3] - r[1] > 40),
|
|
key=lambda r: abs(r[3] - strip[1]), default=None)
|
|
if not plot:
|
|
continue
|
|
line = [x for x in words[i + 1:i + 6] if abs((x[2] + x[4]) / 2 - cy) < 3]
|
|
ch = next((float(x[0]) for x in line if re.fullmatch(r'\d+(\.\d+)?', x[0])), None)
|
|
if ch is None:
|
|
continue
|
|
# offset axis
|
|
ow = [(float(x[0]), (x[1] + x[3]) / 2) for x in words if re.fullmatch(r'-?\d{1,2}', x[0]) and inside(strip, (x[1] + x[3]) / 2, (x[2] + x[4]) / 2)]
|
|
ow = [o for o in ow if o[0] % 5 == 0]
|
|
if len(ow) < 3:
|
|
continue
|
|
kx = H_SCALE # drawing scale is known (1:500) — fit only the position of offset 0
|
|
x0 = float(np.median([x - kx * o for o, x in ow]))
|
|
# level axis: labels snapped to grey grid lines just below them
|
|
grid_y = sorted({round(q[1], 2) for p in paths if grey(p) for q in p['pts'][:1]
|
|
if inside(plot, q[0], q[1]) and len(p['pts']) == 2 and abs(p['pts'][0][1] - p['pts'][1][1]) < 0.01})
|
|
lw = [(float(x[0]), (x[2] + x[4]) / 2) for x in words if re.fullmatch(r'1[34]\d{2}', x[0]) and inside(plot, (x[1] + x[3]) / 2, (x[2] + x[4]) / 2, 1)]
|
|
pairs = []
|
|
for lev, yc in lw:
|
|
g = [gy for gy in grid_y if 0 < gy - yc < 9]
|
|
pairs.append((lev, min(g) if g else yc + 4.0))
|
|
if len(pairs) < 2:
|
|
continue
|
|
ky = -V_SCALE # 1:250, page y grows downwards
|
|
y0 = float(np.median([y - ky * lev for lev, y in pairs]))
|
|
to_m = lambda X, Y: ((X - x0) / kx, (Y - y0) / ky)
|
|
|
|
ngl, des = [], []
|
|
for p in paths:
|
|
if not black(p) or p['w'] not in (2, 3, 6, 8):
|
|
continue
|
|
if not all(inside(plot, x, y, -1) for x, y in p['pts']):
|
|
continue
|
|
pts = [to_m(x, y) for x, y in p['pts']]
|
|
for a, b in zip(pts, pts[1:]):
|
|
if abs(b[0] - a[0]) < 1e-6:
|
|
continue
|
|
(ngl if p['w'] == 3 else des).append((a, b) if a[0] < b[0] else (b, a))
|
|
if not ngl:
|
|
continue
|
|
lo_x = math.floor(min(a[0] for a, _ in ngl)); hi_x = math.ceil(max(b[0] for _, b in ngl))
|
|
xs = np.round(np.arange(lo_x, hi_x + 1e-9, STEP), 3)
|
|
def sample(segs, agg):
|
|
out = np.full(len(xs), np.nan)
|
|
for (ax, ay), (bx_, by_) in segs:
|
|
sel = (xs >= ax) & (xs <= bx_)
|
|
v = ay + (xs[sel] - ax) * (by_ - ay) / (bx_ - ax)
|
|
out[sel] = np.where(np.isnan(out[sel]), v, agg(out[sel], v))
|
|
return out
|
|
g = sample(ngl, np.maximum)
|
|
ok = ~np.isnan(g)
|
|
g = np.interp(xs, xs[ok], g[ok]) # bridge the gaps between dashes
|
|
top = sample(des, np.maximum); bot = sample(des, np.minimum)
|
|
has = ~np.isnan(bot)
|
|
cut = float(np.nansum(np.where(has, np.clip(g - bot, 0, None), 0)) * STEP)
|
|
fill = float(np.nansum(np.where(has, np.clip(bot - g, 0, None), 0)) * STEP)
|
|
sections.append({'ch': ch, 'x': xs, 'ngl': g, 'top': top, 'bot': bot, 'cutArea': cut, 'fillArea': fill,
|
|
'scale': (round(kx, 3), round(ky, 3))})
|
|
|
|
sections.sort(key=lambda s: s['ch'])
|
|
# de-duplicate (same chainage printed twice)
|
|
uniq = []
|
|
for s in sections:
|
|
if not uniq or abs(uniq[-1]['ch'] - s['ch']) > 0.01:
|
|
uniq.append(s)
|
|
sections = uniq
|
|
|
|
al = load_alignment(dxf_path, ox, oy)
|
|
at, ch_end = al.at, al.length
|
|
# validate + choose offset side against the design surface (surface.json, built from the 0.1 m road contours)
|
|
sf = json.load(open(os.path.join(data_dir, 'surface.json')))
|
|
q = np.frombuffer(base64.b64decode(sf['z']), dtype='<i2').astype(float)
|
|
SZ = np.where(q < 0, np.nan, sf['zmin'] + q / 1000).reshape(sf['ny'], sf['nx'])
|
|
def dem(e, n):
|
|
fx = (e - sf['x0']) / sf['cell']; fy = (n - sf['y0']) / sf['cell']
|
|
if not (0 <= fx < sf['nx'] - 1 and 0 <= fy < sf['ny'] - 1): return np.nan
|
|
i, j = int(fx), int(fy); u, v = fx - i, fy - j
|
|
return SZ[j, i] * (1 - u) * (1 - v) + SZ[j, i + 1] * u * (1 - v) + SZ[j + 1, i] * (1 - u) * v + SZ[j + 1, i + 1] * u * v
|
|
def diffs(sign):
|
|
d = []
|
|
for s in sections:
|
|
p, t = at(s['ch']); nrm = np.array([t[1], -t[0]]) * sign
|
|
for k in range(0, len(s['x']), 20):
|
|
if not np.isnan(s['top'][k]) and abs(s['x'][k]) < 12:
|
|
e, n = p + nrm * s['x'][k]
|
|
z = dem(e, n)
|
|
if not np.isnan(z): d.append(s['top'][k] - z)
|
|
return np.array(d)
|
|
d_r, d_l = diffs(1), diffs(-1)
|
|
sign = 1 if np.nanmean(np.abs(d_r)) <= np.nanmean(np.abs(d_l)) else -1
|
|
dd = d_r if sign == 1 else d_l
|
|
|
|
out_secs = []
|
|
for s in sections:
|
|
p, t = at(s['ch'])
|
|
keep = slice(None, None, 5) # 0.5 m spacing for the browser
|
|
out_secs.append({
|
|
'ch': s['ch'], 'e': round(float(p[0]), 3), 'n': round(float(p[1]), 3),
|
|
'dir': [round(float(t[0]), 5), round(float(t[1]), 5)], 'side': sign,
|
|
'x': s['x'][keep].tolist(), 'ngl': np.round(s['ngl'][keep], 3).tolist(),
|
|
'top': [None if np.isnan(v) else round(float(v), 3) for v in s['top'][keep]],
|
|
'bot': [None if np.isnan(v) else round(float(v), 3) for v in s['bot'][keep]],
|
|
'cutArea': round(s['cutArea'], 2), 'fillArea': round(s['fillArea'], 2),
|
|
})
|
|
vols = []
|
|
for a, b in zip(sections, sections[1:]):
|
|
L = b['ch'] - a['ch']
|
|
whole = lambda c: int(c) if float(c).is_integer() else c # keys read earth:240-250
|
|
vols.append({'from': whole(a['ch']), 'to': whole(b['ch']), 'cut': round((a['cutArea'] + b['cutArea']) / 2 * L, 1),
|
|
'fill': round((a['fillArea'] + b['fillArea']) / 2 * L, 1)})
|
|
res = {
|
|
'source': os.path.basename(pdf_path), 'method': 'end-area between sections; natural ground vs design underside (formation)',
|
|
'chainageOffset': 0.0, 'chainageEnd': round(ch_end, 3),
|
|
'sections': out_secs, 'volumes': vols,
|
|
'alignment': [[round(c, 1), *[round(float(v), 3) for v in at(c)[0]]] for c in np.arange(0, ch_end + 0.01, 2.0)],
|
|
'totals': {'cut': round(sum(v['cut'] for v in vols), 1), 'fill': round(sum(v['fill'] for v in vols), 1)},
|
|
'validation': {'samples': int(len(dd)), 'meanAbs': round(float(np.nanmean(np.abs(dd))), 3),
|
|
'median': round(float(np.nanmedian(dd)), 3), 'p90Abs': round(float(np.nanpercentile(np.abs(dd), 90)), 3),
|
|
'note': 'design surface read from cross-sections vs. design surface from 0.1 m road contours (DWG)'},
|
|
}
|
|
json.dump(res, open(os.path.join(data_dir, 'earthworks.json'), 'w'), separators=(',', ':'))
|
|
print(f"sections {len(out_secs)} ch {out_secs[0]['ch']}..{out_secs[-1]['ch']} scales {sorted({s['scale'] for s in sections})[:4]}")
|
|
print(f"alignment length {ch_end:.3f} m (drawing says 674.412)")
|
|
print('totals', res['totals'], 'validation', res['validation'])
|
|
for s in out_secs[:: max(1, len(out_secs) // 8)]:
|
|
print(f" ch {s['ch']:7.1f} cut {s['cutArea']:6.1f} m² fill {s['fillArea']:6.1f} m²")
|
|
|
|
|
|
if __name__ == '__main__':
|
|
main(*sys.argv[1:4])
|