import json
import math
import os
import xml.etree.ElementTree as ET
from dataclasses import dataclass, fields
from pathlib import Path
from typing import Literal
import numpy as np
import shapely
from PIL import Image, ImageDraw
from shapely import constrained_delaunay_triangles
from shapely.geometry import Polygon
from openworld_radio_twin import __version__
from openworld_radio_twin.environment import (
RF_MATERIAL_PROFILE_VERSION,
RF_MATERIAL_PROFILES,
EnvironmentData,
SurfaceClass,
TerrainGrid,
TerrainSource,
)
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.jsonio import dump_json_with_compact_rows
from openworld_radio_twin.models import BuildingFeature
from openworld_radio_twin.rendering import rgba_lookup
from openworld_radio_twin.scene_format import (
SCENE_ENVIRONMENT_ARRAY_PATHS,
SCENE_FORMAT,
SCENE_HEIGHT_MAP_PATH,
SCENE_LAYOUT_PATHS,
SCENE_LAYOUT_VERSION,
SCENE_MEASUREMENT_SURFACES_DIRECTORY,
SCENE_MESH_DIRECTORY,
SCENE_PROVENANCE_PATH,
SCENE_SOURCE_BUILDINGS_PATH,
SCENE_TERRAIN_DIRECTORY,
SCENE_TERRAIN_INFO_PATH,
SCENE_VOXEL_DEPTH_DIRECTORY,
SCENE_VOXEL_SLICES_DIRECTORY,
SCENE_XML_PATH,
)
from openworld_radio_twin.simulation.native_geometry import (
column_intervals,
ground_projection,
has_ground_faces,
horizontal_enu_vertices,
)
Point2D = tuple[float, float]
Point3D = tuple[float, float, float]
Face = tuple[int, int, int]
REFERENCE_VOXEL_PITCH_M = 3.0
MAX_HEIGHT_MAP_SIDE = 4096
MATERIAL_PREVIEW_COLORS = {
"itu_medium_dry_ground": "0.35 0.28 0.20",
"itu_wood": "0.22 0.38 0.20",
"itu_wet_ground": "0.12 0.28 0.36",
"itu_metal": "0.42 0.46 0.48",
"itu_concrete": "0.62 0.64 0.63",
}
BUILDING_FOOTING_MARGIN_M = 0.05
MEASUREMENT_SURFACE_SCHEMA_VERSION = 2
MaterialProfileName = Literal["itu", "uniform"]
[docs]
@dataclass
class Mesh:
"""Triangle mesh as Python lists or as ``(n, 3)`` NumPy arrays; both serialize alike."""
vertices: list[Point3D] | np.ndarray
faces: list[Face] | np.ndarray
[docs]
@classmethod
def empty(cls) -> "Mesh":
return cls(vertices=[], faces=[])
[docs]
@dataclass(frozen=True)
class LocalBuilding:
feature_id: str
exterior: tuple[Point2D, ...]
holes: tuple[tuple[Point2D, ...], ...]
minimum_height_m: float
vertical_extent_m: float
terrain_anchor_m: float = 0.0
anchor_surface_class: SurfaceClass = SurfaceClass.GROUND
native_vertices: tuple[Point3D, ...] = ()
native_triangles: tuple[Face, ...] = ()
native_surface_types: tuple[str, ...] = ()
source_lod: str | None = None
source_vertical_datum: str | None = None
source_ground_elevation_m: float | None = None
parent_reference_elevation_m: float | None = None
parent_building_id: str | None = None
placement_rule: str = "footprint_extrusion"
@property
def roof_height_m(self) -> float:
return self.minimum_height_m + self.vertical_extent_m
@property
def scene_minimum_height_m(self) -> float:
return self.terrain_anchor_m + self.minimum_height_m
@property
def scene_roof_height_m(self) -> float:
return self.terrain_anchor_m + self.roof_height_m
[docs]
def localize_buildings(
buildings: list[BuildingFeature],
origin_longitude: float,
origin_latitude: float,
environment: EnvironmentData | None = None,
) -> tuple[LocalBuilding, ...]:
frame = LocalFrame.at(origin_longitude, origin_latitude)
native = {}
groups = {}
for building in buildings:
if building.geometry is None:
continue
points = horizontal_enu_vertices(building.geometry, frame)
native[building.feature_id] = points
key = (building.source, building.parent_building_id or building.feature_id)
group = groups.setdefault(
key,
{
"reference": building.geometry.ground_elevation_m,
"datum": building.geometry.vertical_datum,
"ground": [],
},
)
if group["datum"] != building.geometry.vertical_datum:
raise ValueError("Parts of one building cannot mix vertical datums")
# A provider's part ordering must not change the common placement.
group["reference"] = min(group["reference"], building.geometry.ground_elevation_m)
group["ground"].append(ground_projection(points, building.geometry))
for group in groups.values():
footprint = shapely.union_all(group["ground"])
group["anchor"] = (
_terrain_footing_elevation(footprint, environment.terrain)
if environment is not None
else 0.0
)
localized = []
for building in buildings:
exterior = _ring_to_local(frame, building.coordinates)
holes = tuple(_ring_to_local(frame, ring) for ring in building.holes)
polygon = Polygon(exterior, holes)
terrain_anchor = (
_terrain_footing_elevation(polygon, environment.terrain)
if environment is not None and not polygon.is_empty and building.geometry is None
else 0.0
)
native_fields = {}
minimum, extent = building.min_height_m, building.height_m
if building.geometry is not None:
geometry = building.geometry
group = groups[(building.source, building.parent_building_id or building.feature_id)]
terrain_anchor = group["anchor"]
points = native[building.feature_id].copy()
points[:, 2] -= group["reference"]
minimum = float(points[:, 2].min())
extent = float(np.ptp(points[:, 2]))
native_fields = {
"native_vertices": tuple(map(tuple, points)),
"native_triangles": tuple(geometry.triangles),
"native_surface_types": tuple(geometry.surface_types),
"source_lod": geometry.lod,
"source_vertical_datum": geometry.vertical_datum,
"source_ground_elevation_m": geometry.ground_elevation_m,
"parent_reference_elevation_m": group["reference"],
"parent_building_id": building.parent_building_id,
"placement_rule": (
"parent_rigid_ground_alignment_not_datum_conversion"
if has_ground_faces(geometry)
else "parent_rigid_ground_alignment_not_datum_conversion"
";shell_projection_no_ground_faces"
),
}
centroid = polygon.centroid
anchor_surface = (
environment.surface_class_at(float(centroid.x), float(centroid.y))
if environment is not None and not polygon.is_empty
else SurfaceClass.GROUND
)
localized.append(
LocalBuilding(
feature_id=building.feature_id,
exterior=exterior,
holes=holes,
minimum_height_m=minimum,
vertical_extent_m=extent,
terrain_anchor_m=terrain_anchor,
anchor_surface_class=anchor_surface,
**native_fields,
)
)
return tuple(localized)
def _ring_to_local(frame: LocalFrame, ring: list[list[float]]) -> tuple[Point2D, ...]:
"""Convert one WGS84 ring to local ENU with a single transformer call."""
coordinates = np.asarray(ring, dtype=np.float64).reshape(-1, 2)
east, north, _ = frame.transformer.transform(
coordinates[:, 0], coordinates[:, 1], np.zeros(len(coordinates))
)
return tuple(zip(np.asarray(east).tolist(), np.asarray(north).tolist(), strict=True))
def _terrain_footing_elevation(polygon: Polygon, terrain: TerrainGrid) -> float:
"""Place a building below every compiled terrain triangle it intersects."""
minimum_east, minimum_north, maximum_east, maximum_north = polygon.bounds
column_start = max(0, int(np.searchsorted(terrain.east_m, minimum_east, side="left") - 1))
column_stop = min(
terrain.east_m.size - 1,
int(np.searchsorted(terrain.east_m, maximum_east, side="right")),
)
row_start = max(0, int(np.searchsorted(terrain.north_m, minimum_north, side="left") - 1))
row_stop = min(
terrain.north_m.size - 1,
int(np.searchsorted(terrain.north_m, maximum_north, side="right")),
)
if row_start < row_stop and column_start < column_stop:
rows, columns = np.meshgrid(
np.arange(row_start, row_stop), np.arange(column_start, column_stop), indexing="ij"
)
rows = rows.reshape(-1)
columns = columns.reshape(-1)
east, north = terrain.east_m, terrain.north_m
south_west = np.column_stack((east[columns], north[rows]))
south_east = np.column_stack((east[columns + 1], north[rows]))
north_east = np.column_stack((east[columns + 1], north[rows + 1]))
north_west = np.column_stack((east[columns], north[rows + 1]))
# The same SW--NE split as ``_regular_grid_mesh``: two triangles per cell.
corners = np.empty((rows.size, 2, 3, 2), dtype=np.float64)
corners[:, 0] = np.stack((south_west, south_east, north_east), axis=1)
corners[:, 1] = np.stack((south_west, north_east, north_west), axis=1)
hits = shapely.intersects(polygon, shapely.polygons(corners.reshape(-1, 3, 2)))
if hits.any():
vertex_rows = np.empty((rows.size, 2, 3), dtype=np.intp)
vertex_columns = np.empty((rows.size, 2, 3), dtype=np.intp)
vertex_rows[:, 0] = np.stack((rows, rows, rows + 1), axis=1)
vertex_columns[:, 0] = np.stack((columns, columns + 1, columns + 1), axis=1)
vertex_rows[:, 1] = np.stack((rows, rows + 1, rows + 1), axis=1)
vertex_columns[:, 1] = np.stack((columns, columns + 1, columns), axis=1)
elevations = terrain.local_elevation_m[
vertex_rows.reshape(-1, 3)[hits], vertex_columns.reshape(-1, 3)[hits]
]
return float(elevations.min()) - BUILDING_FOOTING_MARGIN_M
centroid = polygon.centroid
return terrain.elevation_at(float(centroid.x), float(centroid.y)) - BUILDING_FOOTING_MARGIN_M
[docs]
@dataclass(frozen=True)
class BuildingMeshAssets:
index: int
feature_id: str
rooftop_path: Path
wall_path: Path
rooftop_triangles: int
wall_triangles: int
[docs]
@dataclass(frozen=True)
class SceneAssets:
directory: Path
scene_xml_path: Path
provenance_path: Path
ground_path: Path
terrain_paths: dict[SurfaceClass, Path]
environment_metadata_path: Path | None
height_map_path: Path
building_assets: tuple[BuildingMeshAssets, ...]
local_buildings: tuple[LocalBuilding, ...]
origin_longitude: float
origin_latitude: float
radius_m: float
height_map_pixel_size_m: float
triangle_count: int
environment: EnvironmentData | None
material_profile: MaterialProfileName
@property
def building_count(self) -> int:
return len(self.building_assets)
[docs]
@dataclass(frozen=True)
class VoxelAssets:
directory: Path
voxel_path: Path
slice_path: Path
scene_info_path: Path
depth_array_path: Path
depth_image_path: Path
depth_preview_path: Path
depth_info_path: Path
pitch_m: float
slice_shape: tuple[int, int]
[docs]
@dataclass(frozen=True)
class MeasurementSurfaceAssets:
path: Path
cell_shape: tuple[int, int]
triangle_areas_m2: np.ndarray
cell_centers: np.ndarray
[docs]
def geographic_to_local(
longitude: float,
latitude: float,
origin_longitude: float,
origin_latitude: float,
) -> Point2D:
east, north, _ = LocalFrame.at(origin_longitude, origin_latitude).to_local(longitude, latitude)
return east, north
def _clean_polygon(points: list[Point2D]) -> list[Point2D]:
cleaned: list[Point2D] = []
for point in points:
if not cleaned or math.dist(point, cleaned[-1]) > 0.01:
cleaned.append(point)
if len(cleaned) > 1 and math.dist(cleaned[0], cleaned[-1]) <= 0.01:
cleaned.pop()
return cleaned
BUILDING_MESH_SCHEMA_VERSION = 4
def _oriented_triangles(polygon: Polygon) -> np.ndarray:
"""Tile the complete footprint, including concavities, excluding its holes.
Returns an ``(n, 3, 2)`` array of counter-clockwise triangles in the order the
constrained Delaunay triangulation emits them.
"""
coordinates = shapely.get_coordinates(constrained_delaunay_triangles(polygon))
if coordinates.size == 0:
return np.empty((0, 3, 2), dtype=np.float64)
triangles = coordinates.reshape(-1, 4, 2)[:, :3, :]
# Shapely's ring signed area for a closed triangle ring, evaluated in the same order.
ax, ay = triangles[:, 0, 0], triangles[:, 0, 1]
bx, by = triangles[:, 1, 0], triangles[:, 1, 1]
cx, cy = triangles[:, 2, 0], triangles[:, 2, 1]
signed_area = ((bx * (cy - ay)) + (cx * (ay - by)) + (ax * (by - cy))) / 2.0
clockwise = signed_area < 0
if clockwise.any():
triangles = triangles.copy()
triangles[clockwise] = triangles[clockwise][:, [0, 2, 1], :]
return triangles
def _polygon_triangles(polygon: Polygon):
"""Legacy generator over the oriented triangles as Shapely polygons."""
for triangle in _oriented_triangles(polygon):
yield Polygon(triangle)
[docs]
def triangulate_polygon(points: list[Point2D]) -> tuple[list[Point2D], list[Face]]:
"""Triangulate a simple polygon while retaining the legacy public helper."""
cleaned = _clean_polygon(points)
if len(cleaned) < 3:
return cleaned, []
polygon = Polygon(cleaned)
vertices = cleaned.copy()
faces: list[Face] = []
if polygon.is_empty or not polygon.is_valid:
return vertices, faces
for triangle in _polygon_triangles(polygon):
face = []
for x, y in list(triangle.exterior.coords)[:3]:
point = (float(x), float(y))
try:
face.append(vertices.index(point))
except ValueError:
vertices.append(point)
face.append(len(vertices) - 1)
faces.append(tuple(face))
return vertices, faces
def _append_building(
roof: Mesh,
walls: Mesh,
footprint: list[Point2D],
holes: list[list[Point2D]],
min_height_m: float,
height_m: float,
terrain_anchor_m: float = 0.0,
) -> bool:
footprint = _clean_polygon(footprint)
holes = [_clean_polygon(ring) for ring in holes]
if len(footprint) < 3:
return False
try:
polygon = Polygon(footprint, [ring for ring in holes if len(ring) >= 3])
except (ValueError, TypeError):
# A ring the source left with too few distinct points cannot form a
# polygon; treat it like any other invalid footprint and skip it.
return False
if polygon.is_empty or not polygon.is_valid or polygon.area <= 0:
return False
minimum_scene_height_m = terrain_anchor_m + min_height_m
roof_height_m = minimum_scene_height_m + height_m
triangles = _oriented_triangles(polygon)
roof_count = len(triangles)
if roof_count:
points = triangles.reshape(-1, 2).tolist()
roof_offset = len(roof.vertices)
roof.vertices.extend((x, y, roof_height_m) for x, y in points)
roof.faces.extend(
(roof_offset + 3 * index, roof_offset + 3 * index + 1, roof_offset + 3 * index + 2)
for index in range(roof_count)
)
bottom_offset = len(walls.vertices)
walls.vertices.extend((x, y, minimum_scene_height_m) for x, y in points)
walls.faces.extend(
(
bottom_offset + 3 * index,
bottom_offset + 3 * index + 2,
bottom_offset + 3 * index + 1,
)
for index in range(roof_count)
)
for ring in [footprint, *holes]:
for start, end in zip(ring, ring[1:] + ring[:1], strict=True):
offset = len(walls.vertices)
walls.vertices.extend(
[
(start[0], start[1], minimum_scene_height_m),
(end[0], end[1], minimum_scene_height_m),
(end[0], end[1], roof_height_m),
(start[0], start[1], roof_height_m),
]
)
walls.faces.extend([(offset, offset + 2, offset + 1), (offset, offset + 3, offset + 2)])
return roof_count > 0
[docs]
def write_binary_ply(path: Path, mesh: Mesh) -> None:
header = (
"ply\n"
"format binary_little_endian 1.0\n"
"comment generated by OpenWorld Radio Twin\n"
f"element vertex {len(mesh.vertices)}\n"
"property float x\n"
"property float y\n"
"property float z\n"
f"element face {len(mesh.faces)}\n"
"property list uchar int vertex_indices\n"
"end_header\n"
).encode("ascii")
vertices = np.ascontiguousarray(
np.asarray(mesh.vertices, dtype=np.float64).reshape(-1, 3), dtype="<f4"
)
faces = np.empty(len(mesh.faces), dtype=np.dtype([("count", "u1"), ("indices", "<i4", (3,))]))
faces["count"] = 3
faces["indices"] = np.asarray(mesh.faces, dtype=np.int64).reshape(-1, 3)
with path.open("wb") as handle:
handle.write(header)
handle.write(vertices.tobytes())
handle.write(faces.tobytes())
def _regular_grid_mesh(
east_m: np.ndarray,
north_m: np.ndarray,
elevation_m: np.ndarray,
) -> tuple[Mesh, np.ndarray, np.ndarray]:
rows, columns = elevation_m.shape
east_grid, north_grid = np.meshgrid(
np.asarray(east_m, dtype=np.float64), np.asarray(north_m, dtype=np.float64)
)
vertex_array = np.column_stack(
(
east_grid.reshape(-1),
north_grid.reshape(-1),
np.asarray(elevation_m, dtype=np.float64).reshape(-1),
)
)
south_west = (
np.arange(rows - 1, dtype=np.int64)[:, None] * columns
+ np.arange(columns - 1, dtype=np.int64)[None, :]
).reshape(-1)
south_east = south_west + 1
north_west = south_west + columns
north_east = north_west + 1
faces = np.empty((south_west.size, 2, 3), dtype=np.int64)
faces[:, 0] = np.stack((south_west, south_east, north_east), axis=1)
faces[:, 1] = np.stack((south_west, north_east, north_west), axis=1)
faces = faces.reshape(-1, 3)
triangles = vertex_array[faces]
normals = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
# ``norm`` over the last axis rounds exactly like the per-triangle norm did.
areas = (np.linalg.norm(normals, axis=1) / 2.0).reshape(rows - 1, columns - 1, 2)
centers = ((vertex_array[south_west] + vertex_array[north_east]) / 2.0).reshape(
rows - 1, columns - 1, 3
)
return Mesh(vertex_array, faces), areas, centers
def _mesh_for_cells(mesh: Mesh, cell_mask: np.ndarray) -> Mesh:
"""Keep the two faces of every selected cell and renumber the vertices they use."""
faces = np.asarray(mesh.faces, dtype=np.int64).reshape(-1, 2, 3)
selected = faces[np.asarray(cell_mask, dtype=bool).reshape(-1)].reshape(-1, 3)
if selected.size == 0:
return Mesh.empty()
used, remapped = np.unique(selected, return_inverse=True)
vertices = np.asarray(mesh.vertices, dtype=np.float64).reshape(-1, 3)[used]
return Mesh(vertices=vertices, faces=remapped.reshape(-1, 3))
[docs]
def write_measurement_surface(
path: Path,
assets: SceneAssets,
resolution_m: float | tuple[float, float],
receiver_height_m: float,
) -> MeasurementSurfaceAssets:
east_resolution_m, north_resolution_m = _cell_size(resolution_m)
columns = max(1, math.ceil(2 * assets.radius_m / east_resolution_m))
rows = max(1, math.ceil(2 * assets.radius_m / north_resolution_m))
east_axis = np.linspace(-assets.radius_m, assets.radius_m, columns + 1, dtype=np.float64)
north_axis = np.linspace(-assets.radius_m, assets.radius_m, rows + 1, dtype=np.float64)
east_grid, north_grid = np.meshgrid(east_axis, north_axis)
elevations = (
np.zeros(east_grid.shape, dtype=np.float64)
if assets.environment is None
else _compiled_terrain_elevation(assets.environment.terrain, east_grid, north_grid)
)
elevations += receiver_height_m
mesh, areas, centers = _regular_grid_mesh(east_axis, north_axis, elevations)
write_binary_ply(path, mesh)
_write_measurement_surface_metadata(
path,
assets,
east_axis,
north_axis,
areas,
centers,
receiver_height_m,
)
return MeasurementSurfaceAssets(path, (rows, columns), areas, centers)
def _cell_size(value: float | tuple[float, float]) -> tuple[float, float]:
if isinstance(value, tuple):
if len(value) != 2:
raise ValueError("cell size must contain east and north resolutions")
east_m, north_m = value
else:
east_m = north_m = value
if east_m <= 0 or north_m <= 0:
raise ValueError("cell size must be positive")
return float(east_m), float(north_m)
def _write_measurement_surface_metadata(
path: Path,
assets: SceneAssets,
east_axis: np.ndarray,
north_axis: np.ndarray,
triangle_areas_m2: np.ndarray,
cell_centers: np.ndarray,
receiver_height_m: float,
) -> None:
arrays_path = path.with_suffix(".npz")
np.savez_compressed(
arrays_path,
east_axis_m=east_axis,
north_axis_m=north_axis,
triangle_areas_m2=triangle_areas_m2,
cell_centers_enu_m=cell_centers,
)
metadata = {
"schema_version": MEASUREMENT_SURFACE_SCHEMA_VERSION,
"vertex_elevation_sampling": "compiled_sw_ne_triangles",
"mesh": path.name,
"arrays": arrays_path.name,
"cell_shape": list(triangle_areas_m2.shape[:2]),
"receiver_height_agl_m": receiver_height_m,
"row_order": "south_to_north",
"column_order": "west_to_east",
"faces_per_cell": 2,
"face_order": [
["south_west", "south_east", "north_east"],
["south_west", "north_east", "north_west"],
],
"face_index": "2 * (row * columns + column) + triangle_in_cell",
"cell_reduction": "triangle-area-weighted mean in the linear domain",
"bounds": {
"west_m": -assets.radius_m,
"south_m": -assets.radius_m,
"east_m": assets.radius_m,
"north_m": assets.radius_m,
},
"terrain_resolution_m": (
assets.environment.terrain.spacing_m if assets.environment is not None else None
),
}
temporary = path.with_suffix(".json.tmp")
temporary.write_text(json.dumps(metadata, indent=2), encoding="utf-8")
temporary.replace(path.with_suffix(".json"))
[docs]
def terrain_elevation(assets: SceneAssets, east_m: float, north_m: float) -> float:
"""Return compiled triangle elevation in the scene's local ENU frame."""
if assets.environment is None:
return 0.0
return float(_compiled_terrain_elevation(assets.environment.terrain, east_m, north_m))
def _compiled_terrain_elevation(
terrain: TerrainGrid, east_m: float | np.ndarray, north_m: float | np.ndarray
) -> np.ndarray:
"""Sample the SW--NE triangles emitted by ``_regular_grid_mesh``.
Source-grid bilinear interpolation remains appropriate for DEM processing,
but radio-device placement must reference the actual piecewise-planar mesh.
Array inputs also avoid a Python loop for every receiver-surface vertex.
"""
east = np.clip(east_m, terrain.east_m[0], terrain.east_m[-1])
north = np.clip(north_m, terrain.north_m[0], terrain.north_m[-1])
col = np.clip(np.searchsorted(terrain.east_m, east) - 1, 0, terrain.east_m.size - 2)
row = np.clip(np.searchsorted(terrain.north_m, north) - 1, 0, terrain.north_m.size - 2)
x = (east - terrain.east_m[col]) / (terrain.east_m[col + 1] - terrain.east_m[col])
y = (north - terrain.north_m[row]) / (terrain.north_m[row + 1] - terrain.north_m[row])
z = terrain.local_elevation_m
sw, se, ne, nw = z[row, col], z[row, col + 1], z[row + 1, col + 1], z[row + 1, col]
return np.where(
y <= x,
(1 - x) * sw + (x - y) * se + y * ne,
(1 - y) * sw + x * ne + (y - x) * nw,
)
[docs]
def transmitter_scene_position(
assets: SceneAssets,
longitude: float,
latitude: float,
height_agl_m: float,
) -> Point3D:
"""Convert a geographic TX position with AGL height to Sionna scene coordinates."""
east_m, north_m, _ = LocalFrame.at(
assets.origin_longitude,
assets.origin_latitude,
).to_local(longitude, latitude, 0.0)
return (
east_m,
north_m,
terrain_elevation(assets, east_m, north_m) + height_agl_m,
)
[docs]
def measurement_surface_path(
assets: SceneAssets,
resolution_m: float | tuple[float, float],
receiver_height_m: float,
) -> Path:
"""Return the stable shared path for one receiver-grid definition."""
directory = assets.directory / SCENE_MEASUREMENT_SURFACES_DIRECTORY
directory.mkdir(exist_ok=True)
east_resolution_m, north_resolution_m = _cell_size(resolution_m)
resolution = _number_token(east_resolution_m)
if east_resolution_m != north_resolution_m:
resolution = f"{resolution}x{_number_token(north_resolution_m)}"
receiver_height = _number_token(receiver_height_m)
return directory / f"surface_{resolution}m_rx_{receiver_height}m_agl.ply"
def _write_environment_assets(
directory: Path,
environment: EnvironmentData,
) -> Path:
terrain_directory = directory / SCENE_TERRAIN_DIRECTORY
terrain_directory.mkdir(exist_ok=True)
arrays = {
"source_elevation": environment.terrain.source_elevation_m,
"model_elevation": environment.terrain.model_elevation_m,
"local_elevation": environment.terrain.local_elevation_m,
"surface_classes": environment.surface_cells,
}
for name, array in arrays.items():
np.save(directory / SCENE_ENVIRONMENT_ARRAY_PATHS[name], array, allow_pickle=False)
source = environment.terrain.source
metadata = {
"coordinate_frame": "WGS84 topocentric ENU",
"axis_order": ["south_to_north_row", "west_to_east_column"],
"horizontal_unit": "m",
"elevation_unit": "m",
"arrays": {
name: {
"path": SCENE_ENVIRONMENT_ARRAY_PATHS[name].name,
"shape": list(array.shape),
"dtype": str(array.dtype),
"sampling": "cells" if name == "surface_classes" else "nodes",
"unit": "class_id" if name == "surface_classes" else "m",
}
for name, array in arrays.items()
},
"east_m": environment.terrain.east_m.tolist(),
"north_m": environment.terrain.north_m.tolist(),
"source_origin_elevation_m": environment.terrain.source_origin_elevation_m,
"model_origin_elevation_m": environment.terrain.origin_elevation_m,
"local_elevation_definition": "model_elevation_m - model_elevation_at_ENU_origin",
"model_operations": list(environment.terrain.model_operations),
"terrain_source": {
"provider": source.provider,
"dataset": source.dataset,
"vertical_datum": source.vertical_datum,
"requested_spacing_m": source.requested_spacing_m,
"source_resolution_m": source.source_resolution_m,
"interpolation": source.interpolation,
"source_url": source.source_url,
"license": source.license,
"attributes": source.attributes,
},
"surface_source": environment.surface_provenance,
"overture_release": environment.overture_release,
"surface_classes": {str(int(item)): item.name.lower() for item in SurfaceClass},
"rf_material_profile_version": RF_MATERIAL_PROFILE_VERSION,
"rf_material_profiles": {
item.name.lower(): {
"sionna_material": RF_MATERIAL_PROFILES[item].sionna_material,
"basis": RF_MATERIAL_PROFILES[item].basis,
"confidence": RF_MATERIAL_PROFILES[item].confidence,
"note": RF_MATERIAL_PROFILES[item].note,
}
for item in SurfaceClass
},
}
path = directory / SCENE_TERRAIN_INFO_PATH
path.write_text(json.dumps(metadata, indent=2), encoding="utf-8")
return path
[docs]
def load_environment_assets(directory: Path) -> EnvironmentData:
"""Reconstruct the solver environment from a persisted scene asset directory."""
metadata = json.loads((directory / SCENE_TERRAIN_INFO_PATH).read_text(encoding="utf-8"))
arrays = {
name: np.load(directory / path, allow_pickle=False)
for name, path in SCENE_ENVIRONMENT_ARRAY_PATHS.items()
}
source = TerrainSource(**metadata["terrain_source"])
terrain = TerrainGrid(
east_m=np.asarray(metadata["east_m"], dtype=np.float64),
north_m=np.asarray(metadata["north_m"], dtype=np.float64),
source_elevation_m=arrays["source_elevation"],
model_elevation_m=arrays["model_elevation"],
local_elevation_m=arrays["local_elevation"],
source_origin_elevation_m=float(metadata["source_origin_elevation_m"]),
origin_elevation_m=float(metadata["model_origin_elevation_m"]),
source=source,
model_operations=tuple(metadata["model_operations"]),
)
return EnvironmentData(
terrain=terrain,
surface_features=(),
surface_cells=arrays["surface_classes"],
overture_release=metadata.get("overture_release"),
surface_provenance=metadata["surface_source"],
)
def _footprint_window(
building: LocalBuilding,
radius_m: float,
width: int,
height: int,
) -> tuple[int, int, np.ndarray] | None:
"""Rasterize a footprint inside its pixel bounding box.
Returns ``(row, column, mask)`` where ``mask`` covers the pixels from that
top-left corner, or ``None`` when the footprint lies outside the raster.
Pixel coordinates are computed exactly as for the full raster and shifted by
whole pixels, so the scan conversion is unchanged.
"""
def pixel_ring(ring: tuple[Point2D, ...]) -> list[tuple[float, float]]:
return [
(
(east + radius_m) * width / (2.0 * radius_m),
(radius_m - north) * height / (2.0 * radius_m),
)
for east, north in ring
]
rings = [pixel_ring(building.exterior), *(pixel_ring(hole) for hole in building.holes)]
xs = [x for ring in rings for x, _ in ring]
ys = [y for ring in rings for _, y in ring]
if not xs:
return None
column = max(0, math.floor(min(xs)) - 1)
column_stop = min(width, math.ceil(max(xs)) + 2)
row = max(0, math.floor(min(ys)) - 1)
row_stop = min(height, math.ceil(max(ys)) + 2)
if column >= column_stop or row >= row_stop:
return None
image = Image.new("L", (column_stop - column, row_stop - row), 0)
draw = ImageDraw.Draw(image)
draw.polygon([(x - column, y - row) for x, y in rings[0]], fill=255)
for hole in rings[1:]:
draw.polygon([(x - column, y - row) for x, y in hole], fill=0)
return row, column, np.asarray(image, dtype=np.uint8) > 0
[docs]
def rasterize_receiver_slice(
buildings: tuple[LocalBuilding, ...],
radius_m: float,
receiver_height_m: float,
pitch_m: float = REFERENCE_VOXEL_PITCH_M,
environment: EnvironmentData | None = None,
) -> tuple[np.ndarray, np.ndarray]:
"""Rasterize the north-up receiver plane used by exports and the web API."""
if pitch_m <= 0:
raise ValueError("Voxel pitch must be positive")
size = max(1, math.ceil(2.0 * radius_m / pitch_m))
mask = np.zeros((size, size), dtype=bool)
roof_heights = np.full((size, size), np.nan, dtype=np.float32)
for building in buildings:
if building.native_vertices:
east = -radius_m + (np.arange(size) + 0.5) * (2 * radius_m / size)
north = radius_m - (np.arange(size) + 0.5) * (2 * radius_m / size)
rows, columns, lower, upper = column_intervals(
building.native_vertices, building.native_triangles, east, north
)
ground = (
np.zeros(len(rows))
if environment is None
else _compiled_terrain_elevation(environment.terrain, east[columns], north[rows])
)
lower = lower + building.terrain_anchor_m - ground
upper = upper + building.terrain_anchor_m - ground
included = (lower <= receiver_height_m) & (receiver_height_m <= upper)
rr, cc = rows[included], columns[included]
mask[rr, cc] = True
np.fmax.at(roof_heights, (rr, cc), upper[included])
continue
if not building.minimum_height_m <= receiver_height_m <= building.roof_height_m:
continue
window = _footprint_window(building, radius_m, size, size)
if window is None:
continue
row, column, footprint = window
region = (slice(row, row + footprint.shape[0]), slice(column, column + footprint.shape[1]))
mask[region] |= footprint
heights = roof_heights[region]
current = heights[footprint]
heights[footprint] = np.where(
np.isfinite(current),
np.maximum(current, building.roof_height_m),
building.roof_height_m,
)
return mask, roof_heights
def _write_height_map(
path: Path,
buildings: tuple[LocalBuilding, ...],
radius_m: float,
environment: EnvironmentData | None = None,
) -> float:
size = max(1, min(MAX_HEIGHT_MAP_SIDE, math.ceil(2.0 * radius_m)))
height_map = np.zeros((size, size), dtype=np.float32)
for building in buildings:
if building.native_vertices:
east = -radius_m + (np.arange(size) + 0.5) * (2 * radius_m / size)
north = radius_m - (np.arange(size) + 0.5) * (2 * radius_m / size)
rows, columns, _, upper = column_intervals(
building.native_vertices, building.native_triangles, east, north
)
ground = (
np.zeros(len(rows))
if environment is None
else _compiled_terrain_elevation(environment.terrain, east[columns], north[rows])
)
np.maximum.at(height_map, (rows, columns), upper + building.terrain_anchor_m - ground)
continue
window = _footprint_window(building, radius_m, size, size)
if window is None:
continue
row, column, mask = window
region = height_map[row : row + mask.shape[0], column : column + mask.shape[1]]
region[mask] = np.maximum(region[mask], building.roof_height_m)
np.save(path, height_map, allow_pickle=False)
return 2.0 * radius_m / size
def _add_diffuse_material(scene: ET.Element, material_id: str, rgb: str) -> None:
two_sided = ET.SubElement(scene, "bsdf", type="twosided", id=material_id)
diffuse = ET.SubElement(two_sided, "bsdf", type="diffuse")
ET.SubElement(diffuse, "rgb", name="reflectance", value=rgb)
def _add_shape(scene: ET.Element, shape_id: str, filename: str, material_id: str) -> None:
shape = ET.SubElement(scene, "shape", type="ply", id=shape_id)
ET.SubElement(shape, "string", name="filename", value=filename)
ET.SubElement(shape, "ref", id=material_id, name="bsdf")
ET.SubElement(shape, "boolean", name="face_normals", value="true")
def _write_scene_xml(
path: Path,
assets: tuple[BuildingMeshAssets, ...],
terrain_paths: dict[SurfaceClass, Path],
origin_longitude: float,
origin_latitude: float,
radius_m: float,
maximum_height_m: float,
vertical_datum: str,
material_profile: MaterialProfileName,
) -> None:
frame = LocalFrame.at(origin_longitude, origin_latitude)
corners = frame.corners(-radius_m, -radius_m, radius_m, radius_m)
longitudes = [point[0] for point in corners]
latitudes = [point[1] for point in corners]
scene = ET.Element("scene", version="2.1.0")
defaults = {
"spp": "4096",
"resx": "1024",
"resy": "1024",
"scenegen_version": __version__,
"scenegen_min_lat": f"{min(latitudes):.12f}",
"scenegen_max_lat": f"{max(latitudes):.12f}",
"scenegen_min_lon": f"{min(longitudes):.12f}",
"scenegen_max_lon": f"{max(longitudes):.12f}",
"scenegen_center_lat": f"{origin_latitude:.12f}",
"scenegen_center_lon": f"{origin_longitude:.12f}",
"scenegen_bbox_width": f"{2.0 * radius_m:.6f}",
"scenegen_bbox_length": f"{2.0 * radius_m:.6f}",
"owrt_local_frame": "WGS84 topocentric ENU",
"owrt_horizontal_unit": "m",
"owrt_vertical_datum": vertical_datum,
"owrt_maximum_roof_height_m": f"{maximum_height_m:.6f}",
"owrt_material_profile": material_profile,
"scenegen_ground_material": (
"semantic surface material profiles" if material_profile == "itu" else "itu_concrete"
),
"scenegen_rooftop_material": ("itu_metal" if material_profile == "itu" else "itu_concrete"),
"scenegen_wall_material": "itu_concrete",
}
for name, value in defaults.items():
ET.SubElement(scene, "default", name=name, value=value)
integrator = ET.SubElement(scene, "integrator", type="path")
ET.SubElement(integrator, "integer", name="max_depth", value="12")
material_ids = {profile.sionna_material for profile in RF_MATERIAL_PROFILES.values()} | {
"itu_metal",
"itu_concrete",
}
for material_id in sorted(material_ids):
_add_diffuse_material(scene, material_id, MATERIAL_PREVIEW_COLORS[material_id])
for surface_class, terrain_path in terrain_paths.items():
profile = RF_MATERIAL_PROFILES[surface_class]
_add_shape(
scene,
f"mesh-terrain-{surface_class.name.lower()}",
f"mesh/{terrain_path.name}",
profile.sionna_material if material_profile == "itu" else "itu_concrete",
)
for asset in assets:
_add_shape(
scene,
f"mesh-building_{asset.index}_rooftop",
f"mesh/{asset.rooftop_path.name}",
"itu_metal" if material_profile == "itu" else "itu_concrete",
)
_add_shape(
scene,
f"mesh-building_{asset.index}_wall",
f"mesh/{asset.wall_path.name}",
"itu_concrete",
)
ET.indent(scene, space=" ")
ET.ElementTree(scene).write(path, encoding="utf-8", xml_declaration=True)
def _write_scene_provenance(
path: Path,
origin_longitude: float,
origin_latitude: float,
radius_m: float,
vertical_reference: str,
material_profile: MaterialProfileName,
terrain_enabled: bool,
building_assets: tuple[BuildingMeshAssets, ...],
local_buildings: tuple[LocalBuilding, ...],
) -> None:
"""Persist the geographic-to-local contract shared by every scene consumer."""
frame = LocalFrame.at(origin_longitude, origin_latitude)
payload = {
"building_assets": [
{
"index": asset.index,
"feature_id": asset.feature_id,
"rooftop_triangles": asset.rooftop_triangles,
"wall_triangles": asset.wall_triangles,
}
for asset in building_assets
],
"local_buildings": [_local_building_record(building) for building in local_buildings],
"schema_version": 1,
"scene_layout": {
"format": SCENE_FORMAT,
"version": SCENE_LAYOUT_VERSION,
"paths": dict(SCENE_LAYOUT_PATHS),
},
"coordinate_contract": {
"geographic_crs": "EPSG:4326",
"local_frame": "topocentric_enu",
"origin": {
"longitude": origin_longitude,
"latitude": origin_latitude,
"altitude_m": 0.0,
},
"axis_order": ["east", "north", "up"],
"horizontal_unit": "m",
"vertical_unit": "m",
"geodetic_origin_height_m": 0.0,
"geodetic_origin_vertical_datum": "WGS84 ellipsoidal height",
"terrain_source_vertical_datum": vertical_reference,
"scene_z_reference": (
"modeled terrain elevation at ENU origin" if terrain_enabled else "local z=0 plane"
),
"local_frame_up_component": (
"ellipsoidal up relative to the geodetic origin; not an AGL height"
),
},
"geometry_contract": {
"building_footing": (
"minimum vertex elevation of every intersecting compiled terrain triangle"
if terrain_enabled
else "local z=0 plane"
),
"building_footing_margin_m": (BUILDING_FOOTING_MARGIN_M if terrain_enabled else 0.0),
"boundary_policy": "complete footprints within square ENU bounds",
},
"material_contract": {
"profile": material_profile,
"profile_version": RF_MATERIAL_PROFILE_VERSION,
},
"bounds": {
"shape": "square",
"radius_m": radius_m,
"local_envelope_m": {
"west": -radius_m,
"south": -radius_m,
"east": radius_m,
"north": radius_m,
},
"wgs84_corners": frame.corners(-radius_m, -radius_m, radius_m, radius_m),
},
}
dump_json_with_compact_rows(path, payload, "local_buildings")
def _local_building_record(building: LocalBuilding) -> dict[str, object]:
"""Field mapping of one placed building without the deep copies of ``asdict``."""
return {field.name: getattr(building, field.name) for field in fields(LocalBuilding)}
[docs]
def build_scene_assets(
directory: Path,
buildings: list[BuildingFeature],
origin_longitude: float,
origin_latitude: float,
radius_m: float,
environment: EnvironmentData | None = None,
include_ground: bool = True,
material_profile: MaterialProfileName = "itu",
) -> SceneAssets:
if material_profile not in {"itu", "uniform"}:
raise ValueError(f"Unknown material profile: {material_profile}")
directory.mkdir(parents=True, exist_ok=True)
mesh_directory = directory / SCENE_MESH_DIRECTORY
mesh_directory.mkdir(exist_ok=True)
if environment is None:
ground = (
Mesh(
vertices=[
(-radius_m, -radius_m, 0.0),
(radius_m, -radius_m, 0.0),
(radius_m, radius_m, 0.0),
(-radius_m, radius_m, 0.0),
],
faces=[(0, 1, 2), (0, 2, 3)],
)
if include_ground
else Mesh.empty()
)
ground_path = mesh_directory / "ground.ply"
write_binary_ply(ground_path, ground)
terrain_paths = {SurfaceClass.GROUND: ground_path} if include_ground else {}
environment_metadata_path = None
terrain_triangle_count = len(ground.faces)
else:
terrain_mesh, _, _ = _regular_grid_mesh(
environment.terrain.east_m,
environment.terrain.north_m,
environment.terrain.local_elevation_m,
)
terrain_paths = {}
for surface_class in SurfaceClass:
class_mesh = _mesh_for_cells(
terrain_mesh,
environment.surface_cells == int(surface_class),
)
if len(class_mesh.faces) == 0:
continue
filename = (
"ground.ply"
if surface_class == SurfaceClass.GROUND
else f"terrain_{surface_class.name.lower()}.ply"
)
path = mesh_directory / filename
write_binary_ply(path, class_mesh)
terrain_paths[surface_class] = path
if SurfaceClass.GROUND not in terrain_paths:
ground_path = mesh_directory / "ground.ply"
write_binary_ply(ground_path, Mesh.empty())
else:
ground_path = terrain_paths[SurfaceClass.GROUND]
environment_metadata_path = _write_environment_assets(directory, environment)
terrain_triangle_count = len(terrain_mesh.faces)
if environment is None:
ground_path = mesh_directory / "ground.ply"
building_assets: list[BuildingMeshAssets] = []
local_buildings: list[LocalBuilding] = []
triangle_count = terrain_triangle_count
for local in localize_buildings(
buildings,
origin_longitude,
origin_latitude,
environment,
):
roof = Mesh.empty()
walls = Mesh.empty()
if local.native_vertices:
vertices = np.asarray(local.native_vertices, dtype=np.float64)
vertices = np.column_stack(
(vertices[:, 0], vertices[:, 1], vertices[:, 2] + local.terrain_anchor_m)
)
roof.vertices = vertices
walls.vertices = vertices
native_faces = np.asarray(local.native_triangles, dtype=np.int64).reshape(-1, 3)
is_roof = np.asarray(local.native_surface_types, dtype=object) == "RoofSurface"
roof.faces = native_faces[is_roof]
walls.faces = native_faces[~is_roof]
elif not _append_building(
roof,
walls,
list(local.exterior),
[list(ring) for ring in local.holes],
local.minimum_height_m,
local.vertical_extent_m,
local.terrain_anchor_m,
):
continue
index = len(building_assets)
rooftop_path = mesh_directory / f"building_{index}_rooftop.ply"
wall_path = mesh_directory / f"building_{index}_wall.ply"
write_binary_ply(rooftop_path, roof)
write_binary_ply(wall_path, walls)
building_assets.append(
BuildingMeshAssets(
index=index,
feature_id=local.feature_id,
rooftop_path=rooftop_path,
wall_path=wall_path,
rooftop_triangles=len(roof.faces),
wall_triangles=len(walls.faces),
)
)
local_buildings.append(local)
triangle_count += len(roof.faces) + len(walls.faces)
local_tuple = tuple(local_buildings)
asset_tuple = tuple(building_assets)
height_map_path = directory / SCENE_HEIGHT_MAP_PATH
height_map_pixel_size_m = _write_height_map(height_map_path, local_tuple, radius_m, environment)
provenance_path = directory / SCENE_PROVENANCE_PATH
vertical_reference = (
environment.terrain.source.vertical_datum if environment else "local ground plane"
)
_write_scene_provenance(
provenance_path,
origin_longitude,
origin_latitude,
radius_m,
vertical_reference,
material_profile,
environment is not None,
asset_tuple,
local_tuple,
)
scene_xml_path = directory / SCENE_XML_PATH
maximum_height = max((item.roof_height_m for item in local_tuple), default=0.0)
_write_scene_xml(
scene_xml_path,
asset_tuple,
terrain_paths,
origin_longitude,
origin_latitude,
radius_m,
maximum_height,
vertical_reference,
material_profile,
)
return SceneAssets(
directory=directory,
scene_xml_path=scene_xml_path,
provenance_path=provenance_path,
ground_path=ground_path,
terrain_paths=terrain_paths,
environment_metadata_path=environment_metadata_path,
height_map_path=height_map_path,
building_assets=asset_tuple,
local_buildings=local_tuple,
origin_longitude=origin_longitude,
origin_latitude=origin_latitude,
radius_m=radius_m,
height_map_pixel_size_m=height_map_pixel_size_m,
triangle_count=triangle_count,
environment=environment,
material_profile=material_profile,
)
def _ply_face_count(path: Path) -> int:
"""Read a native PLY header without loading its vertex/face payload."""
count = None
with path.open("rb") as stream:
if stream.readline() != b"ply\n":
raise ValueError(f"Invalid PLY header: {path}")
for _ in range(1024):
line = stream.readline(4096).strip()
if line.startswith(b"element face "):
count = int(line.split()[-1])
if line == b"end_header":
if count is not None and count >= 0:
return count
break
if not line:
break
raise ValueError(f"Missing PLY face count or end_header: {path}")
[docs]
def load_scene_assets(directory: Path) -> SceneAssets:
"""Restore the scene metadata needed by geographic and derived-asset helpers."""
provenance = json.loads((directory / SCENE_PROVENANCE_PATH).read_text(encoding="utf-8"))
coordinate_contract = provenance["coordinate_contract"]
origin = coordinate_contract["origin"]
radius_m = float(provenance["bounds"]["radius_m"])
environment = (
load_environment_assets(directory)
if (directory / SCENE_TERRAIN_INFO_PATH).is_file()
else None
)
if "local_buildings" in provenance:
local_buildings = tuple(
LocalBuilding(
**{
**item,
"exterior": tuple(tuple(point) for point in item["exterior"]),
"holes": tuple(tuple(tuple(point) for point in ring) for ring in item["holes"]),
"anchor_surface_class": SurfaceClass(item["anchor_surface_class"]),
"native_vertices": tuple(
tuple(point) for point in item.get("native_vertices", [])
),
"native_triangles": tuple(
tuple(face) for face in item.get("native_triangles", [])
),
"native_surface_types": tuple(item.get("native_surface_types", [])),
}
)
for item in provenance["local_buildings"]
)
else:
# Legacy caches without placement records: re-derive them from the source records.
source_path = directory / SCENE_SOURCE_BUILDINGS_PATH
buildings: list[BuildingFeature] = []
if source_path.is_file():
source = json.loads(source_path.read_text(encoding="utf-8"))
buildings = [
BuildingFeature.model_validate(item) for item in source.get("buildings", [])
]
if any(building.geometry is not None for building in buildings):
raise ValueError("Native geometry cache lacks placement records; recompile the scene")
# Legacy indices were assigned only after successful mesh construction.
local_buildings = localize_buildings(
buildings, float(origin["longitude"]), float(origin["latitude"]), environment
)
local_buildings = tuple(
local
for local in local_buildings
if _append_building(
Mesh.empty(),
Mesh.empty(),
list(local.exterior),
[list(ring) for ring in local.holes],
local.minimum_height_m,
local.vertical_extent_m,
local.terrain_anchor_m,
)
)
mesh_directory = directory / SCENE_MESH_DIRECTORY
with os.scandir(mesh_directory) as entries:
mesh_files = {entry.name for entry in entries if entry.is_file()}
rooftop_indices = sorted(
int(name.split("_")[1])
for name in mesh_files
if name.startswith("building_") and name.endswith("_rooftop.ply")
)
if len(rooftop_indices) != len(local_buildings):
raise ValueError("Cannot recover building identities from this cache; recompile the scene")
recorded = {item["index"]: item for item in provenance.get("building_assets", [])}
building_assets = []
for index in rooftop_indices:
rooftop_path = mesh_directory / f"building_{index}_rooftop.ply"
wall_path = mesh_directory / f"building_{index}_wall.ply"
if index != len(building_assets):
raise ValueError("Non-contiguous building mesh indices; recompile the scene")
if wall_path.name not in mesh_files:
raise ValueError(f"Missing wall mesh for building {index}; recompile the scene")
local = local_buildings[index]
record = recorded.get(index)
if record is None:
# Legacy caches predate the placement records: count the faces in the meshes.
counts = _ply_face_count(rooftop_path), _ply_face_count(wall_path)
else:
counts = record["rooftop_triangles"], record["wall_triangles"]
building_assets.append(
BuildingMeshAssets(index, local.feature_id, rooftop_path, wall_path, *counts)
)
face_counts = {}
for asset in building_assets:
face_counts[f"mesh/{asset.rooftop_path.name}"] = asset.rooftop_triangles
face_counts[f"mesh/{asset.wall_path.name}"] = asset.wall_triangles
def face_count(filename: str) -> int:
if filename in face_counts:
return face_counts[filename]
return _ply_face_count(directory / filename)
terrain_paths = {
surface_class: path
for surface_class, filename in {
SurfaceClass.GROUND: "ground.ply",
SurfaceClass.VEGETATION: "terrain_vegetation.ply",
SurfaceClass.PAVED: "terrain_paved.ply",
SurfaceClass.WATER: "terrain_water.ply",
}.items()
if (path := mesh_directory / filename).is_file()
}
height_map_path = directory / SCENE_HEIGHT_MAP_PATH
height_map_side = np.load(height_map_path, mmap_mode="r").shape[0]
return SceneAssets(
directory=directory,
scene_xml_path=directory / SCENE_XML_PATH,
provenance_path=directory / SCENE_PROVENANCE_PATH,
ground_path=mesh_directory / "ground.ply",
terrain_paths=terrain_paths,
environment_metadata_path=(
directory / SCENE_TERRAIN_INFO_PATH if environment is not None else None
),
height_map_path=height_map_path,
building_assets=tuple(building_assets),
local_buildings=local_buildings,
origin_longitude=float(origin["longitude"]),
origin_latitude=float(origin["latitude"]),
radius_m=radius_m,
height_map_pixel_size_m=2.0 * radius_m / height_map_side,
triangle_count=sum(
face_count(node.attrib["value"])
for node in ET.parse(directory / SCENE_XML_PATH)
.getroot()
.findall("shape/string[@name='filename']")
),
environment=environment,
material_profile=provenance.get("material_contract", {}).get("profile", "itu"),
)
def _number_token(value: float) -> str:
return str(int(value)) if float(value).is_integer() else f"{value:g}"
def _write_voxel_depth_assets(
root: Path,
indices: np.ndarray,
size: int,
pitch_m: float,
bbox: np.ndarray,
) -> tuple[Path, Path, Path, Path]:
"""Project a sparse voxel volume to its highest occupied z per XY column."""
depth = np.full((size, size), -np.inf, dtype=np.float32)
if indices.size:
np.maximum.at(depth, (indices[:, 0], indices[:, 1]), indices[:, 2] * pitch_m)
depth[~np.isfinite(depth)] = np.nan
directory = root / SCENE_VOXEL_DEPTH_DIRECTORY
directory.mkdir(exist_ok=True)
pitch_token = f"{pitch_m:.1f}"
stem = f"scene_depth_max_z_pitch{pitch_token}"
array_path = directory / f"{stem}.npy"
image_path = directory / f"{stem}.png"
preview_path = directory / f"{stem}_preview.png"
info_path = directory / "scene_info.json"
np.save(array_path, depth, allow_pickle=False)
valid = np.isfinite(depth)
maximum = float(np.max(depth[valid], initial=0.0))
metres_per_code = maximum / 65534.0 if maximum > 0 else 0.0
encoded = np.zeros(depth.shape, dtype=np.uint16)
if maximum > 0:
encoded[valid] = 1 + np.rint(depth[valid] / metres_per_code).astype(np.uint16)
elif valid.any():
encoded[valid] = 1
Image.fromarray(encoded).save(image_path)
preview = np.zeros((*depth.shape, 4), dtype=np.uint8)
if maximum > 0:
normalized = np.zeros(depth.shape, dtype=np.float32)
normalized[valid] = depth[valid] / maximum
lookup = rgba_lookup("viridis")
preview[valid] = lookup[np.rint(normalized[valid] * 255).astype(np.uint8)]
Image.fromarray(preview, mode="RGBA").save(preview_path)
metadata = {
"voxel_source": f"../vox_slices/scene_voxels_pitch{pitch_token}.npz",
"array": array_path.name,
"png": image_path.name,
"preview_png": preview_path.name,
"value": "maximum_occupied_voxel_z_agl_m",
"unit": "m",
"vertical_datum": "local_ground_plane",
"shape": [size, size],
"row_order": "north_to_south",
"column_order": "west_to_east",
"pitch_m": pitch_m,
"bbox": bbox.tolist(),
"array_encoding": {
"dtype": "float32",
"no_data": "NaN",
},
"png_encoding": {
"dtype": "uint16",
"no_data_code": 0,
"valid_code_minimum": 1,
"valid_code_maximum": 65535,
"normalization": "linear",
"maximum_depth_m": maximum,
"metres_per_code": metres_per_code,
"decode": "depth_m = (code - 1) * metres_per_code for code > 0",
},
"preview_encoding": {
"role": "visualization_only",
"colormap": "viridis",
"normalization": "linear_full_depth_range",
"no_data_rendering": "transparent",
},
"spatial_resampling": "none",
}
info_path.write_text(json.dumps(metadata, indent=2), encoding="utf-8")
return array_path, image_path, preview_path, info_path
[docs]
def write_voxel_assets(
assets: SceneAssets,
receiver_height_m: float,
pitch_m: float = REFERENCE_VOXEL_PITCH_M,
) -> VoxelAssets:
"""Write the sparse solid voxels and north-up receiver-height slice."""
directory = assets.directory / SCENE_VOXEL_SLICES_DIRECTORY
directory.mkdir(exist_ok=True)
size = max(1, math.ceil(2.0 * assets.radius_m / pitch_m))
slice_mask, _ = rasterize_receiver_slice(
assets.local_buildings,
assets.radius_m,
receiver_height_m,
pitch_m,
assets.environment,
)
occupied_indices: list[np.ndarray] = []
for building in assets.local_buildings:
if building.native_vertices:
east = -assets.radius_m + np.arange(size) * pitch_m
north = assets.radius_m - np.arange(size) * pitch_m
rows, columns, lower, upper = column_intervals(
building.native_vertices, building.native_triangles, east, north
)
ground = (
np.zeros(len(rows))
if assets.environment is None
else _compiled_terrain_elevation(
assets.environment.terrain, east[columns], north[rows]
)
)
lower = lower + building.terrain_anchor_m - ground
upper = upper + building.terrain_anchor_m - ground
z_start = np.ceil(lower / pitch_m).astype(np.int64)
z_stop = np.floor(upper / pitch_m).astype(np.int64) + 1
counts = np.maximum(z_stop - z_start, 0)
total = int(counts.sum())
if total:
offsets = np.repeat(np.cumsum(counts) - counts, counts)
z_indices = np.repeat(z_start, counts) + (np.arange(total) - offsets)
occupied_indices.append(
np.column_stack(
(
np.repeat(rows, counts).astype(np.int32),
np.repeat(columns, counts).astype(np.int32),
z_indices.astype(np.int32),
)
)
)
continue
window = _footprint_window(building, assets.radius_m, size, size)
if window is None:
continue
window_row, window_column, mask = window
rows, columns = np.nonzero(mask)
rows = rows + window_row
columns = columns + window_column
z_start = math.ceil(building.minimum_height_m / pitch_m)
z_stop = math.floor(building.roof_height_m / pitch_m)
if rows.size == 0 or z_stop < z_start:
continue
z_indices = np.arange(z_start, z_stop + 1, dtype=np.int32)
repeated = np.column_stack(
[
np.repeat(rows, z_indices.size),
np.repeat(columns, z_indices.size),
np.tile(z_indices, rows.size),
]
)
occupied_indices.append(repeated)
if occupied_indices:
indices = np.unique(np.concatenate(occupied_indices, axis=0), axis=0)
points = np.column_stack(
[
-assets.radius_m + indices[:, 1] * pitch_m,
assets.radius_m - indices[:, 0] * pitch_m,
indices[:, 2] * pitch_m,
]
).astype(np.float32)
else:
indices = np.empty((0, 3), dtype=np.int32)
points = np.empty((0, 3), dtype=np.float32)
pitch_token = f"{pitch_m:.1f}"
voxel_path = directory / f"scene_voxels_pitch{pitch_token}.npz"
maximum_height = max((item.roof_height_m for item in assets.local_buildings), default=0.0)
minimum_height = 0.0
if points.size:
minimum_height = min(0.0, float(points[:, 2].min()))
maximum_height = max(maximum_height, float(points[:, 2].max()))
bbox = np.asarray(
[
-assets.radius_m,
-assets.radius_m,
assets.radius_m,
assets.radius_m,
minimum_height,
maximum_height,
],
dtype=np.float32,
)
np.savez_compressed(
voxel_path,
points_xyz=points,
pitch=np.asarray([pitch_m], dtype=np.float32),
bbox=bbox,
ground_z=np.asarray([0.0], dtype=np.float32),
solid=np.asarray([1], dtype=np.uint8),
)
slice_path = directory / (
f"scene_slice_{_number_token(receiver_height_m)}m_pitch{pitch_token}.png"
)
Image.fromarray(slice_mask.astype(np.uint8) * 255, mode="L").save(slice_path)
scene_info_path = directory / "scene_info.json"
scene_info = {
"scene_xml": "../scene.xml",
"mesh_files": [
str(path.relative_to(assets.directory))
for item in assets.building_assets
for path in (item.rooftop_path, item.wall_path)
],
"exclude_names": ["ground.ply", "terrain.ply", "surface.ply"],
"resolution": pitch_m,
"solid": True,
"sionna_bbox": [
-assets.radius_m,
-assets.radius_m,
assets.radius_m,
assets.radius_m,
minimum_height,
maximum_height,
],
"target_heights": [receiver_height_m],
"coordinate_system": "WGS84 topocentric ENU",
"horizontal_unit": "m",
"vertical_unit": "m AGL",
"height_map": {
"path": "../2D_Building_Height_Map.npy",
"dtype": "float32",
"value": "roof_height_agl_m",
"row_order": "north_to_south",
"pixel_size_m": assets.height_map_pixel_size_m,
},
}
scene_info_path.write_text(json.dumps(scene_info, indent=2), encoding="utf-8")
depth_array_path, depth_image_path, depth_preview_path, depth_info_path = (
_write_voxel_depth_assets(
assets.directory,
indices,
size,
pitch_m,
bbox,
)
)
return VoxelAssets(
directory=directory,
voxel_path=voxel_path,
slice_path=slice_path,
scene_info_path=scene_info_path,
depth_array_path=depth_array_path,
depth_image_path=depth_image_path,
depth_preview_path=depth_preview_path,
depth_info_path=depth_info_path,
pitch_m=pitch_m,
slice_shape=slice_mask.shape,
)
