root
384 lines
11 KiB
Python
384 lines
11 KiB
Python
"""
|
|
Geo service for building search zones and querying OpenStreetMap data via Overpass API.
|
|
"""
|
|
import math
|
|
import json
|
|
import hashlib
|
|
from datetime import datetime, timedelta
|
|
from pathlib import Path
|
|
from typing import List, Dict, Optional, Tuple
|
|
import httpx
|
|
from pydantic import BaseModel
|
|
|
|
|
|
class Zone(BaseModel):
|
|
"""Search zone with geographic features."""
|
|
direction: str # N, NE, E, SE, S, SW, W, NW
|
|
distance_km: float
|
|
forest_pct: float
|
|
road_density: float # km of roads per km²
|
|
water_distance_km: Optional[float]
|
|
settlement_distance_km: Optional[float]
|
|
|
|
|
|
# Cache configuration
|
|
CACHE_DIR = Path("/tmp/overpass_cache")
|
|
CACHE_TTL_HOURS = 24
|
|
OVERPASS_URL = "https://overpass-api.de/api/interpreter"
|
|
|
|
# Direction mappings
|
|
SEARCH_DISTANCES = [500, 1000, 2000, 5000]
|
|
|
|
DIRECTIONS = ["N", "NE", "E", "SE", "S", "SW", "W", "NW"]
|
|
DIRECTION_ANGLES = {
|
|
"N": 0,
|
|
"NE": 45,
|
|
"E": 90,
|
|
"SE": 135,
|
|
"S": 180,
|
|
"SW": 225,
|
|
"W": 270,
|
|
"NW": 315
|
|
}
|
|
|
|
# Search distances computed dynamically from max_distance_km
|
|
def _build_search_distances(max_distance_km: float) -> list:
|
|
max_m = int(max_distance_km * 1000)
|
|
raw = [int(max_m * 0.25), int(max_m * 0.50), int(max_m * 0.75), max_m]
|
|
clamped = [max(200, min(5000, d)) for d in raw]
|
|
return sorted(set(clamped))
|
|
|
|
|
|
def haversine(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
|
|
"""
|
|
Calculate distance between two points on Earth using Haversine formula.
|
|
|
|
Args:
|
|
lat1, lon1: First point coordinates
|
|
lat2, lon2: Second point coordinates
|
|
|
|
Returns:
|
|
Distance in kilometers
|
|
"""
|
|
R = 6371 # Earth radius in km
|
|
|
|
lat1_rad = math.radians(lat1)
|
|
lat2_rad = math.radians(lat2)
|
|
dlat = math.radians(lat2 - lat1)
|
|
dlon = math.radians(lon2 - lon1)
|
|
|
|
a = (math.sin(dlat / 2) ** 2 +
|
|
math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2) ** 2)
|
|
c = 2 * math.asin(math.sqrt(a))
|
|
|
|
return R * c
|
|
|
|
|
|
def get_sector_bounds(lat: float, lon: float, direction: str, radius_m: int) -> Tuple[float, float, float, float]:
|
|
"""
|
|
Calculate bounding box for a sector.
|
|
|
|
Args:
|
|
lat, lon: Center point
|
|
direction: Sector direction (N, NE, E, etc.)
|
|
radius_m: Radius in meters
|
|
|
|
Returns:
|
|
(min_lat, min_lon, max_lat, max_lon)
|
|
"""
|
|
# Convert radius to degrees (approximate)
|
|
radius_deg = radius_m / 111320 # 1 degree ≈ 111.32 km at equator
|
|
|
|
angle = DIRECTION_ANGLES[direction]
|
|
angle_rad = math.radians(angle)
|
|
|
|
# Calculate sector boundaries (45° sectors)
|
|
angle_start = angle - 22.5
|
|
angle_end = angle + 22.5
|
|
|
|
# Simple bounding box (can be optimized for actual sector shape)
|
|
lat_offset = radius_deg * math.cos(angle_rad)
|
|
lon_offset = radius_deg * math.sin(angle_rad) / math.cos(math.radians(lat))
|
|
|
|
min_lat = min(lat, lat + lat_offset) - radius_deg * 0.5
|
|
max_lat = max(lat, lat + lat_offset) + radius_deg * 0.5
|
|
min_lon = min(lon, lon + lon_offset) - radius_deg * 0.5
|
|
max_lon = max(lon, lon + lon_offset) + radius_deg * 0.5
|
|
|
|
return (min_lat, min_lon, max_lat, max_lon)
|
|
|
|
|
|
def get_cache_key(query: str) -> str:
|
|
"""Generate cache key from query."""
|
|
return hashlib.md5(query.encode()).hexdigest()
|
|
|
|
|
|
def get_cached_result(cache_key: str) -> Optional[Dict]:
|
|
"""Get cached Overpass API result if not expired."""
|
|
CACHE_DIR.mkdir(exist_ok=True)
|
|
cache_file = CACHE_DIR / f"{cache_key}.json"
|
|
|
|
if not cache_file.exists():
|
|
return None
|
|
|
|
try:
|
|
with open(cache_file, 'r') as f:
|
|
cached = json.load(f)
|
|
|
|
cached_time = datetime.fromisoformat(cached['timestamp'])
|
|
if datetime.now() - cached_time > timedelta(hours=CACHE_TTL_HOURS):
|
|
cache_file.unlink()
|
|
return None
|
|
|
|
return cached['data']
|
|
except Exception:
|
|
return None
|
|
|
|
|
|
def save_to_cache(cache_key: str, data: Dict):
|
|
"""Save Overpass API result to cache."""
|
|
CACHE_DIR.mkdir(exist_ok=True)
|
|
cache_file = CACHE_DIR / f"{cache_key}.json"
|
|
|
|
try:
|
|
with open(cache_file, 'w') as f:
|
|
json.dump({
|
|
'timestamp': datetime.now().isoformat(),
|
|
'data': data
|
|
}, f)
|
|
except Exception:
|
|
pass
|
|
|
|
|
|
async def query_overpass(query: str) -> Dict:
|
|
"""
|
|
Query Overpass API with caching.
|
|
|
|
Args:
|
|
query: Overpass QL query
|
|
|
|
Returns:
|
|
API response as dict
|
|
"""
|
|
cache_key = get_cache_key(query)
|
|
|
|
# Check cache
|
|
cached = get_cached_result(cache_key)
|
|
if cached is not None:
|
|
return cached
|
|
|
|
# Query API
|
|
async with httpx.AsyncClient(timeout=30.0) as client:
|
|
try:
|
|
response = await client.post(
|
|
OVERPASS_URL,
|
|
data={'data': query},
|
|
headers={'Content-Type': 'application/x-www-form-urlencoded'}
|
|
)
|
|
response.raise_for_status()
|
|
data = response.json()
|
|
|
|
# Save to cache
|
|
save_to_cache(cache_key, data)
|
|
|
|
return data
|
|
except Exception as e:
|
|
# Return empty result on error
|
|
return {'elements': []}
|
|
|
|
|
|
def calculate_road_length(elements: List[Dict]) -> float:
|
|
"""
|
|
Calculate total road length from Overpass way elements.
|
|
|
|
Args:
|
|
elements: List of way elements from Overpass
|
|
|
|
Returns:
|
|
Total length in kilometers
|
|
"""
|
|
total_length = 0.0
|
|
|
|
for element in elements:
|
|
if element.get('type') != 'way':
|
|
continue
|
|
|
|
nodes = element.get('geometry', [])
|
|
if len(nodes) < 2:
|
|
continue
|
|
|
|
# Calculate length by summing distances between consecutive nodes
|
|
for i in range(len(nodes) - 1):
|
|
lat1, lon1 = nodes[i]['lat'], nodes[i]['lon']
|
|
lat2, lon2 = nodes[i + 1]['lat'], nodes[i + 1]['lon']
|
|
total_length += haversine(lat1, lon1, lat2, lon2)
|
|
|
|
return total_length
|
|
|
|
|
|
def find_nearest_distance(lat: float, lon: float, elements: List[Dict]) -> Optional[float]:
|
|
"""
|
|
Find distance to nearest element.
|
|
|
|
Args:
|
|
lat, lon: Reference point
|
|
elements: List of node elements from Overpass
|
|
|
|
Returns:
|
|
Distance in kilometers, or None if no elements
|
|
"""
|
|
if not elements:
|
|
return None
|
|
|
|
min_distance = float('inf')
|
|
|
|
for element in elements:
|
|
if element.get('type') != 'node':
|
|
continue
|
|
|
|
elem_lat = element.get('lat')
|
|
elem_lon = element.get('lon')
|
|
|
|
if elem_lat is None or elem_lon is None:
|
|
continue
|
|
|
|
distance = haversine(lat, lon, elem_lat, elem_lon)
|
|
min_distance = min(min_distance, distance)
|
|
|
|
return min_distance if min_distance != float('inf') else None
|
|
|
|
|
|
def calculate_forest_coverage(elements: List[Dict], radius_m: int) -> float:
|
|
"""
|
|
Estimate forest coverage percentage.
|
|
|
|
Args:
|
|
elements: List of way elements from Overpass
|
|
radius_m: Search radius in meters
|
|
|
|
Returns:
|
|
Forest coverage as percentage (0-100)
|
|
"""
|
|
if not elements:
|
|
return 0.0
|
|
|
|
# Approximate: count forest ways and estimate coverage
|
|
# This is a simplified calculation
|
|
forest_ways = len([e for e in elements if e.get('type') == 'way'])
|
|
|
|
# Rough heuristic: each forest way covers ~0.1 km²
|
|
# Total search area = π * r²
|
|
search_area_km2 = math.pi * (radius_m / 1000) ** 2
|
|
estimated_forest_km2 = forest_ways * 0.1
|
|
|
|
coverage_pct = min(100.0, (estimated_forest_km2 / search_area_km2) * 100)
|
|
|
|
return round(coverage_pct, 1)
|
|
|
|
|
|
async def get_zone_features(lat: float, lon: float, direction: str, radius_m: int) -> Dict:
|
|
"""
|
|
Get geographic features for a zone using Overpass API.
|
|
|
|
Args:
|
|
lat, lon: Center point
|
|
direction: Sector direction
|
|
radius_m: Search radius in meters
|
|
|
|
Returns:
|
|
Dict with roads_km, water_distance_km, settlement_distance_km, forest_pct
|
|
"""
|
|
# Query roads
|
|
roads_query = f"""
|
|
[out:json];
|
|
(
|
|
way[highway](around:{radius_m},{lat},{lon});
|
|
);
|
|
out geom;
|
|
"""
|
|
roads_data = await query_overpass(roads_query)
|
|
roads_km = calculate_road_length(roads_data.get('elements', []))
|
|
|
|
# Query water bodies
|
|
water_query = f"""
|
|
[out:json];
|
|
(
|
|
node[natural=water](around:{radius_m},{lat},{lon});
|
|
way[natural=water](around:{radius_m},{lat},{lon});
|
|
);
|
|
out center;
|
|
"""
|
|
water_data = await query_overpass(water_query)
|
|
water_distance = find_nearest_distance(lat, lon, water_data.get('elements', []))
|
|
|
|
# Query settlements
|
|
settlement_query = f"""
|
|
[out:json];
|
|
(
|
|
node[place~"village|town|city"](around:{radius_m},{lat},{lon});
|
|
);
|
|
out;
|
|
"""
|
|
settlement_data = await query_overpass(settlement_query)
|
|
settlement_distance = find_nearest_distance(lat, lon, settlement_data.get('elements', []))
|
|
|
|
# Query forests
|
|
forest_query = f"""
|
|
[out:json];
|
|
(
|
|
way[landuse=forest](around:{radius_m},{lat},{lon});
|
|
way[natural=wood](around:{radius_m},{lat},{lon});
|
|
);
|
|
out geom;
|
|
"""
|
|
forest_data = await query_overpass(forest_query)
|
|
forest_pct = calculate_forest_coverage(forest_data.get('elements', []), radius_m)
|
|
|
|
# Calculate road density (km of roads per km²)
|
|
search_area_km2 = math.pi * (radius_m / 1000) ** 2
|
|
road_density = roads_km / search_area_km2 if search_area_km2 > 0 else 0.0
|
|
|
|
return {
|
|
'roads_km': roads_km,
|
|
'road_density': round(road_density, 2),
|
|
'water_distance_km': water_distance,
|
|
'settlement_distance_km': settlement_distance,
|
|
'forest_pct': forest_pct
|
|
}
|
|
|
|
|
|
async def build_search_zones(lat: float, lon: float, case_data: dict, max_distance_km: float = 3.0) -> List[Zone]:
|
|
"""
|
|
Build search zones around a point.
|
|
|
|
Creates 8 directional sectors (N, NE, E, SE, S, SW, W, NW) at multiple distances
|
|
(500m, 1000m, 2000m, 5000m) and queries geographic features for each.
|
|
|
|
Args:
|
|
lat: Latitude of search origin
|
|
lon: Longitude of search origin
|
|
case_data: Case information (for future enhancements)
|
|
|
|
Returns:
|
|
List of Zone objects with geographic features
|
|
"""
|
|
zones = []
|
|
|
|
for distance_m in _build_search_distances(max_distance_km):
|
|
for direction in DIRECTIONS:
|
|
# Get features for this zone
|
|
features = await get_zone_features(lat, lon, direction, distance_m)
|
|
|
|
zone = Zone(
|
|
direction=direction,
|
|
distance_km=distance_m / 1000,
|
|
forest_pct=features['forest_pct'],
|
|
road_density=features['road_density'],
|
|
water_distance_km=features['water_distance_km'],
|
|
settlement_distance_km=features['settlement_distance_km']
|
|
)
|
|
|
|
zones.append(zone)
|
|
|
|
return zones
|