From 999aa5d32941ccefbf355ae36b4bac670df6dbb4 Mon Sep 17 00:00:00 2001 From: Jonas Rabenstein Date: Wed, 12 Aug 2026 01:41:30 +0200 Subject: [PATCH] split and snap --- src/cli.rs | 36 +++ src/county.rs | 110 ++++++++ src/crossings.rs | 249 ++++++++++++++++++ src/geometry.rs | 87 +++++++ src/gpx_io.rs | 92 +++++++ src/main.rs | 658 ++--------------------------------------------- src/report.rs | 46 ++++ 7 files changed, 636 insertions(+), 642 deletions(-) create mode 100644 src/cli.rs create mode 100644 src/county.rs create mode 100644 src/crossings.rs create mode 100644 src/geometry.rs create mode 100644 src/gpx_io.rs create mode 100644 src/report.rs diff --git a/src/cli.rs b/src/cli.rs new file mode 100644 index 0000000..dfa6bd6 --- /dev/null +++ b/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, + + /// 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, +} + +impl Args { + pub fn parse_args() -> Self { + Self::parse() + } +} diff --git a/src/county.rs b/src/county.rs new file mode 100644 index 0000000..b9b2b56 --- /dev/null +++ b/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, + 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> { + 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 { + 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 = 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 { + match value { + Value::String(value) => Ok(value.clone()), + Value::Number(value) => Ok(value.to_string()), + _ => { + bail!("Expected string or number property, got {}.", value); + } + } +} diff --git a/src/crossings.rs b/src/crossings.rs new file mode 100644 index 0000000..3745407 --- /dev/null +++ b/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, + 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) -> Result> { + 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 = 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, + end: Point, + candidates: &[County], +) -> Vec { + 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, + start: Coord, + end: Coord, + hits: &mut Vec, +) { + 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, + start: Coord, + end: Coord, + hits: &mut Vec, +) { + 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, + end: Point, + candidates: &[County], + hits: &[BoundaryHit], +) -> Vec { + 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, candidates: &[County]) -> Option { + candidates + .iter() + .find(|county| county.geometry.contains(&point)) + .cloned() +} + +fn deduplicate_crossings(mut crossings: Vec) -> Vec { + 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, b: Point) -> bool { + let dx = a.x() - b.x(); + let dy = a.y() - b.y(); + + dx * dx + dy * dy < 1e-20 +} diff --git a/src/geometry.rs b/src/geometry.rs new file mode 100644 index 0000000..4c21143 --- /dev/null +++ b/src/geometry.rs @@ -0,0 +1,87 @@ +use geo::{Coord, Point}; + +use crate::crossings::BoundaryHit; + +pub fn cross(a: Coord, b: Coord) -> f64 { + a.x * b.y - a.y * b.x +} + +pub fn interpolate(start: Point, end: Point, t: f64) -> Point { + Point::new( + start.x() + t * (end.x() - start.x()), + start.y() + t * (end.y() - start.y()), + ) +} + +pub fn segment_intersection( + p: Coord, + p2: Coord, + q: Coord, + q2: Coord, +) -> Option<(f64, Point)> { + 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) -> Vec { + 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, end: Point, position: f64) -> Point { + interpolate(start, end, position.clamp(0.0, 1.0)) +} diff --git a/src/gpx_io.rs b/src/gpx_io.rs new file mode 100644 index 0000000..a6447c1 --- /dev/null +++ b/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 { + 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> { + 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); + } +} diff --git a/src/main.rs b/src/main.rs index cf3f79d..ede4b0d 100644 --- a/src/main.rs +++ b/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, - - /// 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, -} - -#[derive(Debug, Clone)] -struct County { - name: String, - ags: String, - geometry: Geometry, - 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, - from: County, - to: County, - segment_index: usize, - position: f64, -} - -#[derive(Debug, Clone)] -struct BoundaryHit { - point: Point, - 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 { - 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> { - 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> { - 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 { - 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 = 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 { - 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) -> Result> { - 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 = 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, - end: Point, - candidates: &[County], -) -> Vec { - 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, - start: Coord, - end: Coord, - hits: &mut Vec, -) { - 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, - start: Coord, - end: Coord, - hits: &mut Vec, -) { - 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, - p2: Coord, - q: Coord, - q2: Coord, -) -> Option<(f64, Point)> { - 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, b: Coord) -> f64 { - a.x * b.y - a.y * b.x -} - -fn deduplicate_hits(hits: Vec) -> Vec { - 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, - 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, - end: Point, - candidates: &[County], - hits: &[BoundaryHit], -) -> Vec { - 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, end: Point, t: f64) -> Point { - Point::new( - start.x() + t * (end.x() - start.x()), - start.y() + t * (end.y() - start.y()), - ) -} - -fn county_at_point(point: Point, candidates: &[County]) -> Option { - candidates - .iter() - .find(|county| county.geometry.contains(&point)) - .cloned() -} - -// ----------------------------------------------------------------------------- -// Crossing cleanup -// ----------------------------------------------------------------------------- - -fn deduplicate_crossings(crossings: Vec) -> Vec { - 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, b: Point) -> 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 = 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); - } -} diff --git a/src/report.rs b/src/report.rs new file mode 100644 index 0000000..3281ddf --- /dev/null +++ b/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 = 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, + ); + } +}