diff --git a/geo/CHANGES.md b/geo/CHANGES.md index 187385c92..99469f57d 100644 --- a/geo/CHANGES.md +++ b/geo/CHANGES.md @@ -3,6 +3,8 @@ ## Unreleased * Fix a panic when calculating the haversine closest point to a point intersecting the geometry * +* Add `LineStringSegmentizeHaversine` trait as a an alternative to `LineStringSegmentize` for geographic coordinates. + * ## 0.27.0 diff --git a/geo/src/algorithm/linestring_segment.rs b/geo/src/algorithm/linestring_segment.rs index 426732c16..0ddc61a18 100644 --- a/geo/src/algorithm/linestring_segment.rs +++ b/geo/src/algorithm/linestring_segment.rs @@ -1,5 +1,8 @@ use crate::line_interpolate_point::LineInterpolatePoint; -use crate::{Coord, Densify, EuclideanLength, LineString, LinesIter, MultiLineString}; +use crate::{ + Coord, Densify, DensifyHaversine, EuclideanLength, HaversineLength, LineString, LinesIter, + MultiLineString, +}; /// Segments a LineString into `n` equal length LineStrings as a MultiLineString. /// `None` will be returned when `n` is equal to 0 or when a point @@ -8,140 +11,120 @@ use crate::{Coord, Densify, EuclideanLength, LineString, LinesIter, MultiLineStr /// /// # Examples /// ``` -/// use geo::{LineString, MultiLineString, LineStringSegmentize, Coord}; +/// use geo::{LineString, MultiLineString, LineStringSegmentize}; /// // Create a simple line string /// let lns: LineString = vec![[0.0, 0.0], [1.0, 2.0], [3.0, 6.0]].into(); -/// // Segment it into 6 LineStrings inside of a MultiLineString -/// let segmentized = lns.line_segmentize(6).unwrap(); -/// -/// // Recreate the MultiLineString from scratch -/// // this is the inner vector used to create the MultiLineString -/// let all_lines = vec![ -/// LineString::new(vec![Coord { x: 0.0, y: 0.0 }, Coord { x: 0.5, y: 1.0 }]), -/// LineString::new(vec![Coord { x: 0.5, y: 1.0 }, Coord { x: 1.0, y: 2.0 }]), -/// LineString::new(vec![Coord { x: 1.0, y: 2.0 }, Coord { x: 1.5, y: 3.0 }]), -/// LineString::new(vec![Coord { x: 1.5, y: 3.0 }, Coord { x: 2.0, y: 4.0 }]), -/// LineString::new(vec![Coord { x: 2.0, y: 4.0 }, Coord { x: 2.5, y: 5.0 }]), -/// LineString::new(vec![Coord { x: 2.5, y: 5.0 }, Coord { x: 3.0, y: 6.0 }]) -/// ]; -/// -/// // Create the MultiLineString -/// let mlns = MultiLineString::new(all_lines); -/// -/// // Compare the two -/// assert_eq!(mlns, segmentized); +/// // Segment it into n LineStrings inside of a MultiLineString +/// let n = 6; +/// let segmentized = lns.line_segmentize(n).unwrap(); +/// // Compare the number of elements +/// assert_eq!(n, segmentized.0.len()); ///``` pub trait LineStringSegmentize { fn line_segmentize(&self, n: usize) -> Option; } -impl LineStringSegmentize for LineString { - fn line_segmentize(&self, n: usize) -> Option { - // Return None if n is 0 or the maximum usize - if (n == usize::MIN) || (n == usize::MAX) { - return None; - } else if n == 1 { - let mlns = MultiLineString::from(self.clone()); - return Some(mlns); - } - - // Vec to allocate the new LineString segments Coord Vec - // will be iterated over at end to create new vecs - let mut res_coords: Vec> = Vec::with_capacity(n); - - // calculate total length to track cumulative against - let total_length = self.euclidean_length().abs(); - - // tracks total length - let mut cum_length = 0_f64; - - // calculate the target fraction for the first iteration - // fraction will change based on each iteration - let segment_prop = (1_f64) / (n as f64); - let segment_length = total_length * segment_prop; - - // densify the LineString so that each `Line` segment is not longer - // than the segment length ensuring that we will never partition one - // Line more than once. - let densified = self.densify(segment_length); - - // if the densified line is exactly equal to the number of requested - // segments, return early. This will happen when a LineString has - // exactly 2 coordinates - if densified.lines().count() == n { - let linestrings = densified - .lines() - .map(LineString::from) - .collect::>(); - - return Some(MultiLineString::new(linestrings)); - }; - - // count the number of lines that will be iterated through - let n_lines = densified.lines().count(); - - let lns = densified.lines_iter(); - // instantiate the first Vec - let mut ln_vec: Vec = Vec::new(); - - // iterate through each line segment in the LineString - for (i, segment) in lns.enumerate() { - // All iterations only keep track of the second coordinate - // in the Line. We need to push the first coordinate in the - // first line string to ensure the linestring starts at the - // correct place - if i == 0 { - ln_vec.push(segment.start) - } - - let length = segment.euclidean_length().abs(); - - // update cumulative length - cum_length += length; - - if (cum_length >= segment_length) && (i != (n_lines - 1)) { - let remainder = cum_length - segment_length; - // if we get None, we exit the function and return None - let endpoint = segment.line_interpolate_point((length - remainder) / length)?; - - // add final coord to ln_vec - ln_vec.push(endpoint.into()); +/// Segments a LineString into `n` equal length LineStrings as a MultiLineString +/// using Haversine distance calculations. Use this over `LineStringSegmentize` +/// when using data from a geographic coordinate system. +/// `None` will be returned when `n` is equal to 0 or when a point +/// cannot be interpolated on a `Line` segment. +/// +/// +/// # Examples +/// ``` +/// use geo::{LineString, MultiLineString, LineStringSegmentizeHaversine}; +/// // Create a simple line string +/// let lns: LineString = vec![[0.0, 0.0], [1.0, 2.0], [3.0, 6.0]].into(); +/// // Segment it into n LineStrings inside of a MultiLineString +/// let n = 6; +/// let segmentized = lns.line_segmentize_haversine(n).unwrap(); +/// // Compare the number of elements +/// assert_eq!(n, segmentized.0.len()); +///``` +pub trait LineStringSegmentizeHaversine { + fn line_segmentize_haversine(&self, n: usize) -> Option; +} - // now we drain all elements from the vector into an iterator - // this will be collected into a vector to be pushed into the - // results coord vec of vec - let to_push = ln_vec.drain(..); +macro_rules! implement_segmentize { + ($trait_name:ident, $method_name:ident, $distance_method:ident, $densify_method:ident) => { + impl $trait_name for LineString { + fn $method_name(&self, n: usize) -> Option { + if (n == usize::MIN) || (n == usize::MAX) { + return None; + } else if n == 1 { + let mlns = MultiLineString::from(self.clone()); + return Some(mlns); + } - // now collect & push this vector into the results vector - res_coords.push(to_push.collect::>()); + let mut res_coords: Vec> = Vec::with_capacity(n); + let total_length = self.$distance_method().abs(); + let mut cum_length = 0_f64; + let segment_prop = (1_f64) / (n as f64); + let segment_length = total_length * segment_prop; + let densified = self.$densify_method(segment_length - f64::EPSILON); + + if densified.lines().count() == n { + let linestrings = densified + .lines() + .map(LineString::from) + .collect::>(); + return Some(MultiLineString::new(linestrings)); + } - // now add the last endpoint as the first coord - // and the endpoint of the linesegment as well only - if i != n_lines { - ln_vec.push(endpoint.into()); + let n_lines = densified.lines().count(); + let lns = densified.lines_iter(); + let mut ln_vec: Vec = Vec::new(); + + for (i, segment) in lns.enumerate() { + if i == 0 { + ln_vec.push(segment.start) + } + + let length = segment.$distance_method().abs(); + cum_length += length; + + if (cum_length >= segment_length) && (i != (n_lines - 1)) { + let remainder = cum_length - segment_length; + let endpoint = + segment.line_interpolate_point((length - remainder) / length)?; + + ln_vec.push(endpoint.into()); + let to_push = ln_vec.drain(..); + res_coords.push(to_push.collect::>()); + + if i != n_lines { + ln_vec.push(endpoint.into()); + } + cum_length = remainder; + } + ln_vec.push(segment.end); } - cum_length = remainder; + res_coords.push(ln_vec); + let res_lines = res_coords + .into_iter() + .map(LineString::new) + .collect::>(); + Some(MultiLineString::new(res_lines)) } - - // push the end coordinate into the Vec to continue - // building the linestring - ln_vec.push(segment.end); } + }; +} - // push the last linestring vector which isn't done by the for loop - res_coords.push(ln_vec); - - // collect the coords into vectors of LineStrings so we can createa - // a multi linestring - let res_lines = res_coords - .into_iter() - .map(LineString::new) - .collect::>(); +implement_segmentize!( + LineStringSegmentize, + line_segmentize, + euclidean_length, + densify +); - Some(MultiLineString::new(res_lines)) - } -} +implement_segmentize!( + LineStringSegmentizeHaversine, + line_segmentize_haversine, + haversine_length, + densify_haversine +); #[cfg(test)] mod test { @@ -305,4 +288,78 @@ mod test { let segments = linestring.line_segmentize(2).unwrap(); assert_eq!(segments.0.len(), 2) } + + #[test] + fn tiny_distances() { + // this test is to ensure that at super small distances + // the number of units is still the specified one. + let linestring: LineString = vec![ + [-3.19416, 55.95524], + [-3.19352, 55.95535], + [-3.19288, 55.95546], + ] + .into(); + + let n = 8; + let segments = linestring.line_segmentize(n).unwrap(); + assert_eq!(segments.0.len(), n) + } + + #[test] + fn haversine_n_elems() { + let linestring: LineString = vec![ + [-3.19416, 55.95524], + [-3.19352, 55.95535], + [-3.19288, 55.95546], + ] + .into(); + + let n = 8; + + let segments = linestring.line_segmentize_haversine(n).unwrap(); + assert_eq!(n, segments.0.len()); + } + + #[test] + fn haversine_segment_length() { + let linestring: LineString = vec![ + [-3.19416, 55.95524], + [-3.19352, 55.95535], + [-3.19288, 55.95546], + ] + .into(); + + let n = 8; + + let segments = linestring.line_segmentize_haversine(n).unwrap(); + let lens = segments + .0 + .iter() + .map(|li| li.haversine_length()) + .collect::>(); + + let epsilon = 1e-6; // 6th decimal place which is micrometers + assert!(lens.iter().all(|&x| (x - lens[0]).abs() < epsilon)); + } + + #[test] + fn haversine_total_length() { + let linestring: LineString = vec![ + [-3.19416, 55.95524], + [-3.19352, 55.95535], + [-3.19288, 55.95546], + ] + .into(); + + let n = 8; + + let segments = linestring.line_segmentize_haversine(n).unwrap(); + + // different at 12th decimal which is a picometer + assert_relative_eq!( + linestring.haversine_length(), + segments.haversine_length(), + epsilon = 1e-11 + ); + } } diff --git a/geo/src/algorithm/mod.rs b/geo/src/algorithm/mod.rs index 791ad98fe..4eb8b7956 100644 --- a/geo/src/algorithm/mod.rs +++ b/geo/src/algorithm/mod.rs @@ -190,7 +190,7 @@ pub use lines_iter::LinesIter; /// Split a LineString into n segments pub mod linestring_segment; -pub use linestring_segment::LineStringSegmentize; +pub use linestring_segment::{LineStringSegmentize, LineStringSegmentizeHaversine}; /// Apply a function to all `Coord`s of a `Geometry`. pub mod map_coords; diff --git a/geo/src/lib.rs b/geo/src/lib.rs index 6606c4134..f6e9d95ff 100644 --- a/geo/src/lib.rs +++ b/geo/src/lib.rs @@ -166,6 +166,7 @@ //! - **[`RhumbIntermediate`]**: Calculate intermediate points on a sphere along a rhumb line //! - **[`proj`]**: Project geometries with the `proj` crate (requires the `use-proj` feature) //! - **[`LineStringSegmentize`]**: Segment a LineString into `n` segments. +//! - **[`LineStringSegmentizeHaversine`]**: Segment a LineString using Haversine distance. //! - **[`Transform`]**: Transform a geometry using Proj. //! - **[`RemoveRepeatedPoints`]**: Remove repeated points from a geometry. //!