Cross-Country Mobility Surface (Layered Cost Grid)

Status & scope

Module: lens/traversal/surface.py
Milestone: Phase 2

Purpose

A generic raster-based mobility surface where each cell has a composite traversal cost from stacked layers — roads, rivers, slope, vegetation, threat, satellite coverage, flooding. Routes via A* on the 8-connected grid. Same pattern as ESRI cross-country mobility analysis.

Design: Why Raster, Not Edge Annotation

Edge annotation on 267K OSM nodes is slow (Python loops) and fragile (coverage lookup misses). A raster grid has: - Predictable structure (8-connected, numpy-friendly) - Bulk operations (vectorized cost computation) - Always-admissible A* heuristic (grid distance) - Natural cost surface visualization (render as image/heatmap) - Composable layers (stack = element-wise weighted sum)

The 8 Levels

Each level adds a layer to the cost surface. All share one grid, one router, one evidence model.

Level	Layer	Cost Source	Primitive	Grid Effect
0	Base grid	Flat terrain	passthrough	Uniform cost
1	Roads	OSM road network	step_function	Road cells = low cost, off-road = high
2	Road speed	OSM speed tags	inverse_linear	Faster road = lower cost
3	Vehicle constraints	Weight/height/grade limits	threshold_gate	Blocked cells = inf
4	Terrain	DEM slope + land cover	linear_scale + step_function	Steep/dense = high cost
5	Hydrology	Rivers, lakes, flood extent	threshold_gate + decay	Water = inf, flood zone = ramp
6	Threat	Risk overlay (enemy, IED, jamming)	linear_scale	High threat = high cost
7	Satellite exposure	Coverage surface from TLEs	linear_scale	High visibility = high cost

Public API

`MobilitySurface(bbox, resolution_m) → MobilitySurface`

Creates an empty grid over the bounding box.
resolution_m: cell size in meters (default 500m for tactical, 5000m for strategic).
Grid stored as numpy arrays for fast bulk operations.

`surface.add_layer(name, weight, data, cost_fn) → None`

Adds a cost layer to the surface.
data: 2D numpy array (same shape as grid) or callable (lat, lon) → float.
cost_fn: one of the 10 primitives or a custom (value, params) → float.
weight: 0.0–1.0, all layers sum to 1.0.

`surface.add_roads(osm_path, speed_field) → None`

Rasterize OSM roads onto the grid.
Road cells get low base cost. Off-road cells get high base cost.
Optional speed data sets per-cell travel rate.

`surface.add_dem(dem_array, origin_lat, origin_lon, cell_size_m) → None`

Add slope cost from DEM. Uses Tobler's hiking function.
Resamples DEM to grid resolution if needed.

`surface.add_hydrology(water_mask) → None`

Water cells = impassable (cost inf). Flood zone = decay from water edge.

`surface.add_threat(threat_grid, weight) → None`

Overlay threat surface. Values 0.0–1.0 → scaled cost.

`surface.add_satellite_coverage(coverage_surface, weight) → None`

Overlay satellite coverage. High coverage = high cost (avoidance mode).

`surface.compute() → CostGrid`

Computes weighted composite cost per cell.
Returns CostGrid (2D numpy array + metadata).

`surface.route(origin_lat, origin_lon, dest_lat, dest_lon) → Route`

A* on the 8-connected cost grid.
Haversine heuristic (always admissible on lat/lon grid).
Returns Route with per-segment cost breakdown.

`surface.to_geotiff(path) → None`

Export cost surface as GeoTIFF (for GIS tools / Beacon heatmap).

`surface.to_png(path) → None`

Export as colored PNG (quick visualization).

Dataclasses

CostGrid (frozen=True)

Field	Type	Default	Notes
grid	np.ndarray	—	2D float array, each cell = composite cost
rows	int	—	Grid rows
cols	int	—	Grid cols
bbox	tuple[float,float,float,float]	—	(min_lat, min_lon, max_lat, max_lon)
resolution_m	float	—	Cell size in meters
layer_names	tuple[str, ...]	()	Which layers contributed
layer_weights	tuple[float, ...]	()	Weight per layer

MobilityCell (for per-cell detail when needed)

Field	Type	Default	Notes
row	int	—	Grid row
col	int	—	Grid col
lat	float	—	Cell center latitude
lon	float	—	Cell center longitude
cost	float	—	Composite cost
layer_costs	dict	{}	Per-layer cost breakdown
passable	bool	True	False if any layer = inf

Algorithm

Cost Surface Computation

For each cell (r, c):
    composite = 0
    for layer in layers:
        raw = layer.data[r, c]
        cost = layer.cost_fn(raw, layer.params)
        if cost == inf:
            composite = inf
            break
        composite += layer.weight * cost
    grid[r, c] = composite

A* on Grid

Standard A* with:
- Neighbors: 8-connected (cardinal + diagonal)
- Edge cost: avg(grid[current], grid[neighbor]) × distance
- Heuristic: haversine(current, goal) / max_speed
  (admissible because actual cost ≥ straight-line distance / max_speed)
- Diagonal distance: cell_size × sqrt(2)

Road Rasterization

OSM roads → grid cells: 1. For each road segment, Bresenham line rasterization onto grid 2. Road cells get base cost from speed: cell_size_m / (speed_kmh / 3.6) 3. Off-road cells get high base cost: cell_size_m / (walking_speed / 3.6) 4. Highway type sets speed: motorway=110, trunk=90, primary=70, secondary=50, tertiary=40 5. Bridge cells get concealment flag 6. Tunnel cells get full_concealment flag

Test Fixtures

Test	Input	Expected	File
L0: Empty grid routes	20×20 flat grid	Route exists, cost = distance	tests/test_surface.py
L1: Road preference	Grid with road strip	Route follows road	tests/test_surface.py
L2: Speed preference	Two roads, diff speed	Route takes faster road	tests/test_surface.py
L3: Vehicle blocked	Bridge weight < vehicle	Route avoids bridge	tests/test_surface.py
L4: Slope avoidance	Ridge across path	Route goes around ridge	tests/test_surface.py
L5: River blocking	River across path	Route finds bridge/ford	tests/test_surface.py
L6: Threat avoidance	Threat zone on direct path	Route deviates around threat	tests/test_surface.py
L7: Satellite avoidance	High coverage on direct path	Route prefers low coverage	tests/test_surface.py
Composite: All layers	All 8 layers stacked	Route balances all costs	tests/test_surface.py
Export: GeoTIFF	Computed surface	Valid GeoTIFF file	tests/test_surface.py
Export: PNG	Computed surface	Valid PNG file	tests/test_surface.py
Evidence: Block created	Route computed	Frozen evidence with hashes	tests/test_surface.py

File Layout

lens/traversal/
  surface.py           ← MobilitySurface, CostGrid, A* on grid

Integration Points

Cortex: execute_monitor_block with type: cost_surface renders as heatmap
Athena: LinearScale, StepFunction, ThresholdGate ops compute per-cell costs
Beacon: CostGrid → heat_map MapLayerDef, Route → linestring layer
ESRI: to_geotiff() exports for ArcGIS/QGIS interop

DO NOT

Build a separate router — reuse A* pattern from router.py, adapted for grid
Store the full grid in evidence blocks — store hash + metadata, not the raster
Use floating-point grid indices — integer row/col only
Forget diagonal cost — diagonal neighbors are sqrt(2) × cell_size, not 1×
Over-resolve — 500m cells for tactical (600km = 1200×400 = 480K cells), 5km for strategic

Depends on: component.prism.cost-primitives, component.prism.universal-lens-parser

Required by: component.prism.scenario-runner