import numpy as np
from shapely import contains_xy
from shapely.geometry import Polygon
from openworld_radio_twin.environment import EnvironmentData
from openworld_radio_twin.models import BuildingFeature, GridMapping
from openworld_radio_twin.simulation.native_geometry import column_intervals
from openworld_radio_twin.simulation.scene_builder import (
_compiled_terrain_elevation,
localize_buildings,
)
def _rasterize_buildings(
buildings: list[BuildingFeature],
mapping: GridMapping,
width: int,
height: int,
slice_height_agl_m: float | None,
environment: EnvironmentData | None = None,
) -> tuple[np.ndarray, np.ndarray]:
if width <= 0 or height <= 0:
raise ValueError("Coverage dimensions must be positive")
localized = localize_buildings(
buildings,
mapping.origin_longitude,
mapping.origin_latitude,
environment,
)
east = mapping.first_sample_east_m + np.arange(width) * mapping.column_step_m
north = mapping.first_sample_north_m + np.arange(height) * mapping.row_step_m
mask = np.zeros((height, width), dtype=bool)
roofs = np.full((height, width), np.nan, dtype=np.float32)
for building in localized:
if building.native_vertices:
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 = (
np.ones(len(rows), dtype=bool)
if slice_height_agl_m is None
else (lower <= slice_height_agl_m) & (slice_height_agl_m <= upper)
)
rr, cc = rows[included], columns[included]
mask[rr, cc] = True
np.fmax.at(roofs, (rr, cc), upper[included])
continue
if slice_height_agl_m is not None and not (
building.minimum_height_m <= slice_height_agl_m <= building.roof_height_m
):
continue
polygon = Polygon(building.exterior, building.holes)
if polygon.is_empty or not polygon.is_valid or polygon.area <= 0:
continue
min_east, min_north, max_east, max_north = polygon.bounds
columns = np.flatnonzero((east >= min_east) & (east <= max_east))
rows = np.flatnonzero((north >= min_north) & (north <= max_north))
if not columns.size or not rows.size:
continue
sample_east, sample_north = np.meshgrid(east[columns], north[rows])
occupied = contains_xy(polygon, sample_east, sample_north)
region = np.ix_(rows, columns)
region_mask = mask[region]
region_roofs = roofs[region]
region_mask |= occupied
region_roofs[occupied] = np.fmax(region_roofs[occupied], building.roof_height_m)
mask[region] = region_mask
roofs[region] = region_roofs
roofs[~mask] = np.nan
return mask, roofs
[docs]
def rasterize_building_heights(
buildings: list[BuildingFeature],
mapping: GridMapping,
width: int,
height: int,
receiver_height_agl_m: float,
environment: EnvironmentData | None = None,
) -> tuple[np.ndarray, np.ndarray]:
"""Sample building occupancy at one AGL plane on exact grid cell centers."""
return _rasterize_buildings(
buildings,
mapping,
width,
height,
receiver_height_agl_m,
environment,
)
