36da7f9f20
Road areas (5,221 m²) in 20 m sections along the ring roads and the Igwali Street spur, 1.74 km of kerb measured once, stand lines, and the four survey benchmarks with an indicative ground plane. No services or ground levels are placed; the supplied survey and sewer sheets cannot yield them. Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
215 lines
14 KiB
Python
215 lines
14 KiB
Python
"""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])
|