#!/usr/bin/env python3 """Clip CityGML LoD2 buildings per terrain tile and export as OBJ meshes. Inputs - raw_3dgeb_lod2/LoD2___.gml : source building tiles - export_unity/tile_index.csv : manifest produced by export_heightmaps.py Outputs - export_unity/buildings_obj/.obj : clipped building geometries per tile """ from __future__ import annotations import csv import os from pathlib import Path from osgeo import gdal RAW_GML_DIR = Path("raw_3dgeb_lod2") TILE_INDEX = Path("export_unity/tile_index.csv") OUT_DIR = Path("export_unity/buildings_obj") OUTPUT_FORMAT = "OBJ" # glTF driver is not always available; OBJ is widely supported. gdal.UseExceptions() OUT_DIR.mkdir(parents=True, exist_ok=True) def gml_path_for_tile(tile_id: str) -> Path: """Map a DGM tile id (dgm1_utm_easting_northing) to the corresponding CityGML file.""" parts = tile_id.split("_") if len(parts) != 4 or parts[0] != "dgm1": raise ValueError(f"Unexpected tile_id format: {tile_id}") utm_zone, easting, northing = parts[1:] return RAW_GML_DIR / f"LoD2_{utm_zone}_{easting}_{northing}.gml" def open_dataset(path: str, purpose: str): """Open a dataset and fail fast with context.""" try: ds = gdal.OpenEx(path) except RuntimeError as exc: raise SystemExit(f"{purpose}: {exc}") from exc if ds is None: raise SystemExit(f"{purpose}: GDAL returned None for {path}") return ds def main() -> None: if not TILE_INDEX.exists(): raise SystemExit(f"Tile index missing: {TILE_INDEX}. Run export_heightmaps.py first.") written = 0 skipped = 0 with TILE_INDEX.open(newline="", encoding="utf-8") as f: reader = csv.DictReader(f) for row in reader: try: tile_id = row["tile_id"] xmin = float(row["xmin"]) ymin = float(row["ymin"]) xmax = float(row["xmax"]) ymax = float(row["ymax"]) except (KeyError, ValueError) as exc: print(f"Skipping malformed row {row}: {exc}") skipped += 1 continue gml_path = gml_path_for_tile(tile_id) if not gml_path.exists(): print(f"Missing CityGML for {tile_id}: {gml_path}") skipped += 1 continue out_path = OUT_DIR / f"{tile_id}.obj" src_ds = open_dataset(str(gml_path), f"Could not open {gml_path}") opts = gdal.VectorTranslateOptions( format=OUTPUT_FORMAT, spatFilter=(xmin, ymin, xmax, ymax), layerName="Building", # matches GML layer name ) try: gdal.VectorTranslate(destNameOrDestDS=str(out_path), srcDS=src_ds, options=opts) except RuntimeError as exc: print(f"VectorTranslate failed for {tile_id}: {exc}") skipped += 1 continue written += 1 print(f"Wrote {out_path}") print(f"Summary: wrote {written} building tiles; skipped {skipped}.") if skipped: raise SystemExit(1) if __name__ == "__main__": main()