Add navi-offroute service (extraction #8 — final) (#10)

* Add navi-offroute service (extraction #8 — the last one)

Faithful port of recon's /api/offroute (POST) + /api/mvum (GET) and the
runtime offroute modules into a new :8428 service. Closes the loop: after
this, navi-frontend talks only to navi-backend.

Ported: router.py (OffrouteRouter, EntryPointIndex, 4 route strategies,
in-Python MCP_Geometric least-cost path, Valhalla integration, per-request
osmium extract), mvum.py (MVUMReader over navi.db), cost.py, friction.py,
trails.py, and barriers.py (runtime BarrierReader/WildernessReader only).

NOT ported (per Phase A §3/§15): prototype.py (dead at runtime), barriers.py
build_*_raster (offline GDB→raster prep). DEM imported from shared/dem.py
(PR #9), not duplicated.

Behaviour-faithful changes: hardcoded paths/URLs → env vars; the
profile.offroute.* config (osm_pbf_path/postgis_dsn/densify_interval_m) →
dedicated env vars (router drops deployment_config). Both routes public (no
auth, matching recon). PADUS via libpq peer-auth DSN (dbname=padus) — NO
secret. Owns no DB.

15 hermetic tests (offroute validation + mocked-router shape + close-always;
fixture-SQLite MVUM roads/trails/fallback/null; admin auth + no-secrets +
probe shape). Full suite 119 passed / 1 skipped. Adds scikit-image + rasterio.

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

* navi-offroute: PR #10 review cleanups (4 faithful-port deviations)

1. trails.py — drop recon-era "Run the Phase B rasterization script"
   reference from the not-found error (confusing in navi-offroute context).
2. friction.py — add FileNotFoundError-before-rasterio-open check to
   match barriers/trails consistency.
3. mvum.py — remove dead try/except shapely import + warnings.warn at
   2 sites (shapely is a hard pyproject dep; the fallback was unreachable).
4. router.py — declare psutil in pyproject, drop the silent fallback;
   the MEMORY_LIMIT_GB safety check was silently disabled in prod.

Adds test_friction_reader_raises_file_not_found_when_missing (16 navi-offroute
tests; full suite 120 passed / 1 skipped).

Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

---------

Co-authored-by: zvx-echo6 <mj@k7zvx.com>
Co-authored-by: Claude Opus 4.7 (1M context) <noreply@anthropic.com>
This commit is contained in:
malice 2026-05-22 23:30:43 -06:00 committed by GitHub
commit ae82cee46a
17 changed files with 3960 additions and 0 deletions

View file

@ -88,6 +88,24 @@ All `@require_auth`. The per-service admin endpoints stay localhost-only; this i
the single edge-exposed admin surface (needs a Caddy `@authed_api` edit — see
`deploy/caddy/navi-admin.caddy.notes.md`).
## Run (local) — navi-offroute (extraction #8)
```bash
.venv/bin/pytest services/navi_offroute/tests/ -v
# All paths/URLs env-overridable (deploy/env/navi-offroute.env.example).
# No secrets — PADUS via libpq peer-auth (dbname=padus). DEM via shared/dem.py.
# Needs osmium-tool on the host + scikit-image/rasterio in the venv.
.venv/bin/gunicorn 'services.navi_offroute.app:create_app()' \
--bind 127.0.0.1:8428 --workers 2 --timeout 130
```
`navi-offroute` serves `POST /api/offroute` (off-network effort-based routing —
in-Python least-cost path over a DEM/friction/barriers/trails/MVUM cost grid,
stitched to the road network via Valhalla) and `GET /api/mvum` (Motor Vehicle
Use Map road/trail access lookup). Both public. The `^~ /api/offroute` nginx
block needs a long `proxy_read_timeout` (130s); routes can take ~2 min.
## The admin-info convention (§4.5)
Every service exposes `GET /api/admin/<service-name>/info`, gated by `require_auth`,

View file

@ -0,0 +1,30 @@
# navi-offroute — /etc/navi-backend/navi-offroute.env
# Faithful port of recon's /api/offroute + /api/mvum (extraction #8).
# NO SECRETS — PADUS uses libpq peer-auth (dbname=padus, no password); all other
# inputs are read-only paths/URLs. Owns no DB.
# Shared planet-DEM (via shared/dem.py — same file/var as navi-geo).
NAVI_DEM_PMTILES=/mnt/nas/nav/planet-dem.pmtiles
# Offroute-specific read-only data (all stay external per data-ownership rule).
NAVI_OFFROUTE_OSM_PBF=/mnt/nav/sources/idaho-latest.osm.pbf
NAVI_OFFROUTE_NAVI_DB=/mnt/nav/navi.db
NAVI_OFFROUTE_BARRIERS_TIF=/mnt/nav/worldcover/padus_barriers.tif
NAVI_OFFROUTE_WILDERNESS_TIF=/mnt/nav/worldcover/wilderness.tif
NAVI_OFFROUTE_TRAILS_TIF=/mnt/nav/worldcover/trails.tif
NAVI_OFFROUTE_FRICTION_VRT=/mnt/nav/worldcover/friction/friction_conus.vrt
# PADUS PostGIS — peer-auth DSN, NO password. libpq connects via the local
# socket as the service user (systemd User=zvx; verified zvx has padus access).
# Queries the `entry_points` routing index (a different table than navi-landclass).
NAVI_OFFROUTE_POSTGIS_DSN=dbname=padus
# Valhalla — recon-side network router (HTTP), for the on-network leg.
NAVI_OFFROUTE_VALHALLA_URL=http://localhost:8002
# Cost-grid densification interval, metres (was profile.offroute.densify_interval_m).
NAVI_OFFROUTE_DENSIFY_M=100
# HOST DEP: router.py shells out to `osmium extract` per route — the host needs
# osmium-tool installed (`sudo apt-get install osmium-tool`). admin-info reports
# its version (or 'not installed').

View file

@ -0,0 +1,33 @@
# =============================================================================
# navi-offroute — nginx integration for the navi.echo6.co vhost
#
# TWO blocks. Add INSIDE the existing
# server { server_name navi.echo6.co; ... }
# block, BEFORE the existing `location /api/ { ... }` catch-all.
#
# Both routes are PUBLIC (no forward_auth) — @public_api in Caddy already covers
# them, so NO Caddy edit (extraction #8 / Phase A §12).
#
# KEY DELTA from other navi-* blocks: `^~ /api/offroute` needs a LONG
# proxy_read_timeout (130s). Off-network routing runs an in-Python least-cost
# search + Valhalla + a per-request `osmium extract` and can take up to ~2 min;
# the frontend allows 120s (Phase A §8/§15.8). The matching gunicorn
# `--timeout 130` is set in navi-offroute.service. /api/mvum is fast (default).
# -----------------------------------------------------------------------------
location ^~ /api/offroute {
proxy_pass http://127.0.0.1:8428;
proxy_set_header Host $host;
proxy_set_header X-Real-IP $remote_addr;
proxy_set_header X-Authentik-Username $http_x_authentik_username;
proxy_read_timeout 130s; # long-running routing (see above)
add_header X-Cache-Status BYPASS;
}
location ^~ /api/mvum {
proxy_pass http://127.0.0.1:8428;
proxy_set_header Host $host;
proxy_set_header X-Real-IP $remote_addr;
proxy_set_header X-Authentik-Username $http_x_authentik_username;
proxy_read_timeout 15s;
add_header X-Cache-Status BYPASS;
}

View file

@ -0,0 +1,18 @@
[Unit]
Description=navi-offroute — off-network router + MVUM API (Echo6 navi-backend, extraction #8)
After=network-online.target
Wants=network-online.target
[Service]
# User=zvx is REQUIRED: PADUS PostGIS access uses libpq peer-auth via the local
# socket (dbname=padus, no password). zvx is verified to have padus access; a
# different user would fail to connect at runtime.
User=zvx
WorkingDirectory=/home/zvx/projects/repos/navi-backend
EnvironmentFile=/etc/navi-backend/navi-offroute.env
ExecStart=/home/zvx/projects/repos/navi-backend/.venv/bin/gunicorn 'services.navi_offroute.app:create_app()' --bind 127.0.0.1:8428 --workers 2 --timeout 130
Restart=on-failure
RestartSec=2
[Install]
WantedBy=multi-user.target

View file

@ -24,6 +24,10 @@ dependencies = [
"numpy>=1.24", # planet-DEM tile decode
"Pillow>=10", # planet-DEM Terrarium WebP decode
"pmtiles>=3", # planet-DEM PMTiles reader
# navi-offroute (extraction #8): off-network router + readers.
"scikit-image>=0.22", # MCP_Geometric least-cost pathfinding (router.py)
"rasterio>=1.3", # barriers/wilderness/trails/friction raster readers
"psutil>=5.9", # MEMORY_LIMIT_GB enforcement in router.py
]
[tool.setuptools.packages.find]

View file

@ -0,0 +1,135 @@
"""navi-offroute admin-info endpoint (handoff §4.5).
``GET /api/admin/navi-offroute/info`` Authentik-gated, read-only.
Per Phase A §10 this service has NO secrets: PADUS uses libpq peer-auth
(``dbname=padus``, no password); everything else is non-secret paths/URLs.
Probes are cheap (file existence + version/ping) NO COUNT/DISTINCT against
navi.db, which is a fleet-aggregator hot path (see the navi-geo cold-start
lesson).
"""
import os
import subprocess
import time
import psycopg2
import requests
from flask import Blueprint, jsonify, current_app
from shared.auth import require_auth
from shared.admin_info import build_info_response
from shared.dem import dem_path
from . import router as router_mod
from .mvum import navi_db_path
from .barriers import barriers_tif_path, wilderness_tif_path
from .friction import friction_vrt_path
from .trails import trails_tif_path
bp = Blueprint('offroute_admin', __name__)
PORT = 8428
def _valhalla_dependency():
start = time.monotonic()
try:
resp = requests.get(f"{router_mod.VALHALLA_URL}/status", timeout=3)
latency_ms = round((time.monotonic() - start) * 1000, 1)
ok = resp.status_code == 200
r = {'name': 'valhalla', 'status': 'ok' if ok else 'error', 'latency_ms': latency_ms}
if not ok:
r['error'] = f'HTTP {resp.status_code}'
return r
except Exception as e:
latency_ms = round((time.monotonic() - start) * 1000, 1)
return {'name': 'valhalla', 'status': 'error', 'latency_ms': latency_ms, 'error': type(e).__name__}
def _padus_pg_dependency():
"""SELECT 1 over the peer-auth DSN (no password). Cheap liveness check."""
start = time.monotonic()
conn = None
try:
conn = psycopg2.connect(router_mod.POSTGIS_DSN, connect_timeout=3)
with conn.cursor() as cur:
cur.execute("SELECT 1")
cur.fetchone()
latency_ms = round((time.monotonic() - start) * 1000, 1)
return {'name': 'padus-postgis', 'status': 'ok', 'latency_ms': latency_ms}
except Exception as e:
latency_ms = round((time.monotonic() - start) * 1000, 1)
return {'name': 'padus-postgis', 'status': 'error', 'latency_ms': latency_ms, 'error': type(e).__name__}
finally:
if conn is not None:
try:
conn.close()
except Exception:
pass
def _osmium_dependency():
"""osmium-tool version (router shells out to `osmium extract` per route)."""
start = time.monotonic()
try:
out = subprocess.check_output(['osmium', '--version'], stderr=subprocess.STDOUT,
text=True, timeout=3)
latency_ms = round((time.monotonic() - start) * 1000, 1)
version = out.splitlines()[0].strip() if out else 'unknown'
return {'name': 'osmium-tool', 'status': 'ok', 'latency_ms': latency_ms, 'version': version}
except Exception as e:
latency_ms = round((time.monotonic() - start) * 1000, 1)
return {'name': 'osmium-tool', 'status': 'error', 'latency_ms': latency_ms,
'error': 'not installed' if isinstance(e, FileNotFoundError) else type(e).__name__}
def _file_entry(name, path):
"""Cheap existence/readable report — never errors, never reads contents."""
p = str(path)
return {'name': name, 'path': p, 'exists': os.path.exists(p), 'readable': os.access(p, os.R_OK)}
@bp.route('/api/admin/navi-offroute/info')
@require_auth
def navi_offroute_info():
metrics = current_app.config['METRICS']
osm_pbf = str(router_mod.OSM_PBF_PATH)
info = build_info_response(
service='navi-offroute',
version=current_app.config.get('VERSION', 'unknown'),
port=PORT,
config={},
# No secrets (Phase A §10) — peer-auth DSN carries no password.
env=[
{'name': 'NAVI_OFFROUTE_VALHALLA_URL', 'value': router_mod.VALHALLA_URL},
{'name': 'NAVI_OFFROUTE_POSTGIS_DSN', 'value': router_mod.POSTGIS_DSN},
{'name': 'NAVI_OFFROUTE_DENSIFY_M', 'value': str(router_mod.DENSIFY_INTERVAL_M)},
{'name': 'NAVI_OFFROUTE_OSM_PBF', 'value': osm_pbf},
{'name': 'NAVI_OFFROUTE_NAVI_DB', 'value': str(navi_db_path())},
{'name': 'NAVI_DEM_PMTILES', 'value': str(dem_path())},
{'name': 'NAVI_OFFROUTE_BARRIERS_TIF', 'value': str(barriers_tif_path())},
{'name': 'NAVI_OFFROUTE_WILDERNESS_TIF', 'value': str(wilderness_tif_path())},
{'name': 'NAVI_OFFROUTE_TRAILS_TIF', 'value': str(trails_tif_path())},
{'name': 'NAVI_OFFROUTE_FRICTION_VRT', 'value': str(friction_vrt_path())},
],
dependencies=[
_valhalla_dependency(),
_padus_pg_dependency(),
_osmium_dependency(),
],
filesystem=[
_file_entry('dem', dem_path()),
_file_entry('osm_pbf', osm_pbf),
_file_entry('navi_db', navi_db_path()),
_file_entry('barriers_tif', barriers_tif_path()),
_file_entry('wilderness_tif', wilderness_tif_path()),
_file_entry('trails_tif', trails_tif_path()),
_file_entry('friction_vrt', friction_vrt_path()),
],
runtime={
'uptime_s': round(time.time() - metrics['start_time'], 1),
'request_count': metrics['request_count'],
'last_error_at': metrics['last_error_at'],
},
)
return jsonify(info)

View file

@ -0,0 +1,39 @@
"""navi-offroute Flask application factory + gunicorn entry.
Gunicorn entry:
gunicorn 'services.navi_offroute.app:create_app()' --bind 127.0.0.1:8428 --workers 2
"""
import time
from flask import Flask
from shared.git_sha import git_short_sha
from . import offroute_route, admin
def create_app():
app = Flask(__name__)
app.config['VERSION'] = git_short_sha()
app.config['METRICS'] = {
'start_time': time.time(),
'request_count': 0,
'last_error_at': None,
}
@app.before_request
def _count_request():
app.config['METRICS']['request_count'] += 1
@app.after_request
def _track_errors(response):
if response.status_code >= 500:
app.config['METRICS']['last_error_at'] = time.strftime(
'%Y-%m-%dT%H:%M:%SZ', time.gmtime()
)
return response
app.register_blueprint(offroute_route.bp)
app.register_blueprint(admin.bp)
return app

View file

@ -0,0 +1,163 @@
"""
PAD-US barrier and wilderness layers for OFFROUTE.
Provides access to:
1. Barrier raster (Pub_Access = 'XA' - closed/restricted areas)
2. Wilderness raster (Des_Tp = 'WA' - designated wilderness areas)
Runtime readers only. The offline raster-build functions that rasterize the
PAD-US geodatabase via gdal/ogr (recon's build_barriers_raster /
build_wilderness_raster) are NOT part of the service the rasters are a
read-only input produced out of band (Phase A §3/§15.2).
"""
import os
from pathlib import Path
from typing import Tuple, Optional
import numpy as np
try:
import rasterio
from rasterio.windows import from_bounds
from rasterio.enums import Resampling
except ImportError:
raise ImportError("rasterio is required for barriers layer support")
# Default raster paths (single source of truth); env-overridable via the helpers.
DEFAULT_BARRIERS_PATH = Path("/mnt/nav/worldcover/padus_barriers.tif")
DEFAULT_WILDERNESS_PATH = Path("/mnt/nav/worldcover/wilderness.tif")
def barriers_tif_path() -> Path:
"""Barrier raster path, env-overridable via NAVI_OFFROUTE_BARRIERS_TIF."""
return Path(os.environ.get("NAVI_OFFROUTE_BARRIERS_TIF", str(DEFAULT_BARRIERS_PATH)))
def wilderness_tif_path() -> Path:
"""Wilderness raster path, env-overridable via NAVI_OFFROUTE_WILDERNESS_TIF."""
return Path(os.environ.get("NAVI_OFFROUTE_WILDERNESS_TIF", str(DEFAULT_WILDERNESS_PATH)))
class BarrierReader:
"""Reader for PAD-US barrier raster (closed/restricted areas)."""
def __init__(self, barrier_path: Path = None):
self.barrier_path = Path(barrier_path) if barrier_path else barriers_tif_path()
self._dataset = None
def _open(self):
"""Lazy open the dataset."""
if self._dataset is None:
if not self.barrier_path.exists():
raise FileNotFoundError(f"Barrier raster not found at {self.barrier_path}")
self._dataset = rasterio.open(self.barrier_path)
return self._dataset
def get_barrier_grid(
self,
south: float,
north: float,
west: float,
east: float,
target_shape: Tuple[int, int]
) -> np.ndarray:
"""
Get barrier 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 barrier values:
255 = closed/restricted (impassable when respect_boundaries=True)
0 = public/accessible
"""
ds = self._open()
window = from_bounds(west, south, east, north, ds.transform)
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()
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])
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 = None):
self.wilderness_path = Path(wilderness_path) if wilderness_path else wilderness_tif_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}")
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])
def close(self):
"""Close the dataset."""
if self._dataset is not None:
self._dataset.close()
self._dataset = None

View file

@ -0,0 +1,494 @@
"""
Multi-mode travel cost functions for OFFROUTE.
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 dataclasses import dataclass, field
from typing import Optional, Literal, Dict, Callable
# ═══════════════════════════════════════════════════════════════════════════════
# 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
),
}
# Pragmatic mode friction multiplier for private land
PRAGMATIC_BARRIER_MULTIPLIER = 5.0
# ═══════════════════════════════════════════════════════════════════════════════
# COST GRID COMPUTATION
# ═══════════════════════════════════════════════════════════════════════════════
def compute_cost_grid(
elevation: np.ndarray,
cell_size_m: float,
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,
wilderness: Optional[np.ndarray] = None,
mvum: 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.
Args:
elevation: 2D array of elevation values in meters
cell_size_m: Average cell size in meters
cell_size_lat_m: Cell size in latitude direction (optional)
cell_size_lon_m: Cell size in longitude direction (optional)
friction: Optional 2D array of friction multipliers (WorldCover).
Values should be float (1.0 = baseline, 2.0 = 2x slower).
np.inf marks impassable cells.
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, 5 = road, 15 = track, 25 = foot trail
barriers: Optional 2D array of barrier values (uint8).
255 = closed/restricted area (PAD-US Pub_Access = XA).
wilderness: Optional[np.ndarray] of wilderness values (uint8).
255 = designated wilderness area.
mvum: Optional[np.ndarray] of MVUM access values (uint8).
0 = no MVUM data, 1 = open, 255 = closed to this mode.
MVUM closures respond to boundary_mode (strict/pragmatic/emergency).
Foot mode should pass None (MVUM is motor-vehicle specific).
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'")
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
if cell_size_lon_m is None:
cell_size_lon_m = cell_size_m
rows, cols = elevation.shape
# ─── Compute gradients (in-place where possible) ─────────────────────────
# Use float32 to reduce memory footprint
grade = np.zeros(elevation.shape, dtype=np.float32)
# 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
# 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
# grade = sqrt(dx^2 + dy^2)
np.sqrt(dy_contrib, out=grade)
del dy_contrib # Free memory immediately
# ─── Compute speed based on mode ─────────────────────────────────────────
max_grade_val = np.tan(np.radians(profile.max_slope_deg))
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}")
# ─── 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
del speed_kmh
# ─── Max slope limit ──────────────────────────────────────────────────────
cost[grade > max_grade_val] = np.inf
# ─── NaN elevations ──────────────────────────────────────────────────────
cost[np.isnan(elevation)] = np.inf
# ─── 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 mismatch")
np.multiply(cost, friction, out=cost)
# ─── 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 mismatch")
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
# ─── 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 mismatch")
if boundary_mode == "strict":
np.putmask(cost, barriers == 255, np.inf)
elif boundary_mode == "pragmatic":
barrier_mask = barriers == 255
cost[barrier_mask] *= PRAGMATIC_BARRIER_MULTIPLIER
del barrier_mask
# ─── MVUM closures (motor vehicle restrictions) ──────────────────────────
# MVUM only applies to motorized modes, not foot. Foot mode should pass mvum=None.
# MVUM closures respond to the same boundary_mode as PAD-US barriers:
# "strict" = MVUM-closed road/trail is impassable
# "pragmatic" = MVUM-closed road/trail gets 5× friction penalty
# "emergency" = MVUM closures ignored entirely
if mvum is not None and mode != "foot" and boundary_mode != "emergency":
if mvum.shape != elevation.shape:
raise ValueError(f"MVUM shape mismatch")
# Value 255 = road/trail exists but is closed to this mode
mvum_closed_mask = mvum == 255
if boundary_mode == "strict":
np.putmask(cost, mvum_closed_mask, np.inf)
elif boundary_mode == "pragmatic":
cost[mvum_closed_mask] *= PRAGMATIC_BARRIER_MULTIPLIER
del mvum_closed_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("=" * 70)
print("OFFROUTE Multi-Mode Cost Function Tests")
print("=" * 70)
print("\n[1] Speed functions at various grades:")
print(f"{'Grade':<10} {'Foot':<12} {'MTB':<12} {'ATV':<12} {'Vehicle':<12}")
print("-" * 60)
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")
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}")
print("\n[3] Cost grid test (flat terrain, forest):")
elev = np.ones((10, 10), dtype=np.float32) * 1000
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
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))
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.")

View file

@ -0,0 +1,146 @@
"""
Friction layer reader for OFFROUTE.
Reads friction values from the WorldCover friction VRT and resamples
to match the elevation grid dimensions.
"""
import os
from pathlib import Path
from typing import Tuple, Optional
import numpy as np
try:
import rasterio
from rasterio.windows import from_bounds
from rasterio.enums import Resampling
except ImportError:
raise ImportError("rasterio is required for friction layer support")
# Default path to the friction VRT (single source of truth); env-overridable.
DEFAULT_FRICTION_PATH = Path("/mnt/nav/worldcover/friction/friction_conus.vrt")
def friction_vrt_path() -> Path:
"""Friction VRT path, env-overridable via NAVI_OFFROUTE_FRICTION_VRT."""
return Path(os.environ.get("NAVI_OFFROUTE_FRICTION_VRT", str(DEFAULT_FRICTION_PATH)))
class FrictionReader:
"""Reader for WorldCover friction raster."""
def __init__(self, friction_path: Path = None):
self.friction_path = Path(friction_path) if friction_path else friction_vrt_path()
self._dataset = None
def _open(self):
"""Lazy open the dataset."""
if self._dataset is None:
if not self.friction_path.exists():
raise FileNotFoundError(f"Friction VRT not found at {self.friction_path}")
self._dataset = rasterio.open(self.friction_path)
return self._dataset
def get_friction_grid(
self,
south: float,
north: float,
west: float,
east: float,
target_shape: Tuple[int, int]
) -> np.ndarray:
"""
Get friction values for a bounding box, resampled to target shape.
Args:
south, north, west, east: Bounding box coordinates
target_shape: (rows, cols) to resample to (matches elevation grid)
Returns:
np.ndarray of uint8 friction values, same shape as target_shape.
Values: 10-40 = friction multiplier (divide by 10)
255 = impassable
0 = nodata (treat as impassable)
"""
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
# Use nearest neighbor for categorical data
friction = ds.read(
1,
window=window,
out_shape=target_shape,
resampling=Resampling.nearest
)
return friction
def sample_point(self, lat: float, lon: float) -> int:
"""Sample friction 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 = nodata
# Read single pixel
window = rasterio.windows.Window(col, row, 1, 1)
value = ds.read(1, window=window)
return int(value[0, 0])
def close(self):
"""Close the dataset."""
if self._dataset is not None:
self._dataset.close()
self._dataset = None
def friction_to_multiplier(friction: np.ndarray) -> np.ndarray:
"""
Convert friction values to cost multipliers.
Args:
friction: uint8 array of friction values
Returns:
float32 array of multipliers.
Values 10-40 become 1.0-4.0 (divide by 10).
Values 0 or 255 become np.inf (impassable).
"""
multiplier = friction.astype(np.float32) / 10.0
# Mark impassable cells
multiplier[friction == 0] = np.inf # nodata
multiplier[friction == 255] = np.inf # water/impassable
return multiplier
if __name__ == "__main__":
print("Testing FrictionReader...")
reader = FrictionReader()
# Test point sampling - Murtaugh Lake (should be water = 255)
lake_lat, lake_lon = 42.47, -114.15
lake_friction = reader.sample_point(lake_lat, lake_lon)
print(f"Murtaugh Lake ({lake_lat}, {lake_lon}): friction = {lake_friction}")
print(f" Expected: 255 (water/impassable)")
# Test grid read for small bbox
friction = reader.get_friction_grid(
south=42.4, north=42.5, west=-114.2, east=-114.1,
target_shape=(100, 100)
)
print(f"\nGrid test shape: {friction.shape}")
print(f"Unique values: {np.unique(friction)}")
print(f"Water cells (255): {np.sum(friction == 255)}")
reader.close()
print("\nFrictionReader test complete.")

View file

@ -0,0 +1,618 @@
"""
MVUM (Motor Vehicle Use Map) legal access layer for OFFROUTE.
Queries USFS MVUM data from navi.db and provides rasterized access grids
indicating which roads/trails are open or closed to specific vehicle modes.
MVUM is motor-vehicle specific foot mode should skip this layer entirely.
"""
import os
import re
import sqlite3
import warnings
from datetime import datetime
from pathlib import Path
from typing import Dict, List, Optional, Tuple, Literal
import numpy as np
from shapely import wkb
from shapely.geometry import Point
# Path to navi.db (single source of truth); env-overridable.
DEFAULT_NAVI_DB_PATH = Path("/mnt/nav/navi.db")
def navi_db_path() -> Path:
"""navi.db (MVUM tables) path, env-overridable via NAVI_OFFROUTE_NAVI_DB."""
return Path(os.environ.get("NAVI_OFFROUTE_NAVI_DB", str(DEFAULT_NAVI_DB_PATH)))
def parse_date_range(date_str: str) -> List[Tuple[int, int, int, int]]:
"""
Parse MVUM date range strings like "05/01-11/30" or "06/15-10/15,12/01-03/31".
Returns list of (start_month, start_day, end_month, end_day) tuples.
Returns empty list if unparseable.
"""
if not date_str or date_str.strip() == "":
return []
ranges = []
# Split by comma for multi-period strings
for part in date_str.split(","):
part = part.strip()
# Match MM/DD-MM/DD pattern
match = re.match(r"(\d{1,2})/(\d{1,2})-(\d{1,2})/(\d{1,2})", part)
if match:
try:
sm, sd, em, ed = int(match.group(1)), int(match.group(2)), int(match.group(3)), int(match.group(4))
if 1 <= sm <= 12 and 1 <= sd <= 31 and 1 <= em <= 12 and 1 <= ed <= 31:
ranges.append((sm, sd, em, ed))
except ValueError:
pass
return ranges
def is_date_in_range(month: int, day: int, ranges: List[Tuple[int, int, int, int]]) -> bool:
"""
Check if a given month/day falls within any of the date ranges.
Handles ranges that wrap around year end (e.g., 12/01-03/31).
"""
if not ranges:
return True # No ranges = assume open
date_num = month * 100 + day # Simple numeric comparison
for sm, sd, em, ed in ranges:
start_num = sm * 100 + sd
end_num = em * 100 + ed
if start_num <= end_num:
# Normal range (e.g., 05/01-11/30)
if start_num <= date_num <= end_num:
return True
else:
# Wrapping range (e.g., 12/01-03/31)
if date_num >= start_num or date_num <= end_num:
return True
return False
def check_access(
status_field: Optional[str],
dates_field: Optional[str],
seasonal: Optional[str],
check_date: Optional[Tuple[int, int]] = None
) -> Optional[bool]:
"""
Determine if a road/trail is open to a vehicle type.
Args:
status_field: Value of vehicle-class field (e.g., "open", null)
dates_field: Value of *_DATESOPEN field (e.g., "05/01-11/30")
seasonal: Value of SEASONAL field ("yearlong", "seasonal")
check_date: Optional (month, day) tuple to check against date ranges
Returns:
True = open
False = closed
None = no data (field not populated, defer to SYMBOL)
"""
if status_field is None or status_field.strip() == "":
return None # No data
status = status_field.strip().lower()
if status != "open":
return False # Explicitly closed or restricted
# Status is "open" - check seasonal restrictions
if check_date is not None:
month, day = check_date
# Parse date ranges
if dates_field:
ranges = parse_date_range(dates_field)
if ranges:
return is_date_in_range(month, day, ranges)
# No date field but seasonal = "yearlong" means always open
if seasonal and seasonal.strip().lower() == "yearlong":
return True
# Seasonal with no dates - assume open (data quality issue)
if seasonal and seasonal.strip().lower() == "seasonal":
warnings.warn(f"Seasonal road/trail with no DATESOPEN, assuming open")
return True
return True # Open with no date check
def get_mode_field(mode: str) -> Tuple[str, str]:
"""
Get the MVUM field names for a given travel mode.
Returns (status_field, dates_field) tuple.
"""
mode_mapping = {
"atv": ("atv", "atv_datesopen"),
"motorcycle": ("motorcycle", "motorcycle_datesopen"),
"mtb": ("e_bike_class1", "e_bike_class1_dur"), # Closest analog for e-bikes
"vehicle": ("highclearancevehicle", "highclearancevehicle_datesopen"),
"passenger": ("passengervehicle", "passengervehicle_datesopen"),
}
return mode_mapping.get(mode, ("highclearancevehicle", "highclearancevehicle_datesopen"))
def symbol_to_access(symbol: str, mode: str, maint_level: Optional[str] = None) -> Optional[bool]:
"""
Fallback: interpret SYMBOL field when per-vehicle-class fields are null.
MVUM SYMBOL meanings (roads):
1 = Open to all vehicles
2 = Open to highway legal vehicles only
3 = Road closed to motorized
4 = Road open seasonally
11 = Administrative use only
12 = Decommissioned
For trails, similar logic applies based on TRAILCLASS.
"""
if symbol is None:
return None
sym = str(symbol).strip()
# Symbol 1: Open to all
if sym == "1":
return True
# Symbol 2: Highway legal only
if sym == "2":
# ATVs/motorcycles typically not highway legal
if mode in ("atv", "motorcycle"):
return False
return True
# Symbol 3: Closed to motorized
if sym == "3":
return False
# Symbol 4: Seasonally open (assume open if no date check)
if sym == "4":
return True
# Symbol 11/12: Administrative/decommissioned = closed
if sym in ("11", "12"):
return False
# Unknown symbol - defer
return None
class MVUMReader:
"""
Reader for MVUM data from navi.db.
Queries roads and trails by bounding box and returns access grids.
"""
def __init__(self, db_path: Path = None):
self.db_path = Path(db_path) if db_path else navi_db_path()
self._conn = None
def _get_conn(self) -> sqlite3.Connection:
if self._conn is None:
if not self.db_path.exists():
raise FileNotFoundError(f"navi.db not found at {self.db_path}")
self._conn = sqlite3.connect(str(self.db_path))
self._conn.row_factory = sqlite3.Row
# Load Spatialite extension if available
try:
self._conn.enable_load_extension(True)
self._conn.load_extension("mod_spatialite")
except Exception:
pass # Spatialite not available, will use manual bbox queries
return self._conn
def table_exists(self, table_name: str) -> bool:
"""Check if an MVUM table exists."""
conn = self._get_conn()
cur = conn.execute(
"SELECT name FROM sqlite_master WHERE type='table' AND name=?",
(table_name,)
)
return cur.fetchone() is not None
def query_roads_bbox(
self,
south: float, north: float, west: float, east: float,
mode: str = "atv",
check_date: Optional[Tuple[int, int]] = None
) -> List[Dict]:
"""
Query MVUM roads within a bounding box.
Returns list of dicts with access info for the given mode.
"""
if not self.table_exists("mvum_roads"):
return []
conn = self._get_conn()
# Query using bbox on geometry
# Since we don't have spatialite, we'll query all and filter in Python
# For production, consider pre-computing bbox columns
cur = conn.execute("""
SELECT ogc_fid, id, name, symbol, operationalmaintlevel, seasonal,
atv, atv_datesopen, motorcycle, motorcycle_datesopen,
highclearancevehicle, highclearancevehicle_datesopen,
passengervehicle, passengervehicle_datesopen,
e_bike_class1, e_bike_class1_dur,
shape
FROM mvum_roads
""")
status_field, dates_field = get_mode_field(mode)
results = []
for row in cur:
# Parse geometry to check bbox intersection
# The shape is stored as WKB blob
shape = row["shape"]
if shape is None:
continue
# Quick bbox check using geometry extent
# Since we don't have Spatialite functions, we'll include all
# and let the rasterization handle it
access = check_access(
row[status_field] if status_field in row.keys() else None,
row[dates_field] if dates_field in row.keys() else None,
row["seasonal"],
check_date
)
# Fallback to SYMBOL if no per-vehicle data
if access is None:
access = symbol_to_access(row["symbol"], mode, row["operationalmaintlevel"])
if access is not None:
results.append({
"id": row["id"],
"name": row["name"],
"access": access,
"symbol": row["symbol"],
"maint_level": row["operationalmaintlevel"],
"shape": shape,
})
return results
def query_trails_bbox(
self,
south: float, north: float, west: float, east: float,
mode: str = "atv",
check_date: Optional[Tuple[int, int]] = None
) -> List[Dict]:
"""
Query MVUM trails within a bounding box.
"""
if not self.table_exists("mvum_trails"):
return []
conn = self._get_conn()
cur = conn.execute("""
SELECT ogc_fid, id, name, symbol, seasonal, trailclass,
atv, atv_datesopen, motorcycle, motorcycle_datesopen,
highclearancevehicle, highclearancevehicle_datesopen,
passengervehicle, passengervehicle_datesopen,
e_bike_class1, e_bike_class1_dur,
shape
FROM mvum_trails
""")
status_field, dates_field = get_mode_field(mode)
results = []
for row in cur:
shape = row["shape"]
if shape is None:
continue
access = check_access(
row[status_field] if status_field in row.keys() else None,
row[dates_field] if dates_field in row.keys() else None,
row["seasonal"],
check_date
)
if access is None:
access = symbol_to_access(row["symbol"], mode)
if access is not None:
results.append({
"id": row["id"],
"name": row["name"],
"access": access,
"symbol": row["symbol"],
"trail_class": row["trailclass"],
"shape": shape,
})
return results
def query_nearest(
self,
lat: float, lon: float,
radius_m: float = 50,
table: str = "mvum_roads"
) -> Optional[Dict]:
"""
Query the nearest MVUM feature to a point.
Used for the places panel API.
"""
if not self.table_exists(table):
return None
conn = self._get_conn()
# Convert radius to degrees (approximate)
radius_deg = radius_m / 111000
# Query features in bbox around point
if table == "mvum_roads":
cur = conn.execute("""
SELECT ogc_fid, id, name, forestname, districtname, symbol,
operationalmaintlevel, surfacetype, seasonal, jurisdiction,
passengervehicle, passengervehicle_datesopen,
highclearancevehicle, highclearancevehicle_datesopen,
atv, atv_datesopen, motorcycle, motorcycle_datesopen,
fourwd_gt50inches, fourwd_gt50_datesopen,
twowd_gt50inches, twowd_gt50_datesopen,
e_bike_class1, e_bike_class1_dur,
e_bike_class2, e_bike_class2_dur,
e_bike_class3, e_bike_class3_dur,
shape
FROM mvum_roads
LIMIT 1000
""")
else:
cur = conn.execute("""
SELECT ogc_fid, id, name, forestname, districtname, symbol,
seasonal, jurisdiction, trailclass, trailsystem,
passengervehicle, passengervehicle_datesopen,
highclearancevehicle, highclearancevehicle_datesopen,
atv, atv_datesopen, motorcycle, motorcycle_datesopen,
fourwd_gt50inches, fourwd_gt50_datesopen,
twowd_gt50inches, twowd_gt50_datesopen,
e_bike_class1, e_bike_class1_dur,
e_bike_class2, e_bike_class2_dur,
e_bike_class3, e_bike_class3_dur,
shape
FROM mvum_trails
LIMIT 1000
""")
# Find nearest feature
# This is a simplified approach - for production, use spatial index
query_point = Point(lon, lat)
nearest = None
min_dist = float('inf')
for row in cur:
try:
geom = wkb.loads(row["shape"])
dist = query_point.distance(geom)
if dist < min_dist and dist < radius_deg:
min_dist = dist
nearest = dict(row)
nearest["geometry"] = geom
except Exception:
continue
if nearest:
# Convert geometry to GeoJSON
nearest["geojson"] = nearest["geometry"].__geo_interface__
del nearest["geometry"]
del nearest["shape"]
return nearest
return None
def close(self):
if self._conn:
self._conn.close()
self._conn = None
def get_mvum_access_grid(
south: float, north: float, west: float, east: float,
target_shape: Tuple[int, int],
mode: Literal["foot", "mtb", "atv", "vehicle"] = "atv",
check_date: Optional[str] = None,
db_path: Path = None # None → navi_db_path() (env NAVI_OFFROUTE_NAVI_DB)
) -> np.ndarray:
"""
Get MVUM access grid for pathfinding.
Args:
south, north, west, east: Bounding box (WGS84)
target_shape: (rows, cols) to match elevation grid
mode: Travel mode (foot skips MVUM entirely)
check_date: Optional "MM/DD" string for seasonal checking
db_path: Path to navi.db
Returns:
np.ndarray of uint8:
0 = no MVUM data (defer to existing trail/friction logic)
1 = road/trail is OPEN to this vehicle mode
255 = road/trail EXISTS but is CLOSED to this mode
"""
# Foot mode bypasses MVUM entirely
if mode == "foot":
return np.zeros(target_shape, dtype=np.uint8)
# Parse check_date if provided
parsed_date = None
if check_date:
match = re.match(r"(\d{1,2})/(\d{1,2})", check_date)
if match:
parsed_date = (int(match.group(1)), int(match.group(2)))
# Initialize output grid
grid = np.zeros(target_shape, dtype=np.uint8)
rows, cols = target_shape
# Pixel size
pixel_lat = (north - south) / rows
pixel_lon = (east - west) / cols
reader = MVUMReader(db_path)
try:
# Query roads and trails
roads = reader.query_roads_bbox(south, north, west, east, mode, parsed_date)
trails = reader.query_trails_bbox(south, north, west, east, mode, parsed_date)
# Rasterize features
for features in [roads, trails]:
for feat in features:
try:
geom = wkb.loads(feat["shape"])
# Get geometry bounds
minx, miny, maxx, maxy = geom.bounds
# Check if intersects our bbox
if maxx < west or minx > east or maxy < south or miny > north:
continue
# Rasterize line
value = 1 if feat["access"] else 255
# Simple line rasterization
if geom.geom_type in ("LineString", "MultiLineString"):
if geom.geom_type == "MultiLineString":
coords_list = [list(line.coords) for line in geom.geoms]
else:
coords_list = [list(geom.coords)]
for coords in coords_list:
for i in range(len(coords) - 1):
x1, y1 = coords[i]
x2, y2 = coords[i + 1]
# Convert to pixel coordinates
col1 = int((x1 - west) / pixel_lon)
row1 = int((north - y1) / pixel_lat)
col2 = int((x2 - west) / pixel_lon)
row2 = int((north - y2) / pixel_lat)
# Bresenham's line algorithm
_draw_line(grid, row1, col1, row2, col2, value)
except Exception as e:
continue
finally:
reader.close()
return grid
def _draw_line(grid: np.ndarray, r1: int, c1: int, r2: int, c2: int, value: int):
"""Draw a line on the grid using Bresenham's algorithm."""
rows, cols = grid.shape
dr = abs(r2 - r1)
dc = abs(c2 - c1)
sr = 1 if r1 < r2 else -1
sc = 1 if c1 < c2 else -1
err = dr - dc
r, c = r1, c1
while True:
if 0 <= r < rows and 0 <= c < cols:
# Only overwrite if current value is 0 (no data) or we're marking closed
if grid[r, c] == 0 or value == 255:
grid[r, c] = value
if r == r2 and c == c2:
break
e2 = 2 * err
if e2 > -dc:
err -= dc
r += sr
if e2 < dr:
err += dr
c += sc
if __name__ == "__main__":
import sys
print("=" * 60)
print("MVUM Reader Test")
print("=" * 60)
reader = MVUMReader()
if not reader.table_exists("mvum_roads"):
print("ERROR: mvum_roads table not found in navi.db")
sys.exit(1)
# Test bbox query (Sawtooth NF area)
print("\n[1] Testing bbox query (Sawtooth NF area)...")
roads = reader.query_roads_bbox(
south=43.5, north=44.0, west=-115.0, east=-114.0,
mode="atv"
)
print(f" Found {len(roads)} roads")
open_count = sum(1 for r in roads if r["access"])
closed_count = sum(1 for r in roads if not r["access"])
print(f" Open to ATV: {open_count}")
print(f" Closed to ATV: {closed_count}")
# Test with seasonal date
print("\n[2] Testing with date check (July 15)...")
roads_summer = reader.query_roads_bbox(
south=43.5, north=44.0, west=-115.0, east=-114.0,
mode="atv",
check_date=(7, 15)
)
open_summer = sum(1 for r in roads_summer if r["access"])
print(f" Open to ATV on 07/15: {open_summer}")
print("\n[3] Testing with date check (January 15)...")
roads_winter = reader.query_roads_bbox(
south=43.5, north=44.0, west=-115.0, east=-114.0,
mode="atv",
check_date=(1, 15)
)
open_winter = sum(1 for r in roads_winter if r["access"])
print(f" Open to ATV on 01/15: {open_winter}")
# Test grid generation
print("\n[4] Testing grid generation...")
grid = get_mvum_access_grid(
south=43.5, north=44.0, west=-115.0, east=-114.0,
target_shape=(500, 1000),
mode="atv"
)
print(f" Grid shape: {grid.shape}")
print(f" No data (0): {np.sum(grid == 0)}")
print(f" Open (1): {np.sum(grid == 1)}")
print(f" Closed (255): {np.sum(grid == 255)}")
reader.close()
print("\nDone.")

View file

@ -0,0 +1,152 @@
"""navi-offroute API blueprint — faithful port of recon's offroute routes.
POST /api/offroute off-network effort-based routing (OffrouteRouter)
GET /api/mvum MVUM road/trail access lookup (MVUMReader)
Both public (no auth), matching recon. Same request/response shapes and status
codes. Ported from recon's lib/api.py:api_offroute / api_mvum.
"""
import logging
import re
from flask import Blueprint, request, jsonify
from .router import OffrouteRouter
from .mvum import MVUMReader
logger = logging.getLogger('navi_offroute.route')
bp = Blueprint('offroute', __name__)
VALID_MODES = ("auto", "foot", "mtb", "atv", "vehicle")
VALID_BOUNDARY_MODES = ("strict", "pragmatic", "emergency")
@bp.route("/api/offroute", methods=["POST"])
def api_offroute():
"""
Off-network routing from wilderness to destination.
Request body:
{start:[lat,lon], end:[lat,lon],
mode: auto|foot|mtb|atv|vehicle (default foot),
boundary_mode: strict|pragmatic|emergency (default pragmatic)}
Response: {status:"ok", route:<GeoJSON FeatureCollection>, summary:{...}}
or {status:"error", message}. 400 on bad input / router error, 500 on uncaught.
"""
try:
data = request.get_json()
if not data:
return jsonify({"status": "error", "message": "No JSON body provided"}), 400
start = data.get("start")
end = data.get("end")
if not start or not end:
return jsonify({"status": "error", "message": "Missing start or end coordinates"}), 400
if not isinstance(start, (list, tuple)) or len(start) != 2:
return jsonify({"status": "error", "message": "start must be [lat, lon]"}), 400
if not isinstance(end, (list, tuple)) or len(end) != 2:
return jsonify({"status": "error", "message": "end must be [lat, lon]"}), 400
start_lat, start_lon = float(start[0]), float(start[1])
end_lat, end_lon = float(end[0]), float(end[1])
mode = data.get("mode", "foot")
if mode not in VALID_MODES:
return jsonify({"status": "error", "message": "mode must be auto, foot, mtb, atv, or vehicle"}), 400
boundary_mode = data.get("boundary_mode", "pragmatic")
if boundary_mode not in VALID_BOUNDARY_MODES:
return jsonify({"status": "error", "message": "boundary_mode must be strict, pragmatic, or emergency"}), 400
router = OffrouteRouter()
try:
result = router.route(
start_lat=start_lat, start_lon=start_lon,
end_lat=end_lat, end_lon=end_lon,
mode=mode, boundary_mode=boundary_mode,
)
finally:
router.close()
if result.get("status") == "error":
return jsonify(result), 400
return jsonify(result)
except Exception as e:
logger.exception("Offroute error")
return jsonify({"status": "error", "message": str(e)}), 500
@bp.route("/api/mvum", methods=["GET"])
def api_mvum():
"""MVUM (Motor Vehicle Use Map) access near a point. Roads first, then trails.
GET /api/mvum?lat=&lon=&radius= (radius default 50 m)
Returns {status:"ok", feature:{...}|null}. 400 missing coords, 500 on uncaught.
"""
try:
lat = request.args.get("lat", type=float)
lon = request.args.get("lon", type=float)
radius = request.args.get("radius", 50, type=float)
if lat is None or lon is None:
return jsonify({"status": "error", "message": "lat and lon required"}), 400
reader = MVUMReader()
try:
feature = reader.query_nearest(lat, lon, radius, "mvum_roads")
if feature is None:
feature = reader.query_nearest(lat, lon, radius, "mvum_trails")
if feature is None:
return jsonify({"status": "ok", "feature": None})
access = {
"passenger_vehicle": {"status": feature.get("passengervehicle"),
"dates": feature.get("passengervehicle_datesopen")},
"high_clearance": {"status": feature.get("highclearancevehicle"),
"dates": feature.get("highclearancevehicle_datesopen")},
"atv": {"status": feature.get("atv"), "dates": feature.get("atv_datesopen")},
"motorcycle": {"status": feature.get("motorcycle"),
"dates": feature.get("motorcycle_datesopen")},
"4wd_gt50": {"status": feature.get("fourwd_gt50inches"),
"dates": feature.get("fourwd_gt50_datesopen")},
"2wd_gt50": {"status": feature.get("twowd_gt50inches"),
"dates": feature.get("twowd_gt50_datesopen")},
"e_bike_class1": {"status": feature.get("e_bike_class1"),
"dates": feature.get("e_bike_class1_dur")},
"e_bike_class2": {"status": feature.get("e_bike_class2"),
"dates": feature.get("e_bike_class2_dur")},
"e_bike_class3": {"status": feature.get("e_bike_class3"),
"dates": feature.get("e_bike_class3_dur")},
}
maint_level = feature.get("operationalmaintlevel", "")
maint_num = None
if maint_level:
match = re.match(r"(\d+)", maint_level)
if match:
maint_num = int(match.group(1))
result = {
"id": feature.get("id"),
"name": feature.get("name"),
"forest": feature.get("forestname"),
"district": feature.get("districtname"),
"surface": feature.get("surfacetype"),
"maintenance_level": maint_num,
"seasonal": feature.get("seasonal"),
"symbol": feature.get("symbol"),
"trail_class": feature.get("trailclass"),
"trail_system": feature.get("trailsystem"),
"access": access,
"geometry": feature.get("geojson"),
}
return jsonify({"status": "ok", "feature": result})
finally:
reader.close()
except Exception as e:
logger.exception("MVUM query error")
return jsonify({"status": "error", "message": str(e)}), 500

File diff suppressed because it is too large Load diff

View file

@ -0,0 +1,258 @@
"""Hermetic tests for navi-offroute (extraction #8) — service-shape, not routing
correctness. OffrouteRouter is mocked; MVUM uses a tiny fixture SQLite; admin
probes (Valhalla/PG/osmium) are mocked. No live PostGIS/Valhalla/osmium/DEM.
"""
import sqlite3
import pytest
from shapely import wkb
from shapely.geometry import Point
import services.navi_offroute.offroute_route as route_mod
import services.navi_offroute.admin as admin_mod
from services.navi_offroute.app import create_app
AUTH = {'X-Authentik-Username': 'matt'}
@pytest.fixture
def client():
return create_app().test_client()
# ── /api/offroute — mocked router ─────────────────────────────────────────
class FakeRouter:
instances = []
route_result = {'status': 'ok', 'route': {'type': 'FeatureCollection', 'features': []},
'summary': {'total_distance_km': 1.2, 'total_effort_minutes': 30,
'barrier_crossings': 0, 'mvum_closed_crossings': 0}}
raise_on_init = False
raise_on_route = False
def __init__(self):
if FakeRouter.raise_on_init:
raise RuntimeError('router init boom')
self.closed = False
FakeRouter.instances.append(self)
def route(self, **kwargs):
if FakeRouter.raise_on_route:
raise RuntimeError('route boom')
return FakeRouter.route_result
def close(self):
self.closed = True
@pytest.fixture
def fake_router(monkeypatch):
FakeRouter.instances = []
FakeRouter.raise_on_init = False
FakeRouter.raise_on_route = False
FakeRouter.route_result = {'status': 'ok', 'route': {'type': 'FeatureCollection', 'features': []},
'summary': {'total_distance_km': 1.2, 'total_effort_minutes': 30,
'barrier_crossings': 0, 'mvum_closed_crossings': 0}}
monkeypatch.setattr(route_mod, 'OffrouteRouter', FakeRouter)
return FakeRouter
def _post(client, body):
return client.post('/api/offroute', json=body)
def test_offroute_empty_body_400(client, fake_router):
# Body parses to a falsy value (JSON null) → the "No JSON body provided" 400
# branch. (A *malformed* body makes Flask's get_json() raise BadRequest, which
# the outer except turns into 500 — faithful to recon; not this branch.)
r = client.post('/api/offroute', data='null', content_type='application/json')
assert r.status_code == 400 and r.get_json()['message'] == 'No JSON body provided'
def test_offroute_missing_coords_400(client, fake_router):
assert _post(client, {'start': [43.6, -116.2]}).status_code == 400
def test_offroute_bad_start_shape_400(client, fake_router):
assert _post(client, {'start': [1, 2, 3], 'end': [4, 5]}).status_code == 400
def test_offroute_bad_mode_400(client, fake_router):
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3], 'mode': 'spaceship'})
assert r.status_code == 400 and 'mode must be' in r.get_json()['message']
def test_offroute_bad_boundary_mode_400(client, fake_router):
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3], 'boundary_mode': 'yolo'})
assert r.status_code == 400 and 'boundary_mode must be' in r.get_json()['message']
def test_offroute_happy_path_shape(client, fake_router):
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3], 'mode': 'foot',
'boundary_mode': 'strict'})
assert r.status_code == 200
d = r.get_json()
assert d['status'] == 'ok'
assert d['route']['type'] == 'FeatureCollection'
# the summary keys the UI reads (ManeuverList / DirectionsPanel)
assert {'total_distance_km', 'total_effort_minutes', 'barrier_crossings',
'mvum_closed_crossings'} <= set(d['summary'])
assert fake_router.instances[0].closed is True # always closed
def test_offroute_router_status_error_is_400(client, fake_router):
fake_router.route_result = {'status': 'error', 'message': 'no route found'}
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3]})
assert r.status_code == 400 and r.get_json()['message'] == 'no route found'
def test_offroute_router_init_raises_is_500(client, fake_router):
fake_router.raise_on_init = True
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3]})
assert r.status_code == 500 and r.get_json()['status'] == 'error'
def test_offroute_close_called_even_when_route_raises(client, fake_router):
fake_router.raise_on_route = True
r = _post(client, {'start': [43.6, -116.2], 'end': [43.7, -116.3]})
assert r.status_code == 500 # outer except -> 500
assert fake_router.instances[0].closed is True # finally still closed it
# ── /api/mvum — fixture SQLite ─────────────────────────────────────────────
_ROAD_COLS = ['ogc_fid', 'id', 'name', 'forestname', 'districtname', 'symbol',
'operationalmaintlevel', 'surfacetype', 'seasonal', 'jurisdiction',
'passengervehicle', 'passengervehicle_datesopen',
'highclearancevehicle', 'highclearancevehicle_datesopen',
'atv', 'atv_datesopen', 'motorcycle', 'motorcycle_datesopen',
'fourwd_gt50inches', 'fourwd_gt50_datesopen',
'twowd_gt50inches', 'twowd_gt50_datesopen',
'e_bike_class1', 'e_bike_class1_dur', 'e_bike_class2', 'e_bike_class2_dur',
'e_bike_class3', 'e_bike_class3_dur', 'shape']
_TRAIL_COLS = ['ogc_fid', 'id', 'name', 'forestname', 'districtname', 'symbol',
'seasonal', 'jurisdiction', 'trailclass', 'trailsystem',
'passengervehicle', 'passengervehicle_datesopen',
'highclearancevehicle', 'highclearancevehicle_datesopen',
'atv', 'atv_datesopen', 'motorcycle', 'motorcycle_datesopen',
'fourwd_gt50inches', 'fourwd_gt50_datesopen',
'twowd_gt50inches', 'twowd_gt50_datesopen',
'e_bike_class1', 'e_bike_class1_dur', 'e_bike_class2', 'e_bike_class2_dur',
'e_bike_class3', 'e_bike_class3_dur', 'shape']
def _make_table(conn, table, cols, row):
conn.execute(f"CREATE TABLE {table} ({', '.join(cols)})")
if row is not None:
conn.execute(f"INSERT INTO {table} ({', '.join(cols)}) VALUES ({', '.join(['?'] * len(cols))})",
[row.get(c) for c in cols])
conn.commit()
def _shape_at(lat, lon):
return wkb.dumps(Point(lon, lat))
def _mvum_db(tmp_path, monkeypatch, roads=None, trails=None):
db = tmp_path / 'navi.db'
conn = sqlite3.connect(db)
if roads is not None:
_make_table(conn, 'mvum_roads', _ROAD_COLS, roads)
if trails is not None:
_make_table(conn, 'mvum_trails', _TRAIL_COLS, trails)
conn.close()
monkeypatch.setenv('NAVI_OFFROUTE_NAVI_DB', str(db))
return db
def test_mvum_road_happy_path(client, tmp_path, monkeypatch):
_mvum_db(tmp_path, monkeypatch, roads={
'ogc_fid': 1, 'id': 'FR 123', 'name': 'Some Forest Road',
'forestname': 'Sawtooth National Forest', 'districtname': 'Ketchum RD',
'surfacetype': 'NAT', 'operationalmaintlevel': '2 - HIGH CLEARANCE VEHICLES',
'seasonal': 'Seasonal', 'symbol': 2,
'passengervehicle': 'Open', 'passengervehicle_datesopen': '06/15-10/15',
'atv': 'Open', 'shape': _shape_at(43.6150, -116.2023)})
r = client.get('/api/mvum?lat=43.6150&lon=-116.2023&radius=500')
assert r.status_code == 200
d = r.get_json()
assert d['status'] == 'ok'
f = d['feature']
assert f['id'] == 'FR 123' and f['forest'] == 'Sawtooth National Forest'
assert f['maintenance_level'] == 2 # parsed from "2 - HIGH…"
assert f['access']['passenger_vehicle'] == {'status': 'Open', 'dates': '06/15-10/15'}
assert set(f['access']) == {'passenger_vehicle', 'high_clearance', 'atv', 'motorcycle',
'4wd_gt50', '2wd_gt50', 'e_bike_class1', 'e_bike_class2', 'e_bike_class3'}
def test_mvum_falls_back_to_trails(client, tmp_path, monkeypatch):
# No mvum_roads table → roads query returns None → trails consulted.
_mvum_db(tmp_path, monkeypatch, trails={
'ogc_fid': 1, 'id': 'TR 7', 'name': 'Goat Trail', 'forestname': 'Sawtooth NF',
'trailclass': '2', 'trailsystem': 'Alpine', 'atv': 'Open',
'shape': _shape_at(43.6150, -116.2023)})
f = client.get('/api/mvum?lat=43.6150&lon=-116.2023&radius=500').get_json()['feature']
assert f['id'] == 'TR 7' and f['trail_system'] == 'Alpine'
def test_mvum_no_match_returns_null_feature(client, tmp_path, monkeypatch):
# Road exists but far outside the radius → null feature.
_mvum_db(tmp_path, monkeypatch, roads={
'ogc_fid': 1, 'id': 'FR 999', 'name': 'Far Road',
'shape': _shape_at(0.0, 0.0)})
d = client.get('/api/mvum?lat=43.6150&lon=-116.2023&radius=50').get_json()
assert d == {'status': 'ok', 'feature': None}
def test_mvum_missing_coords_400(client):
assert client.get('/api/mvum?lat=43.6').status_code == 400
def test_friction_reader_raises_file_not_found_when_missing(tmp_path):
"""The FileNotFoundError pre-check (review fix #2) fires before rasterio sees
the path consistent with the barriers/trails readers."""
from services.navi_offroute.friction import FrictionReader
reader = FrictionReader(tmp_path / 'does-not-exist.vrt')
with pytest.raises(FileNotFoundError) as exc:
reader._open()
assert 'Friction VRT not found' in str(exc.value)
# ── admin-info — mocked probes ─────────────────────────────────────────────
def _mock_probes_ok(monkeypatch):
class _Resp:
status_code = 200
monkeypatch.setattr(admin_mod.requests, 'get', lambda *a, **k: _Resp())
class _Cur:
def execute(self, *a): pass
def fetchone(self): return (1,)
def __enter__(self): return self
def __exit__(self, *a): return False
class _Conn:
def cursor(self): return _Cur()
def close(self): pass
monkeypatch.setattr(admin_mod.psycopg2, 'connect', lambda *a, **k: _Conn())
monkeypatch.setattr(admin_mod.subprocess, 'check_output', lambda *a, **k: 'osmium version 1.16.0\n')
def test_admin_info_auth_required(client):
assert client.get('/api/admin/navi-offroute/info').status_code == 401
def test_admin_info_no_secrets_and_probes(client, monkeypatch):
_mock_probes_ok(monkeypatch)
d = client.get('/api/admin/navi-offroute/info', headers=AUTH).get_json()
assert d['service'] == 'navi-offroute' and d['port'] == 8428
# No masked secrets anywhere (Phase A §10 — none exist; DSN is peer-auth).
assert all('...' not in str(e['value']) and e['value'] != '****' for e in d['env'])
assert all('password' not in e['name'].lower() for e in d['env'])
names = {dep['name'] for dep in d['dependencies']}
assert names == {'valhalla', 'padus-postgis', 'osmium-tool'}
# cheap file probes only (no row_count/size-of-db enrichment)
fs_names = {f['name'] for f in d['filesystem']}
assert {'dem', 'osm_pbf', 'navi_db', 'barriers_tif', 'wilderness_tif',
'trails_tif', 'friction_vrt'} == fs_names
assert all(set(f) == {'name', 'path', 'exists', 'readable'} for f in d['filesystem'])

View file

@ -0,0 +1,179 @@
"""
Trail corridor reader for OFFROUTE.
Provides access to the OSM-derived trail raster for pathfinding.
Trail values replace WorldCover friction where trails exist.
Raster values:
0 = no trail (use WorldCover friction)
5 = road (0.1× friction)
15 = track (0.3× friction)
25 = foot trail (0.5× friction)
"""
import os
from pathlib import Path
from typing import Tuple, Optional
import numpy as np
try:
import rasterio
from rasterio.windows import from_bounds
from rasterio.enums import Resampling
except ImportError:
raise ImportError("rasterio is required for trails layer support")
# Default path to the trails raster (single source of truth); env-overridable.
DEFAULT_TRAILS_PATH = Path("/mnt/nav/worldcover/trails.tif")
def trails_tif_path() -> Path:
"""Trails raster path, env-overridable via NAVI_OFFROUTE_TRAILS_TIF."""
return Path(os.environ.get("NAVI_OFFROUTE_TRAILS_TIF", str(DEFAULT_TRAILS_PATH)))
# Trail value to friction multiplier mapping
TRAIL_FRICTION_MAP = {
5: 0.1, # road
15: 0.3, # track
25: 0.5, # foot trail
}
class TrailReader:
"""Reader for OSM-derived trail corridor raster."""
def __init__(self, trails_path: Path = None):
self.trails_path = Path(trails_path) if trails_path else trails_tif_path()
self._dataset = None
def _open(self):
"""Lazy open the dataset."""
if self._dataset is None:
if not self.trails_path.exists():
raise FileNotFoundError(f"Trails raster not found at {self.trails_path}")
self._dataset = rasterio.open(self.trails_path)
return self._dataset
def get_trails_grid(
self,
south: float,
north: float,
west: float,
east: float,
target_shape: Tuple[int, int]
) -> np.ndarray:
"""
Get trail 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 trail values:
0 = no trail
5 = road (0.1× friction)
15 = track (0.3× friction)
25 = foot trail (0.5× friction)
"""
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
# Use nearest neighbor to preserve discrete values
trails = ds.read(
1,
window=window,
out_shape=target_shape,
resampling=Resampling.nearest
)
return trails
def sample_point(self, lat: float, lon: float) -> int:
"""Sample trail 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 = no trail
# Read single pixel
window = rasterio.windows.Window(col, row, 1, 1)
value = ds.read(1, window=window)
return int(value[0, 0])
def close(self):
"""Close the dataset."""
if self._dataset is not None:
self._dataset.close()
self._dataset = None
def trails_to_friction(trails: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
"""
Convert trail values to friction multipliers.
Args:
trails: uint8 array of trail values (0, 5, 15, or 25)
Returns:
Tuple of:
- friction: float32 array of friction multipliers
- has_trail: bool array indicating where trails exist
"""
friction = np.ones_like(trails, dtype=np.float32)
has_trail = trails > 0
# Apply friction values where trails exist
friction[trails == 5] = 0.1 # road
friction[trails == 15] = 0.3 # track
friction[trails == 25] = 0.5 # foot trail
return friction, has_trail
if __name__ == "__main__":
print("Testing TrailReader...")
if not DEFAULT_TRAILS_PATH.exists():
print(f"Trails raster not found at {DEFAULT_TRAILS_PATH}")
print("Run Phase B rasterization first.")
exit(1)
reader = TrailReader()
# Test point sampling - Twin Falls downtown (should have roads)
test_lat, test_lon = 42.563, -114.461
trail_value = reader.sample_point(test_lat, test_lon)
print(f"\nTwin Falls ({test_lat}, {test_lon}): trail value = {trail_value}")
label = {0: "no trail", 5: "road", 15: "track", 25: "trail"}.get(trail_value, "unknown")
print(f" Type: {label}")
# Test grid read for test bbox
trails = reader.get_trails_grid(
south=42.21, north=42.60, west=-114.76, east=-113.79,
target_shape=(400, 1000)
)
print(f"\nGrid test shape: {trails.shape}")
unique, counts = np.unique(trails, return_counts=True)
print("Value distribution:")
for v, c in zip(unique, counts):
pct = 100 * c / trails.size
label = {0: "no trail", 5: "road", 15: "track", 25: "trail"}.get(v, f"unknown({v})")
print(f" {label}: {c:,} pixels ({pct:.2f}%)")
# Test conversion to friction
friction, has_trail = trails_to_friction(trails)
print(f"\nTrail coverage: {100 * np.sum(has_trail) / trails.size:.2f}%")
print(f"Friction range (on trails): {friction[has_trail].min():.1f} - {friction[has_trail].max():.1f}")
reader.close()
print("\nTrailReader test complete.")