from __future__ import annotations
import asyncio
from collections import deque
from concurrent.futures import ThreadPoolExecutor
from dataclasses import dataclass
import numpy as np
from shapely import contains_xy, from_wkb, intersects_xy, union_all
from shapely.geometry import LineString, MultiPolygon, Polygon, box
from shapely.ops import substring, transform
from openworld_radio_twin.config import Settings
from openworld_radio_twin.environment import (
FeatureSource,
SurfaceClass,
SurfaceFeature,
TerrainGrid,
)
from openworld_radio_twin.geodesy import LocalFrame
from openworld_radio_twin.providers.overture_reader import (
SEGMENT_COLUMNS,
SURFACE_COLUMNS,
latest_release,
read_batches,
)
[docs]
@dataclass(frozen=True)
class SurfaceQueryResult:
features: tuple[SurfaceFeature, ...]
cells: np.ndarray
release: str
provenance: dict[str, object]
warnings: list[str]
[docs]
class SurfaceProviderError(RuntimeError):
pass
[docs]
@dataclass(frozen=True)
class SurfaceDecision:
surface_class: SurfaceClass
basis: str
precedence: int
[docs]
@dataclass(frozen=True)
class LinearRule:
start: float
stop: float
value: str | float
PAVED_SURFACES = frozenset(
{
"asphalt",
"cobblestone",
"concrete",
"concrete_plates",
"paved",
"paving_stones",
"recreation_paved",
"rubber",
"sett",
"tartan",
}
)
VEGETATION_SURFACES = frozenset({"grass", "recreation_grass"})
GROUND_SURFACES = frozenset(
{
"compacted",
"dirt",
"earth",
"fine_gravel",
"gravel",
"ground",
"pebblestone",
"recreation_sand",
"sand",
"unpaved",
"wood",
"woodchips",
}
)
VEGETATION_LAND_COVER = frozenset(
{"crop", "forest", "grass", "mangrove", "moss", "shrub", "wetland"}
)
GROUND_LAND_COVER = frozenset({"barren", "snow", "urban"})
VEGETATION_LAND_USE_SUBTYPES = frozenset(
{"agriculture", "cemetery", "golf", "grass", "horticulture", "park"}
)
VEGETATION_LAND_USE_CLASSES = frozenset(
{
"allotments",
"cemetery",
"dog_park",
"driving_range",
"fairway",
"farmland",
"flowerbed",
"forest",
"garden",
"golf_course",
"grass",
"green",
"meadow",
"national_park",
"nature_reserve",
"orchard",
"park",
"plant_nursery",
"recreation_ground",
"rough",
"state_park",
"tee",
"village_green",
"vineyard",
"wilderness_area",
}
)
PAVED_ROAD_SURFACES = frozenset({"metal", "paved", "paving_stones"})
GROUND_ROAD_SURFACES = frozenset({"dirt", "gravel", "unpaved"})
GROUND_ROAD_CLASSES = frozenset({"bridleway", "path", "track"})
MARINE_WATER_CLASSES = frozenset({"bay", "ocean", "sea", "strait"})
ROAD_WIDTH_DEFAULTS_M = {
"motorway": 14.0,
"trunk": 12.0,
"primary": 10.0,
"secondary": 9.0,
"tertiary": 8.0,
"residential": 7.0,
"living_street": 6.0,
"unclassified": 6.0,
"service": 5.0,
"pedestrian": 4.0,
"track": 3.0,
"cycleway": 3.0,
"footway": 2.0,
"path": 2.0,
"bridleway": 2.0,
"steps": 2.0,
"other": 6.0,
}
SURFACE_MAPPING_CONTRACT = {
"version": 1,
"precedence": [
"land_use_semantic",
"land_cover_physical",
"explicit_land_use_surface",
"inferred_road_surface",
"explicit_road_surface",
"water_geometry",
],
"surface_material": {
"paved": sorted(PAVED_SURFACES),
"vegetation": sorted(VEGETATION_SURFACES),
"ground": sorted(GROUND_SURFACES),
},
"land_cover": {
"vegetation": sorted(VEGETATION_LAND_COVER),
"ground": sorted(GROUND_LAND_COVER),
},
"land_use": {
"vegetation_subtypes": sorted(VEGETATION_LAND_USE_SUBTYPES),
"vegetation_classes": sorted(VEGETATION_LAND_USE_CLASSES),
},
"road_surface": {
"paved": sorted(PAVED_ROAD_SURFACES),
"ground": sorted(GROUND_ROAD_SURFACES),
"class_inferred_ground": sorted(GROUND_ROAD_CLASSES),
"otherwise_inferred": "paved",
},
}
[docs]
class OvertureSurfaceProvider:
"""Normalize Overture semantic layers into one non-overlapping surface grid."""
themes = ("land_cover", "land_use", "water", "segment")
def __init__(self, settings: Settings) -> None:
self.settings = settings
self._release: str | None = None
self._query_lock = asyncio.Lock()
[docs]
async def surfaces(
self,
longitude: float,
latitude: float,
radius_m: float,
terrain: TerrainGrid,
) -> SurfaceQueryResult:
async with self._query_lock:
try:
return await asyncio.to_thread(
self._query,
longitude,
latitude,
radius_m,
terrain,
)
except Exception as exc:
raise SurfaceProviderError(f"Overture surface query failed: {exc}") from exc
def _query(
self,
longitude: float,
latitude: float,
radius_m: float,
terrain: TerrainGrid,
) -> SurfaceQueryResult:
if self._release is None:
self._release = latest_release()
assert self._release is not None
frame = LocalFrame.at(longitude, 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]
bbox = min(longitudes), min(latitudes), max(longitudes), max(latitudes)
clip = box(-radius_m, -radius_m, radius_m, radius_m)
features: list[SurfaceFeature] = []
counts: dict[str, int] = {}
warnings = []
# The four themes are independent objects on S3; read and decode them at once
# and keep the theme order when assembling the result.
with ThreadPoolExecutor(
max_workers=len(self.themes), thread_name_prefix="owrt-surfaces"
) as pool:
futures = {
theme: pool.submit(self._query_theme, theme, bbox, longitude, latitude, clip)
for theme in self.themes
}
results = {theme: future.result() for theme, future in futures.items()}
for theme in self.themes:
features.extend(results[theme])
counts[theme] = len(results[theme])
cells = classify_surface_cells(terrain, features)
inferred_widths = sum(item.width_source == "class_default" for item in features)
inferred_surfaces = sum(
item.classification_basis == "road_class_default" for item in features
)
classification_counts: dict[str, int] = {}
for feature in features:
classification_counts[feature.classification_basis] = (
classification_counts.get(feature.classification_basis, 0) + 1
)
source_datasets = sorted(
{
(source.dataset, source.license)
for feature in features
for source in feature.sources
if source.dataset is not None
},
key=lambda item: (item[0], item[1] or ""),
)
update_times = sorted(
source.update_time
for feature in features
for source in feature.sources
if source.update_time is not None
)
provenance: dict[str, object] = {
"provider": "Overture Maps",
"release": self._release,
"classification": "cell-center point-in-polygon with evidence precedence",
"mapping_contract": SURFACE_MAPPING_CONTRACT,
"classification_basis_counts": classification_counts,
"feature_count": len(features),
"inferred_road_width_count": inferred_widths,
"inferred_road_surface_count": inferred_surfaces,
"source_datasets": [
{"dataset": dataset, "license": license_name}
for dataset, license_name in source_datasets
],
"source_update_time_range": (
[update_times[0], update_times[-1]] if update_times else None
),
"road_width_defaults_m": ROAD_WIDTH_DEFAULTS_M,
**{f"{theme}_feature_count": count for theme, count in counts.items()},
}
if inferred_widths:
warnings.append(
f"{inferred_widths} road sections lacked source width; deterministic class widths "
"are recorded as assumptions."
)
if inferred_surfaces:
warnings.append(
f"{inferred_surfaces} road sections lacked source surface; deterministic "
"road-class surfaces are recorded as assumptions."
)
return SurfaceQueryResult(
features=tuple(features),
cells=cells,
release=self._release,
provenance=provenance,
warnings=warnings,
)
def _query_theme(
self,
theme: str,
bbox: tuple[float, float, float, float],
longitude: float,
latitude: float,
clip: Polygon,
) -> list[SurfaceFeature]:
batches = read_batches(
theme,
bbox,
self._release,
SEGMENT_COLUMNS if theme == "segment" else SURFACE_COLUMNS,
connect_timeout=self.settings.overture_connect_timeout_seconds,
request_timeout=self.settings.overture_request_timeout_seconds,
)
# Transformers are not shared between threads; each theme builds its own frame.
return self._read_theme(batches, theme, LocalFrame.at(longitude, latitude), clip)
def _read_theme(
self,
reader: object,
theme: str,
frame: LocalFrame,
clip: Polygon,
) -> list[SurfaceFeature]:
features = []
for batch in reader:
for row in batch.to_pylist():
geometry = from_wkb(row["geometry"])
features.extend(self._normalize(row, theme, geometry, frame, clip))
if len(features) >= self.settings.max_surface_features_per_theme:
return features[: self.settings.max_surface_features_per_theme]
return features
@classmethod
def _normalize(
cls,
row: dict,
theme: str,
geometry,
frame: LocalFrame,
clip: Polygon,
) -> list[SurfaceFeature]:
local = transform(frame.transformer.transform, geometry)
if theme == "segment":
if not isinstance(local, LineString):
return []
return cls._normalize_road(row, local, clip)
decision = _surface_decision(theme, row)
if decision is None or not isinstance(local, (Polygon, MultiPolygon)):
return []
clipped = local.intersection(clip)
if clipped.is_empty:
return []
return [cls._feature(row, theme, clipped, decision)]
@classmethod
def _normalize_road(
cls,
row: dict,
geometry: LineString,
clip: Polygon,
) -> list[SurfaceFeature]:
if str(row.get("subtype") or "").lower() != "road":
return []
surface_rules = _linear_rules(row.get("road_surface"), "value", str)
width_rules = _linear_rules(row.get("width_rules"), "value", float)
boundaries = sorted(
{0.0, 1.0}
| {rule.start for rule in surface_rules + width_rules}
| {rule.stop for rule in surface_rules + width_rules}
)
features = []
for section_index, (start, stop) in enumerate(
zip(boundaries, boundaries[1:], strict=False)
):
midpoint = (start + stop) / 2.0
road_surface = _rule_value(surface_rules, midpoint)
decision = _road_decision(row, road_surface)
width = _rule_value(width_rules, midpoint)
width_m = float(width) if width is not None else _default_road_width(row)
width_source = "overture_width_rule" if width is not None else "class_default"
section = substring(geometry, start, stop, normalized=True)
polygon = section.buffer(width_m / 2.0, cap_style="flat", join_style="mitre")
clipped = polygon.intersection(clip)
if clipped.is_empty:
continue
features.append(
cls._feature(
row,
"segment",
clipped,
decision,
section_index=section_index,
source_surface=str(road_surface) if road_surface is not None else None,
width_m=width_m,
width_source=width_source,
)
)
return features
@staticmethod
def _feature(
row: dict,
theme: str,
geometry: Polygon | MultiPolygon,
decision: SurfaceDecision,
*,
section_index: int | None = None,
source_surface: str | None = None,
width_m: float | None = None,
width_source: str | None = None,
) -> SurfaceFeature:
sources = tuple(
FeatureSource(
property=_text(source.get("property")),
dataset=_text(source.get("dataset")),
record_id=_text(source.get("record_id")),
update_time=_text(source.get("update_time")),
license=_text(source.get("license")),
)
for source in row.get("sources") or []
)
record_id = sources[0].record_id if sources else None
feature_id = f"overture:{theme}:{row.get('id', record_id or len(sources))}"
if section_index is not None:
feature_id = f"{feature_id}:{section_index}"
return SurfaceFeature(
feature_id=feature_id,
surface_class=decision.surface_class,
geometry=geometry,
source_theme=theme,
source_subtype=_text(row.get("subtype")),
source_class=_text(row.get("class")),
source_surface=source_surface or _text(row.get("surface")),
sources=sources,
classification_basis=decision.basis,
precedence=decision.precedence,
width_m=width_m,
width_source=width_source,
)
def _text(value) -> str | None:
return str(value) if value not in (None, "") else None
def _surface_decision(theme: str, row: dict) -> SurfaceDecision | None:
subtype = str(row.get("subtype") or "").lower()
feature_class = str(row.get("class") or "").lower()
surface = str(row.get("surface") or "").lower()
if theme == "water":
return SurfaceDecision(SurfaceClass.WATER, "water_geometry", 50)
surface_class = _surface_material_class(surface)
if surface_class is not None:
return SurfaceDecision(surface_class, "land_use_surface", 30)
if theme == "land_cover":
if subtype in VEGETATION_LAND_COVER:
return SurfaceDecision(SurfaceClass.VEGETATION, "land_cover_subtype", 20)
if subtype in GROUND_LAND_COVER:
return SurfaceDecision(SurfaceClass.GROUND, "land_cover_subtype", 20)
if theme == "land_use" and (
subtype in VEGETATION_LAND_USE_SUBTYPES or feature_class in VEGETATION_LAND_USE_CLASSES
):
return SurfaceDecision(SurfaceClass.VEGETATION, "land_use_semantic", 10)
return None
def _surface_material_class(surface: str) -> SurfaceClass | None:
if surface in PAVED_SURFACES:
return SurfaceClass.PAVED
if surface in VEGETATION_SURFACES:
return SurfaceClass.VEGETATION
if surface in GROUND_SURFACES:
return SurfaceClass.GROUND
return None
def _road_decision(row: dict, road_surface: str | float | None) -> SurfaceDecision:
surface = str(road_surface or "").lower()
if surface in GROUND_ROAD_SURFACES:
return SurfaceDecision(SurfaceClass.GROUND, "road_surface", 40)
if surface in PAVED_ROAD_SURFACES:
return SurfaceDecision(SurfaceClass.PAVED, "road_surface", 40)
road_class = str(row.get("class") or "").lower()
inferred = SurfaceClass.GROUND if road_class in GROUND_ROAD_CLASSES else SurfaceClass.PAVED
return SurfaceDecision(inferred, "road_class_default", 35)
def _linear_rules(items, value_key: str, value_type: type) -> list[LinearRule]:
rules = []
for item in items or []:
value = item.get(value_key)
if value is None:
continue
between = item.get("between")
start, stop = (0.0, 1.0) if between is None else map(float, between)
rules.append(LinearRule(start, stop, value_type(value)))
return rules
def _rule_value(rules: list[LinearRule], position: float) -> str | float | None:
matches = [rule for rule in rules if rule.start <= position <= rule.stop]
if not matches:
return None
return min(matches, key=lambda rule: rule.stop - rule.start).value
def _default_road_width(row: dict) -> float:
road_class = str(row.get("class") or "").lower()
return ROAD_WIDTH_DEFAULTS_M.get(road_class, ROAD_WIDTH_DEFAULTS_M["other"])
[docs]
def classify_surface_cells(
terrain: TerrainGrid,
features: list[SurfaceFeature] | tuple[SurfaceFeature, ...],
) -> np.ndarray:
east = (terrain.east_m[:-1] + terrain.east_m[1:]) / 2.0
north = (terrain.north_m[:-1] + terrain.north_m[1:]) / 2.0
cells = np.full(terrain.cell_shape, SurfaceClass.GROUND, dtype=np.uint8)
for feature in sorted(
features,
key=lambda item: (item.precedence, int(item.surface_class), item.feature_id),
):
minimum_east, minimum_north, maximum_east, maximum_north = feature.geometry.bounds
column_start = int(np.searchsorted(east, minimum_east, side="left"))
column_stop = int(np.searchsorted(east, maximum_east, side="right"))
row_start = int(np.searchsorted(north, minimum_north, side="left"))
row_stop = int(np.searchsorted(north, maximum_north, side="right"))
if column_start >= column_stop or row_start >= row_stop:
continue
x, y = np.meshgrid(
east[column_start:column_stop],
north[row_start:row_stop],
)
mask = contains_xy(feature.geometry, x, y)
window = cells[row_start:row_stop, column_start:column_stop]
window[mask] = int(feature.surface_class)
return cells
def _connected_water_components(mask: np.ndarray) -> list[tuple[np.ndarray, np.ndarray]]:
visited = np.zeros(mask.shape, dtype=bool)
components = []
rows, columns = mask.shape
for start_row, start_column in zip(*np.nonzero(mask), strict=True):
if visited[start_row, start_column]:
continue
queue = deque([(int(start_row), int(start_column))])
visited[start_row, start_column] = True
component_rows = []
component_columns = []
while queue:
row, column = queue.popleft()
component_rows.append(row)
component_columns.append(column)
for row_offset in (-1, 0, 1):
for column_offset in (-1, 0, 1):
neighbor_row = row + row_offset
neighbor_column = column + column_offset
if (
0 <= neighbor_row < rows
and 0 <= neighbor_column < columns
and mask[neighbor_row, neighbor_column]
and not visited[neighbor_row, neighbor_column]
):
visited[neighbor_row, neighbor_column] = True
queue.append((neighbor_row, neighbor_column))
components.append(
(
np.asarray(component_rows, dtype=np.int64),
np.asarray(component_columns, dtype=np.int64),
)
)
return components
[docs]
def apply_water_surface_model(
terrain: TerrainGrid,
features: tuple[SurfaceFeature, ...] | list[SurfaceFeature],
surface_cells: np.ndarray,
) -> TerrainGrid:
"""Replace water-tagged terrain with one horizontal RF interface per water body."""
water_cells = surface_cells == int(SurfaceClass.WATER)
if not np.any(water_cells):
return terrain
marine = [
feature.geometry
for feature in features
if feature.source_theme == "water"
and (feature.source_subtype == "ocean" or feature.source_class in MARINE_WATER_CLASSES)
]
east = (terrain.east_m[:-1] + terrain.east_m[1:]) / 2.0
north = (terrain.north_m[:-1] + terrain.north_m[1:]) / 2.0
cell_east, cell_north = np.meshgrid(east, north)
marine_cells = (
intersects_xy(union_all(marine), cell_east, cell_north) & water_cells
if marine
else np.zeros(water_cells.shape, dtype=bool)
)
model = terrain.model_elevation_m.copy()
source = terrain.source_elevation_m
vertex_columns = terrain.shape[1]
component_records = []
for index, (rows, columns) in enumerate(_connected_water_components(water_cells)):
vertex_indices = np.unique(
np.concatenate(
[
rows * vertex_columns + columns,
rows * vertex_columns + columns + 1,
(rows + 1) * vertex_columns + columns,
(rows + 1) * vertex_columns + columns + 1,
]
)
)
is_marine = bool(np.any(marine_cells[rows, columns]))
level = 0.0 if is_marine else float(np.median(source.reshape(-1)[vertex_indices]))
model.reshape(-1)[vertex_indices] = level
component_records.append(
{
"index": index,
"cell_count": int(rows.size),
"vertex_count": int(vertex_indices.size),
"modeled_elevation_m": level,
"level_basis": (
"source_vertical_datum_zero" if is_marine else "source_dem_vertex_median"
),
}
)
water_vertices = np.zeros(terrain.shape, dtype=bool)
water_vertices[:-1, :-1] |= water_cells
water_vertices[:-1, 1:] |= water_cells
water_vertices[1:, :-1] |= water_cells
water_vertices[1:, 1:] |= water_cells
return terrain.with_model_elevation(
model,
{
"operation": "water_surface_hydroflattening",
"selection": (
"8-connected components of the final Overture water-class grid; all four "
"vertices of each selected cell"
),
"vertical_datum": terrain.source.vertical_datum,
"interpretation": (
"marine components use vertical-datum zero; other water components use the "
"median of their source DEM vertices; tides and local tidal-datum offsets are "
"not modeled"
),
"component_connectivity": 8,
"component_count": len(component_records),
"marine_component_count": sum(
item["level_basis"] == "source_vertical_datum_zero"
for item in component_records
),
"dem_estimated_component_count": sum(
item["level_basis"] == "source_dem_vertex_median"
for item in component_records
),
"water_cell_count": int(np.count_nonzero(water_cells)),
"marine_cell_count": int(np.count_nonzero(marine_cells)),
"vertex_count": int(np.count_nonzero(water_vertices)),
"source_minimum_m": float(np.min(source[water_vertices])),
"source_maximum_m": float(np.max(source[water_vertices])),
"components": component_records,
},
)
