How We Build Scenic Bicycle Routes Using Nature Heterogeneity

The Problem with Shortest-Path Routing

Most bicycle route planners give you the shortest or fastest path between two points. On a bike, this is almost never what you want. A 60 km route that hugs a busy national road is “optimal” by distance, but miserable to ride. What cyclists actually want is a route that:

• avoids heavy traffic and urban sprawl
• passes through varied, interesting landscapes
• discovers quiet back roads you would never find on your own
• balances distance with enjoyment

This post explains how we built an algorithm that does exactly that, using landscape heterogeneity — the mathematical idea that a mosaic of different land covers is more interesting than a monoculture.

What Makes a Landscape “Interesting”?

Intuitively, a route that passes through forest, then meadow, then along a river, then through wetlands, is more interesting than 50 km of continuous pine forest. The first landscape has heterogeneity; the second has homogeneity.

Our core insight: we can measure this heterogeneity from OpenStreetMap data and use it as a routing signal.

The Six Nature Categories

From OSM tags, we extract six nature categories that cyclists actually care about:

Each category gets a weight reflecting its scarcity and sensory impact. Water features are weighted highest (2.0×) because they are rare, visually dominant, and strongly correlate with cooler temperatures and tailwinds. Wetlands get 1.8×. Forests get 1.0×. Meadows and scrub get lower but still positive weights.

Measuring Heterogeneity at a Point

For any location on the map, we compute a richness score by looking at nature features within an 800 m radius:

score = Σ category_weight × packing_factor

The packing factor captures how tightly features are clustered around the center point. A feature 50 m away contributes more than one 700 m away. Specifically:

weight = 1.0 - (distance / 800m)
packing_factor = 0.5 + average_weight

This means features at the very edge contribute 0.5× their category weight, while features clustered at the center contribute up to 1.5×. A cell with water at its center and forest at its edge scores higher than the reverse — which matches human intuition about what feels “scenic.”

The Minimum-Count Gate

A category only counts if it has at least 2 features in the radius. This prevents a single stray tree line from registering as “forest presence.” Categories with 4+ features get an extra 0.3 bonus, rewarding genuine biodiversity hotspots.

Water Bonus

If any water feature is present, we add a flat +1.0 to the score. Water is so strongly correlated with scenic quality that it deserves special treatment.

Sparse-Area Penalty

If the total feature count in the radius is below 5, the score is multiplied by 0.4. This suppresses false positives in genuinely empty agricultural areas where a single hedgerow might otherwise register as a mosaic.

From Points to Zones: Clustering Hotspots

Computing richness for every 300 m cell in a 50×50 km area gives us thousands of scores. The next step is finding hotspots — contiguous regions where heterogeneity is consistently high.

We use a two-step process:

Step 1: Grid Scoring

Compute richness for every cell in a 300 m grid. Cells scoring above a threshold (default 2.0) become candidate hotspots.

Step 2: Single-Link Clustering

Candidate hotspots are clustered using single-link agglomerative clustering with a 1 km chain distance. Two hotspots join the same zone if they are within 1 km of each other, or if a chain of hotspots exists connecting them with each link under 1 km.

Why single-link? Because it captures corridor-shaped scenic areas — like a river valley where heterogeneity is high all along the watercourse, even if individual hotspots are spaced a few hundred meters apart. Complete-link or centroid clustering would break these corridors into fragmented blobs.

Zone Properties

Each zone exposes:

• center_lat, center_lon: score-weighted centroid (not geometric center — high-scoring cells pull the centroid toward them)
• score: sum of member cell scores
• density: score / area (km²), so compact, intense zones rank above sprawling, weak ones
• top_categories: which nature categories dominate the zone
• peak_score: highest individual cell score in the zone

Zones are sorted by density descending. Singletons (zones with only 1 member) are filtered out.

The Density Field: From Waypoints to Corridors

Earlier versions of our algorithm treated zones as waypoints — the router would explicitly visit each zone center. This had a problem: if a zone center happened to be on a bad road, the route was forced to visit it anyway. And in sparse rural areas, the closest zone might consume the entire distance budget, leaving no room for exploration.

Our current approach uses a continuous density field instead of discrete waypoints.

Building the Field

We create a coarse 1 km grid over the entire route bbox. For each cell, its value is the weighted sum of nearby zone densities:

cell_value = Σ zone_density × (1.0 - distance/3000m)

Zones within 3 km contribute, with linear falloff. A cell near three moderate-density zones can have a higher value than a cell near one high-density zone — which is correct, because three overlapping scenic corridors are more interesting than one.

Using the Field in Routing

During graph weight preparation, every road edge gets a bonus proportional to the density field value at its midpoint:

bonus = min(density × 0.0125, 0.25) × adventurousness
weight *= 1.0 - bonus

At maximum adventurousness (1.0), edges near dense mosaic zones get up to 25% cheaper — which means A* will naturally route through scenic corridors without being forced to visit specific points. The route flows through interesting areas organically, like water finding a valley.

Candidate Diversity: Swing Waypoints

With the density field alone, every route at a given adventurousness would follow the same corridor. To give riders choices, we generate diverse candidates using geometric swing waypoints:

Each candidate is just A* from start → waypoint(s) → end. The density field ensures the legs between points naturally flow through scenic areas. Candidates are deduplicated by edge-set Jaccard similarity (threshold 0.5) so riders get genuinely different shapes, not minor variations.

Loop Routes: Same Idea, Different Geometry

For round-trip routes (start = end), we use a polar attractor field: 16 sectors × 10 radial rings around the start point. Each attractor is scored by scenic value × distance bonus. The algorithm then builds loop skeletons (ellipse, teardrop, figure-8, triangle) and routes through them with anti-reuse logic to avoid riding the same road twice.

The key insight is the same: density attracts, geometry diversifies.

Validation: Does It Actually Work?

We validate by comparing our adventurous routes against the shortest path on real OSM data. On a 62 km Włocławek → Toruń test route:

• Base (shortest): 60 km, mostly on national road 91, flat, noisy
• Adventurous: 62 km, weaves through 4 nature zones, passes 2 lakes, 1 wetland complex, and 3 villages. Elevation profile has 3× more variation.

The adventurous route is only 3% longer but visits 8 distinct land-cover types versus 2 for the base route. Rider feedback (informal, n=4) consistently prefers the adventurous route for weekend rides, though they note it requires more navigation attention.

Technical Stack

• OSM extraction: osmium-tool extracts bbox from Poland PBF (~2 GB)
• Graph: NetworkX MultiDiGraph with road-type penalties
• Spatial index: Shapely STRtree for fast proximity queries
• Scoring: Vectorized numpy/shapely bulk queries (30× faster than per-cell loops)
• Cache: Per-tile OPL + bbox-level graph+scenic cache, invalidated by content hash
• Elevation: SRTM3 tiles via Mapzen, bilinear interpolation

Future Work

• Industrial proximity penalty: Large factories and logistics parks create visual and air-quality dead zones that current scoring misses
• Seasonal weights: Forest is more valuable in summer (shade) than winter; wetlands are more valuable in autumn (bird migration)
• Surface-aware heterogeneity: A gravel road through meadow+forest might score higher than asphalt through the same landscape, because the surface matches the terrain

Go Rewrite & Performance Optimization

We recently rewrote the pipeline from Python to Go. The rewrite was driven by two goals: single-binary deployment and faster cold-start times. The Go version is now the production engine.

Why Go?

The Python version worked well but had friction: it required a virtualenv, had ~2 second import overhead, and the webapp needed a separate Python process for each job. A single Go binary eliminates all of that.

What we optimized

The initial Go rewrite was correct but slow — ~36 seconds for a cold run, only 3× faster than Python. Profiling revealed the real bottleneck was not A* routing (which takes ~40 ms) but weight preparation: computing scenic bonuses for every edge in the graph.

PrepareWeights iterates over all ~1 million edges and performs 5 spatial queries per edge (nature, scenic points, settlements, cities, busy roads). That’s 5 million R-tree searches per route request.

Fix: parallel edge processing. We split the edge list into 8 chunks and process each chunk in its own goroutine. Since edges are independent, this is embarrassingly parallel. On an Apple M4 Pro, this alone cut weight preparation from ~23 s to ~4 s.

Fix: dense node index. The A* implementation originally used map[int64]float64 for gScore and cameFrom. Maps are fast, but slice indexing is ~5× faster. We added a BuildIndex() method to the Graph that remaps sparse OSM node IDs to dense [0, n) indices. A* now uses []float64 and []int instead of maps.

Fix: cache key bug. When using -pbf input, the cache key was derived from the temp OPL file path (/tmp/cycle-route-*/...), which is different on every run. The cache never hit. We now key the cache by bbox hash, giving a 1.4× warm-run speedup.

Results

Cold run went from 16 s → 37 s because we moved from a tile-level msgpack cache to a bbox-level cache. The first cold run for a new bbox still needs osmium extraction (~29 s). This is expected — the win is on warm runs.

The warm run is now dominated by graph deserialization (~1 s) and weight preparation (~2.5 s). A* routing itself is ~40 ms — negligible.

Test Route Details

All benchmarks use the same real-world test case:

• Start: Włocławek, PL — 52.6482°N, 19.0672°E
• End: Toruń, PL — 53.0138°N, 18.5984°E
• Distance: ~60 km straight-line, ~62 km adventurous route
• OSM extract: Poland PBF (1.9 GB), extracted to 27 MB OPL for the bbox 18.5°E–19.2°E, 52.5°N–53.1°N
• Hardware: Apple M4 Pro, macOS, Go 1.24
• Profile: gravel, adventurousness 0.5

Detailed Timing Breakdown (Go final, warm run)

Weighted graph cache skips PrepareWeights entirely on repeat requests. The 1.3 s msgpack deserialization is now the bottleneck.

The A* router explores ~130k nodes (26.6% of the 490k-node graph) to find the 62 km adventurous path. With scenic weights, edge costs vary by up to 33% from physical length, so the Haversine heuristic is quite optimistic — but the route quality is worth the exploration cost.

Optimization history

Round 1 — correctness + quick wins:

• Dense node index ([]float64 gScore instead of map[int64]float64)
• Parallel PrepareWeights (8 goroutines, 23 s → 4.5 s)
• Precomputed A* heuristic (eliminates 300k+ repeated Haversine calls)
• Cache key fix (bbox hash instead of temp file path)

Round 2 — three targeted optimizations:

1. Grid-based scenic precomputation. Instead of 5 R-tree queries per edge at request time (5M queries total), we precompute scenic density on a 100 m grid once after graph build. PrepareWeights becomes a single bilinear lookup per edge. Route quality impact: negligible — the adventurous route changed by only 40 m (60.16 km vs 60.20 km).

2. Skip scenic weights for base route. When adventurousness=0, we skip PrepareWeights entirely and use PenaltyWeight (highway penalty only). In compare mode this saves the full weight-prep time for the base route.

3. Nearest-node R-tree. nearestNode() was a linear scan over all 490k nodes. We now build a node R-tree in BuildIndex() and do an expanding-radius search. Cuts nearest-node time from ~0.2 s to ~1 ms.

Round 3 — caching layers + visualization:

4. Route result cache. Hash (start_lat, start_lon, end_lat, end_lon, adventurousness, profile) → cache all output files. Repeat queries: <100 ms. Invalidated when PBF or code changes.

5. Weighted graph cache. Cache edge.Weight values after PrepareWeights (~9 MB msgpack). On hit, skip the entire 1.5 s scenic scoring step. Key includes graph hash + scoring params.

6. Bidirectional A.* Search from both ends simultaneously. On a 490k-node graph: 2.2× faster (66 ms → 30 ms per route). Added lazy ReverseEdgeList — only built when bidirectional search is actually used.

7. PNG candidate maps. Each candidate route gets a static PNG map with base + adventurous + candidate overlay. Viewer shows “🗺️ View map” links per candidate.

The Native PBF Parser Experiment (What We Learned)

We tried replacing the osmium subprocess with a native Go PBF parser using github.com/paulmach/osm. The parser worked for basic routing (~14 s, 2× faster than osmium) but failed to extract scenic data:

Why: Osmium has a C++ spatial index that can look up any node by ID in O(1) when resolving way geometries. Go streaming parsers read the file sequentially — to find 100 nodes referenced by a way, they must scan all 235M nodes or do multiple passes. Extracting full scenic data (nature areas, POIs, cities) requires 3+ passes over a 1.9 GB file, which exceeds 5 minutes.

Decision: Osmium stays as default. Native parser is available as -native-pbf for environments without osmium installed. For eliminating the subprocess, pre-extracting OPL once per region is the practical path.

The Profiling Lesson

We added graph pruning (removing disconnected components) and ReverseEdgeList building to BuildIndex(), thinking they were “free” optimizations. Profiling revealed they added 2–3 seconds to every single run:

• Pruning BFS on 490k nodes: ~2 s per run
• ReverseEdgeList scan over 1M edges: ~1 s per run
• Combined: warm run went from 2.9 s → 5.5 s

Rule: Never add O(n) or O(m) work to the hot path without measuring. Both were reverted. Pruning is only useful when building from raw OSM data (where disconnected ways exist), not when loading from a clean osmium extract.

What the Go rewrite replaced

What’s next

• Bidirectional A as default* — 2.2× routing speedup, zero overhead. Wire into CLI.
• Skip non-road ways during parse — filter buildings/barriers during OPL build. Estimated 60% graph shrink.
• zstd compression for msgpack cache — 220 MB → ~80 MB, faster I/O.
• WebSocket/SSE progress — replace polling with server-sent events for real-time job progress.
• Native PBF parser — available as -native-pbf for environments without osmium. Not default.

Summary

The core idea is simple: interesting cycling happens at the boundaries between ecosystems. By measuring nature heterogeneity from OSM data, turning it into a continuous density field, and letting the router flow through it like water through a watershed, we generate routes that are measurably more varied than shortest-path alternatives — without forcing artificial detours.

The algorithm is fully offline (no APIs), runs in ~2.9 seconds on a warm cache (was ~60 s in Python), and produces routes that cyclists actually want to ride.