"""Shared conversion of source-resolved, metric building shells to feature records."""
from __future__ import annotations
import numpy as np
import shapely
from pyproj import Transformer
from shapely.geometry import Polygon
from openworld_radio_twin.models import BuildingFeature, BuildingGeometry
[docs]
def mesh_feature(
*,
identifier: str,
source: str,
parent_id: str,
vertices,
triangles,
surface_types,
crs: str,
vertical_datum: str,
lod: str,
attributes: dict,
) -> BuildingFeature:
"""Preserve the shell; derive only its footprint and scalar summary.
Z is in metres in the named vertical datum, regardless of horizontal CRS
units. GroundSurface vertices supply the source footing, not a DEM guess.
Disconnected projections are rejected, never replaced by a convex hull.
"""
points = np.asarray(vertices, dtype=np.float64)
if points.ndim != 2 or points.shape[1] != 3 or not np.isfinite(points).all():
raise ValueError("Native mesh has invalid vertices")
if not triangles or len(triangles) != len(surface_types):
raise ValueError("Native mesh needs triangles and matching semantics")
used = sorted({index for face in triangles for index in face})
if min(used) < 0 or max(used) >= len(points):
raise ValueError("Native mesh triangle index is out of range")
remap = {old: new for new, old in enumerate(used)}
ground_indices = {
index
for face, kind in zip(triangles, surface_types, strict=True)
if kind == "GroundSurface"
for index in face
}
if not ground_indices:
raise ValueError("Native mesh has no GroundSurface for terrain alignment")
ground = float(points[sorted(ground_indices), 2].min())
projections = shapely.polygons(points[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("Native mesh has a disconnected or invalid footprint")
transform = Transformer.from_crs(crs, "EPSG:4326", always_xy=True)
def ring_coordinates(ring):
x, y = np.asarray(ring.coords).T
lon, lat = transform.transform(x, y)
return [[float(a), float(b)] for a, b in zip(lon, lat, strict=True)]
points = points[used]
minimum, maximum = float(points[:, 2].min()), float(points[:, 2].max())
return BuildingFeature(
feature_id=f"{source}:{identifier}",
source=source,
feature_type="building_part" if identifier != parent_id else "building",
source_dataset_id=identifier,
parent_building_id=parent_id,
source_record_id=parent_id,
height_m=maximum - minimum,
min_height_m=max(0.0, minimum - ground),
height_source=f"{source}_native_geometry",
coordinates=ring_coordinates(footprint.exterior),
holes=[ring_coordinates(ring) for ring in footprint.interiors],
geometry=BuildingGeometry(
crs=crs,
vertical_datum=vertical_datum,
lod=lod,
ground_elevation_m=ground,
vertices=points.tolist(),
triangles=[tuple(remap[i] for i in face) for face in triangles],
surface_types=surface_types,
),
source_attributes=attributes,
)
