feat(offroute): Phase O4 — multi-mode cost functions (foot/mtb/atv/vehicle)

- Add ModeProfile dataclass for data-driven mode configuration
- Implement three speed functions:
  * Tobler off-path hiking (foot)
  * Herzog wheeled-transport polynomial (mtb/atv)
  * Linear speed degradation (vehicle)
- Add WildernessReader for PAD-US Des_Tp=WA wilderness areas
- Mode-specific terrain friction overrides:
  * Forest impassable for ATV/vehicle, high friction for MTB
  * Wetland/mangrove impassable for all wheeled modes
- Trail access rules:
  * Foot trails (value 25) impassable for ATV/vehicle
- Wilderness blocking for mtb/atv/vehicle modes
- Vehicle mode allows flat grassland/cropland traversal
- Memory optimization: limit entry points, constrain bbox size
- Update router to pass mode and wilderness to cost function
- Add vehicle to API mode validation

Validated all four modes with test route:
- foot: 0.46km off-network, 12.11km network, 89% on trail
- mtb: 0.47km off-network, 13.13km network, 90% on trail
- atv: 0.47km off-network, 12.81km network, 90% on trail
- vehicle: 0.46km off-network, 12.81km network, 89% on trail

Co-Authored-By: Claude Opus 4.5 <noreply@anthropic.com>

lib/api.py

@@ -2768,8 +2768,8 @@ def api_offroute():
 
         # Parse options
         mode = data.get("mode", "foot")
-        if mode not in ("foot", "mtb", "atv"):
+        if mode not in ("foot", "mtb", "atv", "vehicle"):
-            return jsonify({"status": "error", "message": "mode must be foot, mtb, or atv"}), 400
+            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"):

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.
+Provides access to:
-Cells with value 255 represent closed/restricted areas (Pub_Access = XA).
+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.
+Build functions rasterize PAD-US geodatabase to aligned GeoTIFFs.
-Runtime functions read the raster and resample to match elevation grids.
+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":
+    if len(sys.argv) > 1:
-        # Build the raster
+        cmd = sys.argv[1]
-        print("="*60)
+
-        print("PAD-US Barriers Raster Build")
+        if cmd == "build":
-        print("="*60)
+            print("=" * 60)
-        build_barriers_raster()
+            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.")

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
+Supports four travel modes: foot, mtb, atv, vehicle.
-hiking function with off-trail penalty. Optionally applies friction
+Each mode has its own speed function, max slope, trail access rules,
-multipliers from land cover data, trail corridors, and barrier grids.
+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)
+                  0 = no trail, 5 = road, 15 = track, 25 = foot trail
-                  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.
         barriers: Optional 2D array of barrier values (uint8).
-                  255 = closed/restricted area (from PAD-US Pub_Access = XA).
+                  255 = closed/restricted area (PAD-US Pub_Access = XA).
-                  0 = accessible.
+        wilderness: Optional[np.ndarray] of wilderness values (uint8).
-                  If None, no barriers are applied.
+                  255 = designated wilderness area.
-        boundary_mode: How to handle private/restricted land barriers:
+        boundary_mode: How to handle barriers ("strict", "pragmatic", "emergency")
-                  "strict"    - cells with barrier=255 become impassable (np.inf)
+        mode: Travel mode ("foot", "mtb", "atv", "vehicle")
-                  "pragmatic" - cells with barrier=255 get 5.0x friction penalty
-                  "emergency" - barriers are ignored entirely
-                  Default: "pragmatic"
 
     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
+    # ─── Compute gradients (in-place where possible) ─────────────────────────
-    dy = np.zeros_like(elevation)
+    # Use float32 to reduce memory footprint
-    dx = np.zeros_like(elevation)
+    grade = np.zeros(elevation.shape, dtype=np.float32)
 
-    # Central differences for interior, forward/backward at edges
+    # Compute dy contribution to grade squared
-    dy[1:-1, :] = (elevation[:-2, :] - elevation[2:, :]) / (2 * cell_size_lat_m)
+    dy_contrib = np.zeros(elevation.shape, dtype=np.float32)
-    dy[0, :] = (elevation[0, :] - elevation[1, :]) / cell_size_lat_m
+    dy_contrib[1:-1, :] = ((elevation[:-2, :] - elevation[2:, :]) / (2 * cell_size_lat_m)) ** 2
-    dy[-1, :] = (elevation[-2, :] - elevation[-1, :]) / cell_size_lat_m
+    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)
+    # Compute dx contribution and add to dy_contrib in-place
-    dx[:, 0] = (elevation[:, 1] - elevation[:, 0]) / cell_size_lon_m
+    dy_contrib[:, 1:-1] += ((elevation[:, 2:] - elevation[:, :-2]) / (2 * cell_size_lon_m)) ** 2
-    dx[:, -1] = (elevation[:, -1] - elevation[:, -2]) / cell_size_lon_m
+    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 = sqrt(dx^2 + dy^2)
-    grade_magnitude = np.sqrt(dx**2 + dy**2)
+    np.sqrt(dy_contrib, out=grade)
+    del dy_contrib  # Free memory immediately
 
-    # Convert to slope angle in degrees
+    # ─── Compute speed based on mode ─────────────────────────────────────────
-    slope_deg = np.degrees(np.arctan(grade_magnitude))
+    max_grade_val = np.tan(np.radians(profile.max_slope_deg))
 
-    # Compute speed for each cell using Tobler function
+    if profile.speed_function == "tobler":
-    speed_kmh = TOBLER_OFF_TRAIL_MULT * TOBLER_BASE_SPEED * np.exp(-3.5 * np.abs(grade_magnitude + 0.05))
+        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
+    # ─── Max slope limit ──────────────────────────────────────────────────────
-    cost[slope_deg > MAX_SLOPE_DEG] = np.inf
+    cost[grade > max_grade_val] = np.inf
 
-    # Handle NaN elevations (no data)
+    # ─── NaN elevations ──────────────────────────────────────────────────────
     cost[np.isnan(elevation)] = np.inf
 
-    # Build effective friction array
+    # ─── Apply friction in-place ─────────────────────────────────────────────
-    # Start with WorldCover friction if provided, else 1.0
+    # 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(
+            raise ValueError(f"Friction shape mismatch")
-                f"Friction shape {friction.shape} does not match elevation shape {elevation.shape}"
+        np.multiply(cost, friction, out=cost)
-            )
-        effective_friction = friction.copy()
-    else:
-        effective_friction = np.ones(elevation.shape, dtype=np.float32)
 
-    # 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(
+            raise ValueError(f"Trails shape mismatch")
-                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
 
-    # Apply friction to cost
+        for trail_value, trail_friction in profile.trail_friction.items():
-    cost = cost * effective_friction
+            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(
+            raise ValueError(f"Barriers shape mismatch")
-                f"Barriers shape {barriers.shape} does not match elevation shape {elevation.shape}"
-            )
-
-        barrier_mask = barriers == 255
 
         if boundary_mode == "strict":
-            # Mark closed/restricted areas as impassable
+            np.putmask(cost, barriers == 255, np.inf)
-            cost[barrier_mask] = np.inf
         elif boundary_mode == "pragmatic":
-            # Apply friction penalty to closed/restricted areas
+            barrier_mask = barriers == 255
-            cost[barrier_mask] = cost[barrier_mask] * PRAGMATIC_BARRIER_MULTIPLIER
+            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:")
+    print("=" * 70)
-    for grade in [-0.3, -0.1, -0.05, 0.0, 0.05, 0.1, 0.3]:
+    print("OFFROUTE Multi-Mode Cost Function Tests")
-        speed = tobler_speed(grade)
+    print("=" * 70)
-        print(f"  Grade {grade:+.2f}: {speed:.2f} km/h")
 
-    print("\nTesting cost grid computation (no friction, no trails):")
+    print("\n[1] Speed functions at various grades:")
-    elev = np.arange(100).reshape(10, 10).astype(np.float32) * 10
+    print(f"{'Grade':<10} {'Foot':<12} {'MTB':<12} {'ATV':<12} {'Vehicle':<12}")
-    cost = compute_cost_grid(elev, cell_size_m=30.0)
+    print("-" * 60)
-    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("\nTesting cost grid with friction and trails:")
+    for grade_val in [-0.3, -0.1, 0.0, 0.1, 0.2, 0.3]:
-    elev = np.ones((10, 10), dtype=np.float32) * 1000  # flat terrain
+        grade_arr = np.array([grade_val])
-    friction = np.ones((10, 10), dtype=np.float32) * 2.0  # 2.0x friction (forest)
+        foot = tobler_off_path_speed(grade_arr, 6.0)[0]
-    trails = np.zeros((10, 10), dtype=np.uint8)
+        mtb = herzog_wheeled_speed(grade_arr, 12.0)[0]
-    trails[5, :] = 5  # road across middle row
+        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)
+    print("\n[2] Mode profiles:")
-    cost_with_trail = compute_cost_grid(elev, cell_size_m=30.0, friction=friction, trails=trails)
+    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("\n[3] Cost grid test (flat terrain, forest):")
-    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):")
     elev = np.ones((10, 10), dtype=np.float32) * 1000
-    barriers = np.zeros((10, 10), dtype=np.uint8)
+    friction = np.ones((10, 10), dtype=np.float32) * 2.0  # Forest friction
-    barriers[3:7, 3:7] = 255
+    friction_raw = np.ones((10, 10), dtype=np.uint8) * 10  # Tree cover class
 
-    for mode in ["strict", "pragmatic", "emergency"]:
+    trails = np.zeros((10, 10), dtype=np.uint8)
-        cost = compute_cost_grid(elev, cell_size_m=30.0, barriers=barriers, boundary_mode=mode)
+    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_name:8s}: off-trail={off_trail_cost:>8.1f}s, road={road_cost:>6.1f}s, impassable={impassable}")
-        print(f"  {mode:10s}: {impassable} impassable, barrier cell cost = {barrier_cost}")
+
+    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.")

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."""
-        Build the entry point index from OSM PBF.
-
-        Extracts endpoints of highway features that connect to the network.
-        """
         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 = [
+            cmd = ["osmium", "tags-filter", str(osm_pbf_path)]
-                "osmium", "tags-filter",
-                str(osm_pbf_path),
-            ]
-            # Add each highway type as a separate filter expression
             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",
+                "lines", "-t_srs", "EPSG:4326"
-                "-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,
+                            "lat": lat, "lon": lon,
-                            "lon": lon,
+                            "highway_class": highway_class, "name": name
-                            "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,
+                lat REAL NOT NULL, lon REAL NOT NULL,
-                lon REAL NOT NULL,
+                highway_class TEXT NOT NULL, name TEXT
-                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("""
+            conn.execute(
-                INSERT INTO trail_entry_points (lat, lon, highway_class, name)
+                "INSERT INTO trail_entry_points (lat, lon, highway_class, name) VALUES (?, ?, ?, ?)",
-                VALUES (?, ?, ?, ?)
+                (point["lat"], point["lon"], point["highway_class"], point["name"])
-            """, (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
+        # Find entry points near start (limit to nearest 10 to control bbox size)
-        entry_points = self.entry_index.query_radius(
+        MAX_ENTRY_POINTS = 10
-            start_lat, start_lon, DEFAULT_SEARCH_RADIUS_KM
+        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
+        # Limit to nearest entry points to prevent huge bounding boxes
-        # Include start, end, and all entry points
+        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"],
+                south=bbox["south"], north=bbox["north"],
-                north=bbox["north"],
+                west=bbox["west"], east=bbox["east"],
-                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"],
+            south=bbox["south"], north=bbox["north"],
-            north=bbox["north"],
+            west=bbox["west"], east=bbox["east"],
-            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"],
+            south=bbox["south"], north=bbox["north"],
-            north=bbox["north"],
+            west=bbox["west"], east=bbox["east"],
-            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"],
+            south=bbox["south"], north=bbox["north"],
-            north=bbox["north"],
+            west=bbox["west"], east=bbox["east"],
-            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({
+                entry_pixels.append({"row": row, "col": col, "entry_point": ep})
-                    "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(
+            resp = requests.post(f"{VALHALLA_URL}/route", json=valhalla_request, timeout=30)
-                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": {
+            "geometry": {"type": "LineString", "coordinates": wilderness_coords}
-                "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": {
+                "geometry": {"type": "LineString", "coordinates": network_segment["coordinates"]}
-                    "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": {
+            "properties": {"segment_type": "combined", "mode": mode, "boundary_mode": boundary_mode},
-                "segment_type": "combined",
+            "geometry": {"type": "LineString", "coordinates": combined_coords}
-                "mode": mode,
-                "boundary_mode": boundary_mode,
-            },
-            "geometry": {
-                "type": "LineString",
-                "coordinates": combined_coords,
-            }
         }
         features.append(combined_feature)
 
-        geojson = {
+        geojson = {"type": "FeatureCollection", "features": features}
-            "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,
+                "lat": entry_lat, "lon": entry_lon,
-                "lon": entry_lon,
+                "highway_class": entry_class, "name": entry_name,
-                "highway_class": entry_class,
-                "name": entry_name,
             },
             "computation_time_s": time.time() - t0,
         }
 
-        result = {
+        result = {"status": "ok", "route": geojson, "summary": summary}
-            "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
+        for mode in ["foot", "mtb", "atv", "vehicle"]:
-        result = router.route(
+            print(f"\n{'='*60}")
-            start_lat=42.35,
+            print(f"Mode: {mode}")
-            start_lon=-114.30,
+            print("="*60)
-            end_lat=42.5629,
+
-            end_lon=-114.4609,
+            result = router.route(
-            mode="foot",
+                start_lat=42.35, start_lon=-114.30,
-            boundary_mode="pragmatic"
+                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")