293 lines
9.8 KiB
Python
293 lines
9.8 KiB
Python
import os
|
|
import struct
|
|
import asyncio
|
|
import gzip
|
|
import zipfile
|
|
import io
|
|
import numpy as np
|
|
import httpx
|
|
from pathlib import Path
|
|
from typing import List, Optional, Tuple
|
|
|
|
|
|
class TerrainService:
|
|
"""
|
|
SRTM elevation data service with local caching.
|
|
- Stores tiles in RFCP_DATA_PATH/terrain/
|
|
- In-memory LRU cache (max 20 tiles)
|
|
- Auto-downloads from S3 mirror
|
|
- Supports both SRTM1 (3601x3601) and SRTM3 (1201x1201)
|
|
"""
|
|
|
|
SRTM_SOURCES = [
|
|
"https://elevation-tiles-prod.s3.amazonaws.com/skadi/{lat_dir}/{tile_name}.hgt.gz",
|
|
"https://s3.amazonaws.com/elevation-tiles-prod/skadi/{lat_dir}/{tile_name}.hgt.gz",
|
|
]
|
|
|
|
def __init__(self):
|
|
self.data_path = Path(os.environ.get('RFCP_DATA_PATH', './data'))
|
|
self.terrain_path = self.data_path / 'terrain'
|
|
self.terrain_path.mkdir(parents=True, exist_ok=True)
|
|
|
|
# In-memory cache for loaded tiles
|
|
self._tile_cache: dict[str, np.ndarray] = {}
|
|
self._max_cache_tiles = 20 # ~500MB max
|
|
|
|
def get_tile_name(self, lat: float, lon: float) -> str:
|
|
"""Convert lat/lon to SRTM tile name (e.g., N48E035)"""
|
|
lat_int = int(lat) if lat >= 0 else int(lat) - 1
|
|
lon_int = int(lon) if lon >= 0 else int(lon) - 1
|
|
|
|
lat_letter = 'N' if lat_int >= 0 else 'S'
|
|
lon_letter = 'E' if lon_int >= 0 else 'W'
|
|
|
|
return f"{lat_letter}{abs(lat_int):02d}{lon_letter}{abs(lon_int):03d}"
|
|
|
|
def get_tile_path(self, tile_name: str) -> Path:
|
|
"""Get local path for tile"""
|
|
return self.terrain_path / f"{tile_name}.hgt"
|
|
|
|
async def download_tile(self, tile_name: str) -> bool:
|
|
"""Download SRTM tile if not cached locally"""
|
|
tile_path = self.get_tile_path(tile_name)
|
|
|
|
if tile_path.exists():
|
|
return True
|
|
|
|
lat_dir = tile_name[:3] # e.g., "N48"
|
|
|
|
async with httpx.AsyncClient(timeout=60.0) as client:
|
|
for source_url in self.SRTM_SOURCES:
|
|
url = source_url.format(lat_dir=lat_dir, tile_name=tile_name)
|
|
try:
|
|
response = await client.get(url)
|
|
|
|
if response.status_code == 200:
|
|
data = response.content
|
|
|
|
if url.endswith('.gz'):
|
|
data = gzip.decompress(data)
|
|
elif url.endswith('.zip'):
|
|
with zipfile.ZipFile(io.BytesIO(data)) as zf:
|
|
for name in zf.namelist():
|
|
if name.endswith('.hgt'):
|
|
data = zf.read(name)
|
|
break
|
|
|
|
tile_path.write_bytes(data)
|
|
print(f"[Terrain] Downloaded {tile_name} ({len(data)} bytes)")
|
|
return True
|
|
|
|
except Exception as e:
|
|
print(f"[Terrain] Failed from {url}: {e}")
|
|
continue
|
|
|
|
print(f"[Terrain] Could not download {tile_name}")
|
|
return False
|
|
|
|
def _load_tile(self, tile_name: str) -> Optional[np.ndarray]:
|
|
"""Load tile from disk into memory cache"""
|
|
# Check memory cache first
|
|
if tile_name in self._tile_cache:
|
|
return self._tile_cache[tile_name]
|
|
|
|
tile_path = self.get_tile_path(tile_name)
|
|
|
|
if not tile_path.exists():
|
|
return None
|
|
|
|
try:
|
|
data = tile_path.read_bytes()
|
|
|
|
# SRTM HGT format: big-endian signed 16-bit integers
|
|
if len(data) == 3601 * 3601 * 2:
|
|
size = 3601 # SRTM1 (30m)
|
|
elif len(data) == 1201 * 1201 * 2:
|
|
size = 1201 # SRTM3 (90m)
|
|
else:
|
|
print(f"[Terrain] Unknown tile size: {len(data)} bytes for {tile_name}")
|
|
return None
|
|
|
|
tile = np.frombuffer(data, dtype='>i2').reshape((size, size))
|
|
|
|
# Manage memory cache with LRU eviction
|
|
if len(self._tile_cache) >= self._max_cache_tiles:
|
|
oldest = next(iter(self._tile_cache))
|
|
del self._tile_cache[oldest]
|
|
|
|
self._tile_cache[tile_name] = tile
|
|
return tile
|
|
|
|
except Exception as e:
|
|
print(f"[Terrain] Failed to load {tile_name}: {e}")
|
|
return None
|
|
|
|
async def load_tile(self, tile_name: str) -> Optional[np.ndarray]:
|
|
"""Load tile into memory, downloading if needed"""
|
|
# Check memory cache
|
|
if tile_name in self._tile_cache:
|
|
return self._tile_cache[tile_name]
|
|
|
|
# Download if missing
|
|
if not self.get_tile_path(tile_name).exists():
|
|
success = await self.download_tile(tile_name)
|
|
if not success:
|
|
return None
|
|
|
|
return self._load_tile(tile_name)
|
|
|
|
async def get_elevation(self, lat: float, lon: float) -> float:
|
|
"""Get elevation at specific coordinate (meters above sea level)"""
|
|
tile_name = self.get_tile_name(lat, lon)
|
|
tile = await self.load_tile(tile_name)
|
|
|
|
if tile is None:
|
|
return 0.0
|
|
|
|
size = tile.shape[0]
|
|
|
|
# Calculate position within tile
|
|
lat_int = int(lat) if lat >= 0 else int(lat) - 1
|
|
lon_int = int(lon) if lon >= 0 else int(lon) - 1
|
|
|
|
lat_frac = lat - lat_int
|
|
lon_frac = lon - lon_int
|
|
|
|
# Row 0 = north edge, last row = south edge
|
|
row = int((1 - lat_frac) * (size - 1))
|
|
col = int(lon_frac * (size - 1))
|
|
|
|
row = max(0, min(row, size - 1))
|
|
col = max(0, min(col, size - 1))
|
|
|
|
elevation = tile[row, col]
|
|
|
|
# -32768 = void/no data
|
|
if elevation == -32768:
|
|
return 0.0
|
|
|
|
return float(elevation)
|
|
|
|
def get_elevation_sync(self, lat: float, lon: float) -> float:
|
|
"""Sync elevation lookup from memory cache. Returns 0.0 if tile not loaded."""
|
|
tile_name = self.get_tile_name(lat, lon)
|
|
tile = self._tile_cache.get(tile_name)
|
|
if tile is None:
|
|
return 0.0
|
|
|
|
size = tile.shape[0]
|
|
lat_int = int(lat) if lat >= 0 else int(lat) - 1
|
|
lon_int = int(lon) if lon >= 0 else int(lon) - 1
|
|
|
|
row = int((1 - (lat - lat_int)) * (size - 1))
|
|
col = int((lon - lon_int) * (size - 1))
|
|
row = max(0, min(row, size - 1))
|
|
col = max(0, min(col, size - 1))
|
|
|
|
elevation = tile[row, col]
|
|
return 0.0 if elevation == -32768 else float(elevation)
|
|
|
|
async def get_elevation_profile(
|
|
self,
|
|
lat1: float, lon1: float,
|
|
lat2: float, lon2: float,
|
|
num_points: int = 100
|
|
) -> List[dict]:
|
|
"""Get elevation profile between two points"""
|
|
lats = np.linspace(lat1, lat2, num_points)
|
|
lons = np.linspace(lon1, lon2, num_points)
|
|
|
|
total_distance = self.haversine_distance(lat1, lon1, lat2, lon2)
|
|
distances = np.linspace(0, total_distance, num_points)
|
|
|
|
profile = []
|
|
for i, (lat, lon, dist) in enumerate(zip(lats, lons, distances)):
|
|
elev = await self.get_elevation(lat, lon)
|
|
profile.append({
|
|
"lat": float(lat),
|
|
"lon": float(lon),
|
|
"elevation": elev,
|
|
"distance": float(dist)
|
|
})
|
|
|
|
return profile
|
|
|
|
def get_elevation_profile_sync(
|
|
self,
|
|
lat1: float, lon1: float,
|
|
lat2: float, lon2: float,
|
|
num_points: int = 50
|
|
) -> List[dict]:
|
|
"""Sync elevation profile - tiles must be pre-loaded into memory cache."""
|
|
lats = np.linspace(lat1, lat2, num_points)
|
|
lons = np.linspace(lon1, lon2, num_points)
|
|
|
|
total_distance = self.haversine_distance(lat1, lon1, lat2, lon2)
|
|
distances = np.linspace(0, total_distance, num_points)
|
|
|
|
profile = []
|
|
for i in range(num_points):
|
|
profile.append({
|
|
"lat": float(lats[i]),
|
|
"lon": float(lons[i]),
|
|
"elevation": self.get_elevation_sync(float(lats[i]), float(lons[i])),
|
|
"distance": float(distances[i])
|
|
})
|
|
|
|
return profile
|
|
|
|
async def ensure_tiles_for_bbox(
|
|
self,
|
|
min_lat: float, min_lon: float,
|
|
max_lat: float, max_lon: float
|
|
) -> list[str]:
|
|
"""Pre-download all tiles needed for a bounding box"""
|
|
tiles_needed = []
|
|
|
|
for lat in range(int(min_lat), int(max_lat) + 1):
|
|
for lon in range(int(min_lon), int(max_lon) + 1):
|
|
tile_name = self.get_tile_name(lat, lon)
|
|
tiles_needed.append(tile_name)
|
|
|
|
# Download in parallel (batches of 5 to avoid overload)
|
|
downloaded = []
|
|
batch_size = 5
|
|
for i in range(0, len(tiles_needed), batch_size):
|
|
batch = tiles_needed[i:i + batch_size]
|
|
results = await asyncio.gather(*[
|
|
self.download_tile(tile) for tile in batch
|
|
])
|
|
for tile, ok in zip(batch, results):
|
|
if ok:
|
|
downloaded.append(tile)
|
|
|
|
return downloaded
|
|
|
|
def get_cached_tiles(self) -> list[str]:
|
|
"""List all locally cached tile names"""
|
|
return [f.stem for f in self.terrain_path.glob("*.hgt")]
|
|
|
|
def get_cache_size_mb(self) -> float:
|
|
"""Get total terrain cache size in MB"""
|
|
total = sum(f.stat().st_size for f in self.terrain_path.glob("*.hgt"))
|
|
return total / (1024 * 1024)
|
|
|
|
@staticmethod
|
|
def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
|
|
"""Calculate distance between two points in meters"""
|
|
EARTH_RADIUS = 6371000
|
|
|
|
lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])
|
|
|
|
dlat = lat2 - lat1
|
|
dlon = lon2 - lon1
|
|
|
|
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
|
|
c = 2 * np.arcsin(np.sqrt(a))
|
|
|
|
return EARTH_RADIUS * c
|
|
|
|
|
|
# Singleton instance
|
|
terrain_service = TerrainService()
|