from __future__ import annotations
import asyncio
import gzip
import io
import json
import math
from dataclasses import dataclass
import httpx
import numpy as np
from PIL import Image
from pyproj.enums import TransformDirection
from openworld_radio_twin.environment import TerrainGrid, TerrainSource
from openworld_radio_twin.geodesy import LocalFrame
USGS_3DEP_URL = "https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer"
MAPZEN_SKADI_URL = (
"https://s3.amazonaws.com/elevation-tiles-prod/skadi/{latitude_band}/{tile}.hgt.gz"
)
SKADI_VOID = -32768
MIN_TERRESTRIAL_ELEVATION_M = -500.0
[docs]
class TerrainProviderError(RuntimeError):
pass
[docs]
@dataclass(frozen=True)
class TerrainQuery:
longitude: float
latitude: float
radius_m: float
spacing_m: float
@property
def axis(self) -> np.ndarray:
count = max(2, math.ceil(2 * self.radius_m / self.spacing_m) + 1)
return np.linspace(-self.radius_m, self.radius_m, count, dtype=np.float64)
def _local_coordinates(
query: TerrainQuery,
) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
axis = query.axis
east, north = np.meshgrid(axis, axis)
frame = LocalFrame.at(query.longitude, query.latitude)
longitude, latitude, _ = frame.transformer.transform(
east,
north,
np.zeros_like(east),
direction=TransformDirection.INVERSE,
)
return axis, east, np.asarray(longitude), np.asarray(latitude)
FLAT_TERRAIN_PROVIDER = "flat plane"
[docs]
def flat_terrain(query: TerrainQuery) -> TerrainGrid:
"""Return a gridded flat ground at local z=0 so surface classes can be assigned per cell.
No elevation source is queried. The grid uses the same axes as a DEM-backed
terrain at the requested spacing, so semantic surface classification and the
per-class terrain meshes work unchanged on a planar ground.
"""
source = TerrainSource(
provider=FLAT_TERRAIN_PROVIDER,
dataset="planar ground at local z=0",
vertical_datum="local ground plane",
requested_spacing_m=float(query.spacing_m),
source_resolution_m=None,
interpolation="none",
source_url="",
license="not applicable",
attributes={"terrain_mode": "flat", "elevation_source_queried": False},
)
return _normalize(query, np.zeros((query.axis.size, query.axis.size), dtype=np.float32), source)
def _normalize(
query: TerrainQuery,
elevations: np.ndarray,
source: TerrainSource,
) -> TerrainGrid:
axis = query.axis
model = np.asarray(elevations, dtype=np.float32)
center = _sample_regular_grid(model, axis, 0.0, 0.0)
# Subtract before rounding to local float32: the interpolated origin can
# lie between representable absolute elevations, particularly at altitude.
local = (model.astype(np.float64) - center).astype(np.float32)
return TerrainGrid(
east_m=axis,
north_m=axis.copy(),
source_elevation_m=model,
model_elevation_m=model.copy(),
local_elevation_m=local,
source_origin_elevation_m=float(center),
origin_elevation_m=float(center),
source=source,
)
def _sample_regular_grid(
values: np.ndarray,
axis: np.ndarray,
east_m: float,
north_m: float,
) -> float:
column = int(np.clip(np.searchsorted(axis, east_m) - 1, 0, axis.size - 2))
row = int(np.clip(np.searchsorted(axis, north_m) - 1, 0, axis.size - 2))
tx = (east_m - axis[column]) / (axis[column + 1] - axis[column])
ty = (north_m - axis[row]) / (axis[row + 1] - axis[row])
south = values[row, column] * (1 - tx) + values[row, column + 1] * tx
north = values[row + 1, column] * (1 - tx) + values[row + 1, column + 1] * tx
return float(south * (1 - ty) + north * ty)
def _bilinear_rectilinear(
values_north_up: np.ndarray,
west: float,
south: float,
east: float,
north: float,
longitude: np.ndarray,
latitude: np.ndarray,
) -> np.ndarray:
rows, columns = values_north_up.shape
x = (longitude - west) / (east - west) * columns - 0.5
y = (north - latitude) / (north - south) * rows - 0.5
x = np.clip(x, 0, columns - 1)
y = np.clip(y, 0, rows - 1)
x0 = np.floor(x).astype(np.int64)
y0 = np.floor(y).astype(np.int64)
x1 = np.minimum(x0 + 1, columns - 1)
y1 = np.minimum(y0 + 1, rows - 1)
tx = x - x0
ty = y - y0
north_values = values_north_up[y0, x0] * (1 - tx) + values_north_up[y0, x1] * tx
south_values = values_north_up[y1, x0] * (1 - tx) + values_north_up[y1, x1] * tx
return north_values * (1 - ty) + south_values * ty
[docs]
class USGS3DEPTerrainProvider:
def __init__(self, client: httpx.AsyncClient) -> None:
self.client = client
[docs]
@staticmethod
def covers(longitude: float, latitude: float) -> bool:
mainland_and_alaska = -180 <= longitude <= -60 and 18 <= latitude <= 72
pacific_territories = 130 <= longitude <= 180 and 0 <= latitude <= 25
return mainland_and_alaska or pacific_territories
[docs]
async def terrain(self, query: TerrainQuery) -> TerrainGrid:
axis, _, longitudes, latitudes = _local_coordinates(query)
margin = max(query.spacing_m, 2.0) / 111_000.0
west = float(np.min(longitudes) - margin)
east = float(np.max(longitudes) + margin)
south = float(np.min(latitudes) - margin)
north = float(np.max(latitudes) + margin)
size = axis.size + 2
params = {
"bbox": f"{west},{south},{east},{north}",
"bboxSR": 4326,
"imageSR": 4326,
"size": f"{size},{size}",
"format": "tiff",
"pixelType": "F32",
"interpolation": "RSP_BilinearInterpolation",
"f": "image",
}
# The DEM export and the catalog lookup are independent requests.
response, identified = await asyncio.gather(
self.client.get(f"{USGS_3DEP_URL}/exportImage", params=params),
self._identify(query.longitude, query.latitude),
)
response.raise_for_status()
try:
with Image.open(io.BytesIO(response.content)) as image:
source_values = np.asarray(image, dtype=np.float32)
except Exception as exc:
raise TerrainProviderError("USGS 3DEP did not return a readable F32 DEM") from exc
if source_values.shape != (size, size):
raise TerrainProviderError(f"USGS 3DEP returned unexpected shape {source_values.shape}")
if not np.all(np.isfinite(source_values)) or np.min(source_values) < -10_000:
raise TerrainProviderError(
"USGS 3DEP has no complete elevation coverage for this scene"
)
elevations = _bilinear_rectilinear(
source_values,
west,
south,
east,
north,
longitudes,
latitudes,
)
attributes = {
**identified,
"owrt_metadata_sample_longitude": query.longitude,
"owrt_metadata_sample_latitude": query.latitude,
}
source_resolution = _numeric_attribute(
attributes,
"Resolution",
"Resolution_X",
"Pixel Size",
"LowPS",
)
source = TerrainSource(
provider="USGS 3DEP",
dataset="3DEP Bare Earth DEM Dynamic",
vertical_datum=str(
_attribute(attributes, "VerticalDatum", "Vertical Datum") or "source metadata"
),
requested_spacing_m=float(2 * query.radius_m / (axis.size - 1)),
source_resolution_m=source_resolution,
interpolation="USGS bilinear export followed by WGS84-to-ENU bilinear sampling",
source_url=USGS_3DEP_URL,
license="United States public domain",
attributes=attributes,
)
return _normalize(query, elevations, source)
async def _identify(self, longitude: float, latitude: float) -> dict:
params = {
"geometry": json.dumps(
{
"x": longitude,
"y": latitude,
"spatialReference": {"wkid": 4326},
}
),
"geometryType": "esriGeometryPoint",
"returnGeometry": "false",
"returnCatalogItems": "true",
"f": "json",
}
try:
response = await self.client.get(f"{USGS_3DEP_URL}/identify", params=params)
response.raise_for_status()
except httpx.HTTPError as exc:
status = exc.response.status_code if isinstance(exc, httpx.HTTPStatusError) else None
return {
"owrt_catalog_metadata_status": "unavailable",
"owrt_catalog_metadata_http_status": status,
}
payload = response.json()
features = payload.get("catalogItems", {}).get("features", [])
candidates = [
dict(feature.get("attributes", {}))
for feature in features
if feature.get("attributes", {}).get("Category") == 1
]
return {
**min(candidates, key=_pixel_size, default={}),
"owrt_catalog_metadata_status": "available",
}
def _attribute(attributes: dict, *names: str):
lowered = {str(key).lower(): value for key, value in attributes.items()}
for name in names:
value = lowered.get(name.lower())
if value not in (None, ""):
return value
return None
def _numeric_attribute(attributes: dict, *names: str) -> float | None:
value = _attribute(attributes, *names)
if value is None:
return None
try:
return float(value)
except (TypeError, ValueError):
return None
def _pixel_size(attributes: dict) -> float:
value = _numeric_attribute(attributes, "LowPS", "Resolution_X", "MinPS")
return value if value is not None and value > 0 else math.inf
[docs]
class MapzenTerrainProvider:
"""Global raw HGT fallback with explicit void handling and source headers."""
def __init__(self, client: httpx.AsyncClient) -> None:
self.client = client
[docs]
async def terrain(self, query: TerrainQuery) -> TerrainGrid:
axis, _, longitudes, latitudes = _local_coordinates(query)
normalized_longitudes = (longitudes + 180.0) % 360.0 - 180.0
tile_keys = sorted(
{
(math.floor(float(latitude)), math.floor(float(longitude)))
for latitude, longitude in zip(
latitudes.reshape(-1), normalized_longitudes.reshape(-1), strict=True
)
}
)
fetched = await asyncio.gather(*(self._tile(*key) for key in tile_keys))
tiles = {key: values for key, values, _ in fetched}
imagery_sources = sorted(
{source for _, _, sources in fetched for source in sources.split(",") if source}
)
elevations = np.empty_like(latitudes, dtype=np.float64)
source_spacings = []
for latitude_floor, longitude_floor in tile_keys:
tile = tiles[(latitude_floor, longitude_floor)]
source_spacings.append(1.0 / (tile.shape[0] - 1))
mask = (np.floor(latitudes) == latitude_floor) & (
np.floor(normalized_longitudes) == longitude_floor
)
elevations[mask] = _sample_hgt(
tile,
normalized_longitudes[mask] - longitude_floor,
latitude_floor + 1.0 - latitudes[mask],
)
if np.min(elevations) < -12_000 or np.max(elevations) > 9_000:
raise TerrainProviderError(
"Mapzen returned elevations outside the physical Earth range"
)
arc_spacing = max(source_spacings)
source_resolution_m = max(
111_132.0 * arc_spacing,
111_320.0 * math.cos(math.radians(query.latitude)) * arc_spacing,
)
source = TerrainSource(
provider="Mapzen Terrain Tiles",
dataset="Skadi global raw elevation mosaic",
vertical_datum="EGM96 orthometric height",
requested_spacing_m=float(2 * query.radius_m / (axis.size - 1)),
source_resolution_m=source_resolution_m,
interpolation="big-endian int16 HGT decode followed by geographic bilinear sampling",
source_url="https://registry.opendata.aws/terrain-tiles/",
license="Mapzen/AWS Terrain Tiles attribution terms; source datasets vary",
attributes={
"format": "Skadi HGT.gz",
"tile_count": len(tile_keys),
"tiles": ",".join(_skadi_name(*key) for key in tile_keys),
"imagery_sources": ",".join(imagery_sources),
"void_policy": "fail if any -32768 sample participates in interpolation",
"source_range_policy": "fail outside -12000..9000 m",
"rf_surface_policy": (
"reject modeled values below -500 m after explicit marine hydroflattening"
),
"quality_note": (
"Multi-source mosaic; source accuracy and datum consistency vary by region"
),
},
)
return _normalize(query, elevations, source)
async def _tile(self, latitude_floor: int, longitude_floor: int):
tile_name = _skadi_name(latitude_floor, longitude_floor)
url = MAPZEN_SKADI_URL.format(latitude_band=tile_name[:3], tile=tile_name)
response = await self.client.get(url)
response.raise_for_status()
try:
raw = gzip.decompress(response.content)
side = math.isqrt(len(raw) // 2)
if side * side * 2 != len(raw) or side < 2:
raise ValueError("payload is not a square int16 grid")
values = np.frombuffer(raw, dtype=">i2").reshape(side, side)
except (gzip.BadGzipFile, ValueError) as exc:
raise TerrainProviderError(f"Mapzen {tile_name} is not a valid HGT tile") from exc
return (
(latitude_floor, longitude_floor),
values,
response.headers.get("x-imagery-sources", ""),
)
def _skadi_name(latitude_floor: int, longitude_floor: int) -> str:
latitude = f"N{latitude_floor:02d}" if latitude_floor >= 0 else f"S{-latitude_floor:02d}"
longitude = f"E{longitude_floor:03d}" if longitude_floor >= 0 else f"W{-longitude_floor:03d}"
return f"{latitude}{longitude}"
def _sample_hgt(values: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
scale = values.shape[0] - 1
x = np.clip(x * scale, 0, scale)
y = np.clip(y * scale, 0, scale)
x0 = np.floor(x).astype(np.int64)
y0 = np.floor(y).astype(np.int64)
x1 = np.minimum(x0 + 1, scale)
y1 = np.minimum(y0 + 1, scale)
samples = np.stack((values[y0, x0], values[y0, x1], values[y1, x0], values[y1, x1]))
if np.any(samples == SKADI_VOID):
raise TerrainProviderError("Mapzen Skadi has a void in the requested scene")
tx = x - x0
ty = y - y0
top = samples[0] * (1 - tx) + samples[1] * tx
bottom = samples[2] * (1 - tx) + samples[3] * tx
return top * (1 - ty) + bottom * ty
[docs]
class TerrainProvider:
def __init__(self, client: httpx.AsyncClient) -> None:
self.usgs = USGS3DEPTerrainProvider(client)
self.global_fallback = MapzenTerrainProvider(client)
[docs]
async def terrain(
self,
query: TerrainQuery,
allow_bathymetry: bool = False,
) -> tuple[TerrainGrid, list[str]]:
warnings = []
if self.usgs.covers(query.longitude, query.latitude):
try:
terrain = await self.usgs.terrain(query)
if (
terrain.source.attributes.get("owrt_catalog_metadata_status")
== "unavailable"
):
warnings.append(
"USGS catalog metadata was unavailable; DEM samples are retained and "
"source resolution is reported as unknown."
)
return terrain, warnings
except (httpx.HTTPError, TerrainProviderError) as exc:
warnings.append(f"USGS 3DEP unavailable for this scene; using global DEM: {exc}")
try:
terrain = await self.global_fallback.terrain(query)
if not allow_bathymetry:
validate_terrestrial_surface(terrain)
warnings.append(
"Global terrain uses a multi-source Mapzen mosaic; regional source accuracy, "
"datum consistency, and attribution must be reviewed in exported metadata."
)
return terrain, warnings
except (httpx.HTTPError, TerrainProviderError) as exc:
raise TerrainProviderError(f"Global terrain query failed: {exc}") from exc
[docs]
def validate_terrestrial_surface(terrain: TerrainGrid) -> None:
if np.min(terrain.model_elevation_m) < MIN_TERRESTRIAL_ELEVATION_M:
raise TerrainProviderError(
f"{terrain.source.provider} sampled bathymetry in the scene; enable semantic "
"surface data so marine cells can be modeled at mean sea level"
)
