from dataclasses import dataclass
import numpy as np
from openworld_radio_twin.environment import EnvironmentData
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.models import (
BuildingFeature,
CoverageGrid,
GridMapping,
RadioLayer,
SimulationRequest,
)
from openworld_radio_twin.simulation.building_grid import rasterize_building_heights
[docs]
@dataclass(frozen=True)
class MetricDefinition:
label: str
unit: str
logarithmic_offset_db: float
METRICS = {
"path_gain": MetricDefinition("Path gain", "dB", 0.0),
"rss": MetricDefinition("Received power", "dBm", 30.0),
"sinr": MetricDefinition("SINR", "dB", 0.0),
}
[docs]
@dataclass(frozen=True)
class CoverageProducts:
coverage: CoverageGrid
building_mask: np.ndarray
building_heights_m: np.ndarray
[docs]
def finite_max(values: np.ndarray, axis: int) -> np.ndarray:
"""Reduce finite values while preserving an all-missing cell as NaN."""
finite = np.isfinite(values)
reduced = np.max(np.where(finite, values, -np.inf), axis=axis)
return np.where(np.any(finite, axis=axis), reduced, np.nan)
[docs]
@dataclass(frozen=True)
class RadioMapData:
"""Canonical solver products using the reference package's array contract.
All metric arrays are linear float32 values with shape ``(n_tx, rows, columns)``.
Row zero is the south edge so the arrays can be plotted with ``origin="lower"``,
exactly as in ``reference/package``. Coordinates are local ENU meters.
"""
path_gain: np.ndarray
rss: np.ndarray
sinr: np.ndarray
association: np.ndarray
cell_centers: np.ndarray
sample_layout: str = "cell_average"
def __post_init__(self) -> None:
arrays = {
"path_gain": self.path_gain,
"rss": self.rss,
"sinr": self.sinr,
}
shape: tuple[int, int, int] | None = None
for name, value in arrays.items():
if value.ndim != 3:
raise ValueError(f"{name} must have shape (n_tx, rows, columns)")
if shape is None:
shape = value.shape
elif value.shape != shape:
raise ValueError(f"{name} shape {value.shape} does not match {shape}")
if np.any(value < 0):
raise ValueError(f"{name} contains negative linear values")
assert shape is not None
if self.association.shape != shape[1:]:
raise ValueError("association must have shape (rows, columns)")
if self.cell_centers.shape != (*shape[1:], 3):
raise ValueError("cell_centers must have shape (rows, columns, 3)")
if shape[1] > 1 and self.cell_centers[0, 0, 1] > self.cell_centers[-1, 0, 1]:
raise ValueError("canonical arrays must be ordered south-to-north")
@property
def shape(self) -> tuple[int, int, int]:
return self.path_gain.shape
@property
def transmitter_count(self) -> int:
return self.shape[0]
@property
def rows(self) -> int:
return self.shape[1]
@property
def columns(self) -> int:
return self.shape[2]
[docs]
def linear(self, metric: str) -> np.ndarray:
if metric not in METRICS:
raise ValueError(f"Unknown radio metric: {metric}")
return getattr(self, metric)
[docs]
def best_db(self, metric: str, *, zero_floor: bool = False) -> np.ndarray:
definition = METRICS[metric]
linear = self.linear(metric).astype(np.float64, copy=False)
if zero_floor:
with np.errstate(divide="ignore", invalid="ignore"):
db = 10.0 * np.log10(np.maximum(linear, 1e-30))
else:
with np.errstate(divide="ignore", invalid="ignore"):
db = np.where(linear > 0, 10.0 * np.log10(linear), np.nan)
return finite_max(db, axis=0) + definition.logarithmic_offset_db
[docs]
def normalize_solver_arrays(
path_gain: np.ndarray,
rss: np.ndarray,
sinr: np.ndarray,
association: np.ndarray,
cell_centers: np.ndarray,
) -> RadioMapData:
"""Normalize a solver result to float32 and south-to-north row order."""
metrics = [np.asarray(value, dtype=np.float32) for value in (path_gain, rss, sinr)]
metrics = [value[np.newaxis, ...] if value.ndim == 2 else value for value in metrics]
centers = np.asarray(cell_centers, dtype=np.float64)
assoc = np.asarray(association, dtype=np.int32).reshape(metrics[0].shape[1:])
if centers.shape != (*metrics[0].shape[1:], 3):
centers = centers.reshape(*metrics[0].shape[1:], 3)
if metrics[0].shape[1] > 1 and centers[0, 0, 1] > centers[-1, 0, 1]:
metrics = [np.flip(value, axis=1).copy() for value in metrics]
assoc = np.flipud(assoc).copy()
centers = np.flipud(centers).copy()
return RadioMapData(*metrics, assoc, centers)
def _display_limits(values: np.ndarray) -> tuple[float, float, float, float]:
finite = values[np.isfinite(values)]
if finite.size == 0:
raise RuntimeError("Radio map contains no finite samples")
minimum = float(np.min(finite))
maximum = float(np.max(finite))
if finite.size >= 20:
display_minimum, display_maximum = np.percentile(finite, [2, 98]).tolist()
else:
display_minimum, display_maximum = minimum, maximum
if display_maximum <= display_minimum:
display_maximum = display_minimum + 1.0
return minimum, maximum, float(display_minimum), float(display_maximum)
[docs]
def build_coverage_grid(
request: SimulationRequest,
data: RadioMapData,
buildings: list[BuildingFeature],
environment: EnvironmentData | None = None,
) -> CoverageGrid:
return build_coverage_products(request, data, buildings, environment).coverage
[docs]
def build_coverage_products(
request: SimulationRequest,
data: RadioMapData,
buildings: list[BuildingFeature],
environment: EnvironmentData | None = None,
) -> CoverageProducts:
"""Create the north-up JSON view from the same canonical arrays used for export."""
reference = request.scene_origin
frame = LocalFrame.at(reference.longitude, reference.latitude)
centers = np.flipud(data.cell_centers)
column_step = (
float(centers[0, 1, 0] - centers[0, 0, 0])
if data.columns > 1
else float(request.resolution_m)
)
row_step = (
float(centers[1, 0, 1] - centers[0, 0, 1])
if data.rows > 1
else -float(request.resolution_m)
)
west_m = float(np.min(centers[..., 0]) - abs(column_step) / 2)
east_m = float(np.max(centers[..., 0]) + abs(column_step) / 2)
south_m = float(np.min(centers[..., 1]) - abs(row_step) / 2)
north_m = float(np.max(centers[..., 1]) + abs(row_step) / 2)
corners = frame.corners(west_m, south_m, east_m, north_m)
longitudes = [point[0] for point in corners]
latitudes = [point[1] for point in corners]
mapping = GridMapping(
origin_longitude=reference.longitude,
origin_latitude=reference.latitude,
first_sample_east_m=float(centers[0, 0, 0]),
first_sample_north_m=float(centers[0, 0, 1]),
column_step_m=column_step,
row_step_m=row_step,
sample_layout=data.sample_layout,
corner_coordinates=corners,
)
layers: dict[str, RadioLayer] = {}
for metric, definition in METRICS.items():
values = np.flipud(data.best_db(metric))
minimum, maximum, display_minimum, display_maximum = _display_limits(values)
layers[metric] = RadioLayer(
label=definition.label,
unit=definition.unit,
minimum=minimum,
maximum=maximum,
display_minimum=display_minimum,
display_maximum=display_maximum,
)
building_mask, building_heights_m = rasterize_building_heights(
buildings,
mapping,
data.columns,
data.rows,
request.receiver_height_m,
environment,
)
coverage = CoverageGrid(
west=min(longitudes),
south=min(latitudes),
east=max(longitudes),
north=max(latitudes),
width=data.columns,
height=data.rows,
receiver_height_agl_m=request.receiver_height_m,
grid_mapping=mapping,
layers=layers,
building_cell_count=int(building_mask.sum()),
)
return CoverageProducts(coverage, building_mask, building_heights_m)
