Skip to main content

trailgen_core/
enrich.rs

1use crate::geo::{Coord, LineString};
2use crate::hiking::HikingModel;
3use crate::model::{Edge, GradeDistribution, Provenance, Terrain, TerrainEvidence, WalkGraph};
4use crate::{Result, TrailgenError};
5use serde::{Deserialize, Serialize};
6
7#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
8pub struct EnrichmentConfig {
9    pub sample_spacing_m: f64,
10    pub steep_grade_threshold: f64,
11}
12
13impl Default for EnrichmentConfig {
14    fn default() -> Self {
15        Self {
16            sample_spacing_m: 25.0,
17            steep_grade_threshold: 0.15,
18        }
19    }
20}
21
22#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
23pub struct ElevationSample {
24    pub ele_m: f64,
25    pub confidence: f64,
26    pub provenance: Provenance,
27}
28
29pub trait ElevationSampler {
30    fn sample(&self, coord: Coord) -> Option<ElevationSample>;
31}
32
33#[derive(Clone, Debug, Eq, PartialEq, Serialize, Deserialize)]
34pub struct ElevationMosaic<S> {
35    samplers: Vec<S>,
36}
37
38impl<S> ElevationMosaic<S> {
39    pub fn new(samplers: Vec<S>) -> Result<Self> {
40        if samplers.is_empty() {
41            return Err(TrailgenError::InvalidData(
42                "elevation mosaic requires at least one sampler".to_owned(),
43            ));
44        }
45        Ok(Self { samplers })
46    }
47}
48
49impl<S: ElevationSampler> ElevationSampler for ElevationMosaic<S> {
50    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
51        self.samplers
52            .iter()
53            .find_map(|sampler| sampler.sample(coord))
54    }
55}
56
57#[derive(Clone, Copy, Debug, Default, Eq, PartialEq)]
58pub struct EmbeddedElevation;
59
60impl ElevationSampler for EmbeddedElevation {
61    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
62        Some(ElevationSample {
63            ele_m: coord.ele?,
64            confidence: 0.85,
65            provenance: Provenance {
66                source: "embedded-geometry-elevation".to_owned(),
67                layer: None,
68                source_id: None,
69                license: None,
70            },
71        })
72    }
73}
74
75#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
76pub struct PlaneElevation {
77    pub origin: Coord,
78    pub origin_ele_m: f64,
79    pub east_gain_m_per_degree: f64,
80    pub north_gain_m_per_degree: f64,
81    pub confidence: f64,
82}
83
84impl ElevationSampler for PlaneElevation {
85    fn sample(&self, coord: Coord) -> Option<ElevationSample> {
86        Some(ElevationSample {
87            ele_m: self.north_gain_m_per_degree.mul_add(
88                coord.lat - self.origin.lat,
89                self.east_gain_m_per_degree
90                    .mul_add(coord.lon - self.origin.lon, self.origin_ele_m),
91            ),
92            confidence: self.confidence,
93            provenance: Provenance {
94                source: "synthetic-plane-elevation".to_owned(),
95                layer: Some("fixture".to_owned()),
96                source_id: None,
97                license: Some("CC0-fixture".to_owned()),
98            },
99        })
100    }
101}
102
103pub fn enrich_graph<S: ElevationSampler>(
104    graph: &mut WalkGraph,
105    sampler: &S,
106    config: EnrichmentConfig,
107) -> Result<()> {
108    if config.sample_spacing_m <= 0.0 {
109        return Err(TrailgenError::InvalidData(
110            "enrichment sample spacing must be positive".to_owned(),
111        ));
112    }
113    for edge in &mut graph.edges {
114        enrich_edge(edge, sampler, config)?;
115    }
116    Ok(())
117}
118
119fn enrich_edge<S: ElevationSampler>(
120    edge: &mut Edge,
121    sampler: &S,
122    config: EnrichmentConfig,
123) -> Result<()> {
124    let (sampled_line, elevation_provenance, elevation_confidence) =
125        densify_and_sample(&edge.geometry, sampler, config.sample_spacing_m)?;
126    let profile = grade_profile(&sampled_line, config.steep_grade_threshold);
127    edge.geometry = sampled_line;
128    edge.attr.length_m = edge.geometry.length_m();
129    edge.attr.ascent_m = profile.ascent_m;
130    edge.attr.descent_m = profile.descent_m;
131    edge.attr.grade_abs_mean = profile.grade_abs_mean;
132    edge.attr.grade_abs_max = profile.grade_abs_max;
133    edge.attr.sustained_steep_m = profile.sustained_steep_m;
134    edge.attr.grade_distribution = profile.grade_distribution;
135    edge.attr.hill_slope_deg = mean_hill_slope_deg(&edge.geometry, sampler);
136    edge.attr.elevation_provenance = elevation_provenance;
137    if elevation_confidence > 0.0 {
138        edge.attr.confidence = edge.attr.confidence.min(elevation_confidence);
139    }
140    infer_terrain(edge, &profile);
141    HikingModel.apply(edge);
142    Ok(())
143}
144
145fn mean_hill_slope_deg<S: ElevationSampler>(line: &LineString, sampler: &S) -> Option<f64> {
146    const RADIUS_M: f64 = 10.0;
147    const METERS_PER_LATITUDE_DEGREE: f64 = 111_195.0;
148    let mut weighted_slope = 0.0;
149    let mut measured_m = 0.0;
150    for segment in line.points.windows(2) {
151        let distance_m = segment[0].haversine_m(segment[1]);
152        if distance_m <= f64::EPSILON {
153            continue;
154        }
155        let center = segment[0].lerp(segment[1], 0.5);
156        let latitude_delta = RADIUS_M / METERS_PER_LATITUDE_DEGREE;
157        let longitude_delta =
158            RADIUS_M / (METERS_PER_LATITUDE_DEGREE * center.lat.to_radians().cos().max(0.01));
159        let [Some(west), Some(east), Some(south), Some(north)] = [
160            sampler.sample(Coord::new(center.lon - longitude_delta, center.lat)),
161            sampler.sample(Coord::new(center.lon + longitude_delta, center.lat)),
162            sampler.sample(Coord::new(center.lon, center.lat - latitude_delta)),
163            sampler.sample(Coord::new(center.lon, center.lat + latitude_delta)),
164        ] else {
165            continue;
166        };
167        let east_grade = (east.ele_m - west.ele_m) / (2.0 * RADIUS_M);
168        let north_grade = (north.ele_m - south.ele_m) / (2.0 * RADIUS_M);
169        let slope_deg = east_grade.hypot(north_grade).atan().to_degrees();
170        weighted_slope = slope_deg.mul_add(distance_m, weighted_slope);
171        measured_m += distance_m;
172    }
173    (measured_m > 0.0).then_some(weighted_slope / measured_m)
174}
175
176fn densify_and_sample<S: ElevationSampler>(
177    line: &LineString,
178    sampler: &S,
179    spacing_m: f64,
180) -> Result<(LineString, Vec<Provenance>, f64)> {
181    let mut points = Vec::new();
182    let mut provenance = Vec::new();
183    let mut confidence = 1.0;
184    let mut sample_count = 0u32;
185    let mut sample_attempts = 0u32;
186    for segment in line.points.windows(2) {
187        let a = segment[0];
188        let b = segment[1];
189        if points.is_empty() {
190            points.push(sample_coord(
191                a,
192                sampler,
193                &mut provenance,
194                &mut confidence,
195                &mut sample_count,
196                &mut sample_attempts,
197            ));
198        }
199        let length_m = a.haversine_m(b);
200        for distance_m in std::iter::successors(Some(spacing_m), move |d| Some(d + spacing_m))
201            .take_while(|d| *d < length_m)
202        {
203            points.push(sample_coord(
204                a.lerp(b, distance_m / length_m),
205                sampler,
206                &mut provenance,
207                &mut confidence,
208                &mut sample_count,
209                &mut sample_attempts,
210            ));
211        }
212        points.push(sample_coord(
213            b,
214            sampler,
215            &mut provenance,
216            &mut confidence,
217            &mut sample_count,
218            &mut sample_attempts,
219        ));
220    }
221    if sample_count == 0 {
222        return Ok((line.clone(), Vec::new(), 0.0));
223    }
224    if sample_count > 0 && sample_count < sample_attempts {
225        confidence = confidence.min(f64::from(sample_count) / f64::from(sample_attempts));
226    }
227    Ok((LineString::new(points)?, provenance, confidence))
228}
229
230fn sample_coord<S: ElevationSampler>(
231    coord: Coord,
232    sampler: &S,
233    provenance: &mut Vec<Provenance>,
234    confidence: &mut f64,
235    sample_count: &mut u32,
236    sample_attempts: &mut u32,
237) -> Coord {
238    *sample_attempts += 1;
239    sampler
240        .sample(coord)
241        .map_or(Coord { ele: None, ..coord }, |sample| {
242            *sample_count += 1;
243            *confidence = confidence.min(sample.confidence);
244            if !provenance.contains(&sample.provenance) {
245                provenance.push(sample.provenance);
246            }
247            Coord {
248                ele: Some(sample.ele_m),
249                ..coord
250            }
251        })
252}
253
254#[derive(Clone, Copy)]
255struct GradeProfile {
256    ascent_m: f64,
257    descent_m: f64,
258    grade_abs_mean: f64,
259    grade_abs_max: f64,
260    sustained_steep_m: f64,
261    grade_distribution: GradeDistribution,
262}
263
264fn grade_profile(line: &LineString, steep_grade_threshold: f64) -> GradeProfile {
265    let mut ascent_m = 0.0;
266    let mut descent_m = 0.0;
267    let mut graded_m = 0.0;
268    let mut weighted_abs_grade = 0.0;
269    let mut grade_abs_max = 0.0;
270    let mut sustained_steep_m = 0.0;
271    let mut grade_distribution = GradeDistribution::default();
272
273    for segment in line.points.windows(2) {
274        let a = segment[0];
275        let b = segment[1];
276        let distance = a.haversine_m(b);
277        let Some(ele_a) = a.ele else {
278            continue;
279        };
280        let Some(ele_b) = b.ele else {
281            continue;
282        };
283        let rise = ele_b - ele_a;
284        if rise >= 0.0 {
285            ascent_m += rise;
286        } else {
287            descent_m -= rise;
288        }
289        if distance > 0.0 {
290            let abs_grade = (rise / distance).abs();
291            graded_m += distance;
292            weighted_abs_grade = abs_grade.mul_add(distance, weighted_abs_grade);
293            if abs_grade > grade_abs_max {
294                grade_abs_max = abs_grade;
295            }
296            if abs_grade >= steep_grade_threshold {
297                sustained_steep_m += distance;
298            }
299            grade_distribution = grade_distribution.add_segment(distance, abs_grade);
300        }
301    }
302
303    GradeProfile {
304        ascent_m,
305        descent_m,
306        grade_abs_mean: weighted_abs_grade / graded_m.max(1.0),
307        grade_abs_max,
308        sustained_steep_m,
309        grade_distribution,
310    }
311}
312
313struct TerrainInference {
314    terrain: Terrain,
315    confidence: f64,
316    rationale: String,
317    provenance: Option<Provenance>,
318}
319
320fn infer_terrain(edge: &mut Edge, profile: &GradeProfile) {
321    let evicted_enrichment_evidence = evict_enrichment_terrain_evidence(edge);
322    if evicted_enrichment_evidence
323        && edge.attr.terrain_confidence < 0.85
324        && !edge
325            .attr
326            .terrain_evidence
327            .iter()
328            .any(|e| e.terrain == edge.attr.terrain)
329    {
330        edge.attr.terrain = Terrain::Unknown;
331        edge.attr.terrain_confidence = 0.0;
332    }
333
334    let terrain = edge.attr.terrain;
335    let explicit_confidence = edge.attr.terrain_confidence;
336    if terrain != Terrain::Unknown && explicit_confidence >= 0.85 {
337        if !edge
338            .attr
339            .terrain_evidence
340            .iter()
341            .any(|e| e.terrain == terrain)
342        {
343            upsert_terrain_evidence(
344                edge,
345                TerrainEvidence {
346                    terrain,
347                    confidence: explicit_confidence,
348                    rationale: "explicit source terrain tag".to_owned(),
349                    provenance: edge.attr.provenance.first().cloned(),
350                },
351            );
352        }
353        edge.attr.terrain_confidence = edge.attr.terrain_confidence.max(explicit_confidence);
354        return;
355    }
356
357    let inference = terrain_inference(edge, profile);
358    if inference.confidence < edge.attr.terrain_confidence {
359        return;
360    }
361    edge.attr.terrain = inference.terrain;
362    edge.attr.terrain_confidence = inference.confidence;
363    edge.attr.confidence = edge
364        .attr
365        .confidence
366        .min(inference.confidence.mul_add(0.35, 0.65));
367    upsert_terrain_evidence(
368        edge,
369        TerrainEvidence {
370            terrain: inference.terrain,
371            confidence: inference.confidence,
372            rationale: inference.rationale,
373            provenance: inference.provenance,
374        },
375    );
376}
377
378fn terrain_inference(edge: &Edge, profile: &GradeProfile) -> TerrainInference {
379    let grade_basis = grade_basis(edge, profile);
380    let edge_provenance = edge.attr.provenance.first().cloned();
381    if edge.attr.road_exposure >= 0.75 {
382        TerrainInference {
383            terrain: Terrain::Road,
384            confidence: 0.70,
385            rationale: format!(
386                "inferred from road exposure {:.0}% with {grade_basis}",
387                edge.attr.road_exposure * 100.0
388            ),
389            provenance: edge
390                .attr
391                .crossings
392                .iter()
393                .find(|x| x.kind == crate::model::CrossingKind::Road)
394                .map(|x| x.provenance.clone())
395                .or(edge_provenance),
396        }
397    } else if profile.grade_abs_max >= 0.32 {
398        TerrainInference {
399            terrain: Terrain::Scramble,
400            confidence: 0.52,
401            rationale: format!("inferred from savage sampled grade: {grade_basis}"),
402            provenance: edge
403                .attr
404                .elevation_provenance
405                .first()
406                .cloned()
407                .or(edge_provenance),
408        }
409    } else if profile.grade_abs_max >= 0.22 {
410        TerrainInference {
411            terrain: Terrain::Talus,
412            confidence: 0.45,
413            rationale: format!("inferred from steep sampled grade: {grade_basis}"),
414            provenance: edge
415                .attr
416                .elevation_provenance
417                .first()
418                .cloned()
419                .or(edge_provenance),
420        }
421    } else if edge.attr.way_kind.pathless() {
422        TerrainInference {
423            terrain: Terrain::Unknown,
424            confidence: 0.0,
425            rationale: format!("no surrounding-terrain evidence for pathless route: {grade_basis}"),
426            provenance: edge_provenance,
427        }
428    } else {
429        TerrainInference {
430            terrain: Terrain::Trail,
431            confidence: 0.35,
432            rationale: format!("inferred default hiking surface: {grade_basis}"),
433            provenance: edge_provenance,
434        }
435    }
436}
437
438fn grade_basis(edge: &Edge, profile: &GradeProfile) -> String {
439    if edge.attr.elevation_provenance.is_empty() {
440        format!(
441            "no sampled elevation source; road exposure {:.0}%",
442            edge.attr.road_exposure * 100.0
443        )
444    } else {
445        let bins = profile.grade_distribution;
446        let total = bins.total_m().max(1.0);
447        format!(
448            "max grade {:.1}%, mean grade {:.1}%, steep {:.0}%, savage {:.0}%",
449            profile.grade_abs_max * 100.0,
450            profile.grade_abs_mean * 100.0,
451            bins.steep_m / total * 100.0,
452            bins.savage_m / total * 100.0
453        )
454    }
455}
456
457fn evict_enrichment_terrain_evidence(edge: &mut Edge) -> bool {
458    let old_len = edge.attr.terrain_evidence.len();
459    edge.attr.terrain_evidence.retain(|e| {
460        !(e.rationale.starts_with("inferred ")
461            || matches!(
462                e.rationale.as_str(),
463                "high road-exposure fraction"
464                    | "very steep sampled grade"
465                    | "steep sampled grade"
466                    | "default low-confidence hiking surface"
467            ))
468    });
469    edge.attr.terrain_evidence.len() != old_len
470}
471
472fn upsert_terrain_evidence(edge: &mut Edge, evidence: TerrainEvidence) {
473    if let Some(existing) = edge.attr.terrain_evidence.iter_mut().find(|e| {
474        e.terrain == evidence.terrain
475            && e.rationale == evidence.rationale
476            && e.provenance == evidence.provenance
477    }) {
478        existing.confidence = existing.confidence.max(evidence.confidence);
479    } else {
480        edge.attr.terrain_evidence.push(evidence);
481    }
482}