import io
import json
import math
import tempfile
import xml.etree.ElementTree as ET
import zipfile
from pathlib import Path
import numpy as np
from PIL import Image
from openworld_radio_twin.environment import EnvironmentData
from openworld_radio_twin.models import SimulationRequest, SimulationResponse
from openworld_radio_twin.rendering import rgba_lookup
from openworld_radio_twin.scene_format import (
SCENE_ENVIRONMENT_ARRAY_PATHS,
SCENE_HEIGHT_MAP_PATH,
SCENE_MESH_DIRECTORY,
SCENE_PROVENANCE_PATH,
SCENE_TERRAIN_INFO_PATH,
SCENE_VOXEL_DEPTH_DIRECTORY,
SCENE_VOXEL_SLICES_DIRECTORY,
SCENE_XML_PATH,
)
from openworld_radio_twin.simulation.radio_map import RadioMapData
from openworld_radio_twin.simulation.scene_builder import (
SceneAssets,
build_scene_assets,
measurement_surface_path,
transmitter_scene_position,
write_measurement_surface,
write_voxel_assets,
)
REQUIRED_SINGLE_SCENE_PATHS = {
SCENE_HEIGHT_MAP_PATH.as_posix(),
"arrays/path_gain/path_gain.npz",
"arrays/rss/rss.npz",
"arrays/sinr/sinr.npz",
"images/association.png",
"images/path_gain_db.png",
"images/path_gain_db_masked.png",
"images/rss_db.png",
"images/rss_db_masked.png",
"images/sinr.png",
"images/sinr_masked.png",
(SCENE_MESH_DIRECTORY / "ground.ply").as_posix(),
SCENE_PROVENANCE_PATH.as_posix(),
SCENE_XML_PATH.as_posix(),
"tx_info.json",
(SCENE_VOXEL_SLICES_DIRECTORY / "scene_info.json").as_posix(),
(SCENE_VOXEL_DEPTH_DIRECTORY / "scene_depth_max_z_pitch3.0.npy").as_posix(),
(SCENE_VOXEL_DEPTH_DIRECTORY / "scene_depth_max_z_pitch3.0.png").as_posix(),
(SCENE_VOXEL_DEPTH_DIRECTORY / "scene_depth_max_z_pitch3.0_preview.png").as_posix(),
(SCENE_VOXEL_DEPTH_DIRECTORY / "scene_info.json").as_posix(),
}
MIN_RENDERED_IMAGE_SIDE = 1800
ASSOCIATION_COLORS = np.asarray(
[
[40, 199, 183],
[244, 162, 97],
[77, 121, 255],
[230, 92, 118],
[162, 112, 214],
[118, 164, 75],
[235, 198, 68],
[73, 172, 220],
],
dtype=np.uint8,
)
def _display_range(values: np.ndarray) -> tuple[float, float]:
finite = values[np.isfinite(values)]
if not finite.size:
return 0.0, 1.0
if finite.size >= 20:
minimum, maximum = np.percentile(finite, [2, 98]).tolist()
else:
minimum, maximum = float(np.min(finite)), float(np.max(finite))
if maximum <= minimum:
maximum = minimum + 1.0
return float(minimum), float(maximum)
def _save_clean_png(
values: np.ndarray,
path: Path,
*,
cmap: str = "viridis",
building_mask: np.ndarray | None = None,
) -> tuple[float, float]:
"""Render exact solver cells by integer pixel replication, without interpolation."""
if values.ndim != 2:
raise ValueError("PNG source must be a two-dimensional solver grid")
if building_mask is not None and building_mask.shape != values.shape:
raise ValueError("Building mask and PNG source grid must have the same shape")
minimum, maximum = _display_range(values)
normalized = np.clip((values - minimum) / max(1e-30, maximum - minimum), 0.0, 1.0)
lookup = rgba_lookup(cmap)
indices = np.nan_to_num(normalized, nan=0.0)
rgba = lookup[np.rint(indices * 255).astype(np.uint8)]
invalid = ~np.isfinite(values)
rgba[invalid] = (0, 0, 0, 0)
if building_mask is not None:
rgba[building_mask] = (0, 0, 0, 255)
rows, columns = values.shape
pixels_per_cell = max(1, math.ceil(MIN_RENDERED_IMAGE_SIDE / max(rows, columns)))
north_up = np.flipud(rgba)
rendered = np.repeat(
np.repeat(north_up, pixels_per_cell, axis=0),
pixels_per_cell,
axis=1,
)
Image.fromarray(rendered, mode="RGBA").save(path)
return minimum, maximum
[docs]
def render_radio_map_png(
data: RadioMapData,
metric: str,
building_mask: np.ndarray,
*,
mask_buildings: bool,
feather_edges: bool,
) -> bytes:
"""Render one pixel per solver cell without spatial resampling."""
if metric == "association":
values = data.association
valid = values >= 0
rgba = np.zeros((*values.shape, 4), dtype=np.uint8)
rgba[valid, :3] = ASSOCIATION_COLORS[values[valid] % len(ASSOCIATION_COLORS)]
else:
values = data.best_db(metric)
valid = np.isfinite(values)
minimum, maximum = _display_range(values)
normalized = np.clip((values - minimum) / (maximum - minimum), 0.0, 1.0)
lookup = rgba_lookup({"path_gain": "viridis", "rss": "inferno", "sinr": "jet"}[metric])
indices = np.rint(np.nan_to_num(normalized, nan=0.0) * 255).astype(np.uint8)
rgba = lookup[indices]
rgba[~valid] = (0, 0, 0, 0)
rgba[valid, 3] = 188
if feather_edges:
rows, columns = values.shape
x = np.linspace(-1.0, 1.0, columns) if columns > 1 else np.zeros(1)
y = np.linspace(-1.0, 1.0, rows) if rows > 1 else np.zeros(1)
edge = np.clip((1.0 - np.hypot(y[:, None], x[None, :])) / 0.2, 0.0, 1.0)
feather = edge * edge * (3.0 - 2.0 * edge)
rgba[:, :, 3] = np.rint(rgba[:, :, 3] * feather).astype(np.uint8)
if mask_buildings:
rgba[np.flipud(building_mask)] = (0, 0, 0, 230)
output = io.BytesIO()
Image.fromarray(np.flipud(rgba), mode="RGBA").save(output, format="PNG")
return output.getvalue()
def _look_at(
request: SimulationRequest,
tx_index: int,
position: tuple[float, float, float],
) -> list[float]:
tx = request.transmitters[tx_index]
azimuth = math.radians(tx.azimuth_deg)
downtilt = math.radians(tx.downtilt_deg)
horizontal = 100.0 * math.cos(downtilt)
return [
position[0] + horizontal * math.sin(azimuth),
position[1] + horizontal * math.cos(azimuth),
position[2] - 100.0 * math.sin(downtilt),
]
def _tx_info(
request: SimulationRequest,
result: SimulationResponse,
data: RadioMapData,
maximum_roof_height_m: float,
slice_path: Path,
slice_shape: tuple[int, int],
slice_pitch_m: float,
height_map_pixel_size_m: float,
image_products: dict[str, object],
scene: SceneAssets,
) -> dict[str, object]:
tx_list = []
for index, tx in enumerate(request.transmitters):
position = transmitter_scene_position(
scene,
tx.position.longitude,
tx.position.latitude,
tx.position.altitude_m,
)
tx_list.append(
{
"name": f"tx-{index}",
"position": list(position),
"look_at": _look_at(request, index, position),
"power_dbm": tx.power_dbm,
"frequency_ghz": tx.frequency_ghz,
"azimuth_deg": tx.azimuth_deg,
"downtilt_deg": tx.downtilt_deg,
"antenna_pattern": tx.antenna_pattern,
"geographic_position": {
"longitude": tx.position.longitude,
"latitude": tx.position.latitude,
"height_agl_m": tx.position.altitude_m,
},
}
)
sector_array = request.reference_transmitter.antenna_pattern == "sector"
mask_applied = result.coverage.building_cell_count > 0
return {
"simulation_id": result.simulation_id,
"request": request.model_dump(mode="json"),
"scene_name": "scene",
"scene_folder": ".",
"scene_xml": SCENE_XML_PATH.as_posix(),
"bbox_min_xyz": [-request.radius_m, -request.radius_m, 0.0],
"bbox_max_xyz": [
request.radius_m,
request.radius_m,
maximum_roof_height_m,
],
"rm_center": [0.0, 0.0, request.receiver_height_m],
"rm_size": [2.0 * request.radius_m, 2.0 * request.radius_m],
"rm_orientation": [0.0, 0.0, 0.0],
"params": {
"max_depth": request.max_depth,
"samples_per_tx": request.samples_per_tx,
"cell_size": [
abs(result.coverage.grid_mapping.column_step_m),
abs(result.coverage.grid_mapping.row_step_m),
],
"rm_height_above_ground": request.receiver_height_m,
"assoc_metric": request.association_metric,
"seed": request.seed,
"n_tx": len(request.transmitters),
"engine": result.engine,
},
"tx_array": {
"num_rows": 8 if sector_array else 1,
"num_cols": 2 if sector_array else 1,
"pattern": "tr38901" if sector_array else "iso",
"polarization": "V",
},
"rx_array": {
"num_rows": 1,
"num_cols": 1,
"pattern": "iso",
"polarization": "V",
},
"tx_list": tx_list,
"postprocess": {
"enabled": True,
"rm_height_above_ground": request.receiver_height_m,
"slice_png": str(slice_path),
"slice_shape": list(slice_shape),
"slice_pitch_m": slice_pitch_m,
"slice_row_order": "north_to_south",
"coverage_mask_source": "building_footprints_sampled_at_solver_cell_centers",
"mask_applied": mask_applied,
},
"array_contract": {
"shape": [data.transmitter_count, data.rows, data.columns],
"axis_order": ["transmitter", "south_to_north_row", "west_to_east_column"],
"path_gain": {
"path": "arrays/path_gain/path_gain.npz",
"key": "path_gain",
"unit": "linear ratio",
},
"rss": {"path": "arrays/rss/rss.npz", "key": "rss", "unit": "W"},
"sinr": {"path": "arrays/sinr/sinr.npz", "key": "sinr", "unit": "linear ratio"},
"image_logarithm": "10*log10(linear); RSS display images use dBm (+30 dB)",
"height_map": {
"path": SCENE_HEIGHT_MAP_PATH.as_posix(),
"dtype": "float32",
"value": "roof_height_agl_m",
"row_order": "north_to_south",
"pixel_size_m": height_map_pixel_size_m,
},
},
"image_products": image_products,
"grid_mapping": result.coverage.grid_mapping.model_dump(mode="json"),
"coordinate_system": {
"geographic_crs": "EPSG:4326",
"local_frame": "WGS84 topocentric ENU",
"origin_longitude": scene.origin_longitude,
"origin_latitude": scene.origin_latitude,
"horizontal_unit": "m",
"vertical_datum": "local ENU origin",
"vertical_unit": "m",
"transmitter_input_height_reference": "terrain AGL",
"transmitter_scene_z": "terrain_elevation_m + height_agl_m",
},
"engine_details": result.engine_details,
"data_provenance": result.data_provenance,
"warnings": result.warnings,
}
[docs]
def write_radio_products(
root: Path,
result: SimulationResponse,
data: RadioMapData,
building_mask: np.ndarray,
array_metrics: list[str] | tuple[str, ...] = ("path_gain", "rss", "sinr"),
image_metrics: list[str] | tuple[str, ...] = ("path_gain", "rss", "sinr"),
association_image: bool = True,
) -> dict[str, object]:
arrays = root / "arrays"
images = root / "images"
for metric in array_metrics:
(arrays / metric).mkdir(parents=True, exist_ok=True)
np.savez_compressed(arrays / metric / f"{metric}.npz", **{metric: data.linear(metric)})
if image_metrics or association_image:
images.mkdir(exist_ok=True)
building_mask = np.flipud(building_mask)
metric_specs = {
"path_gain": ("path_gain_db.png", "path_gain_db_masked.png", "viridis", "dB"),
"rss": ("rss_db.png", "rss_db_masked.png", "inferno", "dBm"),
"sinr": ("sinr.png", "sinr_masked.png", "jet", "dB"),
}
metric_products = {}
for metric in image_metrics:
plain_name, masked_name, cmap, unit = metric_specs[metric]
values = data.best_db(metric)
display_minimum, display_maximum = _save_clean_png(
values,
images / plain_name,
cmap=cmap,
)
_save_clean_png(
values,
images / masked_name,
cmap=cmap,
building_mask=building_mask,
)
metric_products[metric] = {
"unit": unit,
"colormap": cmap,
"display_minimum": display_minimum,
"display_maximum": display_maximum,
}
if association_image:
association = np.where(
data.association >= 0,
data.association.astype(np.float32),
np.nan,
)
_save_clean_png(
association,
images / "association.png",
cmap="viridis",
)
pixels_per_cell = max(
1,
math.ceil(MIN_RENDERED_IMAGE_SIDE / max(data.rows, data.columns)),
)
return {
"source_grid_shape": [data.rows, data.columns],
"rendered_size_px": [
data.columns * pixels_per_cell,
data.rows * pixels_per_cell,
],
"pixels_per_solver_cell": pixels_per_cell,
"row_order": "north_to_south",
"spatial_resampling": "none_integer_cell_replication",
"no_data_rendering": "transparent",
"masked_building_rendering": "opaque_black",
"value_source": "best_transmitter_from_in_memory_linear_solver_output",
"array_metrics": list(array_metrics),
"image_metrics": list(image_metrics),
"association_image": association_image,
"metrics": metric_products,
}
[docs]
def validate_artifact_tree(root: Path, data: RadioMapData) -> None:
"""Fail export when the portable tree diverges from the declared contract."""
paths = {path.relative_to(root).as_posix() for path in root.rglob("*") if path.is_file()}
missing = REQUIRED_SINGLE_SCENE_PATHS - paths
if missing:
raise RuntimeError(f"Artifact tree is missing: {', '.join(sorted(missing))}")
expected_shape = data.shape
for metric in ("path_gain", "rss", "sinr"):
path = root / "arrays" / metric / f"{metric}.npz"
with np.load(path) as archive:
if archive.files != [metric]:
raise RuntimeError(f"{path.name} must contain only the '{metric}' key")
values = archive[metric]
if values.shape != expected_shape or values.dtype != np.float32:
raise RuntimeError(
f"{path.name} has {values.dtype} {values.shape}; "
f"expected float32 {expected_shape}"
)
pixels_per_cell = max(
1,
math.ceil(MIN_RENDERED_IMAGE_SIDE / max(data.rows, data.columns)),
)
expected_image_size = (
data.columns * pixels_per_cell,
data.rows * pixels_per_cell,
)
for image_path in (root / "images").glob("*.png"):
with Image.open(image_path) as image:
if image.size != expected_image_size or image.mode != "RGBA":
raise RuntimeError(f"{image_path.name} must be an {expected_image_size} RGBA image")
depth = np.load(root / SCENE_VOXEL_DEPTH_DIRECTORY / "scene_depth_max_z_pitch3.0.npy")
with Image.open(root / SCENE_VOXEL_DEPTH_DIRECTORY / "scene_depth_max_z_pitch3.0.png") as image:
if image.size != (depth.shape[1], depth.shape[0]) or image.mode not in {"I", "I;16"}:
raise RuntimeError("Voxel depth PNG must match its numeric depth grid")
with Image.open(next((root / SCENE_VOXEL_SLICES_DIRECTORY).glob("scene_slice_*.png"))) as image:
if image.size != (depth.shape[1], depth.shape[0]):
raise RuntimeError("Voxel depth and receiver slice must use the same XY grid")
for shape in ET.parse(root / SCENE_XML_PATH).getroot().findall("shape"):
filename = shape.find("string[@name='filename']")
if filename is None or not (root / filename.attrib["value"]).is_file():
raise RuntimeError(f"scene.xml shape {shape.attrib.get('id')} has no usable mesh")
metadata = json.loads((root / "tx_info.json").read_text(encoding="utf-8"))
if metadata["array_contract"]["shape"] != list(expected_shape):
raise RuntimeError("tx_info.json array shape does not match the NPZ products")
[docs]
def build_artifact_bundle(
request: SimulationRequest,
result: SimulationResponse,
data: RadioMapData,
building_mask: np.ndarray,
environment: EnvironmentData | None = None,
) -> bytes:
"""Build a portable single-scene version of the reference/package layout."""
with tempfile.TemporaryDirectory(prefix="owrt-artifact-") as temp_directory:
root = Path(temp_directory)
reference = request.scene_origin
scene = build_scene_assets(
root,
result.buildings,
reference.longitude,
reference.latitude,
request.radius_m,
environment,
True,
request.material_profile,
)
measurement = None
if environment is not None:
measurement = write_measurement_surface(
measurement_surface_path(
scene,
request.resolution_m,
request.receiver_height_m,
),
scene,
request.resolution_m,
request.receiver_height_m,
)
voxels = write_voxel_assets(scene, request.receiver_height_m)
image_products = write_radio_products(root, result, data, building_mask)
maximum_roof_height = max(
(building.roof_height_m for building in scene.local_buildings), default=0.0
)
info = _tx_info(
request,
result,
data,
maximum_roof_height,
voxels.slice_path.relative_to(root),
voxels.slice_shape,
voxels.pitch_m,
scene.height_map_pixel_size_m,
image_products,
scene,
)
info["postprocess"]["voxel_depth"] = {
"array": str(voxels.depth_array_path.relative_to(root)),
"png": str(voxels.depth_image_path.relative_to(root)),
"preview_png": str(voxels.depth_preview_path.relative_to(root)),
"shape": list(voxels.slice_shape),
"value": "maximum_occupied_voxel_z_agl_m",
"pitch_m": voxels.pitch_m,
}
if environment is not None:
assert measurement is not None
info["measurement_surface"] = {
"path": str(measurement.path.relative_to(root)),
"cell_shape": list(measurement.cell_shape),
"receiver_height_agl_m": request.receiver_height_m,
"reduction": "two triangles per cell; surface-area-weighted linear mean",
}
info["environment"] = {
"metadata": SCENE_TERRAIN_INFO_PATH.as_posix(),
**{name: path.as_posix() for name, path in SCENE_ENVIRONMENT_ARRAY_PATHS.items()},
}
(root / "tx_info.json").write_text(json.dumps(info, indent=2), encoding="utf-8")
validate_artifact_tree(root, data)
output = io.BytesIO()
with zipfile.ZipFile(output, "w", compression=zipfile.ZIP_DEFLATED) as archive:
for path in sorted(root.rglob("*")):
if path.is_file():
archive.write(path, path.relative_to(root).as_posix())
return output.getvalue()
