diff --git a/lib/api.py b/lib/api.py index cd32e33..a127866 100644 --- a/lib/api.py +++ b/lib/api.py @@ -2768,8 +2768,8 @@ def api_offroute(): # Parse options mode = data.get("mode", "foot") - if mode not in ("foot", "mtb", "atv"): - return jsonify({"status": "error", "message": "mode must be foot, mtb, or atv"}), 400 + if mode not in ("foot", "mtb", "atv", "vehicle"): + return jsonify({"status": "error", "message": "mode must be foot, mtb, atv, or vehicle"}), 400 boundary_mode = data.get("boundary_mode", "pragmatic") if boundary_mode not in ("strict", "pragmatic", "emergency"): diff --git a/lib/offroute/barriers.py b/lib/offroute/barriers.py index 7fcad75..f68e892 100644 --- a/lib/offroute/barriers.py +++ b/lib/offroute/barriers.py @@ -1,11 +1,12 @@ """ -PAD-US barrier layer for OFFROUTE. +PAD-US barrier and wilderness layers for OFFROUTE. -Provides access to the PAD-US land ownership raster for routing decisions. -Cells with value 255 represent closed/restricted areas (Pub_Access = XA). +Provides access to: +1. Barrier raster (Pub_Access = 'XA' - closed/restricted areas) +2. Wilderness raster (Des_Tp = 'WA' - designated wilderness areas) -Build function rasterizes PAD-US geodatabase to aligned GeoTIFF. -Runtime functions read the raster and resample to match elevation grids. +Build functions rasterize PAD-US geodatabase to aligned GeoTIFFs. +Runtime functions read the rasters and resample to match elevation grids. """ import numpy as np from pathlib import Path @@ -23,6 +24,7 @@ except ImportError: # Paths DEFAULT_BARRIERS_PATH = Path("/mnt/nav/worldcover/padus_barriers.tif") +DEFAULT_WILDERNESS_PATH = Path("/mnt/nav/worldcover/wilderness.tif") PADUS_GDB_PATH = Path("/mnt/nav/padus/PADUS4_0_Geodatabase.gdb") PADUS_LAYER = "PADUS4_0Combined_Proclamation_Marine_Fee_Designation_Easement" @@ -39,7 +41,7 @@ PIXEL_SIZE = 0.0003 # ~33m class BarrierReader: - """Reader for PAD-US barrier raster.""" + """Reader for PAD-US barrier raster (closed/restricted areas).""" def __init__(self, barrier_path: Path = DEFAULT_BARRIERS_PATH): self.barrier_path = barrier_path @@ -77,32 +79,86 @@ class BarrierReader: 0 = public/accessible """ ds = self._open() - - # Create a window from the bounding box window = from_bounds(west, south, east, north, ds.transform) - - # Read with resampling to target shape barriers = ds.read( 1, window=window, out_shape=target_shape, resampling=Resampling.nearest ) - return barriers def sample_point(self, lat: float, lon: float) -> int: """Sample barrier value at a single point.""" ds = self._open() - - # Get pixel coordinates row, col = ds.index(lon, lat) - - # Check bounds if row < 0 or row >= ds.height or col < 0 or col >= ds.width: - return 0 # Out of bounds = accessible + return 0 + window = rasterio.windows.Window(col, row, 1, 1) + value = ds.read(1, window=window) + return int(value[0, 0]) - # Read single pixel + def close(self): + """Close the dataset.""" + if self._dataset is not None: + self._dataset.close() + self._dataset = None + + +class WildernessReader: + """Reader for PAD-US wilderness raster (designated wilderness areas).""" + + def __init__(self, wilderness_path: Path = DEFAULT_WILDERNESS_PATH): + self.wilderness_path = wilderness_path + self._dataset = None + + def _open(self): + """Lazy open the dataset.""" + if self._dataset is None: + if not self.wilderness_path.exists(): + raise FileNotFoundError( + f"Wilderness raster not found at {self.wilderness_path}. " + f"Run build_wilderness_raster() first." + ) + self._dataset = rasterio.open(self.wilderness_path) + return self._dataset + + def get_wilderness_grid( + self, + south: float, + north: float, + west: float, + east: float, + target_shape: Tuple[int, int] + ) -> np.ndarray: + """ + Get wilderness values for a bounding box, resampled to target shape. + + Args: + south, north, west, east: Bounding box coordinates (WGS84) + target_shape: (rows, cols) to resample to (matches elevation grid) + + Returns: + np.ndarray of uint8 wilderness values: + 255 = designated wilderness area + 0 = not wilderness + """ + ds = self._open() + window = from_bounds(west, south, east, north, ds.transform) + wilderness = ds.read( + 1, + window=window, + out_shape=target_shape, + resampling=Resampling.nearest + ) + return wilderness + + def sample_point(self, lat: float, lon: float) -> int: + """Sample wilderness value at a single point.""" + ds = self._open() + row, col = ds.index(lon, lat) + if row < 0 or row >= ds.height or col < 0 or col >= ds.width: + return 0 window = rasterio.windows.Window(col, row, 1, 1) value = ds.read(1, window=window) return int(value[0, 0]) @@ -124,22 +180,12 @@ def build_barriers_raster( Build the PAD-US barriers raster from the source geodatabase. Extracts polygons where Pub_Access = 'XA' (Closed) and rasterizes them. - - Args: - output_path: Output GeoTIFF path - gdb_path: Path to PAD-US geodatabase - pixel_size: Pixel size in degrees - bounds: CONUS bounding box - - Returns: - Path to the created raster """ import shutil if not gdb_path.exists(): raise FileNotFoundError(f"PAD-US geodatabase not found at {gdb_path}") - # Check for required tools if not shutil.which('ogr2ogr'): raise RuntimeError("ogr2ogr not found. Install GDAL.") if not shutil.which('gdal_rasterize'): @@ -154,7 +200,6 @@ def build_barriers_raster( print(f" Bounds: {bounds}") with tempfile.TemporaryDirectory() as tmpdir: - # Step 1: Extract closed areas and reproject to WGS84 closed_gpkg = Path(tmpdir) / "closed_areas.gpkg" print(f"\n[1/3] Extracting closed areas (Pub_Access = 'XA')...") @@ -176,28 +221,23 @@ def build_barriers_raster( print(f"STDERR: {result.stderr}") raise RuntimeError(f"ogr2ogr failed: {result.stderr}") - # Check feature count info_cmd = ["ogrinfo", "-so", str(closed_gpkg), "closed_areas"] info_result = subprocess.run(info_cmd, capture_output=True, text=True) print(f" Extraction result:\n{info_result.stdout}") - # Step 2: Create empty raster print(f"\n[2/3] Creating raster grid...") width = int((bounds['east'] - bounds['west']) / pixel_size) height = int((bounds['north'] - bounds['south']) / pixel_size) - print(f" Grid size: {width} x {height} pixels") - print(f" Memory estimate: {width * height / 1e6:.1f} MB") - # Step 3: Rasterize print(f"\n[3/3] Rasterizing closed areas...") rasterize_cmd = [ "gdal_rasterize", "-burn", "255", "-init", "0", - "-a_nodata", "0", # No nodata - 0 means accessible + "-a_nodata", "0", "-te", str(bounds['west']), str(bounds['south']), str(bounds['east']), str(bounds['north']), "-tr", str(pixel_size), str(pixel_size), @@ -214,14 +254,10 @@ def build_barriers_raster( print(f"STDERR: {result.stderr}") raise RuntimeError(f"gdal_rasterize failed: {result.stderr}") - # Verify output print(f"\n[Done] Verifying output...") with rasterio.open(output_path) as ds: print(f" Size: {ds.width} x {ds.height}") print(f" CRS: {ds.crs}") - print(f" Bounds: {ds.bounds}") - - # Sample a few tiles to check sample = ds.read(1, window=rasterio.windows.Window(0, 0, 1000, 1000)) closed_count = np.sum(sample == 255) print(f" Sample (1000x1000): {closed_count} closed cells") @@ -232,17 +268,140 @@ def build_barriers_raster( return output_path +def build_wilderness_raster( + output_path: Path = DEFAULT_WILDERNESS_PATH, + gdb_path: Path = PADUS_GDB_PATH, + pixel_size: float = PIXEL_SIZE, + bounds: dict = CONUS_BOUNDS, +) -> Path: + """ + Build the PAD-US wilderness raster from the source geodatabase. + + Extracts polygons where Des_Tp = 'WA' (Wilderness Area) and rasterizes them. + """ + import shutil + + if not gdb_path.exists(): + raise FileNotFoundError(f"PAD-US geodatabase not found at {gdb_path}") + + if not shutil.which('ogr2ogr'): + raise RuntimeError("ogr2ogr not found. Install GDAL.") + if not shutil.which('gdal_rasterize'): + raise RuntimeError("gdal_rasterize not found. Install GDAL.") + + output_path.parent.mkdir(parents=True, exist_ok=True) + + print(f"Building PAD-US wilderness raster...") + print(f" Source: {gdb_path}") + print(f" Output: {output_path}") + print(f" Pixel size: {pixel_size} degrees (~{pixel_size * 111000:.0f}m)") + print(f" Bounds: {bounds}") + + with tempfile.TemporaryDirectory() as tmpdir: + wilderness_gpkg = Path(tmpdir) / "wilderness_areas.gpkg" + + print(f"\n[1/3] Extracting wilderness areas (Des_Tp = 'WA')...") + + ogr_cmd = [ + "ogr2ogr", + "-f", "GPKG", + str(wilderness_gpkg), + str(gdb_path), + PADUS_LAYER, + "-where", "Des_Tp = 'WA'", + "-t_srs", "EPSG:4326", + "-nlt", "MULTIPOLYGON", + "-nln", "wilderness_areas", + ] + + result = subprocess.run(ogr_cmd, capture_output=True, text=True) + if result.returncode != 0: + print(f"STDERR: {result.stderr}") + raise RuntimeError(f"ogr2ogr failed: {result.stderr}") + + info_cmd = ["ogrinfo", "-so", str(wilderness_gpkg), "wilderness_areas"] + info_result = subprocess.run(info_cmd, capture_output=True, text=True) + print(f" Extraction result:\n{info_result.stdout}") + + print(f"\n[2/3] Creating raster grid...") + + width = int((bounds['east'] - bounds['west']) / pixel_size) + height = int((bounds['north'] - bounds['south']) / pixel_size) + print(f" Grid size: {width} x {height} pixels") + + print(f"\n[3/3] Rasterizing wilderness areas...") + + rasterize_cmd = [ + "gdal_rasterize", + "-burn", "255", + "-init", "0", + "-a_nodata", "0", + "-te", str(bounds['west']), str(bounds['south']), + str(bounds['east']), str(bounds['north']), + "-tr", str(pixel_size), str(pixel_size), + "-ot", "Byte", + "-co", "COMPRESS=LZW", + "-co", "TILED=YES", + "-l", "wilderness_areas", + str(wilderness_gpkg), + str(output_path), + ] + + result = subprocess.run(rasterize_cmd, capture_output=True, text=True) + if result.returncode != 0: + print(f"STDERR: {result.stderr}") + raise RuntimeError(f"gdal_rasterize failed: {result.stderr}") + + print(f"\n[Done] Verifying output...") + with rasterio.open(output_path) as ds: + print(f" Size: {ds.width} x {ds.height}") + print(f" CRS: {ds.crs}") + sample = ds.read(1, window=rasterio.windows.Window(0, 0, 1000, 1000)) + wilderness_count = np.sum(sample == 255) + print(f" Sample (1000x1000): {wilderness_count} wilderness cells") + + file_size = output_path.stat().st_size / (1024**2) + print(f" File size: {file_size:.1f} MB") + + return output_path + + if __name__ == "__main__": import sys - if len(sys.argv) > 1 and sys.argv[1] == "build": - # Build the raster - print("="*60) - print("PAD-US Barriers Raster Build") - print("="*60) - build_barriers_raster() + if len(sys.argv) > 1: + cmd = sys.argv[1] + + if cmd == "build": + print("=" * 60) + print("PAD-US Barriers Raster Build") + print("=" * 60) + build_barriers_raster() + + elif cmd == "build-wilderness": + print("=" * 60) + print("PAD-US Wilderness Raster Build") + print("=" * 60) + build_wilderness_raster() + + elif cmd == "build-all": + print("=" * 60) + print("Building all PAD-US rasters") + print("=" * 60) + build_barriers_raster() + print("\n") + build_wilderness_raster() + + else: + print(f"Unknown command: {cmd}") + print("Usage:") + print(" python barriers.py build # Build barriers raster") + print(" python barriers.py build-wilderness # Build wilderness raster") + print(" python barriers.py build-all # Build both rasters") + sys.exit(1) + else: - # Test the reader + # Test readers print("Testing BarrierReader...") if not DEFAULT_BARRIERS_PATH.exists(): @@ -251,16 +410,31 @@ if __name__ == "__main__": sys.exit(1) reader = BarrierReader() - - # Test grid read for Idaho area barriers = reader.get_barrier_grid( south=42.2, north=42.6, west=-114.8, east=-113.8, target_shape=(400, 1000) ) - print(f"\nGrid test shape: {barriers.shape}") + print(f"\nBarrier grid shape: {barriers.shape}") print(f"Unique values: {np.unique(barriers)}") closed_cells = np.sum(barriers == 255) print(f"Closed cells: {closed_cells} ({100*closed_cells/barriers.size:.2f}%)") - reader.close() - print("\nBarrierReader test complete.") + + print("\nTesting WildernessReader...") + + if not DEFAULT_WILDERNESS_PATH.exists(): + print(f"Wilderness raster not found at {DEFAULT_WILDERNESS_PATH}") + print(f"Run: python barriers.py build-wilderness") + else: + wilderness_reader = WildernessReader() + wilderness = wilderness_reader.get_wilderness_grid( + south=42.2, north=42.6, west=-114.8, east=-113.8, + target_shape=(400, 1000) + ) + print(f"Wilderness grid shape: {wilderness.shape}") + print(f"Unique values: {np.unique(wilderness)}") + wilderness_cells = np.sum(wilderness == 255) + print(f"Wilderness cells: {wilderness_cells} ({100*wilderness_cells/wilderness.size:.2f}%)") + wilderness_reader.close() + + print("\nDone.") diff --git a/lib/offroute/cost.py b/lib/offroute/cost.py index 5f3618c..c3b6a5a 100644 --- a/lib/offroute/cost.py +++ b/lib/offroute/cost.py @@ -1,43 +1,207 @@ """ -Tobler off-path hiking cost function for OFFROUTE. +Multi-mode travel cost functions for OFFROUTE. -Computes travel time cost based on terrain slope using Tobler's -hiking function with off-trail penalty. Optionally applies friction -multipliers from land cover data, trail corridors, and barrier grids. +Supports four travel modes: foot, mtb, atv, vehicle. +Each mode has its own speed function, max slope, trail access rules, +and terrain friction overrides. + +Mode profiles are data-driven — adding a new mode means adding a profile entry. """ import math import numpy as np -from typing import Optional, Literal +from dataclasses import dataclass, field +from typing import Optional, Literal, Dict, Callable -# Maximum passable slope in degrees -MAX_SLOPE_DEG = 40.0 +# ═══════════════════════════════════════════════════════════════════════════════ +# SPEED FUNCTIONS +# ═══════════════════════════════════════════════════════════════════════════════ + +def tobler_off_path_speed(grade: np.ndarray, base_speed: float = 6.0) -> np.ndarray: + """ + Tobler off-path hiking function. + + W = 0.6 * base_speed * exp(-3.5 * |S + 0.05|) + + Peak ~3.6 km/h at grade = -0.05 (slight downhill). + The 0.6 multiplier is the off-trail penalty. + """ + return 0.6 * base_speed * np.exp(-3.5 * np.abs(grade + 0.05)) + + +def herzog_wheeled_speed(grade: np.ndarray, base_speed: float = 12.0) -> np.ndarray: + """ + Herzog wheeled-transport polynomial. + + Relative speed factor: + 1 / (1337.8·S^6 + 278.19·S^5 − 517.39·S^4 − 78.199·S^3 + 93.419·S^2 + 19.825·|S| + 1.64) + + Multiply by base_speed to get km/h. + """ + S = grade + S_abs = np.abs(S) + + # Herzog polynomial (returns relative speed factor 0-1) + denom = (1337.8 * S**6 + 278.19 * S**5 - 517.39 * S**4 + - 78.199 * S**3 + 93.419 * S**2 + 19.825 * S_abs + 1.64) + + # Avoid division by zero and negative speeds + denom = np.maximum(denom, 0.1) + rel_speed = 1.0 / denom + + # Clamp relative speed to reasonable bounds (0.05 to 1.5) + rel_speed = np.clip(rel_speed, 0.05, 1.5) + + return base_speed * rel_speed + + +def linear_degrade_speed(grade: np.ndarray, base_speed: float = 40.0, max_grade: float = 0.364) -> np.ndarray: + """ + Linear speed degradation with slope. + + speed = base_speed * max(0, 1 - |grade| / max_grade) + + max_grade = tan(20°) ≈ 0.364 for 20° max slope. + """ + speed = base_speed * np.maximum(0, 1.0 - np.abs(grade) / max_grade) + return np.maximum(speed, 0.1) # Minimum crawl speed + + +# ═══════════════════════════════════════════════════════════════════════════════ +# MODE PROFILES (Data-driven configuration) +# ═══════════════════════════════════════════════════════════════════════════════ + +@dataclass +class ModeProfile: + """Configuration for a travel mode.""" + + name: str + description: str + + # Speed function parameters + speed_function: str # "tobler", "herzog", "linear" + base_speed_kmh: float + max_slope_deg: float + + # Trail access: trail_value -> friction multiplier (None = impassable) + # Trail values: 5=road, 15=track, 25=foot trail + trail_friction: Dict[int, Optional[float]] = field(default_factory=dict) + + # Off-trail terrain friction overrides (by WorldCover class) + # These MULTIPLY the base WorldCover friction + # None = use default, np.inf = impassable + # WorldCover values: 10=tree, 20=shrub, 30=grass, 40=crop, 50=urban, + # 60=bare, 80=water, 90=wetland, 95=mangrove, 100=moss + terrain_friction_override: Dict[int, Optional[float]] = field(default_factory=dict) + + # Should wilderness areas be impassable? + wilderness_impassable: bool = False + + # For vehicle mode: can traverse off-trail flat terrain? + off_trail_flat_threshold_deg: float = 0.0 # 0 = no off-trail allowed + off_trail_flat_friction: float = np.inf # friction if allowed + + +# Define all mode profiles +MODE_PROFILES: Dict[str, ModeProfile] = { + "foot": ModeProfile( + name="foot", + description="Hiking on foot (Tobler off-path model)", + speed_function="tobler", + base_speed_kmh=6.0, + max_slope_deg=40.0, + trail_friction={ + 5: 0.1, # road + 15: 0.3, # track + 25: 0.5, # foot trail + }, + terrain_friction_override={ + # Use default WorldCover friction for foot mode + }, + wilderness_impassable=False, + ), + + "mtb": ModeProfile( + name="mtb", + description="Mountain bike / dirt bike (Herzog wheeled model)", + speed_function="herzog", + base_speed_kmh=12.0, + max_slope_deg=25.0, + trail_friction={ + 5: 0.1, # road + 15: 0.2, # track + 25: 0.5, # foot trail (rideable but slow) + }, + terrain_friction_override={ + 30: 2.0, # Grassland: rideable but slow + 20: 4.0, # Shrubland: barely rideable + 10: 8.0, # Tree cover/forest: effectively impassable + 60: 3.0, # Bare/rocky + 90: np.inf, # Wetland: impassable + 95: np.inf, # Mangrove: impassable + 80: np.inf, # Water: impassable + }, + wilderness_impassable=True, + ), + + "atv": ModeProfile( + name="atv", + description="ATV / side-by-side (Herzog wheeled model, higher base speed)", + speed_function="herzog", + base_speed_kmh=25.0, + max_slope_deg=30.0, + trail_friction={ + 5: 0.1, # road + 15: 0.3, # track + 25: None, # foot trail: impassable (too narrow) + }, + terrain_friction_override={ + 30: 1.5, # Grassland: passable + 20: 3.0, # Shrubland: rough + 10: np.inf, # Forest: impassable + 60: 2.0, # Bare/rocky + 90: np.inf, # Wetland: impassable + 95: np.inf, # Mangrove: impassable + 80: np.inf, # Water: impassable + }, + wilderness_impassable=True, + ), + + "vehicle": ModeProfile( + name="vehicle", + description="4x4 truck / jeep (linear speed degradation)", + speed_function="linear", + base_speed_kmh=40.0, + max_slope_deg=20.0, + trail_friction={ + 5: 0.1, # road + 15: 0.5, # track (rough but passable) + 25: None, # foot trail: impassable + }, + terrain_friction_override={ + # All off-trail terrain is impassable by default + 10: np.inf, # Forest + 20: np.inf, # Shrubland + 30: np.inf, # Grassland (except flat - see below) + 40: np.inf, # Cropland (except flat - see below) + 60: np.inf, # Bare + 90: np.inf, # Wetland + 95: np.inf, # Mangrove + 80: np.inf, # Water + }, + wilderness_impassable=True, + off_trail_flat_threshold_deg=5.0, # Can drive on flat fields + off_trail_flat_friction=5.0, # But very slow + ), +} -# Tobler off-path parameters -TOBLER_BASE_SPEED = 6.0 -TOBLER_OFF_TRAIL_MULT = 0.6 # Pragmatic mode friction multiplier for private land PRAGMATIC_BARRIER_MULTIPLIER = 5.0 -# Trail value to friction multiplier mapping -# Trail friction REPLACES land cover friction (a road through forest is still easy) -TRAIL_FRICTION_MAP = { - 5: 0.1, # road - 15: 0.3, # track - 25: 0.5, # foot trail -} - - -def tobler_speed(grade: float) -> float: - """ - Calculate hiking speed using Tobler's off-path function. - - speed_kmh = 0.6 * 6.0 * exp(-3.5 * |grade + 0.05|) - - Peak speed is ~3.6 km/h at grade = -0.05 (slight downhill). - """ - return TOBLER_OFF_TRAIL_MULT * TOBLER_BASE_SPEED * math.exp(-3.5 * abs(grade + 0.05)) +# ═══════════════════════════════════════════════════════════════════════════════ +# COST GRID COMPUTATION +# ═══════════════════════════════════════════════════════════════════════════════ def compute_cost_grid( elevation: np.ndarray, @@ -45,16 +209,16 @@ def compute_cost_grid( cell_size_lat_m: float = None, cell_size_lon_m: float = None, friction: Optional[np.ndarray] = None, + friction_raw: Optional[np.ndarray] = None, trails: Optional[np.ndarray] = None, barriers: Optional[np.ndarray] = None, - boundary_mode: Literal["strict", "pragmatic", "emergency"] = "pragmatic" + wilderness: Optional[np.ndarray] = None, + boundary_mode: Literal["strict", "pragmatic", "emergency"] = "pragmatic", + mode: Literal["foot", "mtb", "atv", "vehicle"] = "foot" ) -> np.ndarray: """ Compute isotropic travel cost grid from elevation data. - Each cell's cost represents the time (in seconds) to traverse that cell, - based on the average slope from neighboring cells. - Args: elevation: 2D array of elevation values in meters cell_size_m: Average cell size in meters @@ -63,30 +227,29 @@ def compute_cost_grid( friction: Optional 2D array of friction multipliers (WorldCover). Values should be float (1.0 = baseline, 2.0 = 2x slower). np.inf marks impassable cells. - If None, no friction is applied (backward compatible). + friction_raw: Optional 2D array of raw WorldCover class values (uint8). + Used for mode-specific terrain overrides. + Values: 10=tree, 20=shrub, 30=grass, etc. trails: Optional 2D array of trail values (uint8). - 0 = no trail (use friction) - 5 = road (0.1× friction, replaces WorldCover) - 15 = track (0.3× friction, replaces WorldCover) - 25 = foot trail (0.5× friction, replaces WorldCover) - Trail friction REPLACES land cover friction where trails exist. - If None, no trail burn-in is applied. + 0 = no trail, 5 = road, 15 = track, 25 = foot trail barriers: Optional 2D array of barrier values (uint8). - 255 = closed/restricted area (from PAD-US Pub_Access = XA). - 0 = accessible. - If None, no barriers are applied. - boundary_mode: How to handle private/restricted land barriers: - "strict" - cells with barrier=255 become impassable (np.inf) - "pragmatic" - cells with barrier=255 get 5.0x friction penalty - "emergency" - barriers are ignored entirely - Default: "pragmatic" + 255 = closed/restricted area (PAD-US Pub_Access = XA). + wilderness: Optional[np.ndarray] of wilderness values (uint8). + 255 = designated wilderness area. + boundary_mode: How to handle barriers ("strict", "pragmatic", "emergency") + mode: Travel mode ("foot", "mtb", "atv", "vehicle") Returns: 2D array of travel cost in seconds per cell. np.inf for impassable cells. """ if boundary_mode not in ("strict", "pragmatic", "emergency"): - raise ValueError(f"boundary_mode must be 'strict', 'pragmatic', or 'emergency', got '{boundary_mode}'") + raise ValueError(f"boundary_mode must be 'strict', 'pragmatic', or 'emergency'") + + if mode not in MODE_PROFILES: + raise ValueError(f"mode must be one of {list(MODE_PROFILES.keys())}") + + profile = MODE_PROFILES[mode] if cell_size_lat_m is None: cell_size_lat_m = cell_size_m @@ -95,120 +258,212 @@ def compute_cost_grid( rows, cols = elevation.shape - # Compute gradients in both directions - dy = np.zeros_like(elevation) - dx = np.zeros_like(elevation) + # ─── Compute gradients (in-place where possible) ───────────────────────── + # Use float32 to reduce memory footprint + grade = np.zeros(elevation.shape, dtype=np.float32) - # Central differences for interior, forward/backward at edges - dy[1:-1, :] = (elevation[:-2, :] - elevation[2:, :]) / (2 * cell_size_lat_m) - dy[0, :] = (elevation[0, :] - elevation[1, :]) / cell_size_lat_m - dy[-1, :] = (elevation[-2, :] - elevation[-1, :]) / cell_size_lat_m + # Compute dy contribution to grade squared + dy_contrib = np.zeros(elevation.shape, dtype=np.float32) + dy_contrib[1:-1, :] = ((elevation[:-2, :] - elevation[2:, :]) / (2 * cell_size_lat_m)) ** 2 + dy_contrib[0, :] = ((elevation[0, :] - elevation[1, :]) / cell_size_lat_m) ** 2 + dy_contrib[-1, :] = ((elevation[-2, :] - elevation[-1, :]) / cell_size_lat_m) ** 2 - dx[:, 1:-1] = (elevation[:, 2:] - elevation[:, :-2]) / (2 * cell_size_lon_m) - dx[:, 0] = (elevation[:, 1] - elevation[:, 0]) / cell_size_lon_m - dx[:, -1] = (elevation[:, -1] - elevation[:, -2]) / cell_size_lon_m + # Compute dx contribution and add to dy_contrib in-place + dy_contrib[:, 1:-1] += ((elevation[:, 2:] - elevation[:, :-2]) / (2 * cell_size_lon_m)) ** 2 + dy_contrib[:, 0] += ((elevation[:, 1] - elevation[:, 0]) / cell_size_lon_m) ** 2 + dy_contrib[:, -1] += ((elevation[:, -1] - elevation[:, -2]) / cell_size_lon_m) ** 2 - # Compute slope magnitude (grade = rise/run) - grade_magnitude = np.sqrt(dx**2 + dy**2) + # grade = sqrt(dx^2 + dy^2) + np.sqrt(dy_contrib, out=grade) + del dy_contrib # Free memory immediately - # Convert to slope angle in degrees - slope_deg = np.degrees(np.arctan(grade_magnitude)) + # ─── Compute speed based on mode ───────────────────────────────────────── + max_grade_val = np.tan(np.radians(profile.max_slope_deg)) - # Compute speed for each cell using Tobler function - speed_kmh = TOBLER_OFF_TRAIL_MULT * TOBLER_BASE_SPEED * np.exp(-3.5 * np.abs(grade_magnitude + 0.05)) + if profile.speed_function == "tobler": + speed_kmh = tobler_off_path_speed(grade, profile.base_speed_kmh) + elif profile.speed_function == "herzog": + speed_kmh = herzog_wheeled_speed(grade, profile.base_speed_kmh) + elif profile.speed_function == "linear": + speed_kmh = linear_degrade_speed(grade, profile.base_speed_kmh, max_grade_val) + else: + raise ValueError(f"Unknown speed function: {profile.speed_function}") - # Convert speed to time cost (seconds to traverse one cell) + # ─── Base cost (seconds per cell) ───────────────────────────────────────── avg_cell_size = (cell_size_lat_m + cell_size_lon_m) / 2 - cost = avg_cell_size * 3.6 / speed_kmh + cost = (avg_cell_size * 3.6) / speed_kmh + del speed_kmh - # Set impassable cells (slope > MAX_SLOPE_DEG) to infinity - cost[slope_deg > MAX_SLOPE_DEG] = np.inf + # ─── Max slope limit ────────────────────────────────────────────────────── + cost[grade > max_grade_val] = np.inf - # Handle NaN elevations (no data) + # ─── NaN elevations ────────────────────────────────────────────────────── cost[np.isnan(elevation)] = np.inf - # Build effective friction array - # Start with WorldCover friction if provided, else 1.0 + # ─── Apply friction in-place ───────────────────────────────────────────── + # Instead of creating effective_friction copy, apply directly to cost + + # Start with base friction if friction is not None: if friction.shape != elevation.shape: - raise ValueError( - f"Friction shape {friction.shape} does not match elevation shape {elevation.shape}" - ) - effective_friction = friction.copy() - else: - effective_friction = np.ones(elevation.shape, dtype=np.float32) + raise ValueError(f"Friction shape mismatch") + np.multiply(cost, friction, out=cost) - # Apply trail burn-in: trails REPLACE land cover friction + # ─── Mode-specific terrain friction overrides (memory-efficient) ───────── + if friction_raw is not None and profile.terrain_friction_override: + if friction_raw.shape != elevation.shape: + raise ValueError(f"Friction_raw shape mismatch") + + # Process all overrides without creating large intermediate masks + for wc_class, override in profile.terrain_friction_override.items(): + if override is not None: + if override == np.inf: + # Use np.where for in-place-like behavior + np.putmask(cost, friction_raw == wc_class, np.inf) + else: + # Multiply cost where friction_raw matches + # Using a loop with putmask is more memory efficient + mask = friction_raw == wc_class + cost[mask] *= override + del mask + + # ─── Vehicle mode: allow flat grassland/cropland ───────────────────────── + if mode == "vehicle" and profile.off_trail_flat_threshold_deg > 0: + if friction_raw is not None: + # Compute slope in degrees for flat terrain check + slope_deg = np.degrees(np.arctan(grade)) + # Flat grassland or cropland - recompute cost for these cells + flat_field_mask = ( + (slope_deg <= profile.off_trail_flat_threshold_deg) & + ((friction_raw == 30) | (friction_raw == 40)) + ) + del slope_deg + # Recalculate cost for these cells with flat field friction + if np.any(flat_field_mask): + base_time = avg_cell_size * 3.6 / linear_degrade_speed( + grade[flat_field_mask], profile.base_speed_kmh, max_grade_val + ) + cost[flat_field_mask] = base_time * profile.off_trail_flat_friction + del base_time + del flat_field_mask + + # ─── Trail friction (mode-specific) ────────────────────────────────────── if trails is not None: if trails.shape != elevation.shape: - raise ValueError( - f"Trails shape {trails.shape} does not match elevation shape {elevation.shape}" - ) - # Replace friction where trails exist - for trail_value, trail_friction in TRAIL_FRICTION_MAP.items(): - trail_mask = trails == trail_value - effective_friction[trail_mask] = trail_friction + raise ValueError(f"Trails shape mismatch") - # Apply friction to cost - cost = cost * effective_friction + for trail_value, trail_friction in profile.trail_friction.items(): + if trail_friction is None: + # Impassable for this mode + np.putmask(cost, trails == trail_value, np.inf) + else: + # Trail friction REPLACES terrain friction + # Recalculate cost = base_time * trail_friction + trail_mask = trails == trail_value + if np.any(trail_mask): + # Get base travel time (without friction) + if profile.speed_function == "tobler": + trail_speed = tobler_off_path_speed(grade[trail_mask], profile.base_speed_kmh) + elif profile.speed_function == "herzog": + trail_speed = herzog_wheeled_speed(grade[trail_mask], profile.base_speed_kmh) + else: + trail_speed = linear_degrade_speed( + grade[trail_mask], profile.base_speed_kmh, max_grade_val + ) + cost[trail_mask] = (avg_cell_size * 3.6 / trail_speed) * trail_friction + del trail_speed + del trail_mask - # Apply barriers based on boundary_mode + # ─── Wilderness areas (mode-specific) ──────────────────────────────────── + if wilderness is not None and profile.wilderness_impassable: + if wilderness.shape != elevation.shape: + raise ValueError(f"Wilderness shape mismatch") + np.putmask(cost, wilderness == 255, np.inf) + + # ─── Barriers (private land) ───────────────────────────────────────────── if barriers is not None and boundary_mode != "emergency": if barriers.shape != elevation.shape: - raise ValueError( - f"Barriers shape {barriers.shape} does not match elevation shape {elevation.shape}" - ) - - barrier_mask = barriers == 255 + raise ValueError(f"Barriers shape mismatch") if boundary_mode == "strict": - # Mark closed/restricted areas as impassable - cost[barrier_mask] = np.inf + np.putmask(cost, barriers == 255, np.inf) elif boundary_mode == "pragmatic": - # Apply friction penalty to closed/restricted areas - cost[barrier_mask] = cost[barrier_mask] * PRAGMATIC_BARRIER_MULTIPLIER + barrier_mask = barriers == 255 + cost[barrier_mask] *= PRAGMATIC_BARRIER_MULTIPLIER + del barrier_mask return cost +# ═══════════════════════════════════════════════════════════════════════════════ +# LEGACY API (backward compatibility) +# ═══════════════════════════════════════════════════════════════════════════════ + +def tobler_speed(grade: float) -> float: + """Legacy single-value Tobler speed function.""" + return 0.6 * 6.0 * math.exp(-3.5 * abs(grade + 0.05)) + + +# ═══════════════════════════════════════════════════════════════════════════════ +# TESTING +# ═══════════════════════════════════════════════════════════════════════════════ + if __name__ == "__main__": - print("Testing Tobler speed function:") - for grade in [-0.3, -0.1, -0.05, 0.0, 0.05, 0.1, 0.3]: - speed = tobler_speed(grade) - print(f" Grade {grade:+.2f}: {speed:.2f} km/h") + print("=" * 70) + print("OFFROUTE Multi-Mode Cost Function Tests") + print("=" * 70) - print("\nTesting cost grid computation (no friction, no trails):") - elev = np.arange(100).reshape(10, 10).astype(np.float32) * 10 - cost = compute_cost_grid(elev, cell_size_m=30.0) - print(f" Elevation range: {elev.min():.0f} - {elev.max():.0f} m") - finite = cost[~np.isinf(cost)] - if len(finite) > 0: - print(f" Cost range: {finite.min():.1f} - {finite.max():.1f} s") - else: - print(f" All cells impassable (test data too steep)") + print("\n[1] Speed functions at various grades:") + print(f"{'Grade':<10} {'Foot':<12} {'MTB':<12} {'ATV':<12} {'Vehicle':<12}") + print("-" * 60) - print("\nTesting cost grid with friction and trails:") - elev = np.ones((10, 10), dtype=np.float32) * 1000 # flat terrain - friction = np.ones((10, 10), dtype=np.float32) * 2.0 # 2.0x friction (forest) - trails = np.zeros((10, 10), dtype=np.uint8) - trails[5, :] = 5 # road across middle row + for grade_val in [-0.3, -0.1, 0.0, 0.1, 0.2, 0.3]: + grade_arr = np.array([grade_val]) + foot = tobler_off_path_speed(grade_arr, 6.0)[0] + mtb = herzog_wheeled_speed(grade_arr, 12.0)[0] + atv = herzog_wheeled_speed(grade_arr, 25.0)[0] + veh = linear_degrade_speed(grade_arr, 40.0, np.tan(np.radians(20)))[0] + print(f"{grade_val:+.2f} {foot:>6.2f} km/h {mtb:>6.2f} km/h {atv:>6.2f} km/h {veh:>6.2f} km/h") - cost_no_trail = compute_cost_grid(elev, cell_size_m=30.0, friction=friction) - cost_with_trail = compute_cost_grid(elev, cell_size_m=30.0, friction=friction, trails=trails) + print("\n[2] Mode profiles:") + for name, profile in MODE_PROFILES.items(): + print(f"\n {name.upper()}: {profile.description}") + print(f" Max slope: {profile.max_slope_deg}°") + print(f" Trail access: {profile.trail_friction}") + print(f" Wilderness blocked: {profile.wilderness_impassable}") - base_cost = 30 * 3.6 / (0.6 * 6.0 * np.exp(-3.5 * 0.05)) - print(f" Base cost (flat, 30m cell): {base_cost:.1f} s") - print(f" Forest cell (2.0x friction): {cost_no_trail[0, 0]:.1f} s") - print(f" Road cell (0.1x friction, replaces forest): {cost_with_trail[5, 0]:.1f} s") - print(f" Road friction advantage: {cost_no_trail[0, 0] / cost_with_trail[5, 0]:.1f}x faster") - - print("\nTesting cost grid with barriers (three modes):") + print("\n[3] Cost grid test (flat terrain, forest):") elev = np.ones((10, 10), dtype=np.float32) * 1000 - barriers = np.zeros((10, 10), dtype=np.uint8) - barriers[3:7, 3:7] = 255 + friction = np.ones((10, 10), dtype=np.float32) * 2.0 # Forest friction + friction_raw = np.ones((10, 10), dtype=np.uint8) * 10 # Tree cover class - for mode in ["strict", "pragmatic", "emergency"]: - cost = compute_cost_grid(elev, cell_size_m=30.0, barriers=barriers, boundary_mode=mode) + trails = np.zeros((10, 10), dtype=np.uint8) + trails[5, :] = 5 # Road across middle + + for mode_name in ["foot", "mtb", "atv", "vehicle"]: + cost = compute_cost_grid( + elev, cell_size_m=30.0, + friction=friction, + friction_raw=friction_raw, + trails=trails, + mode=mode_name + ) + off_trail_cost = cost[0, 0] + road_cost = cost[5, 0] impassable = np.sum(np.isinf(cost)) - barrier_cost = cost[5, 5] if not np.isinf(cost[5, 5]) else "inf" - print(f" {mode:10s}: {impassable} impassable, barrier cell cost = {barrier_cost}") + print(f" {mode_name:8s}: off-trail={off_trail_cost:>8.1f}s, road={road_cost:>6.1f}s, impassable={impassable}") + + print("\n[4] Wilderness blocking test:") + wilderness = np.zeros((10, 10), dtype=np.uint8) + wilderness[3:7, 3:7] = 255 + + for mode_name in ["foot", "mtb", "atv", "vehicle"]: + cost = compute_cost_grid( + elev, cell_size_m=30.0, + wilderness=wilderness, + mode=mode_name + ) + wilderness_impassable = np.sum(np.isinf(cost[3:7, 3:7])) + print(f" {mode_name:8s}: wilderness cells impassable = {wilderness_impassable}/16") + + print("\nDone.") diff --git a/lib/offroute/router.py b/lib/offroute/router.py index 57d6ce5..d44bbb4 100644 --- a/lib/offroute/router.py +++ b/lib/offroute/router.py @@ -6,6 +6,8 @@ Connects the raster pathfinder (wilderness segment) to Valhalla (on-network segm Entry points are extracted from OSM highways and stored in /mnt/nav/navi.db. The pathfinder routes from a wilderness start to the nearest entry point, then Valhalla completes the route to the destination. + +Supports four travel modes: foot, mtb, atv, vehicle. """ import json import math @@ -23,7 +25,7 @@ from skimage.graph import MCP_Geometric from .dem import DEMReader from .cost import compute_cost_grid from .friction import FrictionReader, friction_to_multiplier -from .barriers import BarrierReader +from .barriers import BarrierReader, WildernessReader, DEFAULT_WILDERNESS_PATH from .trails import TrailReader # Paths @@ -45,6 +47,7 @@ MODE_TO_COSTING = { "foot": "pedestrian", "mtb": "bicycle", "atv": "auto", + "vehicle": "auto", } @@ -120,7 +123,6 @@ class EntryPointIndex: def query_radius(self, lat: float, lon: float, radius_km: float) -> List[Dict]: """Query entry points within radius of a point.""" - # Approximate bbox for the radius lat_delta = radius_km / 111.0 lon_delta = radius_km / (111.0 * math.cos(math.radians(lat))) @@ -129,7 +131,6 @@ class EntryPointIndex: lon - lon_delta, lon + lon_delta ) - # Filter by actual distance and add distance field result = [] for p in points: dist = haversine_distance(lat, lon, p['lat'], p['lon']) @@ -140,17 +141,12 @@ class EntryPointIndex: return sorted(result, key=lambda x: x['distance_m']) def build_index(self, osm_pbf_path: Path = OSM_PBF_PATH) -> Dict: - """ - Build the entry point index from OSM PBF. - - Extracts endpoints of highway features that connect to the network. - """ + """Build the entry point index from OSM PBF.""" if not osm_pbf_path.exists(): raise FileNotFoundError(f"OSM PBF not found: {osm_pbf_path}") print(f"Building trail entry point index from {osm_pbf_path}...") - # Highway types to extract (routable network entry points) highway_types = [ "primary", "secondary", "tertiary", "unclassified", "residential", "service", "track", "path", "footway", "bridleway" @@ -159,42 +155,29 @@ class EntryPointIndex: stats = {"total": 0, "by_class": {}} with tempfile.TemporaryDirectory() as tmpdir: - # Extract highways to GeoJSON geojson_path = Path(tmpdir) / "highways.geojson" - # Build osmium tags-filter expressions (one per highway type) print(f" Extracting highways with osmium...") - cmd = [ - "osmium", "tags-filter", - str(osm_pbf_path), - ] - # Add each highway type as a separate filter expression + cmd = ["osmium", "tags-filter", str(osm_pbf_path)] for ht in highway_types: cmd.append(f"w/highway={ht}") cmd.extend(["-o", str(Path(tmpdir) / "filtered.osm.pbf"), "--overwrite"]) - subprocess.run(cmd, check=True, capture_output=True) - # Convert to GeoJSON print(f" Converting to GeoJSON with ogr2ogr...") cmd = [ "ogr2ogr", "-f", "GeoJSON", str(geojson_path), str(Path(tmpdir) / "filtered.osm.pbf"), - "lines", - "-t_srs", "EPSG:4326" + "lines", "-t_srs", "EPSG:4326" ] subprocess.run(cmd, check=True, capture_output=True) - # Parse GeoJSON and extract endpoints print(f" Extracting entry points...") with open(geojson_path) as f: data = json.load(f) - # Collect unique points (endpoints) - # Key: (lat, lon) rounded to 5 decimal places (~1m precision) points = {} - for feature in data.get("features", []): props = feature.get("properties", {}) geom = feature.get("geometry", {}) @@ -209,62 +192,48 @@ class EntryPointIndex: highway_class = props.get("highway", "unknown") name = props.get("name", "") - # Extract endpoints for coord in [coords[0], coords[-1]]: lon, lat = coord[0], coord[1] key = (round(lat, 5), round(lon, 5)) if key not in points: points[key] = { - "lat": lat, - "lon": lon, - "highway_class": highway_class, - "name": name + "lat": lat, "lon": lon, + "highway_class": highway_class, "name": name } else: - # Keep the "best" highway class (roads > tracks > paths) existing = points[key] if self._highway_priority(highway_class) < self._highway_priority(existing["highway_class"]): points[key]["highway_class"] = highway_class if name and not existing["name"]: points[key]["name"] = name - # Create/update database print(f" Writing {len(points)} entry points to {self.db_path}...") self.db_path.parent.mkdir(parents=True, exist_ok=True) conn = self._get_conn() - # Create table conn.execute(""" CREATE TABLE IF NOT EXISTS trail_entry_points ( id INTEGER PRIMARY KEY AUTOINCREMENT, - lat REAL NOT NULL, - lon REAL NOT NULL, - highway_class TEXT NOT NULL, - name TEXT + lat REAL NOT NULL, lon REAL NOT NULL, + highway_class TEXT NOT NULL, name TEXT ) """) - - # Clear existing data conn.execute("DELETE FROM trail_entry_points") - # Insert new points for point in points.values(): - conn.execute(""" - INSERT INTO trail_entry_points (lat, lon, highway_class, name) - VALUES (?, ?, ?, ?) - """, (point["lat"], point["lon"], point["highway_class"], point["name"])) - + conn.execute( + "INSERT INTO trail_entry_points (lat, lon, highway_class, name) VALUES (?, ?, ?, ?)", + (point["lat"], point["lon"], point["highway_class"], point["name"]) + ) stats["total"] += 1 hc = point["highway_class"] stats["by_class"][hc] = stats["by_class"].get(hc, 0) + 1 - # Create spatial index conn.execute("CREATE INDEX IF NOT EXISTS idx_entry_lat ON trail_entry_points(lat)") conn.execute("CREATE INDEX IF NOT EXISTS idx_entry_lon ON trail_entry_points(lon)") conn.execute("CREATE INDEX IF NOT EXISTS idx_entry_latlon ON trail_entry_points(lat, lon)") - conn.commit() print(f" Done. Total: {stats['total']} entry points") @@ -291,12 +260,15 @@ class EntryPointIndex: class OffrouteRouter: """ OFFROUTE Router — orchestrates wilderness pathfinding and Valhalla stitching. + + Supports modes: foot, mtb, atv, vehicle """ def __init__(self): self.dem_reader = None self.friction_reader = None self.barrier_reader = None + self.wilderness_reader = None self.trail_reader = None self.entry_index = EntryPointIndex() @@ -308,6 +280,8 @@ class OffrouteRouter: self.friction_reader = FrictionReader() if self.barrier_reader is None: self.barrier_reader = BarrierReader() + if self.wilderness_reader is None and DEFAULT_WILDERNESS_PATH.exists(): + self.wilderness_reader = WildernessReader() if self.trail_reader is None: self.trail_reader = TrailReader() @@ -317,16 +291,25 @@ class OffrouteRouter: start_lon: float, end_lat: float, end_lon: float, - mode: Literal["foot", "mtb", "atv"] = "foot", + mode: Literal["foot", "mtb", "atv", "vehicle"] = "foot", boundary_mode: Literal["strict", "pragmatic", "emergency"] = "pragmatic" ) -> Dict: """ Route from a wilderness start point to a destination. + Args: + start_lat, start_lon: Starting coordinates (wilderness) + end_lat, end_lon: Destination coordinates + mode: Travel mode (foot, mtb, atv, vehicle) + boundary_mode: How to handle private land (strict, pragmatic, emergency) + Returns a GeoJSON FeatureCollection with wilderness and network segments. """ t0 = time.time() + if mode not in MODE_TO_COSTING: + return {"status": "error", "message": f"Unknown mode: {mode}"} + # Ensure entry point index exists if not self.entry_index.table_exists() or self.entry_index.get_entry_point_count() == 0: return { @@ -334,28 +317,27 @@ class OffrouteRouter: "message": "Trail entry point index not built. Run build_entry_index() first." } - # Find entry points near start - entry_points = self.entry_index.query_radius( - start_lat, start_lon, DEFAULT_SEARCH_RADIUS_KM - ) + # Find entry points near start (limit to nearest 10 to control bbox size) + MAX_ENTRY_POINTS = 10 + entry_points = self.entry_index.query_radius(start_lat, start_lon, DEFAULT_SEARCH_RADIUS_KM) if not entry_points: - # Try expanded radius - entry_points = self.entry_index.query_radius( - start_lat, start_lon, EXPANDED_SEARCH_RADIUS_KM - ) + entry_points = self.entry_index.query_radius(start_lat, start_lon, EXPANDED_SEARCH_RADIUS_KM) if not entry_points: return { "status": "error", "message": f"No trail entry points found within {EXPANDED_SEARCH_RADIUS_KM}km of start" } - # Build bbox for pathfinding grid - # Include start, end, and all entry points + # Limit to nearest entry points to prevent huge bounding boxes + entry_points = entry_points[:MAX_ENTRY_POINTS] + + # Build bbox with max size limit (prevent OOM on large areas) + MAX_BBOX_DEGREES = 0.5 # ~55km at mid-latitudes all_lats = [start_lat, end_lat] + [p["lat"] for p in entry_points] all_lons = [start_lon, end_lon] + [p["lon"] for p in entry_points] - padding = 0.05 # ~5km padding + padding = 0.05 bbox = { "south": min(all_lats) - padding, "north": max(all_lats) + padding, @@ -363,16 +345,28 @@ class OffrouteRouter: "east": max(all_lons) + padding, } + # Clamp bbox size to prevent memory exhaustion + lat_span = bbox["north"] - bbox["south"] + lon_span = bbox["east"] - bbox["west"] + if lat_span > MAX_BBOX_DEGREES or lon_span > MAX_BBOX_DEGREES: + center_lat = (bbox["south"] + bbox["north"]) / 2 + center_lon = (bbox["west"] + bbox["east"]) / 2 + half_span = MAX_BBOX_DEGREES / 2 + bbox = { + "south": center_lat - half_span, + "north": center_lat + half_span, + "west": center_lon - half_span, + "east": center_lon + half_span, + } + # Initialize readers self._init_readers() # Load elevation try: elevation, meta = self.dem_reader.get_elevation_grid( - south=bbox["south"], - north=bbox["north"], - west=bbox["west"], - east=bbox["east"], + south=bbox["south"], north=bbox["north"], + west=bbox["west"], east=bbox["east"], ) except Exception as e: return {"status": "error", "message": f"Failed to load elevation: {e}"} @@ -382,62 +376,69 @@ class OffrouteRouter: if mem > MEMORY_LIMIT_GB: return {"status": "error", "message": f"Memory limit exceeded: {mem:.1f}GB > {MEMORY_LIMIT_GB}GB"} - # Load friction + # Load friction (both processed and raw for mode-specific overrides) friction_raw = self.friction_reader.get_friction_grid( - south=bbox["south"], - north=bbox["north"], - west=bbox["west"], - east=bbox["east"], + south=bbox["south"], north=bbox["north"], + west=bbox["west"], east=bbox["east"], target_shape=elevation.shape ) friction_mult = friction_to_multiplier(friction_raw) # Load barriers barriers = self.barrier_reader.get_barrier_grid( - south=bbox["south"], - north=bbox["north"], - west=bbox["west"], - east=bbox["east"], + south=bbox["south"], north=bbox["north"], + west=bbox["west"], east=bbox["east"], target_shape=elevation.shape ) + # Load wilderness (if available and mode requires it) + wilderness = None + if self.wilderness_reader is not None and mode in ("mtb", "atv", "vehicle"): + wilderness = self.wilderness_reader.get_wilderness_grid( + south=bbox["south"], north=bbox["north"], + west=bbox["west"], east=bbox["east"], + target_shape=elevation.shape + ) + # Load trails trails = self.trail_reader.get_trails_grid( - south=bbox["south"], - north=bbox["north"], - west=bbox["west"], - east=bbox["east"], + south=bbox["south"], north=bbox["north"], + west=bbox["west"], east=bbox["east"], target_shape=elevation.shape ) - # Compute cost grid + # Compute cost grid with mode-specific parameters cost = compute_cost_grid( elevation, cell_size_m=meta["cell_size_m"], friction=friction_mult, + friction_raw=friction_raw, trails=trails, barriers=barriers, + wilderness=wilderness, boundary_mode=boundary_mode, + mode=mode, ) + # Free intermediate arrays to reduce memory before MCP + # Note: Keep trails and barriers - needed for path statistics + del friction_mult, friction_raw, wilderness + import gc + gc.collect() + # Convert start to pixel coordinates start_row, start_col = self.dem_reader.latlon_to_pixel(start_lat, start_lon, meta) - # Validate start is in bounds rows, cols = elevation.shape if not (0 <= start_row < rows and 0 <= start_col < cols): return {"status": "error", "message": "Start point outside grid bounds"} - # Mark entry points on the grid + # Mark entry points on grid entry_pixels = [] for ep in entry_points: row, col = self.dem_reader.latlon_to_pixel(ep["lat"], ep["lon"], meta) if 0 <= row < rows and 0 <= col < cols: - entry_pixels.append({ - "row": row, - "col": col, - "entry_point": ep - }) + entry_pixels.append({"row": row, "col": col, "entry_point": ep}) if not entry_pixels: return {"status": "error", "message": "No entry points map to grid bounds"} @@ -465,18 +466,21 @@ class OffrouteRouter: # Traceback wilderness path path_indices = mcp.traceback((best_entry["row"], best_entry["col"])) - # Convert to coordinates and collect stats + # Convert to coordinates wilderness_coords = [] elevations = [] trail_values = [] + barrier_crossings = 0 for row, col in path_indices: lat, lon = self.dem_reader.pixel_to_latlon(row, col, meta) wilderness_coords.append([lon, lat]) elevations.append(elevation[row, col]) trail_values.append(trails[row, col]) + if barriers[row, col] == 255: + barrier_crossings += 1 - # Calculate wilderness segment stats + # Calculate stats wilderness_distance_m = 0 for i in range(1, len(wilderness_coords)): lon1, lat1 = wilderness_coords[i-1] @@ -493,13 +497,16 @@ class OffrouteRouter: total_cells = len(trail_arr) on_trail_pct = float(100 * on_trail_cells / total_cells) if total_cells > 0 else 0 - # Entry point reached + # Free trails and barriers now that path stats are computed + del trails, barriers + + # Entry point entry_lat = best_entry["entry_point"]["lat"] entry_lon = best_entry["entry_point"]["lon"] entry_class = best_entry["entry_point"]["highway_class"] entry_name = best_entry["entry_point"].get("name", "") - # Call Valhalla for on-network segment + # Call Valhalla valhalla_costing = MODE_TO_COSTING.get(mode, "pedestrian") valhalla_request = { @@ -515,11 +522,7 @@ class OffrouteRouter: valhalla_error = None try: - resp = requests.post( - f"{VALHALLA_URL}/route", - json=valhalla_request, - timeout=30 - ) + resp = requests.post(f"{VALHALLA_URL}/route", json=valhalla_request, timeout=30) if resp.status_code == 200: valhalla_data = resp.json() @@ -529,11 +532,8 @@ class OffrouteRouter: if legs: leg = legs[0] shape = leg.get("shape", "") - - # Decode polyline6 network_coords = self._decode_polyline(shape) - # Extract maneuvers maneuvers = [] for m in leg.get("maneuvers", []): maneuvers.append({ @@ -560,7 +560,6 @@ class OffrouteRouter: # Build response features = [] - # Feature 1: Wilderness segment wilderness_feature = { "type": "Feature", "properties": { @@ -572,15 +571,13 @@ class OffrouteRouter: "boundary_mode": boundary_mode, "on_trail_pct": on_trail_pct, "cell_count": total_cells, + "barrier_crossings": barrier_crossings, + "mode": mode, }, - "geometry": { - "type": "LineString", - "coordinates": wilderness_coords, - } + "geometry": {"type": "LineString", "coordinates": wilderness_coords} } features.append(wilderness_feature) - # Feature 2: Network segment (if available) if network_segment: network_feature = { "type": "Feature", @@ -590,40 +587,23 @@ class OffrouteRouter: "duration_minutes": network_segment["duration_minutes"], "maneuvers": network_segment["maneuvers"], }, - "geometry": { - "type": "LineString", - "coordinates": network_segment["coordinates"], - } + "geometry": {"type": "LineString", "coordinates": network_segment["coordinates"]} } features.append(network_feature) - # Build combined route coordinates combined_coords = wilderness_coords.copy() if network_segment: - # Skip first point of network segment (it's the same as last wilderness point) combined_coords.extend(network_segment["coordinates"][1:]) - # Feature 3: Combined route combined_feature = { "type": "Feature", - "properties": { - "segment_type": "combined", - "mode": mode, - "boundary_mode": boundary_mode, - }, - "geometry": { - "type": "LineString", - "coordinates": combined_coords, - } + "properties": {"segment_type": "combined", "mode": mode, "boundary_mode": boundary_mode}, + "geometry": {"type": "LineString", "coordinates": combined_coords} } features.append(combined_feature) - geojson = { - "type": "FeatureCollection", - "features": features, - } + geojson = {"type": "FeatureCollection", "features": features} - # Build summary total_distance_km = wilderness_distance_m / 1000 total_effort_minutes = best_cost / 60 @@ -639,22 +619,17 @@ class OffrouteRouter: "network_distance_km": float(network_segment["distance_km"]) if network_segment else 0, "network_duration_minutes": float(network_segment["duration_minutes"]) if network_segment else 0, "on_trail_pct": on_trail_pct, + "barrier_crossings": barrier_crossings, "boundary_mode": boundary_mode, "mode": mode, "entry_point": { - "lat": entry_lat, - "lon": entry_lon, - "highway_class": entry_class, - "name": entry_name, + "lat": entry_lat, "lon": entry_lon, + "highway_class": entry_class, "name": entry_name, }, "computation_time_s": time.time() - t0, } - result = { - "status": "ok", - "route": geojson, - "summary": summary, - } + result = {"status": "ok", "route": geojson, "summary": summary} if valhalla_error: result["warning"] = f"Network segment incomplete: {valhalla_error}" @@ -669,7 +644,6 @@ class OffrouteRouter: lon = 0 while index < len(encoded): - # Latitude shift = 0 result = 0 while True: @@ -682,7 +656,6 @@ class OffrouteRouter: dlat = ~(result >> 1) if result & 1 else result >> 1 lat += dlat - # Longitude shift = 0 result = 0 while True: @@ -707,6 +680,8 @@ class OffrouteRouter: self.friction_reader.close() if self.barrier_reader: self.barrier_reader.close() + if self.wilderness_reader: + self.wilderness_reader.close() if self.trail_reader: self.trail_reader.close() self.entry_index.close() @@ -729,24 +704,33 @@ if __name__ == "__main__": print(f"\nDone. Total entry points: {stats['total']}") elif len(sys.argv) > 1 and sys.argv[1] == "test": - print("Testing router...") + print("Testing router (all modes)...") router = OffrouteRouter() - # Test route: wilderness to Twin Falls - result = router.route( - start_lat=42.35, - start_lon=-114.30, - end_lat=42.5629, - end_lon=-114.4609, - mode="foot", - boundary_mode="pragmatic" - ) + for mode in ["foot", "mtb", "atv", "vehicle"]: + print(f"\n{'='*60}") + print(f"Mode: {mode}") + print("="*60) + + result = router.route( + start_lat=42.35, start_lon=-114.30, + end_lat=42.5629, end_lon=-114.4609, + mode=mode, boundary_mode="pragmatic" + ) + + if result["status"] == "ok": + s = result["summary"] + print(f" Wilderness: {s['wilderness_distance_km']:.2f} km, {s['wilderness_effort_minutes']:.1f} min") + print(f" Network: {s['network_distance_km']:.2f} km, {s['network_duration_minutes']:.1f} min") + print(f" On-trail: {s['on_trail_pct']:.1f}%") + print(f" Entry: {s['entry_point']['highway_class']}") + else: + print(f" ERROR: {result['message']}") - print(json.dumps(result, indent=2, default=str)) router.close() else: print("Usage:") print(" python router.py build # Build entry point index") - print(" python router.py test # Test route") + print(" python router.py test # Test all modes")