Geographic scenes with the native Sionna RT workflow

This tutorial loads a WGS84 region as a native sionna.rt.Scene, previews it with Sionna’s interactive WebGL viewer, and compares radio maps produced from two environment configurations. OWRT only extends scene loading, geographic transmitter placement, and terrain-following measurement surfaces; arrays, radio devices, solvers, radio maps, and previews remain Sionna objects.

First use downloads public geodata and compiles a cached scene. Radio-map cells execute only when requested. A CUDA-capable environment is recommended.

Install

Activate the owrt environment from the repository root, install the ray-tracing and tutorial extras into that environment, then open this notebook.

conda activate owrt
pip install -e ".[rt,tutorial]"
jupyter lab examples/tutorial_sionna_scene_comparison.ipynb

Scene controls

Control	Choices	Effect
buildings	"auto", "none"	Load source-derived buildings or omit them
terrain	"elevation", "flat", "none"	Use a DEM, a planar ground, or no ground geometry
material_profile	"itu", "uniform"	Map semantic surfaces to RF proxies or use one material
terrain_resolution_m	positive metres	Set DEM sampling and terrain-mesh spacing
device	"auto", "cuda", "cpu"	Select the Sionna/Mitsuba execution backend

The first four choices identify compiled geometry in the cache. Frequency, arrays, transmitters, ray budget, depth, and seed remain ordinary Sionna settings. The quality profile below changes geometry sampling, receiver-grid density, propagation sampling, and render resolution together. It cannot add facade or roof detail absent from the source building geometry.

import json
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np

from openworld_radio_twin import rt as owrt

# Select the backend before importing sionna.rt.
DEVICE = "auto"  # "auto" | "cuda" | "cpu"
CACHE_DIR = Path(".cache/owrt-tutorial")
LOCATION = dict(latitude=52.3762, longitude=4.8993, radius_m=250)
TX = dict(latitude=52.3770, longitude=4.9003, height_agl=25.0)
FREQUENCY_HZ = 3.5e9
TX_POWER_DBM = 30.0
RX_HEIGHT_M = 1.5
SEED = 2026

QUALITY = "high"  # "draft" | "high" | "ultra"
QUALITY_PRESETS = {
    "draft": dict(
        terrain_resolution_m=5.0,
        cell_size=(5.0, 5.0),
        samples_per_tx=1_000_000,
        max_depth=5,
        preview_resolution=(800, 520),
        render_resolution=(800, 520),
        render_samples=128,
    ),
    "high": dict(
        terrain_resolution_m=2.0,
        cell_size=(2.0, 2.0),
        samples_per_tx=20_000_000,
        max_depth=8,
        preview_resolution=(1100, 700),
        render_resolution=(1200, 760),
        render_samples=512,
    ),
    "ultra": dict(
        terrain_resolution_m=1.0,
        cell_size=(1.0, 1.0),
        samples_per_tx=100_000_000,
        max_depth=8,
        preview_resolution=(1400, 900),
        render_resolution=(1600, 1000),
        render_samples=1024,
    ),
}
quality = QUALITY_PRESETS[QUALITY]
TERRAIN_RESOLUTION_M = quality["terrain_resolution_m"]
CELL_SIZE = quality["cell_size"]
SAMPLES_PER_TX = quality["samples_per_tx"]
MAX_DEPTH = quality["max_depth"]
PREVIEW_RESOLUTION = quality["preview_resolution"]
RENDER_RESOLUTION = quality["render_resolution"]
RENDER_SAMPLES = quality["render_samples"]

receiver_cells = round(
    (2 * LOCATION["radius_m"] / CELL_SIZE[0]) * (2 * LOCATION["radius_m"] / CELL_SIZE[1])
)
print(
    f"Quality: {QUALITY}; receiver cells: {receiver_cells:,}; "
    f"rays/Tx: {SAMPLES_PER_TX:,}; depth: {MAX_DEPTH}"
)

SCENARIOS = {
    "Full context": dict(
        buildings="auto",
        terrain="elevation",
        material_profile="itu",
        terrain_resolution_m=TERRAIN_RESOLUTION_M,
    ),
    "Flat terrain": dict(
        buildings="auto",
        terrain="flat",
        material_profile="itu",
        terrain_resolution_m=TERRAIN_RESOLUTION_M,
    ),
    "No buildings": dict(
        buildings="none",
        terrain="elevation",
        material_profile="itu",
        terrain_resolution_m=TERRAIN_RESOLUTION_M,
    ),
    "Uniform materials": dict(
        buildings="auto",
        terrain="elevation",
        material_profile="uniform",
        terrain_resolution_m=TERRAIN_RESOLUTION_M,
    ),
    "Free space": dict(
        buildings="none",
        terrain="none",
        material_profile="uniform",
        terrain_resolution_m=TERRAIN_RESOLUTION_M,
    ),
}

SCENARIOS

Load a native scene

owrt.load_scene(...) selects the execution backend, compiles or reuses the geographic scene, and delegates to Sionna’s loader. The returned object is not wrapped or subclassed.

reference_scene = owrt.load_scene(
    **LOCATION,
    **SCENARIOS["Full context"],
    cache_dir=CACHE_DIR,
    device=DEVICE,
)

# Sionna is imported after OWRT has selected the Mitsuba backend.
import sionna.rt as rt  # noqa: E402

print(type(reference_scene))
print(owrt.scene_path(reference_scene))

_scenes = {}
_radio_maps = {}


def _configure_scene(scene):
    scene.frequency = FREQUENCY_HZ
    array = rt.PlanarArray(num_rows=1, num_cols=1, pattern="iso", polarization="V")
    scene.tx_array = array
    scene.rx_array = array
    scene.add(
        rt.Transmitter(
            name="tutorial-tx",
            position=owrt.position(scene, **TX),
            power_dbm=TX_POWER_DBM,
        )
    )
    return scene


_scenes["Full context"] = _configure_scene(reference_scene)


def get_scene(name):
    if name not in SCENARIOS:
        raise KeyError(f"Unknown scenario: {name}")
    if name not in _scenes:
        scene = owrt.load_scene(**LOCATION, **SCENARIOS[name], cache_dir=CACHE_DIR, device=DEVICE)
        _scenes[name] = _configure_scene(scene)
    return _scenes[name]


def solve(name):
    if name not in _radio_maps:
        scene = get_scene(name)
        surface = owrt.measurement_surface(scene, cell_size=CELL_SIZE, height=RX_HEIGHT_M)
        _radio_maps[name] = rt.RadioMapSolver()(
            scene,
            measurement_surface=surface,
            samples_per_tx=SAMPLES_PER_TX,
            max_depth=MAX_DEPTH,
            seed=SEED,
        )
    return get_scene(name), _radio_maps[name]

Optional native Sionna preview

Set the three variables below and run the cell. The returned viewer is Sionna’s native interactive preview: drag to orbit, scroll to zoom, and use its controls to inspect geometry and coverage. Each scenario is compiled and solved at most once per kernel session.

scene.preview(...) requires a working Jupyter widget frontend. If VS Code reports a jupyter-ipywidget-renderer or jupyter-threejs error, use the Plotly 3D cell below or open this notebook with JupyterLab from the same owrt environment.

PREVIEW_SCENARIO = "Full context"  # any key in SCENARIOS
PREVIEW_METRIC = "rss"  # "rss" | "path_gain" | "sinr"
OVERLAY_RADIO_MAP = True
RUN_NATIVE_PREVIEW = False  # set True in JupyterLab when jupyter-threejs works

if RUN_NATIVE_PREVIEW and OVERLAY_RADIO_MAP:
    preview_scene, preview_radio_map = solve(PREVIEW_SCENARIO)
    preview_scene.preview(
        radio_map=preview_radio_map,
        rm_metric=PREVIEW_METRIC,
        resolution=PREVIEW_RESOLUTION,
        rm_cmap="viridis",
    )
elif RUN_NATIVE_PREVIEW:
    get_scene(PREVIEW_SCENARIO).preview(resolution=PREVIEW_RESOLUTION)
else:
    print("Native widget skipped; run the Plotly 3D cell below.")

VS Code-compatible interactive 3D preview

This Plotly viewer reads the exact compiled OWRT meshes and overlays the native Sionna radio map. It supports orbit, pan, zoom, hover values, legend-based layer toggles, and camera export without ipywidgets or jupyter-threejs. The solver result is unchanged; max_mesh_faces only bounds browser rendering cost.

from openworld_radio_twin.plotly_preview import interactive_preview

plotly_scene, plotly_radio_map = solve(PREVIEW_SCENARIO)
plotly_figure = interactive_preview(
    plotly_scene,
    radio_map=plotly_radio_map,
    metric=PREVIEW_METRIC,
    cell_size=CELL_SIZE,
    receiver_height_m=RX_HEIGHT_M,
    title=f"{PREVIEW_SCENARIO} · {PREVIEW_METRIC.upper()} · {QUALITY}",
    max_mesh_faces=600_000,
)
plotly_figure.show(config={"scrollZoom": True, "displaylogo": False})

Controlled radio-map comparison

The comparison below holds the WGS84 extent, transmitter, frequency, receiver height, grid, ray budget, depth, and seed fixed. It reduces Sionna’s two triangular mesh bins per logical receiver cell with the same area-weighted linear mean used by the dataset exporter. The first two panels share a color scale; the third reports right - left.

METRIC_LABELS = {
    "rss": ("Received power", "dBm", 30.0),
    "path_gain": ("Path gain", "dB", 0.0),
    "sinr": ("SINR", "dB", 0.0),
}


def _measurement_record(scene):
    info_path = owrt.scene_path(scene) / "scene_info.json"
    info = json.loads(info_path.read_text(encoding="utf-8"))
    records = info["shared_assets"]["measurement_surfaces"]
    for record in reversed(records):
        same_cell = np.allclose(record["cell_size_m"], CELL_SIZE)
        same_height = np.isclose(record["receiver_height_agl_m"], RX_HEIGHT_M)
        if same_cell and same_height:
            return record
    raise RuntimeError("The requested measurement-surface metadata was not found")


def metric_grid(scene, radio_map, metric):
    record = _measurement_record(scene)
    with np.load(owrt.scene_path(scene) / record["arrays"]) as arrays:
        weights = arrays["triangle_areas_m2"].astype(np.float64)
    rows, columns = weights.shape[:2]

    values = np.asarray(getattr(radio_map, metric).numpy(), dtype=np.float64)
    if values.ndim == 1:
        values = values[None, :]
    triangles = values.reshape(values.shape[0], rows, columns, 2)
    linear = np.sum(triangles * weights[None, ...], axis=-1)
    linear /= np.sum(weights, axis=-1)[None, ...]
    strongest = np.max(linear, axis=0)

    offset = METRIC_LABELS[metric][2]
    with np.errstate(divide="ignore", invalid="ignore"):
        return np.where(strongest > 0, 10.0 * np.log10(strongest) + offset, np.nan)


def compare(left_name, right_name, metric="rss"):
    left_scene, left_rm = solve(left_name)
    right_scene, right_rm = solve(right_name)
    left = metric_grid(left_scene, left_rm, metric)
    right = metric_grid(right_scene, right_rm, metric)
    common = np.isfinite(left) & np.isfinite(right)
    if not np.any(common):
        raise RuntimeError("The selected radio maps have no common finite cells")
    delta = np.where(common, right - left, np.nan)

    pooled = np.concatenate([left[np.isfinite(left)], right[np.isfinite(right)]])
    vmin, vmax = np.nanpercentile(pooled, [2, 98])
    delta_limit = max(float(np.nanpercentile(np.abs(delta[common]), 98)), 0.5)
    extent = (
        -LOCATION["radius_m"],
        LOCATION["radius_m"],
        -LOCATION["radius_m"],
        LOCATION["radius_m"],
    )

    field_cmap = plt.colormaps["viridis"].with_extremes(bad="#ECEFF1")
    delta_cmap = plt.colormaps["RdBu_r"].with_extremes(bad="#ECEFF1")
    fig, axes = plt.subplots(1, 3, figsize=(13.2, 4.1), constrained_layout=True)
    first = axes[0].imshow(
        left, origin="lower", extent=extent, cmap=field_cmap, vmin=vmin, vmax=vmax
    )
    axes[1].imshow(right, origin="lower", extent=extent, cmap=field_cmap, vmin=vmin, vmax=vmax)
    change = axes[2].imshow(
        delta, origin="lower", extent=extent, cmap=delta_cmap, vmin=-delta_limit, vmax=delta_limit
    )
    axes[0].set_title(left_name)
    axes[1].set_title(right_name)
    axes[2].set_title(f"Change: {right_name} - {left_name}")
    for axis in axes:
        axis.set_xlabel("East (m)")
        axis.set_ylabel("North (m)")
        axis.set_aspect("equal")
    label, unit, _ = METRIC_LABELS[metric]
    fig.colorbar(first, ax=axes[:2], shrink=0.88, label=f"{label} ({unit})")
    fig.colorbar(change, ax=axes[2], shrink=0.88, label=f"Difference ({unit})")
    fig.suptitle("Environment-controlled Sionna RT comparison", fontweight="bold")

    absolute = np.abs(delta[common])
    print(f"Common finite cells: {common.sum():,} / {common.size:,}")
    print(f"Median |difference|: {np.median(absolute):.3f} {unit}")
    print(f"95th percentile |difference|: {np.percentile(absolute, 95):.3f} {unit}")
    return fig

LEFT_SCENARIO = "Full context"
RIGHT_SCENARIO = "Flat terrain"
COMPARISON_METRIC = "rss"

if LEFT_SCENARIO == RIGHT_SCENARIO:
    raise ValueError("Choose two different configurations")
compare(LEFT_SCENARIO, RIGHT_SCENARIO, COMPARISON_METRIC)

Interpreting the comparison

Use one controlled change at a time: compare Full context with Flat terrain for terrain geometry, No buildings for building geometry, or Uniform materials for semantic material assignment. The tutorial statistics describe this run only. Increase SAMPLES_PER_TX, repeat seeds, and report uncertainty before treating differences as experimental evidence.

Compiled scenes persist under CACHE_DIR; rerunning the notebook reuses matching geometry. Changing frequency, transmitter power, arrays, solver depth, or ray budget does not recompile the scene.