"""Berlin's official tiled CityGML LoD2 buildings, retaining source geometry."""
from __future__ import annotations
import asyncio
import hashlib
import io
import math
import re
import zipfile
from urllib.parse import urlsplit
from xml.etree.ElementTree import ParseError
import httpx
import numpy as np
from defusedxml import ElementTree as ET
from defusedxml.common import DefusedXmlException
from pyproj import Transformer
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.providers.base import BuildingProviderError, BuildingQueryResult
from openworld_radio_twin.providers.cityjson import triangulate_surface
from openworld_radio_twin.providers.native_mesh import mesh_feature
INDEX_URL = "https://gdi.berlin.de/data/a_lod2/atom/0.atom"
ATTRIBUTION = "Geoportal Berlin / 3D-Gebäudemodelle LoD2"
LICENSE_URL = "https://www.govdata.de/dl-de/zero-2-0"
CRS_NAME = "urn:adv:crs:ETRS89_UTM33*DE_DHHN2016_NH"
GML = "{http://www.opengis.net/gml}"
ATOM = "{http://www.w3.org/2005/Atom}"
WGS84_TO_UTM = Transformer.from_crs(4326, 25833, always_xy=True)
MAX_TILE_BYTES = 128 * 1024 * 1024
MAX_XML_BYTES = 512 * 1024 * 1024
def _name(element):
return element.tag.rsplit("}", 1)[-1]
[docs]
def tile_index(content: bytes) -> tuple[dict[tuple[int, int], str], str]:
"""Use only tile links actually advertised by the official Atom service."""
root = ET.fromstring(content)
tiles = {}
for link in root.iter(f"{ATOM}link"):
url = link.get("href", "")
parsed = urlsplit(url)
match = re.fullmatch(r"/data/a_lod2/atom/LoD2_(\d+)_(\d+)\.zip", parsed.path)
if link.get("rel") != "section":
continue
if not match or parsed.scheme != "https" or parsed.netloc != "gdi.berlin.de":
raise ValueError("Berlin tile link left the official download service")
tiles[tuple(map(int, match.groups()))] = url
if not tiles:
raise ValueError("Berlin Atom index contains no LoD2 tiles")
return tiles, root.findtext(f"{ATOM}updated", default="unknown")
def _attributes(node):
attributes = {}
for child in node:
if child.get("name"):
attributes[child.get("name")] = next(iter(child), child).text
elif _name(child) in {"creationDate", "function", "roofType", "measuredHeight"}:
attributes[_name(child)] = child.text
return attributes
[docs]
def decode_building(node, *, parent_id=None, inherited=None):
"""Decode owned semantic surfaces once, not duplicate lod2Solid XLinks."""
identifier = node.get(f"{GML}id")
if not identifier:
raise ValueError("CityGML building has no gml:id")
parent_id = parent_id or identifier
attributes = {**(inherited or {}), **_attributes(node)}
vertices, lookup, faces, kinds = [], {}, [], []
zero_area, source_count = 0, 0
for boundary in node:
if _name(boundary) != "boundedBy":
continue
for surface in boundary:
kind = _name(surface)
for polygon in surface.iter(f"{GML}Polygon"):
rings = []
for role in ("exterior", "interior"):
for ring in polygon.findall(f"{GML}{role}/{GML}LinearRing"):
coordinates = ring.find(f"{GML}posList")
if coordinates is None or coordinates.get("srsDimension", "3") != "3":
raise ValueError("CityGML ring requires a 3D posList")
values = np.asarray((coordinates.text or "").split(), dtype=float)
if len(values) % 3 or not np.isfinite(values).all():
raise ValueError("CityGML ring has invalid coordinates")
points = values.reshape(-1, 3)
if len(points) < 4 or not np.array_equal(points[0], points[-1]):
raise ValueError("CityGML LinearRing must be closed")
indices = []
for point in points[:-1]:
key = tuple(point)
if key not in lookup:
lookup[key] = len(vertices)
vertices.append(key)
indices.append(lookup[key])
rings.append(indices)
triangles = triangulate_surface(np.asarray(vertices), rings)
source_count += 1
zero_area += not triangles
faces.extend(triangles)
kinds.extend([kind] * len(triangles))
result = []
if faces:
result.append(
mesh_feature(
identifier=identifier,
source="berlin",
parent_id=parent_id,
vertices=vertices,
triangles=faces,
surface_types=kinds,
crs="EPSG:25833",
vertical_datum="DHHN2016",
lod="2",
attributes={
**attributes,
"attribution": ATTRIBUTION,
"license": LICENSE_URL,
"source_crs": CRS_NAME,
"source_surface_count": source_count,
"omitted_zero_area_surface_count": zero_area,
},
)
)
for container in node:
if _name(container) == "consistsOfBuildingPart":
for part in container:
result.extend(decode_building(part, parent_id=parent_id, inherited=attributes))
if not result:
raise ValueError(f"{identifier} has no owned LoD2 semantic surfaces")
return result
[docs]
def decode_tile(content: bytes, bounds: tuple[float, float, float, float]):
"""Stream one CityGML tile; reject/report whole invalid parent records."""
buildings, warnings = [], []
with zipfile.ZipFile(io.BytesIO(content)) as archive:
members = [i for i in archive.infolist() if i.filename.lower().endswith((".gml", ".xml"))]
if not members or sum(i.file_size for i in members) > MAX_XML_BYTES:
raise ValueError("Berlin archive has no CityGML or exceeds the expanded size limit")
for member in members:
with archive.open(member) as stream:
crs_verified = False
for _, element in ET.iterparse(stream, events=("end",)):
if element.tag == f"{GML}Envelope":
if element.get("srsName") != CRS_NAME:
raise ValueError(
"Berlin CityGML has an unsupported coordinate reference"
)
crs_verified = True
if _name(element) != "cityObjectMember":
continue
if not crs_verified:
raise ValueError("Berlin CityGML is missing its coordinate reference")
for node in element:
if _name(node) != "Building":
continue
# Bound filtering before triangulation avoids meshing a whole 1 km tile.
positions = [
np.fromstring(p.text or "", sep=" ") for p in node.iter(f"{GML}posList")
]
try:
points = np.concatenate(positions).reshape(-1, 3)
if (
points[:, 0].max() < bounds[0]
or points[:, 0].min() > bounds[2]
or points[:, 1].max() < bounds[1]
or points[:, 1].min() > bounds[3]
):
continue
buildings.extend(decode_building(node))
except (ValueError, TypeError, KeyError) as exc:
warnings.append(
f"Skipped Berlin building {node.get(f'{GML}id')}: {exc}"
)
element.clear()
return buildings, warnings
[docs]
class BerlinProvider:
"""Acquire intersecting tiles only, with bounded I/O and per-instance reuse."""
def __init__(self, client: httpx.AsyncClient):
self.client = client
self._lock = asyncio.Lock()
self._index = None
async def _download(self, url, maximum):
chunks, size = [], 0
async with self.client.stream("GET", url, timeout=120) as response:
response.raise_for_status()
async for chunk in response.aiter_bytes():
size += len(chunk)
if size > maximum:
raise ValueError("Berlin download exceeds size limit")
chunks.append(chunk)
return b"".join(chunks)
[docs]
async def buildings(self, latitude, longitude, radius_m, limit):
if limit < 1 or not math.isfinite(radius_m) or radius_m <= 0:
raise ValueError("Building limit and radius must be positive")
if not (-90 <= latitude <= 90 and -180 <= longitude <= 180):
raise ValueError("Invalid query latitude/longitude")
corners = LocalFrame.at(longitude, latitude).corners(
-radius_m, -radius_m, radius_m, radius_m
)
x, y = WGS84_TO_UTM.transform(*zip(*corners, strict=True))
bounds = (min(x), min(y), max(x), max(y))
if not all(math.isfinite(value) for value in bounds):
raise ValueError("Query cannot be projected into Berlin's coordinate system")
low_x, low_y, high_x, high_y = [math.floor(value / 1000) for value in bounds]
if (high_x - low_x + 1) * (high_y - low_y + 1) > 64:
raise ValueError("Berlin query exceeds 64 tiles; use smaller scenes")
buildings, warnings, seen = [], [], set()
async with self._lock:
try:
if self._index is None:
self._index = tile_index(await self._download(INDEX_URL, 4 * 1024 * 1024))
index, release = self._index
keys = [
(east, north)
for east in range(low_x, high_x + 1)
for north in range(low_y, high_y + 1)
]
urls = [index[key] for key in keys if key in index]
if not urls:
raise ValueError("Query does not intersect Berlin's published LoD2 tiles")
if len(urls) < len(keys):
warnings.append("Query extends beyond Berlin's published tile coverage")
if len(urls) > 64:
raise ValueError(
"Berlin query exceeds 64 tiles; split the scene into smaller areas"
)
for url in urls:
content = await self._download(url, MAX_TILE_BYTES)
records, rejected = await asyncio.to_thread(decode_tile, content, bounds)
warnings.extend(rejected)
digest = hashlib.sha256(content).hexdigest()
for record in records:
if record.feature_id in seen:
continue
seen.add(record.feature_id)
record.source_release = release
record.source_attributes.update(source_url=url, source_sha256=digest)
buildings.append(record)
if len(buildings) >= limit:
warnings.append(f"Berlin building limit ({limit}) reached")
break
if len(buildings) >= limit:
break
except (
httpx.HTTPError,
ValueError,
KeyError,
zipfile.BadZipFile,
ParseError,
DefusedXmlException,
) as exc:
raise BuildingProviderError(f"Berlin LoD2 acquisition failed: {exc}") from exc
return BuildingQueryResult(
buildings=buildings, provider="Berlin LoD2", release=release, warnings=warnings
)
