split and snap

src/cli.rs

@@ -0,0 +1,36 @@
+use clap::Parser;
+use std::path::PathBuf;
+
+#[derive(Parser, Debug)]
+#[command(
+    name = "borderpoi-rs",
+    version,
+    about = "Find administrative border crossings along a GPX track"
+)]
+pub struct Args {
+    /// Input GPX track. Reads from stdin when omitted or set to '-'.
+    #[arg(short, long)]
+    pub track: Option<PathBuf>,
+
+    /// Boundary GeoJSON in WGS84 / EPSG:4326.
+    #[arg(short, long)]
+    pub boundaries: PathBuf,
+
+    /// GeoJSON property containing the county name.
+    #[arg(long, default_value = "GEN")]
+    pub name_field: String,
+
+    /// GeoJSON property containing the administrative identifier.
+    #[arg(long, default_value = "ARS")]
+    pub id_field: String,
+
+    /// Output GPX. Writes to stdout when omitted or set to '-'.
+    #[arg(short, long)]
+    pub output: Option<PathBuf>,
+}
+
+impl Args {
+    pub fn parse_args() -> Self {
+        Self::parse()
+    }
+}

src/county.rs

@@ -0,0 +1,110 @@
+use anyhow::{bail, Context, Result};
+use geo::{algorithm::bounding_rect::BoundingRect, Geometry};
+use geojson::{Feature, GeoJson};
+use rstar::{RTreeObject, AABB};
+use serde_json::Value;
+use std::{fs::File, io::BufReader, path::Path};
+
+#[derive(Debug, Clone)]
+pub struct County {
+    pub name: String,
+    pub id: String,
+    pub geometry: Geometry<f64>,
+    envelope: AABB<[f64; 2]>,
+}
+
+impl RTreeObject for County {
+    type Envelope = AABB<[f64; 2]>;
+
+    fn envelope(&self) -> Self::Envelope {
+        self.envelope
+    }
+}
+
+pub fn load_counties(path: &Path, name_field: &str, id_field: &str) -> Result<Vec<County>> {
+    let file =
+        File::open(path).with_context(|| format!("Failed to open GeoJSON: {}", path.display()))?;
+
+    let geojson: GeoJson = serde_json::from_reader(BufReader::new(file))
+        .context("Failed to parse boundary GeoJSON")?;
+
+    let features = match geojson {
+        GeoJson::FeatureCollection(collection) => collection.features,
+        GeoJson::Feature(feature) => {
+            vec![feature]
+        }
+        GeoJson::Geometry(_) => {
+            bail!("Boundary GeoJSON must be a FeatureCollection or Feature.");
+        }
+    };
+
+    let mut counties = Vec::with_capacity(features.len());
+
+    for feature in features {
+        counties.push(feature_to_county(&feature, name_field, id_field)?);
+    }
+
+    if counties.is_empty() {
+        bail!("No boundary features were found.");
+    }
+
+    Ok(counties)
+}
+
+fn feature_to_county(feature: &Feature, name_field: &str, id_field: &str) -> Result<County> {
+    let properties = feature
+        .properties
+        .as_ref()
+        .context("Boundary feature has no properties.")?;
+
+    let name = property_as_string(
+        properties
+            .get(name_field)
+            .with_context(|| format!("Missing name field '{}'.", name_field))?,
+    )?;
+
+    let id = property_as_string(
+        properties
+            .get(id_field)
+            .with_context(|| format!("Missing ID field '{}'.", id_field))?,
+    )?;
+
+    let geojson_geometry = feature
+        .geometry
+        .as_ref()
+        .context("Boundary feature has no geometry.")?;
+
+    let geometry: Geometry<f64> = geojson_geometry
+        .try_into()
+        .context("Failed to convert GeoJSON geometry.")?;
+
+    match &geometry {
+        Geometry::Polygon(_) | Geometry::MultiPolygon(_) => {}
+        _ => {
+            bail!("Feature '{}' is not a Polygon or MultiPolygon.", name);
+        }
+    }
+
+    let bbox = geometry
+        .bounding_rect()
+        .with_context(|| format!("Feature '{}' has no bounding box.", name))?;
+
+    let envelope = AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]);
+
+    Ok(County {
+        name,
+        id,
+        geometry,
+        envelope,
+    })
+}
+
+fn property_as_string(value: &Value) -> Result<String> {
+    match value {
+        Value::String(value) => Ok(value.clone()),
+        Value::Number(value) => Ok(value.to_string()),
+        _ => {
+            bail!("Expected string or number property, got {}.", value);
+        }
+    }
+}

src/crossings.rs

@@ -0,0 +1,249 @@
+use anyhow::{Context, Result};
+use geo::{algorithm::bounding_rect::BoundingRect, Contains, Coord, Geometry, LineString, Point};
+use rstar::{RTree, AABB};
+use std::cmp::Ordering;
+
+use crate::{
+    county::County,
+    geometry::{deduplicate_hits, interpolate, segment_intersection, snap_to_segment},
+};
+
+#[derive(Debug, Clone)]
+pub struct BoundaryHit {
+    pub position: f64,
+}
+
+#[derive(Debug, Clone)]
+struct Transition {
+    position: f64,
+    from: County,
+    to: County,
+}
+
+#[derive(Debug, Clone)]
+pub struct Crossing {
+    pub point: Point<f64>,
+    pub from: County,
+    pub to: County,
+    pub segment_index: usize,
+    pub position: f64,
+}
+
+#[derive(Debug, Clone, Copy)]
+pub struct TrackPoint {
+    pub lon: f64,
+    pub lat: f64,
+}
+
+/// Find all administrative boundary crossings along a GPX track.
+///
+/// The original track is never modified. Crossing points are reconstructed
+/// from the original track segment and the intersection parameter, ensuring
+/// that the resulting waypoint lies exactly on the original segment.
+pub fn find_crossings(track: &[TrackPoint], tree: &RTree<County>) -> Result<Vec<Crossing>> {
+    let mut crossings = Vec::new();
+
+    for (segment_index, pair) in track.windows(2).enumerate() {
+        let start = Point::new(pair[0].lon, pair[0].lat);
+        let end = Point::new(pair[1].lon, pair[1].lat);
+
+        let segment = LineString::from(vec![
+            Coord {
+                x: start.x(),
+                y: start.y(),
+            },
+            Coord {
+                x: end.x(),
+                y: end.y(),
+            },
+        ]);
+
+        let bbox = segment
+            .bounding_rect()
+            .context("Track segment has no bounding box.")?;
+
+        let envelope =
+            AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]);
+
+        let candidates: Vec<County> = tree
+            .locate_in_envelope_intersecting(envelope)
+            .cloned()
+            .collect();
+
+        if candidates.is_empty() {
+            continue;
+        }
+
+        let hits = collect_segment_hits(start, end, &candidates);
+
+        if hits.is_empty() {
+            continue;
+        }
+
+        let transitions = reconstruct_transitions(start, end, &candidates, &hits);
+
+        for transition in transitions {
+            // Reconstruct the point from the original GPX segment.
+            // This is the actual snapping step.
+            let snapped_point = snap_to_segment(start, end, transition.position);
+
+            crossings.push(Crossing {
+                point: snapped_point,
+                from: transition.from,
+                to: transition.to,
+                segment_index,
+                position: transition.position,
+            });
+        }
+    }
+
+    Ok(deduplicate_crossings(crossings))
+}
+
+fn collect_segment_hits(
+    start: Point<f64>,
+    end: Point<f64>,
+    candidates: &[County],
+) -> Vec<BoundaryHit> {
+    let start_coord = Coord {
+        x: start.x(),
+        y: start.y(),
+    };
+
+    let end_coord = Coord {
+        x: end.x(),
+        y: end.y(),
+    };
+
+    let mut hits = Vec::new();
+
+    for county in candidates {
+        collect_geometry_intersections(&county.geometry, start_coord, end_coord, &mut hits);
+    }
+
+    deduplicate_hits(hits)
+}
+
+fn collect_geometry_intersections(
+    geometry: &Geometry<f64>,
+    start: Coord<f64>,
+    end: Coord<f64>,
+    hits: &mut Vec<BoundaryHit>,
+) {
+    match geometry {
+        Geometry::Polygon(polygon) => {
+            collect_ring_intersections(polygon.exterior(), start, end, hits);
+
+            for interior in polygon.interiors() {
+                collect_ring_intersections(interior, start, end, hits);
+            }
+        }
+
+        Geometry::MultiPolygon(multipolygon) => {
+            for polygon in &multipolygon.0 {
+                collect_ring_intersections(polygon.exterior(), start, end, hits);
+
+                for interior in polygon.interiors() {
+                    collect_ring_intersections(interior, start, end, hits);
+                }
+            }
+        }
+
+        _ => {}
+    }
+}
+
+fn collect_ring_intersections(
+    ring: &LineString<f64>,
+    start: Coord<f64>,
+    end: Coord<f64>,
+    hits: &mut Vec<BoundaryHit>,
+) {
+    for edge in ring.lines() {
+        if let Some((position, _point)) = segment_intersection(start, end, edge.start, edge.end) {
+            hits.push(BoundaryHit { position });
+        }
+    }
+}
+
+fn reconstruct_transitions(
+    start: Point<f64>,
+    end: Point<f64>,
+    candidates: &[County],
+    hits: &[BoundaryHit],
+) -> Vec<Transition> {
+    let mut transitions = Vec::new();
+
+    for hit in hits {
+        // Sample slightly before and after the intersection.
+        //
+        // This lets us determine whether the track actually changes
+        // administrative area rather than merely touching a boundary.
+        let before_t = (hit.position - 1e-8).max(0.0);
+        let after_t = (hit.position + 1e-8).min(1.0);
+
+        let before = interpolate(start, end, before_t);
+        let after = interpolate(start, end, after_t);
+
+        let from = county_at_point(before, candidates);
+        let to = county_at_point(after, candidates);
+
+        let (Some(from), Some(to)) = (from, to) else {
+            continue;
+        };
+
+        // Ignore boundary touches where the track remains in the same county.
+        if from.id == to.id {
+            continue;
+        }
+
+        transitions.push(Transition {
+            position: hit.position,
+            from,
+            to,
+        });
+    }
+
+    transitions
+}
+
+fn county_at_point(point: Point<f64>, candidates: &[County]) -> Option<County> {
+    candidates
+        .iter()
+        .find(|county| county.geometry.contains(&point))
+        .cloned()
+}
+
+fn deduplicate_crossings(mut crossings: Vec<Crossing>) -> Vec<Crossing> {
+    crossings.sort_by(|a, b| {
+        a.segment_index.cmp(&b.segment_index).then_with(|| {
+            a.position
+                .partial_cmp(&b.position)
+                .unwrap_or(Ordering::Equal)
+        })
+    });
+
+    let mut result = Vec::new();
+
+    for crossing in crossings {
+        let duplicate = result.iter().any(|existing: &Crossing| {
+            existing.segment_index == crossing.segment_index
+                && existing.from.id == crossing.from.id
+                && existing.to.id == crossing.to.id
+                && same_point(existing.point, crossing.point)
+        });
+
+        if !duplicate {
+            result.push(crossing);
+        }
+    }
+
+    result
+}
+
+fn same_point(a: Point<f64>, b: Point<f64>) -> bool {
+    let dx = a.x() - b.x();
+    let dy = a.y() - b.y();
+
+    dx * dx + dy * dy < 1e-20
+}

src/geometry.rs

@@ -0,0 +1,87 @@
+use geo::{Coord, Point};
+
+use crate::crossings::BoundaryHit;
+
+pub fn cross(a: Coord<f64>, b: Coord<f64>) -> f64 {
+    a.x * b.y - a.y * b.x
+}
+
+pub fn interpolate(start: Point<f64>, end: Point<f64>, t: f64) -> Point<f64> {
+    Point::new(
+        start.x() + t * (end.x() - start.x()),
+        start.y() + t * (end.y() - start.y()),
+    )
+}
+
+pub fn segment_intersection(
+    p: Coord<f64>,
+    p2: Coord<f64>,
+    q: Coord<f64>,
+    q2: Coord<f64>,
+) -> Option<(f64, Point<f64>)> {
+    let r = Coord {
+        x: p2.x - p.x,
+        y: p2.y - p.y,
+    };
+
+    let s = Coord {
+        x: q2.x - q.x,
+        y: q2.y - q.y,
+    };
+
+    let denominator = cross(r, s);
+
+    if denominator.abs() < 1e-14 {
+        return None;
+    }
+
+    let qp = Coord {
+        x: q.x - p.x,
+        y: q.y - p.y,
+    };
+
+    let t = cross(qp, s) / denominator;
+    let u = cross(qp, r) / denominator;
+
+    const EPSILON: f64 = 1e-10;
+
+    if !(-EPSILON..=1.0 + EPSILON).contains(&t) || !(-EPSILON..=1.0 + EPSILON).contains(&u) {
+        return None;
+    }
+
+    let t = t.clamp(0.0, 1.0);
+
+    let point = Point::new(p.x + t * r.x, p.y + t * r.y);
+
+    Some((t, point))
+}
+
+pub fn deduplicate_hits(mut hits: Vec<BoundaryHit>) -> Vec<BoundaryHit> {
+    hits.sort_by(|a, b| {
+        a.position
+            .partial_cmp(&b.position)
+            .unwrap_or(std::cmp::Ordering::Equal)
+    });
+
+    let mut result = Vec::new();
+
+    for hit in hits {
+        let duplicate = result
+            .iter()
+            .any(|existing: &BoundaryHit| (existing.position - hit.position).abs() < 1e-9);
+
+        if !duplicate {
+            result.push(hit);
+        }
+    }
+
+    result
+}
+
+/// Recalculate the point directly from the original track segment.
+///
+/// This guarantees that the returned point lies on the segment rather than
+/// relying on the coordinates calculated during the intersection operation.
+pub fn snap_to_segment(start: Point<f64>, end: Point<f64>, position: f64) -> Point<f64> {
+    interpolate(start, end, position.clamp(0.0, 1.0))
+}

src/gpx_io.rs

@@ -0,0 +1,92 @@
+use anyhow::{bail, Context, Result};
+use gpx::{read, write, Gpx, Waypoint};
+use std::{
+    fs::File,
+    io::{self, BufReader, BufWriter},
+    path::Path,
+};
+
+use crate::crossings::{Crossing, TrackPoint};
+
+pub fn read_gpx(path: Option<&Path>) -> Result<Gpx> {
+    match path {
+        Some(path) if path.as_os_str() != "-" => {
+            let file = File::open(path)
+                .with_context(|| format!("Failed to open GPX: {}", path.display()))?;
+
+            read(BufReader::new(file)).context("Failed to parse GPX")
+        }
+
+        _ => {
+            let stdin = io::stdin();
+
+            read(stdin.lock()).context("Failed to parse GPX from stdin")
+        }
+    }
+}
+
+pub fn write_gpx(gpx: &Gpx, path: Option<&Path>) -> Result<()> {
+    match path {
+        Some(path) if path.as_os_str() != "-" => {
+            let file = File::create(path)
+                .with_context(|| format!("Failed to create output GPX: {}", path.display()))?;
+
+            write(gpx, BufWriter::new(file)).context("Failed to write GPX")?;
+        }
+
+        _ => {
+            let stdout = io::stdout();
+
+            write(gpx, BufWriter::new(stdout.lock())).context("Failed to write GPX to stdout")?;
+        }
+    }
+
+    Ok(())
+}
+
+pub fn extract_track_points(gpx: &Gpx) -> Result<Vec<TrackPoint>> {
+    let mut points = Vec::new();
+
+    for track in &gpx.tracks {
+        for segment in &track.segments {
+            for waypoint in &segment.points {
+                let point = waypoint.point();
+
+                points.push(TrackPoint {
+                    lon: point.x(),
+                    lat: point.y(),
+                });
+            }
+        }
+    }
+
+    if points.len() < 2 {
+        bail!("The GPX contains fewer than two track points.");
+    }
+
+    Ok(points)
+}
+
+pub fn append_crossing_waypoints(gpx: &mut Gpx, crossings: &[Crossing]) {
+    for (index, crossing) in crossings.iter().enumerate() {
+        let mut waypoint = Waypoint::new(crossing.point);
+
+        waypoint.name = Some(format!(
+            "Border {:02}: {} -> {}",
+            index + 1,
+            crossing.from.name,
+            crossing.to.name,
+        ));
+
+        waypoint.description = Some(format!(
+            "Administrative border crossing: {} [{}] -> {} [{}]",
+            crossing.from.name, crossing.from.id, crossing.to.name, crossing.to.id,
+        ));
+
+        waypoint.comment = Some(format!("{} -> {}", crossing.from.name, crossing.to.name,));
+
+        waypoint.type_ = Some("administrative_boundary".to_string());
+
+        gpx.waypoints.push(waypoint);
+    }
+}

src/main.rs

@@ -1,103 +1,24 @@
-use anyhow::{bail, Context, Result};
-use clap::Parser;
+mod cli;
+mod county;
+mod crossings;
+mod geometry;
+mod gpx_io;
+mod report;
 
-use geo::{algorithm::bounding_rect::BoundingRect, Contains, Coord, Geometry, LineString, Point};
+use anyhow::Result;
 
-use geojson::{Feature, GeoJson};
-
-use gpx::{read, write, Gpx, Waypoint};
-
-use rstar::{RTree, RTreeObject, AABB};
-
-use serde_json::Value;
-
-use std::{
-    cmp::Ordering,
-    collections::HashSet,
-    fs::File,
-    io::{self, BufReader, BufWriter},
-    path::{Path, PathBuf},
-};
-
-#[derive(Parser, Debug)]
-#[command(
-    name = "borderpoi-rs",
-    version,
-    about = "Find administrative border crossings along a GPX track"
-)]
-struct Args {
-    /// Input GPX track. Reads from stdin when omitted or set to '-'.
-    #[arg(short, long)]
-    track: Option<PathBuf>,
-
-    /// Boundary GeoJSON in WGS84 / EPSG:4326.
-    #[arg(short, long)]
-    boundaries: PathBuf,
-
-    /// GeoJSON property containing the county name.
-    #[arg(long, default_value = "GEN")]
-    name_field: String,
-
-    /// GeoJSON property containing the administrative identifier.
-    #[arg(long, default_value = "AGS")]
-    id_field: String,
-
-    /// Output GPX. Writes to stdout when omitted or set to '-'.
-    #[arg(short, long)]
-    output: Option<PathBuf>,
-}
-
-#[derive(Debug, Clone)]
-struct County {
-    name: String,
-    ags: String,
-    geometry: Geometry<f64>,
-    envelope: AABB<[f64; 2]>,
-}
-
-impl RTreeObject for County {
-    type Envelope = AABB<[f64; 2]>;
-
-    fn envelope(&self) -> Self::Envelope {
-        self.envelope
-    }
-}
-
-#[derive(Debug, Clone, Copy)]
-struct TrackPoint {
-    lon: f64,
-    lat: f64,
-}
-
-#[derive(Debug, Clone)]
-struct Crossing {
-    point: Point<f64>,
-    from: County,
-    to: County,
-    segment_index: usize,
-    position: f64,
-}
-
-#[derive(Debug, Clone)]
-struct BoundaryHit {
-    point: Point<f64>,
-    position: f64,
-}
-
-// -----------------------------------------------------------------------------
-// Main
-// -----------------------------------------------------------------------------
+use cli::Args;
+use county::load_counties;
+use crossings::find_crossings;
+use gpx_io::{append_crossing_waypoints, read_gpx, write_gpx};
+use report::print_report;
 
 fn main() -> Result<()> {
-    let args = Args::parse();
-
-    eprintln!("borderpoi-rs");
-    eprintln!("============");
-    eprintln!();
+    let args = Args::parse_args();
 
     let gpx = read_gpx(args.track.as_deref())?;
 
-    let track = extract_track_points(&gpx)?;
+    let track = gpx_io::extract_track_points(&gpx)?;
 
     eprintln!("Track points: {}", track.len());
 
@@ -105,9 +26,9 @@ fn main() -> Result<()> {
 
     eprintln!("Loaded counties: {}", counties.len());
 
-    let tree = RTree::bulk_load(counties);
+    let tree = rstar::RTree::bulk_load(counties);
 
-    eprintln!("Spatial index: R*-tree");
+    eprintln!("Spatial index: R-tree");
 
     let crossings = find_crossings(&track, &tree)?;
 
@@ -121,550 +42,3 @@ fn main() -> Result<()> {
 
     Ok(())
 }
-
-// -----------------------------------------------------------------------------
-// GPX I/O
-// -----------------------------------------------------------------------------
-
-fn read_gpx(path: Option<&Path>) -> Result<Gpx> {
-    match path {
-        Some(path) if path.as_os_str() != "-" => {
-            let file = File::open(path)
-                .with_context(|| format!("Failed to open GPX: {}", path.display()))?;
-
-            read(BufReader::new(file)).context("Failed to parse GPX")
-        }
-
-        _ => {
-            let stdin = io::stdin();
-            read(stdin.lock()).context("Failed to parse GPX from stdin")
-        }
-    }
-}
-
-fn write_gpx(gpx: &Gpx, path: Option<&Path>) -> Result<()> {
-    match path {
-        Some(path) if path.as_os_str() != "-" => {
-            let file = File::create(path)
-                .with_context(|| format!("Failed to create output GPX: {}", path.display()))?;
-
-            write(gpx, BufWriter::new(file)).context("Failed to write GPX")?;
-        }
-
-        _ => {
-            let stdout = io::stdout();
-
-            write(gpx, BufWriter::new(stdout.lock())).context("Failed to write GPX to stdout")?;
-        }
-    }
-
-    Ok(())
-}
-
-fn extract_track_points(gpx: &Gpx) -> Result<Vec<TrackPoint>> {
-    let mut points = Vec::new();
-
-    for track in &gpx.tracks {
-        for segment in &track.segments {
-            for waypoint in &segment.points {
-                let point = waypoint.point();
-
-                points.push(TrackPoint {
-                    lon: point.x(),
-                    lat: point.y(),
-                });
-            }
-        }
-    }
-
-    if points.len() < 2 {
-        bail!("The GPX contains fewer than two track points.");
-    }
-
-    Ok(points)
-}
-
-// -----------------------------------------------------------------------------
-// GeoJSON
-// -----------------------------------------------------------------------------
-
-fn load_counties(path: &Path, name_field: &str, id_field: &str) -> Result<Vec<County>> {
-    let file =
-        File::open(path).with_context(|| format!("Failed to open GeoJSON: {}", path.display()))?;
-
-    let geojson: GeoJson = serde_json::from_reader(BufReader::new(file))
-        .context("Failed to parse boundary GeoJSON")?;
-
-    let features = match geojson {
-        GeoJson::FeatureCollection(collection) => collection.features,
-
-        GeoJson::Feature(feature) => {
-            vec![feature]
-        }
-
-        GeoJson::Geometry(_) => {
-            bail!("Boundary GeoJSON must be a FeatureCollection or Feature.");
-        }
-    };
-
-    let mut counties = Vec::with_capacity(features.len());
-
-    for feature in features {
-        counties.push(feature_to_county(&feature, name_field, id_field)?);
-    }
-
-    if counties.is_empty() {
-        bail!("No boundary features were found.");
-    }
-
-    Ok(counties)
-}
-
-fn feature_to_county(feature: &Feature, name_field: &str, id_field: &str) -> Result<County> {
-    let properties = feature
-        .properties
-        .as_ref()
-        .context("Boundary feature has no properties.")?;
-
-    let name = property_as_string(
-        properties
-            .get(name_field)
-            .with_context(|| format!("Missing name field '{}'.", name_field))?,
-    )?;
-
-    let ags = property_as_string(
-        properties
-            .get(id_field)
-            .with_context(|| format!("Missing ID field '{}'.", id_field))?,
-    )?;
-
-    let geojson_geometry = feature
-        .geometry
-        .as_ref()
-        .context("Boundary feature has no geometry.")?;
-
-    let geometry: Geometry<f64> = geojson_geometry
-        .try_into()
-        .context("Failed to convert GeoJSON geometry.")?;
-
-    match geometry {
-        Geometry::Polygon(_) | Geometry::MultiPolygon(_) => {}
-
-        _ => {
-            bail!("Feature '{}' is not a Polygon or MultiPolygon.", name);
-        }
-    }
-
-    let bbox = geometry
-        .bounding_rect()
-        .with_context(|| format!("Feature '{}' has no bounding box.", name))?;
-
-    let envelope = AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]);
-
-    Ok(County {
-        name,
-        ags,
-        geometry,
-        envelope,
-    })
-}
-
-fn property_as_string(value: &Value) -> Result<String> {
-    match value {
-        Value::String(value) => Ok(value.clone()),
-
-        Value::Number(value) => Ok(value.to_string()),
-
-        _ => {
-            bail!("Expected string or number property, got {}.", value);
-        }
-    }
-}
-
-// -----------------------------------------------------------------------------
-// Track processing
-// -----------------------------------------------------------------------------
-
-/// Find all county transitions along the complete track.
-///
-/// Every track segment is processed independently. A segment may cross
-/// multiple administrative boundaries. All intersections are therefore
-/// collected, sorted along the segment, and converted into transitions.
-fn find_crossings(track: &[TrackPoint], tree: &RTree<County>) -> Result<Vec<Crossing>> {
-    let mut crossings = Vec::new();
-
-    for (segment_index, pair) in track.windows(2).enumerate() {
-        let start = Point::new(pair[0].lon, pair[0].lat);
-
-        let end = Point::new(pair[1].lon, pair[1].lat);
-
-        let segment = LineString::from(vec![
-            Coord {
-                x: start.x(),
-                y: start.y(),
-            },
-            Coord {
-                x: end.x(),
-                y: end.y(),
-            },
-        ]);
-
-        let bbox = segment
-            .bounding_rect()
-            .context("Track segment has no bounding box.")?;
-
-        let envelope =
-            AABB::from_corners([bbox.min().x, bbox.min().y], [bbox.max().x, bbox.max().y]);
-
-        let candidates: Vec<County> = tree
-            .locate_in_envelope_intersecting(envelope)
-            .cloned()
-            .collect();
-
-        if candidates.is_empty() {
-            continue;
-        }
-
-        let hits = collect_segment_hits(start, end, &candidates);
-
-        if hits.is_empty() {
-            continue;
-        }
-
-        let transitions = reconstruct_transitions(start, end, &candidates, &hits);
-
-        for transition in transitions {
-            crossings.push(Crossing {
-                point: transition.point,
-                from: transition.from,
-                to: transition.to,
-                segment_index,
-                position: transition.position,
-            });
-        }
-    }
-
-    Ok(deduplicate_crossings(crossings))
-}
-
-// -----------------------------------------------------------------------------
-// Boundary intersections
-// -----------------------------------------------------------------------------
-
-fn collect_segment_hits(
-    start: Point<f64>,
-    end: Point<f64>,
-    candidates: &[County],
-) -> Vec<BoundaryHit> {
-    let start_coord = Coord {
-        x: start.x(),
-        y: start.y(),
-    };
-
-    let end_coord = Coord {
-        x: end.x(),
-        y: end.y(),
-    };
-
-    let mut hits = Vec::new();
-
-    for county in candidates {
-        collect_geometry_intersections(&county.geometry, start_coord, end_coord, &mut hits);
-    }
-
-    hits.sort_by(|a, b| {
-        a.position
-            .partial_cmp(&b.position)
-            .unwrap_or(Ordering::Equal)
-    });
-
-    deduplicate_hits(hits)
-}
-
-fn collect_geometry_intersections(
-    geometry: &Geometry<f64>,
-    start: Coord<f64>,
-    end: Coord<f64>,
-    hits: &mut Vec<BoundaryHit>,
-) {
-    match geometry {
-        Geometry::Polygon(polygon) => {
-            collect_ring_intersections(polygon.exterior(), start, end, hits);
-
-            for interior in polygon.interiors() {
-                collect_ring_intersections(interior, start, end, hits);
-            }
-        }
-
-        Geometry::MultiPolygon(multipolygon) => {
-            for polygon in &multipolygon.0 {
-                collect_ring_intersections(polygon.exterior(), start, end, hits);
-
-                for interior in polygon.interiors() {
-                    collect_ring_intersections(interior, start, end, hits);
-                }
-            }
-        }
-
-        _ => {}
-    }
-}
-
-fn collect_ring_intersections(
-    ring: &LineString<f64>,
-    start: Coord<f64>,
-    end: Coord<f64>,
-    hits: &mut Vec<BoundaryHit>,
-) {
-    for edge in ring.lines() {
-        if let Some((position, point)) = segment_intersection(start, end, edge.start, edge.end) {
-            hits.push(BoundaryHit { point, position });
-        }
-    }
-}
-
-fn segment_intersection(
-    p: Coord<f64>,
-    p2: Coord<f64>,
-    q: Coord<f64>,
-    q2: Coord<f64>,
-) -> Option<(f64, Point<f64>)> {
-    let r = Coord {
-        x: p2.x - p.x,
-        y: p2.y - p.y,
-    };
-
-    let s = Coord {
-        x: q2.x - q.x,
-        y: q2.y - q.y,
-    };
-
-    let denominator = cross(r, s);
-
-    if denominator.abs() < 1e-14 {
-        return None;
-    }
-
-    let qp = Coord {
-        x: q.x - p.x,
-        y: q.y - p.y,
-    };
-
-    let t = cross(qp, s) / denominator;
-    let u = cross(qp, r) / denominator;
-
-    const EPSILON: f64 = 1e-10;
-
-    if !(-EPSILON..=1.0 + EPSILON).contains(&t) || !(-EPSILON..=1.0 + EPSILON).contains(&u) {
-        return None;
-    }
-
-    let point = Point::new(p.x + t * r.x, p.y + t * r.y);
-
-    Some((t.clamp(0.0, 1.0), point))
-}
-
-fn cross(a: Coord<f64>, b: Coord<f64>) -> f64 {
-    a.x * b.y - a.y * b.x
-}
-
-fn deduplicate_hits(hits: Vec<BoundaryHit>) -> Vec<BoundaryHit> {
-    let mut result = Vec::new();
-
-    for hit in hits {
-        let duplicate = result
-            .iter()
-            .any(|existing: &BoundaryHit| (existing.position - hit.position).abs() < 1e-9);
-
-        if !duplicate {
-            result.push(hit);
-        }
-    }
-
-    result
-}
-
-// -----------------------------------------------------------------------------
-// County reconstruction
-// -----------------------------------------------------------------------------
-
-#[derive(Debug, Clone)]
-struct Transition {
-    point: Point<f64>,
-    position: f64,
-    from: County,
-    to: County,
-}
-
-/// Reconstruct the county sequence between all boundary intersections.
-///
-/// The county immediately before an intersection is sampled slightly before
-/// the intersection; the county immediately after is sampled slightly after
-/// it. This avoids ambiguity when the intersection lies exactly on a shared
-/// polygon boundary.
-fn reconstruct_transitions(
-    start: Point<f64>,
-    end: Point<f64>,
-    candidates: &[County],
-    hits: &[BoundaryHit],
-) -> Vec<Transition> {
-    let mut transitions = Vec::new();
-
-    for hit in hits {
-        let before_t = (hit.position - 1e-8).max(0.0);
-
-        let after_t = (hit.position + 1e-8).min(1.0);
-
-        let before = interpolate(start, end, before_t);
-
-        let after = interpolate(start, end, after_t);
-
-        let from = county_at_point(before, candidates);
-
-        let to = county_at_point(after, candidates);
-
-        let (Some(from), Some(to)) = (from, to) else {
-            continue;
-        };
-
-        if from.ags == to.ags {
-            continue;
-        }
-
-        transitions.push(Transition {
-            point: hit.point,
-            position: hit.position,
-            from,
-            to,
-        });
-    }
-
-    transitions
-}
-
-fn interpolate(start: Point<f64>, end: Point<f64>, t: f64) -> Point<f64> {
-    Point::new(
-        start.x() + t * (end.x() - start.x()),
-        start.y() + t * (end.y() - start.y()),
-    )
-}
-
-fn county_at_point(point: Point<f64>, candidates: &[County]) -> Option<County> {
-    candidates
-        .iter()
-        .find(|county| county.geometry.contains(&point))
-        .cloned()
-}
-
-// -----------------------------------------------------------------------------
-// Crossing cleanup
-// -----------------------------------------------------------------------------
-
-fn deduplicate_crossings(crossings: Vec<Crossing>) -> Vec<Crossing> {
-    let mut result = Vec::new();
-
-    for crossing in crossings {
-        let duplicate = result.iter().any(|existing: &Crossing| {
-            existing.segment_index == crossing.segment_index
-                && existing.from.ags == crossing.from.ags
-                && existing.to.ags == crossing.to.ags
-                && same_point(existing.point, crossing.point)
-        });
-
-        if !duplicate {
-            result.push(crossing);
-        }
-    }
-
-    result.sort_by(|a, b| {
-        a.segment_index.cmp(&b.segment_index).then_with(|| {
-            a.position
-                .partial_cmp(&b.position)
-                .unwrap_or(Ordering::Equal)
-        })
-    });
-
-    result
-}
-
-fn same_point(a: Point<f64>, b: Point<f64>) -> bool {
-    let dx = a.x() - b.x();
-    let dy = a.y() - b.y();
-
-    dx * dx + dy * dy < 1e-10
-}
-
-// -----------------------------------------------------------------------------
-// Report
-// -----------------------------------------------------------------------------
-
-fn print_report(crossings: &[Crossing]) {
-    eprintln!();
-    eprintln!("County visits");
-    eprintln!("=============");
-
-    if crossings.is_empty() {
-        eprintln!("No county crossings found.");
-        return;
-    }
-
-    let mut visits = Vec::new();
-
-    visits.push(crossings[0].from.clone());
-
-    for crossing in crossings {
-        visits.push(crossing.to.clone());
-    }
-
-    for (index, county) in visits.iter().enumerate() {
-        eprintln!("{:3}  {} [{}]", index + 1, county.name, county.ags);
-    }
-
-    let unique: HashSet<String> = visits.iter().map(|county| county.ags.clone()).collect();
-
-    eprintln!();
-    eprintln!("Visits:          {}", visits.len());
-
-    eprintln!("Unique counties: {}", unique.len());
-
-    eprintln!();
-    eprintln!("Border crossings");
-    eprintln!("=================");
-
-    for (index, crossing) in crossings.iter().enumerate() {
-        eprintln!(
-            "{:3}  {} -> {}",
-            index + 1,
-            crossing.from.name,
-            crossing.to.name
-        );
-    }
-}
-
-// -----------------------------------------------------------------------------
-// GPX POIs
-// -----------------------------------------------------------------------------
-
-fn append_crossing_waypoints(gpx: &mut Gpx, crossings: &[Crossing]) {
-    for (index, crossing) in crossings.iter().enumerate() {
-        let mut waypoint = Waypoint::new(crossing.point);
-
-        waypoint.name = Some(format!(
-            "Border {:02}: {} -> {}",
-            index + 1,
-            crossing.from.name,
-            crossing.to.name
-        ));
-
-        waypoint.description = Some(format!(
-            "Administrative border crossing: {} [{}] -> {} [{}]",
-            crossing.from.name, crossing.from.ags, crossing.to.name, crossing.to.ags
-        ));
-
-        waypoint.comment = Some(format!("{} -> {}", crossing.from.name, crossing.to.name));
-
-        waypoint.type_ = Some("administrative_boundary".to_string());
-
-        gpx.waypoints.push(waypoint);
-    }
-}

src/report.rs

@@ -0,0 +1,46 @@
+use std::collections::HashSet;
+
+use crate::crossings::Crossing;
+
+pub fn print_report(crossings: &[Crossing]) {
+    eprintln!();
+    eprintln!("County visits");
+    eprintln!("=============");
+
+    if crossings.is_empty() {
+        eprintln!("No county crossings found.");
+        return;
+    }
+
+    let mut visits = Vec::new();
+
+    visits.push(crossings[0].from.clone());
+
+    for crossing in crossings {
+        visits.push(crossing.to.clone());
+    }
+
+    for (index, county) in visits.iter().enumerate() {
+        eprintln!("{:3}  {} [{}]", index + 1, county.name, county.id,);
+    }
+
+    let unique: HashSet<String> = visits.iter().map(|county| county.id.clone()).collect();
+
+    eprintln!();
+    eprintln!("Visits:          {}", visits.len());
+
+    eprintln!("Unique counties: {}", unique.len());
+
+    eprintln!();
+    eprintln!("Border crossings");
+    eprintln!("=================");
+
+    for (index, crossing) in crossings.iter().enumerate() {
+        eprintln!(
+            "{:3}  {} -> {}",
+            index + 1,
+            crossing.from.name,
+            crossing.to.name,
+        );
+    }
+}