323 lines
18 KiB
Python
323 lines
18 KiB
Python
"""Build Cedar Gate's Map 2 model from the C3277 road/services DXFs and Chris's ground contours.
|
|
|
|
Requires numpy, scipy, shapely>=2.1 and ezdxf. See README.md for conversion/build commands.
|
|
All coordinates use local east/north metres; CAD road x=-Lo Y, y=-Lo X.
|
|
Chris's ground drawing stores positive Lo Y/X and is reversed on both axes.
|
|
"""
|
|
import argparse
|
|
import base64
|
|
from collections import defaultdict
|
|
import hashlib
|
|
import io
|
|
import json
|
|
import math
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
from scipy.interpolate import LinearNDInterpolator
|
|
from scipy.spatial import cKDTree
|
|
from shapely import make_valid, voronoi_polygons
|
|
from shapely.geometry import LineString, MultiPoint, Point, Polygon, box
|
|
from shapely.ops import unary_union
|
|
|
|
|
|
def load_dxf(path):
|
|
"""Open a LibreDWG DXF after repairing its multiline text and unsupported annotation payloads, not geometry."""
|
|
import ezdxf
|
|
lines = Path(path).read_text(errors="replace").splitlines()
|
|
pairs = []
|
|
i = 0
|
|
while i < len(lines):
|
|
try:
|
|
code = int(lines[i])
|
|
except ValueError:
|
|
if not pairs or pairs[-1][0] not in (1, 3):
|
|
raise ValueError(f"Unexpected DXF group at line {i+1}")
|
|
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))
|
|
return ezdxf.read(io.StringIO("\n".join(clean) + "\n"))
|
|
|
|
|
|
def read_dxf(path):
|
|
from ezdxf.path import make_path, from_hatch
|
|
doc = load_dxf(path)
|
|
omitted = []
|
|
def flatten(entities):
|
|
for entity in entities:
|
|
if entity.dxftype() == "INSERT":
|
|
if entity.dxf.name not in doc.blocks:
|
|
omitted.append(entity.dxf.name)
|
|
continue
|
|
yield from flatten(entity.virtual_entities())
|
|
else:
|
|
yield entity
|
|
layers = defaultdict(list)
|
|
for entity in flatten(doc.modelspace()):
|
|
layer, kind = entity.dxf.layer, entity.dxftype()
|
|
if not layer.startswith(("RDS_LO$", "SER$", "FRL_CONTOURS$", "ARCH$0$ERF", "ARCH$0$BLDG_OUTLINE", "BOUNDARY_NEW", "Chris - Cedar Gate")):
|
|
continue
|
|
row = {"type": kind}
|
|
if kind in ("LINE", "ARC", "LWPOLYLINE", "POLYLINE", "CIRCLE"):
|
|
row["pts"] = [list(p) for p in make_path(entity).flattening(.025)]
|
|
elif kind == "HATCH":
|
|
row["rings"] = [[list(p) for p in path.flattening(.025)] for path in from_hatch(entity)]
|
|
elif kind in ("TEXT", "MTEXT"):
|
|
row["text"] = entity.plain_text() if kind == "MTEXT" else entity.dxf.text
|
|
row["pos"] = list(entity.dxf.insert)
|
|
else:
|
|
continue
|
|
layers[layer].append(row)
|
|
return dict(layers), sorted(set(omitted))
|
|
|
|
|
|
def polygons(geometry):
|
|
if geometry.is_empty:
|
|
return []
|
|
if geometry.geom_type == "Polygon":
|
|
return [geometry]
|
|
if hasattr(geometry, "geoms"):
|
|
return [p for g in geometry.geoms for p in polygons(g)]
|
|
return []
|
|
|
|
|
|
def hatch_geometry(rows):
|
|
"""Even/odd nesting handles holes and disjoint outer hatch boundaries."""
|
|
shapes = []
|
|
for row in rows:
|
|
rings = [make_valid(Polygon([p[:2] for p in ring])) for ring in row.get("rings", []) if len(ring) >= 3]
|
|
shape = Polygon()
|
|
for ring in rings:
|
|
shape = shape.symmetric_difference(ring)
|
|
shapes.append(shape)
|
|
return unary_union(shapes)
|
|
|
|
|
|
def alignment(rows, start, expected):
|
|
"""Trim the PI tangents at the actual circular-arc endpoints, then order the pieces."""
|
|
lines = [x["pts"] for x in rows if x["type"] == "LINE"]
|
|
arcs = [x["pts"] for x in rows if x["type"] == "ARC"]
|
|
segments = [[line[0][:2], line[-1][:2]] for line in lines]
|
|
for arc in arcs:
|
|
matched = []
|
|
for endpoint in (arc[0][:2], arc[-1][:2]):
|
|
matches = [i for i, line in enumerate(lines) if LineString(line).distance(Point(endpoint)) < .002]
|
|
if len(matches) != 1:
|
|
raise ValueError("Arc endpoint does not identify one PI tangent")
|
|
matched.append((matches[0], endpoint))
|
|
(a, pa), (b, pb) = matched
|
|
distance, ia, ib = min((math.dist(lines[a][i][:2], lines[b][j][:2]), i, j) for i in (0, -1) for j in (0, -1))
|
|
if distance > .01:
|
|
raise ValueError("Curve PI tangents do not meet")
|
|
segments[a][ia] = pa
|
|
segments[b][ib] = pb
|
|
segments.extend([[p[:2] for p in arc] for arc in arcs])
|
|
current, points = start[:2], [start[:2]]
|
|
while segments:
|
|
distance, index, reverse = min((Point(current).distance(Point(p)), i, rev) for i, segment in enumerate(segments) for rev, p in ((False, segment[0]), (True, segment[-1])))
|
|
if distance > .01:
|
|
raise ValueError(f"Disconnected centreline: {distance:.3f} m")
|
|
segment = segments.pop(index)
|
|
if reverse:
|
|
segment.reverse()
|
|
points.extend(segment[1:]); current = segment[-1]
|
|
result = LineString(points)
|
|
if abs(result.length - expected) > .12:
|
|
raise ValueError(f"Alignment {result.length:.3f} m differs from drawing {expected} m")
|
|
return result
|
|
|
|
|
|
def surface(contours, origin, source, cell=1):
|
|
samples = []
|
|
for contour in contours:
|
|
line = LineString(contour["pts"])
|
|
for distance in np.linspace(0, line.length, max(2, math.ceil(line.length/.75)+1)):
|
|
p = line.interpolate(distance)
|
|
samples.append((p.x-origin[0], p.y-origin[1], contour["z"]))
|
|
samples = np.array(samples)
|
|
# Identical XY vertices are averaged to avoid triangulation ambiguity.
|
|
xy, inverse = np.unique(samples[:, :2].round(4), axis=0, return_inverse=True)
|
|
z = np.bincount(inverse, weights=samples[:, 2])/np.bincount(inverse)
|
|
lower = np.floor(xy.min(axis=0))-2; upper = np.ceil(xy.max(axis=0))+2
|
|
nx, ny = (np.ceil((upper-lower)/cell).astype(int)+1).tolist()
|
|
xx, yy = np.meshgrid(lower[0]+np.arange(nx)*cell, lower[1]+np.arange(ny)*cell)
|
|
grid = LinearNDInterpolator(xy, z)(xx, yy)
|
|
distance, _ = cKDTree(xy).query(np.column_stack([xx.ravel(), yy.ravel()]))
|
|
grid[distance.reshape(grid.shape) > 15] = np.nan # no long-range extrapolation
|
|
zmin = math.floor(float(np.nanmin(grid))*1000)/1000
|
|
if (float(np.nanmax(grid))-zmin)*1000 >= 32767:
|
|
raise ValueError("Surface exceeds the renderer's int16 elevation range")
|
|
encoded = np.full(grid.shape, -1, dtype="<i2")
|
|
valid = np.isfinite(grid)
|
|
encoded[valid] = np.round((grid[valid]-zmin)*1000).astype("<i2")
|
|
return {"x0": float(lower[0]), "y0": float(lower[1]), "cell": cell, "nx": nx, "ny": ny, "zmin": zmin,
|
|
"z": base64.b64encode(encoded.tobytes()).decode(), "encoding": "int16 little-endian, mm above zmin, -1 = outside supplied contour coverage, row-major",
|
|
"source": source, "check": {"pegs": 0, "medianAbs": 0, "p90Abs": 0}}, grid
|
|
|
|
|
|
def layer(name, thickness, unit, color):
|
|
return {"id": name.lower().replace(" ", "-"), "name": name, "thk": thickness, "unit": unit, "color": color}
|
|
|
|
|
|
def build(raw, output):
|
|
road, services = raw["road"], raw["services"]
|
|
prefix = "RDS_LO$0$CC_RDS_NEW_"
|
|
internal = road[prefix+"INT_CL"]
|
|
external = road[prefix+"EXT_CL"]
|
|
ext_lines = [x for x in external if x["type"] == "LINE"]
|
|
ext_arcs = [x for x in external if x["type"] == "ARC"]
|
|
routes = [
|
|
("cedar-gate", "Cedar Gate Road", alignment(internal, internal[0]["pts"][0], 385.41)),
|
|
("libertas", "Libertas Avenue", alignment(ext_lines[:5]+[ext_arcs[1]], ext_lines[0]["pts"][0], 173.12)),
|
|
("furrow", "Furrow Road", alignment(ext_lines[5:9]+[ext_arcs[0]], ext_lines[5]["pts"][0], 237.64)),
|
|
("access", "Access road", alignment(ext_lines[9:], ext_lines[9]["pts"][0], 36.75)),
|
|
]
|
|
origin = (-66720., -2850540.)
|
|
local = lambda pts: [[round(p[0]-origin[0], 4), round(p[1]-origin[1], 4)] for p in pts]
|
|
geometries = {
|
|
"internal_paving": hatch_geometry(road[prefix+"INT_HATCH"]),
|
|
"asphalt": hatch_geometry(road[prefix+"EXT_HATCH_ASPHALT"]),
|
|
"external_paving": hatch_geometry(road[prefix+"EXT_HATCH_PAVING"]),
|
|
"walkway": hatch_geometry(road[prefix+"EXT_WW_HATCH"]),
|
|
}
|
|
# Partition by nearest centreline samples (1 m spacing). A union of Voronoi cells
|
|
# gives a complete, disjoint allocation, including at intersections and the loop closure.
|
|
samples, owners, sections = [], [], []
|
|
for route_id, name, line in routes:
|
|
for lo in range(0, math.ceil(line.length), 20):
|
|
hi = min(lo+20, line.length)
|
|
key = f"{route_id}-ch{lo}-{min(lo+20, math.ceil(line.length))}"
|
|
sections.append({"key": key, "name": name, "from": lo, "to": round(hi, 2)})
|
|
for ch in np.arange(lo+.5, hi, 1):
|
|
point = line.interpolate(ch); samples.append((point.x, point.y)); owners.append(key)
|
|
extent = unary_union(list(geometries.values())).envelope.buffer(25)
|
|
cells = voronoi_polygons(MultiPoint(samples), extend_to=extent, ordered=True)
|
|
grouped = defaultdict(list)
|
|
for key, cell in zip(owners, cells.geoms):
|
|
grouped[key].append(cell)
|
|
bands = {key: unary_union(value) for key, value in grouped.items()}
|
|
areas = []
|
|
for kind, geometry in geometries.items():
|
|
for key, band in bands.items():
|
|
for polygon in polygons(geometry.intersection(band)):
|
|
if polygon.area <= 1e-7:
|
|
continue
|
|
rings = [polygon.exterior, *polygon.interiors]
|
|
areas.append({"kind": kind, "section": key, "area": polygon.area, "rings": [local(r.coords) for r in rings]})
|
|
allocated = sum(x["area"] for x in areas if x["kind"] == kind)
|
|
if abs(allocated-geometry.area) > .01:
|
|
raise ValueError(f"Area allocation mismatch: {kind} {allocated} vs {geometry.area}")
|
|
# Kerbs are single centre/face lines in this drawing; do not apply Bushveld's half-linework rule.
|
|
kerbs = []
|
|
for suffix in ("INT_KERBS-300", "EXT_KERBS-500", "EXT_KERBS_FIG.7_F"):
|
|
for row in road[prefix+suffix]:
|
|
if "pts" not in row:
|
|
continue
|
|
line = LineString([p[:2] for p in row["pts"]])
|
|
for key, band in bands.items():
|
|
clipped = line.intersection(band)
|
|
parts = list(clipped.geoms) if hasattr(clipped, "geoms") else [clipped]
|
|
for part in parts:
|
|
if part.geom_type == "LineString" and part.length > .01:
|
|
kerbs.append({"section": key, "length": part.length, "pts": local(part.coords)})
|
|
context = []
|
|
for suffix, kind in [("INT_EDGE", "tar_edge"), ("EXT_EDGE", "tar_edge"), ("INT_RM", "marking"), ("EXT_RM", "marking"), ("EXT_RES", "reserve")]:
|
|
for row in road.get(prefix+suffix, []):
|
|
if "pts" in row:
|
|
context.append({"kind": kind, "pts": local(row["pts"])})
|
|
for name in ("ARCH$0$ERF", "ARCH$0$BLDG_OUTLINE", "BOUNDARY_NEW"):
|
|
for row in road.get(name, []):
|
|
if "pts" in row:
|
|
context.append({"kind": "reserve" if name=="BOUNDARY_NEW" else "other_road", "pts": local(row["pts"])})
|
|
contours = []
|
|
for name, rows in road.items():
|
|
if not name.startswith("FRL_CONTOURS$"):
|
|
continue
|
|
for row in rows:
|
|
if "pts" in row and len(row["pts"]) > 1:
|
|
contours.append({"z": row["pts"][0][2], "major": "MCONTOUR" in name, "pts": [p[:2] for p in row["pts"]]})
|
|
ground = []
|
|
for rows in raw["ground"].values():
|
|
for row in rows:
|
|
if "pts" in row and len(row["pts"]) > 1:
|
|
ground.append({"z": row["pts"][0][2], "major": False, "pts": [[-p[0], -p[1]] for p in row["pts"]]})
|
|
sf, grid = surface(contours, origin, "C3277-910-001 Rev A: FRL_CONTOURS 0.1 m finished-level contours; triangulated at 1 m. Not an as-built survey.")
|
|
ground_surface, _ = surface(ground, origin, "Cedar Gate Contours Chris.dwg: ground-level contours; positive Lo Y/X reversed to CAD east/north. Survey date/datum confirmation outstanding.")
|
|
builds = {
|
|
"internal_paving": {"name": "Internal road paving", "color": 0xb9b2a6, "layers": [
|
|
layer("In-situ roadbed G9 150mm", .15, "m²", 0x7d6446), layer("Selected G7 150mm", .15, "m³", 0xae9466),
|
|
layer("C4 subbase G5 cement 150mm", .15, "m³", 0xa6a49a), layer("Bedding sand 20mm", .02, "m³", 0xd8c48f), layer("80mm interlocking pavers", .08, "m²", 0xb9b2a6)]},
|
|
"asphalt": {"name": "External asphalt surface", "color": 0x34373b, "layers": [layer("Asphalt 30mm", .03, "m²", 0x34373b)]},
|
|
"external_paving": {"name": "Access road paving surface", "color": 0xb9b2a6, "layers": [layer("Bedding sand 20mm", .02, "m³", 0xd8c48f), layer("80mm interlocking pavers", .08, "m²", 0xb9b2a6)]},
|
|
"walkway": {"name": "Pedestrian walkway surface", "color": 0xdeb261, "layers": [layer("Walkway area", 0, "m²", 0xdeb261)]},
|
|
}
|
|
notes = [
|
|
"Cedar Gate = Equestria Ext. 289, C3277. Nellmapius Ext. 25 drawings are a separate project and are excluded.",
|
|
"20 m planning sections follow the four drawn centrelines. Nearest-centreline allocation sampled at 1 m; hatch areas are conserved without overlaps.",
|
|
"Internal road layers follow C3277-913-001 Rev A, typical section 1b/3/4/5: pavers, bedding, C4 subbase, selected G7 and roadbed. Confirm adopted contract detail.",
|
|
"External asphalt, access paving and walkway quantities cover the shown surface only; unverified lower layers are excluded.",
|
|
"Kerbs are measured once on the 300 mm and 500 mm lines and Fig. 7 front face; Fig. 7 back faces/channels/edge beams are excluded to avoid duplicate length.",
|
|
"Service routes show plan positions draped over the finished-level surface. Pipe inverts, burial depths, diameters and construction progress are not inferred.",
|
|
"Chris ground contours have separate coverage and are available as Natural ground. They are not substituted for the finished-level contours.",
|
|
"Bulk cut/fill quantities are unavailable until the road cross-sections and formation levels are reconciled; no zero-volume claim is made.",
|
|
"No setting-out peg schedule was supplied; surface-to-peg validation is unavailable. Benchmarks are not treated as design pegs.",
|
|
]
|
|
model_road = {"project": "Cedar Gate — Equestria Ext. 289 (C3277)", "origin": {"cadX": origin[0], "cadY": origin[1]},
|
|
"builds": builds, "sections": sections, "areas": areas, "kerbs": kerbs, "kerbQuantityFactor": 1,
|
|
"context": context, "contours": [{**c, "pts": local(c["pts"])} for c in contours], "pegs": [], "water": [], "buffer": [], "notes": notes}
|
|
plan_routes = []
|
|
for service, suffixes in [("stormwater", ["SW_NEW_INT", "SW_NEW_EXT"]), ("sewer", ["SR_NEW_INT", "SR_NEW_EXT"]), ("water", ["WR_NEW_INT", "WR_NEW_EXT"])]:
|
|
for suffix in suffixes:
|
|
for i, row in enumerate(services.get("SER$0$CC_"+suffix, [])):
|
|
if "pts" in row and len(row["pts"]) > 1:
|
|
pts = row["pts"]
|
|
plan_routes.append({"id": f"{suffix.lower()}-{i+1}", "service": service, "pts": local(pts),
|
|
"length": LineString(pts).length, "source": "C3277 services Rev B — SER$0$CC_"+suffix})
|
|
sv = {"stormwater": {"structures": [], "pipes": [], "notes": [notes[5]]}, "sewer": {"manholes": [], "pipes": [], "notes": [notes[5]]},
|
|
"water": {"points": [], "pipeLabels": {}}, "holds": [], "planRoutes": plan_routes}
|
|
earth = {"source": notes[7], "method": "Not quantified", "chainageOffset": 0, "chainageEnd": round(routes[0][2].length, 2),
|
|
"sections": [], "volumes": [], "alignment": [[round(ch, 3), *local([routes[0][2].interpolate(ch).coords[0]])[0]] for ch in np.arange(0, routes[0][2].length, 2)],
|
|
"totals": {"cut": 0, "fill": 0}, "validation": {"samples": 0, "meanAbs": 0, "median": 0, "p90Abs": 0, "note": notes[7]}, "groundSurface": ground_surface}
|
|
schedules = {"signs": [], "signCounts": {}, "markingLabels": [], "markingCounts": {}, "studs": [], "studCounts": {},
|
|
"markings": {}, "paintAreas": {}, "waterFittings": {}, "notes": ["Signs, studs and fittings are not quantified in this initial model."]}
|
|
output.mkdir(parents=True, exist_ok=True)
|
|
for name, data in [("bushveld", model_road), ("services", sv), ("earthworks", earth), ("schedules", schedules), ("surface", sf)]:
|
|
(output/f"{name}.json").write_text(json.dumps(data, separators=(",", ":"), allow_nan=False))
|
|
report = {"project": model_road["project"], "alignments": {name: {"length_m": round(line.length, 3)} for _, name, line in routes},
|
|
"sections": len(sections), "areas_m2": {k: round(v.area, 3) for k, v in geometries.items()}, "kerbs_m": round(sum(x["length"] for x in kerbs), 3),
|
|
"design_contours": len(contours), "ground_contours": len(ground), "service_plan_routes": len(plan_routes),
|
|
"service_plan_lengths_m": {s: round(sum(x["length"] for x in plan_routes if x["service"]==s), 3) for s in ("stormwater", "sewer", "water")},
|
|
"design_surface": {"nx": sf["nx"], "ny": sf["ny"], "valid_cells": int(np.isfinite(grid).sum())}, "notes": notes,
|
|
"sha256": {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in output.glob("*.json") if p.name != "validation.json"}}
|
|
(output/"validation.json").write_text(json.dumps(report, indent=2))
|
|
print(json.dumps({k:v for k,v in report.items() if k not in ("notes", "sha256")}, indent=2))
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description=__doc__)
|
|
parser.add_argument("dxf_folder", type=Path)
|
|
parser.add_argument("output", type=Path)
|
|
parser.add_argument("--raw-cache", type=Path, help="Optional already-extracted CAD geometry JSON for local rebuilds")
|
|
args = parser.parse_args()
|
|
if args.raw_cache:
|
|
raw = json.loads(args.raw_cache.read_text())
|
|
else:
|
|
raw = {}
|
|
for key, pattern in [("road", "*910-001*.dxf"), ("services", "*920to940*REVB*.dxf"), ("ground", "*Cedar Gate*Chris.dxf")]:
|
|
raw[key], omitted = read_dxf(next(args.dxf_folder.glob(pattern)))
|
|
if omitted:
|
|
print(f"{key}: unresolved CAD blocks (excluded): {omitted}")
|
|
build(raw, args.output)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|