"""3DBAG CityJSON LoD2.2 acquisition and source-preserving geometry decoding."""
from __future__ import annotations
import asyncio
from collections import deque
from urllib.parse import parse_qs, urljoin, urlsplit
import httpx
import numpy as np
import shapely
from pyproj import Transformer
from shapely.geometry import Polygon
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.models import BuildingFeature, BuildingGeometry
from openworld_radio_twin.providers.base import BuildingProviderError, BuildingQueryResult
from openworld_radio_twin.providers.cityjson import geometry_surfaces, triangulate_surface
API_URL = "https://api.3dbag.nl/collections/pand/items"
PAGE_SIZE = 100 # The service caps ``limit`` at 100 items per page.
PAGE_CONCURRENCY = 4
ATTRIBUTION = "© 3DBAG by tudelft3d and 3DGI"
LICENSE_URL = "https://creativecommons.org/licenses/by/4.0/"
COPYRIGHT_URL = "https://docs.3dbag.nl/en/copyright/"
RD_TO_WGS84 = Transformer.from_crs("EPSG:28992", "EPSG:4326", always_xy=True)
WGS84_TO_RD = Transformer.from_crs("EPSG:4326", "EPSG:28992", always_xy=True)
[docs]
def decode_feature(feature: dict, metadata: dict) -> list[BuildingFeature]:
"""Decode only LoD2.2, retaining source Z in NAP and semantic triangles.
The parent LoD0 footprint is not emitted as another building. Unsupported
or invalid geometry raises instead of silently becoming a LoD1 extrusion.
"""
reference = metadata.get("metadata", {}).get("referenceSystem", "")
if not reference.endswith("/7415") and reference != "EPSG:7415":
raise ValueError(f"3DBAG response has unsupported CRS: {reference}")
transform = metadata["transform"]
vertices = np.asarray(feature["vertices"], dtype=np.float64)
vertices = vertices * np.asarray(transform["scale"]) + np.asarray(transform["translate"])
if vertices.ndim != 2 or vertices.shape[1] != 3 or not np.isfinite(vertices).all():
raise ValueError("3DBAG response has invalid vertices")
objects = feature["CityObjects"]
result = []
for object_id, obj in objects.items():
candidates = [g for g in obj.get("geometry", []) if str(g.get("lod")) == "2.2"]
if not candidates:
continue
parents = obj.get("parents", [])
parent_id = parents[0] if parents else object_id
attributes = {
**objects.get(parent_id, {}).get("attributes", {}),
**obj.get("attributes", {}),
}
triangles, surface_types = [], []
zero_area_faces = 0
source_faces = 0
for geometry in candidates:
for rings, surface_type in geometry_surfaces(geometry):
source_faces += 1
faces = triangulate_surface(vertices, rings)
zero_area_faces += not faces
triangles.extend(faces)
surface_types.extend([surface_type] * len(faces))
if not triangles:
raise ValueError(f"{object_id} has no LoD2.2 faces")
# Project the complete shell, including overhangs, for source bounds.
projections = shapely.polygons(vertices[np.asarray(triangles, dtype=np.intp)][:, :, :2])
footprint = shapely.union_all(projections[shapely.area(projections) > 1e-9])
if not isinstance(footprint, Polygon) or not footprint.is_valid:
raise ValueError(f"{object_id} has a disconnected or invalid footprint")
used = sorted({index for face in triangles for index in face})
remap = {old: new for new, old in enumerate(used)}
points = vertices[used]
ground = attributes.get("b3_h_maaiveld")
if ground is None:
raise ValueError(f"{object_id} has no source ground elevation")
ground = float(ground)
minimum = float(points[:, 2].min())
maximum = float(points[:, 2].max())
def geographic_ring(ring):
x, y = np.asarray(ring.coords).T
lon, lat = RD_TO_WGS84.transform(x, y)
return [[float(a), float(b)] for a, b in zip(lon, lat, strict=True)]
result.append(
BuildingFeature(
feature_id=f"3dbag:{object_id}",
source="3DBAG",
feature_type=("building_part" if obj.get("type") == "BuildingPart" else "building"),
source_dataset_id=object_id,
parent_building_id=parent_id,
source_record_id=parent_id,
height_m=maximum - minimum,
min_height_m=max(0.0, minimum - ground),
height_source="3dbag_lod22_geometry",
coordinates=geographic_ring(footprint.exterior),
holes=[geographic_ring(ring) for ring in footprint.interiors],
geometry=BuildingGeometry(
crs="EPSG:7415",
vertical_datum="NAP",
lod="2.2",
ground_elevation_m=ground,
vertices=points.tolist(),
triangles=[tuple(remap[index] for index in face) for face in triangles],
surface_types=surface_types,
),
source_attributes={
**attributes,
"attribution": ATTRIBUTION,
"license": LICENSE_URL,
"copyright_url": COPYRIGHT_URL,
"source_url": f"{API_URL}/{parent_id}",
"source_vertical_datum": "NAP",
"source_surface_count": source_faces,
"omitted_zero_area_surface_count": zero_area_faces,
},
)
)
if not result:
raise ValueError("3DBAG feature has no LoD2.2 geometry")
return result
def _next_link(url: str, page: dict) -> str:
links = [link["href"] for link in page.get("links", []) if link.get("rel") == "next"]
return urljoin(url, links[0]) if links else ""
def _offset_url(url: str, offset: int) -> str:
parsed = urlsplit(url)
params = {key: values[-1] for key, values in parse_qs(parsed.query).items()}
params["offset"] = str(offset)
return str(httpx.URL(f"{parsed.scheme}://{parsed.netloc}{parsed.path}", params=params))
[docs]
class ThreeDBagProvider:
"""Read paginated official API responses with bounded request concurrency."""
def __init__(self, client: httpx.AsyncClient) -> None:
self.client = client
self._lock = asyncio.Lock()
[docs]
async def buildings(
self,
latitude: float,
longitude: float,
radius_m: float,
limit: int,
) -> BuildingQueryResult:
if limit < 1:
raise ValueError("Building limit must be positive")
frame = LocalFrame.at(longitude, latitude)
corners = frame.corners(-radius_m, -radius_m, radius_m, radius_m)
x, y = WGS84_TO_RD.transform(*zip(*corners, strict=True))
bbox = ",".join(str(value) for value in (min(x), min(y), max(x), max(y)))
page_size = min(PAGE_SIZE, limit)
url = str(httpx.URL(API_URL, params={"bbox": bbox, "limit": page_size}))
buildings, warnings, seen, visited = [], [], set(), set()
def collect(page: dict) -> bool:
"""Decode one page; return whether the building limit has been reached."""
for feature in page["features"]:
feature_id = feature.get("id")
if feature_id in seen:
continue
seen.add(feature_id)
try:
decoded = decode_feature(feature, page["metadata"])
except (KeyError, TypeError, ValueError) as exc:
warnings.append(f"Skipped 3DBAG feature {feature_id}: {exc}")
continue
buildings.extend(decoded[: limit - len(buildings)])
if len(buildings) >= limit:
warnings.append(f"3DBAG building limit ({limit}) reached.")
return True
return False
async with self._lock:
try:
page = await self._page(url, visited)
done = await asyncio.to_thread(collect, page)
next_url = _next_link(url, page)
matched = page.get("numberMatched")
if not done and next_url and isinstance(matched, int):
# The service reports the total and accepts ``offset``: keep a few
# page requests in flight, at most ``limit`` items since a building
# needs at least one, and decode the pages in offset order. The
# last page's link continues the sequential walk.
window: deque[tuple[int, asyncio.Task[dict]]] = deque()
async def decode_oldest() -> None:
nonlocal done, next_url
offset, task = window.popleft()
page = await task
done = await asyncio.to_thread(collect, page)
next_url = _next_link(_offset_url(url, offset), page)
try:
for offset in range(page_size, min(matched, limit), page_size):
if done:
break
request = self._page(_offset_url(url, offset), visited)
window.append((offset, asyncio.create_task(request)))
if len(window) == PAGE_CONCURRENCY:
await decode_oldest()
while window and not done:
await decode_oldest()
finally:
for _, task in window:
task.cancel()
while not done and next_url:
page = await self._page(next_url, visited)
done = await asyncio.to_thread(collect, page)
next_url = _next_link(next_url, page)
except (httpx.HTTPError, KeyError, TypeError, ValueError) as exc:
raise BuildingProviderError(f"3DBAG query failed: {exc}") from exc
if seen and not buildings:
raise BuildingProviderError(
"3DBAG returned no usable LoD2.2 buildings: " + "; ".join(warnings[:3])
)
zero_area_count = sum(
int(building.source_attributes["omitted_zero_area_surface_count"])
for building in buildings
)
if zero_area_count:
warnings.append(
f"Omitted {zero_area_count} collinear zero-area source faces; "
"all nondegenerate source faces of those parts were retained."
)
return BuildingQueryResult(buildings=buildings, provider="3DBAG", warnings=warnings)
async def _page(self, url: str, visited: set[str]) -> dict:
parsed = urlsplit(url)
if (parsed.scheme, parsed.netloc, parsed.path) != (
"https",
"api.3dbag.nl",
"/collections/pand/items",
):
raise ValueError("3DBAG pagination left the official endpoint")
if url in visited:
raise ValueError("3DBAG pagination contains a cycle")
visited.add(url)
response = await self.client.get(url, timeout=60)
response.raise_for_status()
return response.json()
