192 lines
6.3 KiB
Python
192 lines
6.3 KiB
Python
import struct
|
|
import asyncio
|
|
import aiofiles
|
|
import httpx
|
|
from pathlib import Path
|
|
from typing import List, Optional, Tuple
|
|
import numpy as np
|
|
|
|
|
|
class TerrainService:
|
|
"""
|
|
SRTM elevation data service
|
|
- Downloads and caches .hgt tiles
|
|
- Provides elevation lookups
|
|
- Generates elevation profiles
|
|
"""
|
|
|
|
# SRTM tile dimensions (1 arc-second = 3601x3601, 3 arc-second = 1201x1201)
|
|
TILE_SIZE = 3601 # 1 arc-second (30m resolution)
|
|
|
|
# Mirror URLs for SRTM data (USGS requires login, use mirrors)
|
|
SRTM_MIRRORS = [
|
|
"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, cache_dir: str = "/opt/rfcp/backend/data/srtm"):
|
|
self.cache_dir = Path(cache_dir)
|
|
self.cache_dir.mkdir(exist_ok=True, parents=True)
|
|
self._tile_cache: dict[str, np.ndarray] = {} # In-memory cache
|
|
self._max_cached_tiles = 10 # Limit memory usage
|
|
|
|
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.cache_dir / f"{tile_name}.hgt"
|
|
|
|
async def download_tile(self, tile_name: str) -> bool:
|
|
"""Download SRTM tile from mirror"""
|
|
import gzip
|
|
|
|
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 mirror in self.SRTM_MIRRORS:
|
|
url = mirror.format(lat_dir=lat_dir, tile_name=tile_name)
|
|
try:
|
|
response = await client.get(url)
|
|
if response.status_code == 200:
|
|
# Decompress gzip
|
|
decompressed = gzip.decompress(response.content)
|
|
|
|
async with aiofiles.open(tile_path, 'wb') as f:
|
|
await f.write(decompressed)
|
|
|
|
print(f"Downloaded {tile_name} from {mirror}")
|
|
return True
|
|
except Exception as e:
|
|
print(f"Failed to download from {mirror}: {e}")
|
|
continue
|
|
|
|
print(f"Failed to download tile {tile_name}")
|
|
return False
|
|
|
|
async def load_tile(self, tile_name: str) -> Optional[np.ndarray]:
|
|
"""Load tile into memory (with caching)"""
|
|
# Check memory cache
|
|
if tile_name in self._tile_cache:
|
|
return self._tile_cache[tile_name]
|
|
|
|
tile_path = self.get_tile_path(tile_name)
|
|
|
|
# Download if missing
|
|
if not tile_path.exists():
|
|
success = await self.download_tile(tile_name)
|
|
if not success:
|
|
return None
|
|
|
|
# Read HGT file (big-endian signed 16-bit integers)
|
|
try:
|
|
async with aiofiles.open(tile_path, 'rb') as f:
|
|
data = await f.read()
|
|
|
|
# Parse as numpy array
|
|
arr = np.frombuffer(data, dtype='>i2').reshape(self.TILE_SIZE, self.TILE_SIZE)
|
|
|
|
# Manage cache size
|
|
if len(self._tile_cache) >= self._max_cached_tiles:
|
|
# Remove oldest entry
|
|
oldest = next(iter(self._tile_cache))
|
|
del self._tile_cache[oldest]
|
|
|
|
self._tile_cache[tile_name] = arr
|
|
return arr
|
|
|
|
except Exception as e:
|
|
print(f"Error loading tile {tile_name}: {e}")
|
|
return None
|
|
|
|
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 # No data, assume sea level
|
|
|
|
# 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, row 3600 = south edge
|
|
row = int((1 - lat_frac) * (self.TILE_SIZE - 1))
|
|
col = int(lon_frac * (self.TILE_SIZE - 1))
|
|
|
|
# Clamp to valid range
|
|
row = max(0, min(row, self.TILE_SIZE - 1))
|
|
col = max(0, min(col, self.TILE_SIZE - 1))
|
|
|
|
elevation = tile[row, col]
|
|
|
|
# -32768 = void/no data
|
|
if elevation == -32768:
|
|
return 0.0
|
|
|
|
return 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
|
|
|
|
Returns list of {lat, lon, elevation, distance} dicts
|
|
"""
|
|
lats = np.linspace(lat1, lat2, num_points)
|
|
lons = np.linspace(lon1, lon2, num_points)
|
|
|
|
# Calculate cumulative distances
|
|
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
|
|
|
|
@staticmethod
|
|
def haversine_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
|
|
"""Calculate distance between two points in meters"""
|
|
EARTH_RADIUS = 6371000 # meters
|
|
|
|
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()
|