Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add DensifyGeodesic trait #1208

Closed
wants to merge 12 commits into from
2 changes: 2 additions & 0 deletions geo/CHANGES.md
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
* <https://github.com/georust/geo/pull/1201>
* Add `StitchTriangles` trait which implements a new kind of combining algorithm for `Triangle`s
* <https://github.com/georust/geo/pull/1087>
* Add `DensifyGeodesic` trait which implements densification by interpolating points on a spheroid
* <https://github.com/georust/geo/pull/1208>

## 0.28.0

Expand Down
268 changes: 268 additions & 0 deletions geo/src/algorithm/densify_geodesic.rs
Original file line number Diff line number Diff line change
@@ -0,0 +1,268 @@
use num_traits::ToPrimitive;

use crate::{
CoordFloat, CoordsIter, Line, LineString, MultiLineString, MultiPolygon, Point, Polygon, Rect,
Triangle,
};

use crate::{GeodesicIntermediate, GeodesicLength};

/// Returns a new geometry on a spheroid containing both existing and new interpolated coordinates with
/// a maximum distance of `max_meters` between them.
///
/// Note: `max_meters` must be greater than 0.
///
/// ## Units
///
/// `max_meters`: meters
///
/// # Examples
/// ```
/// use approx::assert_relative_eq;
///
/// use geo::{coord, GeodesicLength, Line, LineString};
/// use geo::DensifyGeodesic;
///
/// let line = Line::new(coord! {x: 10.0, y: 20.0}, coord! { x: 125.0, y: 25.00 });
/// // known output
/// let output: LineString = vec![[10.0, 20.0], [65.879360, 37.722253], [125.0, 25.00]].into();
/// // densify
/// let dense = line.densify_geodesic(5703861.471800622);
///
/// assert_relative_eq!(dense, output, epsilon = 1.0e-6);
///```
pub trait DensifyGeodesic<F: CoordFloat> {
type Output;

fn densify_geodesic(&self, max_meters: F) -> Self::Output;
}

// Helper for densification trait
fn densify_line(line: Line<f64>, container: &mut Vec<Point<f64>>, max_meters: f64) {
assert!(max_meters > 0.0);

container.push(line.start_point());

let num_segments = (line.geodesic_length() / max_meters)
.ceil()
.to_u64()
.unwrap();
// distance "unit" for this line segment
let frac = 1.0 / num_segments as f64;

let start = line.start;
let end = line.end;

for segment_idx in 1..num_segments {
let ratio = frac * segment_idx as f64;
let interpolated_point = Point::from(start).geodesic_intermediate(&Point::from(end), ratio);
container.push(interpolated_point);
}
}

impl DensifyGeodesic<f64> for MultiPolygon<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = MultiPolygon<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
MultiPolygon::new(
self.iter()
.map(|polygon| polygon.densify_geodesic(max_meters))
.collect(),
)
}
}

impl DensifyGeodesic<f64> for Polygon<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = Polygon<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
let densified_exterior = self.exterior().densify_geodesic(max_meters);
let densified_interiors = self
.interiors()
.iter()
.map(|ring| ring.densify_geodesic(max_meters))
.collect();
Polygon::new(densified_exterior, densified_interiors)
}
}

impl DensifyGeodesic<f64> for MultiLineString<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = MultiLineString<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
MultiLineString::new(
self.iter()
.map(|linestring| linestring.densify_geodesic(max_meters))
.collect(),
)
}
}

impl DensifyGeodesic<f64> for LineString<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = LineString<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
if self.coords_count() == 0 {
return LineString::new(vec![]);
}

let mut new_line = vec![];
self.lines()
.for_each(|line| densify_line(line, &mut new_line, max_meters));
// we're done, push the last coordinate on to finish
new_line.push(self.points().last().unwrap());
LineString::from(new_line)
}
}

impl DensifyGeodesic<f64> for Line<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = LineString<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
let mut new_line = vec![];
densify_line(*self, &mut new_line, max_meters);
// we're done, push the last coordinate on to finish
new_line.push(self.end_point());
LineString::from(new_line)
}
}

impl DensifyGeodesic<f64> for Triangle<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = Polygon<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
self.to_polygon().densify_geodesic(max_meters)
}
}

impl DensifyGeodesic<f64> for Rect<f64>
where
Line<f64>: GeodesicLength<f64>,
LineString<f64>: GeodesicLength<f64>,
{
type Output = Polygon<f64>;

fn densify_geodesic(&self, max_meters: f64) -> Self::Output {
self.to_polygon().densify_geodesic(max_meters)
}
}

#[cfg(test)]
mod tests {
use super::*;
use crate::coord;

grevgeny marked this conversation as resolved.
Show resolved Hide resolved
#[test]
fn test_polygon_densify() {
let exterior: LineString = vec![
[4.925, 45.804],
[4.732, 45.941],
[4.935, 46.513],
[5.821, 46.103],
[5.627, 45.611],
[5.355, 45.883],
[4.925, 45.804],
]
.into();

let polygon = Polygon::new(exterior, vec![]);

let output_exterior: LineString = vec![
[4.925, 45.804],
[4.732, 45.941],
[4.832972865149862, 46.22705224065524],
[4.935, 46.513],
[5.379653814979939, 46.30886184400083],
[5.821, 46.103],
[5.723572275808633, 45.85704648840237],
[5.627, 45.611],
[5.355, 45.883],
[4.925, 45.804],
]
.into();

let dense = polygon.densify_geodesic(50000.0);
assert_relative_eq!(dense.exterior(), &output_exterior, epsilon = 1.0e-6);
}

#[test]
fn test_linestring_densify() {
let linestring: LineString = vec![
[-3.202, 55.9471],
[-3.2012, 55.9476],
[-3.1994, 55.9476],
[-3.1977, 55.9481],
[-3.196, 55.9483],
[-3.1947, 55.9487],
[-3.1944, 55.9488],
[-3.1944, 55.949],
]
.into();

let output: LineString = vec![
[-3.202, 55.9471],
[-3.2012, 55.9476],
[-3.2002999999999995, 55.94760000327935],
[-3.1994, 55.9476],
[-3.1985500054877773, 55.94785000292509],
[-3.1977, 55.9481],
[-3.196, 55.9483],
[-3.1947, 55.9487],
[-3.1944, 55.9488],
[-3.1944, 55.949],
]
.into();

let dense = linestring.densify_geodesic(110.0);
assert_relative_eq!(dense, output, epsilon = 1.0e-6);
}

#[test]
fn test_line_densify() {
let output: LineString = vec![[10.0, 20.0], [65.879360, 37.722253], [125.0, 25.00]].into();
let line = Line::new(coord! {x: 10.0, y: 20.0}, coord! { x: 125.0, y: 25.00 });
let dense = line.densify_geodesic(5703861.471800622);
assert_relative_eq!(dense, output, epsilon = 1.0e-6);
}

#[test]
fn test_empty_linestring() {
let linestring: LineString<f64> = LineString::new(vec![]);
let dense = linestring.densify_geodesic(10.0);
assert_eq!(0, dense.coords_count());
}

#[test]
fn test_no_op_densify() {
let linestring: LineString = vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0]].into();
let output: LineString = vec![[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [3.0, 0.0]].into();
// Use a very large max_meters to ensure no points are added.
let dense = linestring.densify_geodesic(1000000.0);
// Check that the densified linestring is identical to the original.
assert_relative_eq!(dense, output, epsilon = 1.0e-6);
}
}
4 changes: 4 additions & 0 deletions geo/src/algorithm/mod.rs
Original file line number Diff line number Diff line change
Expand Up @@ -84,6 +84,10 @@ pub use densify::Densify;
pub mod densify_haversine;
pub use densify_haversine::DensifyHaversine;

/// Densify geometry components on a geodesic.
pub mod densify_geodesic;
pub use densify_geodesic::DensifyGeodesic;

/// Dimensionality of a geometry and its boundary, based on OGC-SFA.
pub mod dimensions;
pub use dimensions::HasDimensions;
Expand Down
1 change: 1 addition & 0 deletions geo/src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -158,6 +158,7 @@
//! - **[`ChaikinSmoothing`]**: Smoothen `LineString`, `Polygon`, `MultiLineString` and `MultiPolygon` using Chaikin's algorithm.
//! - **[`Densify`]**: Densify linear geometry components by interpolating points
//! - **[`DensifyHaversine`]**: Densify spherical geometry by interpolating points on a sphere
//! - **[`DensifyGeodesic`]**: Densify geometry by interpolating points on a spheroid
//! - **[`GeodesicDestination`]**: Given a start point, bearing, and distance, calculate the destination point on a [geodesic](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid)
//! - **[`GeodesicIntermediate`]**: Calculate intermediate points on a [geodesic](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid)
//! - **[`HaversineDestination`]**: Given a start point, bearing, and distance, calculate the destination point on a sphere assuming travel on a great circle
Expand Down
Loading