geo 0.29.3

Geospatial primitives and algorithms
Documentation
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
#![doc(html_logo_url = "https://raw.githubusercontent.com/georust/meta/master/logo/logo.png")]

//! The `geo` crate provides geospatial primitive types and algorithms.
//!
//! # Types
//!
//! - **[`Coord`]**: A two-dimensional coordinate. All geometry types are composed of [`Coord`]s, though [`Coord`] itself is not a [`Geometry`] type.
//! - **[`Point`]**: A single point represented by one [`Coord`]
//! - **[`MultiPoint`]**: A collection of [`Point`]s
//! - **[`Line`]**: A line segment represented by two [`Coord`]s
//! - **[`LineString`]**: A series of contiguous line segments represented by two or more
//!   [`Coord`]s
//! - **[`MultiLineString`]**: A collection of [`LineString`]s
//! - **[`Polygon`]**: A bounded area represented by one [`LineString`] exterior ring, and zero or
//!   more [`LineString`] interior rings
//! - **[`MultiPolygon`]**: A collection of [`Polygon`]s
//! - **[`Rect`]**: An axis-aligned bounded rectangle represented by minimum and maximum
//!   [`Coord`]s
//! - **[`Triangle`]**: A bounded area represented by three [`Coord`] vertices
//! - **[`GeometryCollection`]**: A collection of [`Geometry`]s
//! - **[`Geometry`]**: An enumeration of all geometry types, excluding [`Coord`]
//!
//! The preceding types are reexported from the [`geo-types`] crate. Consider using that crate
//! if you only need access to these types and no other `geo` functionality.
//!
//! ## Semantics
//!
//! The geospatial types provided here aim to adhere to the [OpenGIS Simple feature access][OGC-SFA]
//! standards. Thus, the types here are inter-operable with other implementations of the standards:
//! [JTS], [GEOS], etc.
//!
//! # Algorithms
//!
//! ## Measures
//!
//! Algorithms for measures along a line, and how a line is measured.
//!
//! ### Metric Spaces
//!
//! - **[`Euclidean`]**: The [Euclidean plane] measures distance with the pythagorean formula. Not suitable for lon/lat geometries.
//! - **[`Haversine`]**: The [Haversine Formula] measures distance on a sphere. Only suitable for lon/lat geometries.
//! - **[`Geodesic`]**: Geodesic methods based on [Karney (2013)] more accurately reflect the shape of the Earth, but are slower than Haversine. Only suitable for lon/lat geometries.
//! - **[`Rhumb`]**: [Rhumb line] (a.k.a. loxodrome) measures can be useful for navigation applications where maintaining a constant bearing or direction is important. Only suitable for lon/lat geometries.
//!
//! ### Operations on Metric Spaces
//!
//! - **[`Distance`]**: Calculate the minimum distance between two geometries.
//! - **[`Length`]**: Calculate the length of a `Line`, `LineString`, or `MultiLineString`.
//! - **[`Bearing`]**: Calculate the bearing between two points.
//!
//! - **[`Destination`]**: Calculate the destination point from an origin point, given a bearing and a distance.
//! - **[`InterpolatePoint`]**: Interpolate points along a line.
//! - **[`Densify`]**: Insert points into a geometry so there is never more than `max_segment_length` between points.
//!
//! ### Misc measures
//!
//! - **[`HausdorffDistance`]**: Calculate "the maximum of the distances from a point in any of the sets to the nearest point in the other set." (Rote, 1991)
//! - **[`VincentyDistance`]**: Calculate the minimum geodesic distance between geometries using Vincenty’s formula
//! - **[`VincentyLength`]**: Calculate the geodesic length of a geometry using Vincenty’s formula
//! - **[`FrechetDistance`]**: Calculate the similarity between [`LineString`]s using the Fréchet distance
//!
//! ## Area
//!
//! - **[`Area`]**: Calculate the planar area of a geometry
//! - **[`ChamberlainDuquetteArea`]**: Calculate the geodesic area of a geometry on a sphere using the algorithm presented in _Some Algorithms for Polygons on a Sphere_ by Chamberlain and Duquette (2007)
//! - **[`GeodesicArea`]**: Calculate the geodesic area and perimeter of a geometry on an ellipsoid using the algorithm presented in _Algorithms for geodesics_ by Charles Karney (2013)
//!
//! ## Boolean Operations
//!
//! - **[`BooleanOps`]**: combine or split (Multi)Polygons using intersecton, union, xor, or difference operations.
//!
//! ## Outlier Detection
//!
//! - **[`OutlierDetection`]**: Detect outliers in a group of points using [LOF](https://en.wikipedia.org/wiki/Local_outlier_factor)
//!
//! ## Simplification
//!
//! - **[`Simplify`]**: Simplify a geometry using the Ramer–Douglas–Peucker algorithm
//! - **[`SimplifyIdx`]**: Calculate a simplified geometry using the Ramer–Douglas–Peucker algorithm, returning coordinate indices
//! - **[`SimplifyVw`]**: Simplify a geometry using the Visvalingam-Whyatt algorithm
//! - **[`SimplifyVwPreserve`]**: Simplify a geometry using a topology-preserving variant of the Visvalingam-Whyatt algorithm
//! - **[`SimplifyVwIdx`]**: Calculate a simplified geometry using the Visvalingam-Whyatt algorithm, returning coordinate indices
//!
//! ## Query
//!
//! - **[`ClosestPoint`]**: Find the point on a geometry
//!   closest to a given point
//! - **[`HaversineClosestPoint`]**: Find the point on a geometry
//!   closest to a given point on a sphere using spherical coordinates and lines being great arcs.
//! - **[`IsConvex`]**: Calculate the convexity of a
//!   [`LineString`]
//! - **[`LineInterpolatePoint`]**:
//!   Generates a point that lies a given fraction along the line
//! - **[`LineLocatePoint`]**: Calculate the
//!   fraction of a line’s total length representing the location of the closest point on the
//!   line to the given point
//!
//! ## Topology
//!
//! - **[`Contains`]**: Calculate if a geometry contains another
//!   geometry
//! - **[`CoordinatePosition`]**: Calculate
//!   the position of a coordinate relative to a geometry
//! - **[`HasDimensions`]**: Determine the dimensions of a geometry
//! - **[`Intersects`]**: Calculate if a geometry intersects
//!   another geometry
//! - **[`line_intersection`]**: Calculates the
//!   intersection, if any, between two lines.
//! - **[`Relate`]**: Topologically relate two geometries based on
//!   [DE-9IM](https://en.wikipedia.org/wiki/DE-9IM) semantics.
//! - **[`Within`]**: Calculate if a geometry lies completely within another geometry.
//!
//! ## Triangulation
//!
//! - **[`TriangulateEarcut`](triangulate_earcut)**: Triangulate polygons using the earcut algorithm. Requires the `"earcutr"` feature, which is enabled by default.
//!
//! ## Winding
//!
//! - **[`Orient`]**: Apply a specified winding [`Direction`](orient::Direction) to a [`Polygon`]’s interior and exterior rings
//! - **[`Winding`]**: Calculate and manipulate the [`WindingOrder`](winding_order::WindingOrder) of a [`LineString`]
//!
//! ## Iteration
//!
//! - **[`CoordsIter`]**: Iterate over the coordinates of a geometry
//! - **[`MapCoords`]**: Map a function over all the coordinates
//!   in a geometry, returning a new geometry
//! - **[`MapCoordsInPlace`]**: Map a function over all the
//!   coordinates in a geometry in-place
//! - **[`LinesIter`]**: Iterate over lines of a geometry
//!
//! ## Boundary
//!
//! - **[`BoundingRect`]**: Calculate the axis-aligned
//!   bounding rectangle of a geometry
//! - **[`MinimumRotatedRect`]**: Calculate the
//!   minimum bounding box of a geometry
//! - **[`ConcaveHull`]**: Calculate the concave hull of a
//!   geometry
//! - **[`ConvexHull`]**: Calculate the convex hull of a
//!   geometry
//! - **[`Extremes`]**: Calculate the extreme coordinates and
//!   indices of a geometry
//!
//! ## Affine transformations
//!
//! - **[`Rotate`]**: Rotate a geometry around its centroid
//! - **[`Scale`]**: Scale a geometry up or down by a factor
//! - **[`Skew`]**: Skew a geometry by shearing angles along the `x` and `y` dimension
//! - **[`Translate`]**: Translate a geometry along its axis
//! - **[`AffineOps`]**: generalised composable affine operations
//!
//! ## Conversion
//!
//! - **[`Convert`]**: Convert (infalliby) the type of a geometry’s coordinate value
//! - **[`TryConvert`]**: Convert (falliby) the type of a geometry’s coordinate value
//! - **[`ToDegrees`]**: Radians to degrees coordinate transforms for a given geometry.
//! - **[`ToRadians`]**: Degrees to radians coordinate transforms for a given geometry.
//!
//! ## Miscellaneous
//!
//! - **[`Centroid`]**: Calculate the centroid of a geometry
//! - **[`ChaikinSmoothing`]**: Smoothen `LineString`, `Polygon`, `MultiLineString` and `MultiPolygon` using Chaikin's algorithm.
//! - **[`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.
//!
//! # Spatial Indexing
//!
//! `geo` geometries (`Point`, `Line`, `LineString`, `Polygon`) can be used with the [rstar](https://docs.rs/rstar/0.12.0/rstar/struct.RTree.html#usage)
//! R*-tree crate for fast distance and nearest-neighbour queries. Multi- geometries can be added to the tree by iterating over
//! their members and adding them. Note in particular the availability of the [`bulk_load`](https://docs.rs/rstar/0.12.0/rstar/struct.RTree.html#method.bulk_load)
//! method and [`GeomWithData`](https://docs.rs/rstar/0.12.0/rstar/primitives/struct.GeomWithData.html) struct.
//!
//! # Features
//!
//! The following optional [Cargo features] are available:
//!
//! - `earcutr`:
//!     - Enables the `earcutr` crate, which provides triangulation of polygons using the earcut algorithm.
//!     - ☑ Enabled by default.
//! - `proj-network`:
//!     - Enables [network grid] support for the [`proj` crate].
//!     - After enabling this feature, [further configuration][proj crate file download] is required to use the network grid.
//!     - ☐ Disabled by default.
//! - `use-proj`:
//!     - Enables coordinate conversion and transformation of `Point` geometries using the [`proj` crate]
//!     - ☐ Disabled by default.
//! - `use-serde`:
//!     - Allows geometry types to be serialized and deserialized with [Serde].
//!     - ☐ Disabled by default.
//! - `multithreading`:
//!     - Enables multithreading support for the `i_overlay` crate (via Rayon), and activates the `multithreading` flag
//!       in `geo-types`, enabling multi-threaded iteration over `Multi*` geometries.
//!     - ☑ Enabled by default.
//!
//! # Ecosystem
//!
//! There’s a wide variety of `geo`-compatible crates in the ecosystem that offer functionality not
//! included in the `geo` crate, including:
//!
//! * Reading and writing file formats (e.g. [GeoJSON][geojson crate], [WKT][wkt crate],
//!   [shapefile][shapefile crate])
//! * [Latitude and longitude parsing][latlng crate]
//! * [Label placement][polylabel crate]
//! * [Geocoding][geocoding crate]
//! * [and much more...][georust website]
//!
//! [Euclidean plane]: https://en.wikipedia.org/wiki/Euclidean_plane
//! [`geo-types`]: https://crates.io/crates/geo-types
//! [haversine formula]: https://en.wikipedia.org/wiki/Haversine_formula//
//! [`proj` crate]: https://github.com/georust/proj
//! [geojson crate]: https://crates.io/crates/geojson
//! [Karney (2013)]:  https://arxiv.org/pdf/1109.4448.pdf
//! [wkt crate]: https://crates.io/crates/wkt
//! [shapefile crate]: https://crates.io/crates/shapefile
//! [latlng crate]: https://crates.io/crates/latlon
//! [polylabel crate]: https://crates.io/crates/polylabel
//! [geocoding crate]: https://crates.io/crates/geocoding
//! [georust website]: https://georust.org
//! [Cargo features]: https://doc.rust-lang.org/cargo/reference/features.html
//! [GEOS]: https://trac.osgeo.org/geos
//! [JTS]: https://github.com/locationtech/jts
//! [network grid]: https://proj.org/usage/network.html
//! [OGC-SFA]: https://www.ogc.org/standards/sfa
//! [proj crate file download]: https://docs.rs/proj/*/proj/#grid-file-download
//! [rhumb line]: https://en.wikipedia.org/wiki/Rhumb_line
//! [Serde]: https://serde.rs/

#[cfg(feature = "use-serde")]
#[macro_use]
extern crate serde;

pub use crate::algorithm::*;
pub use crate::types::Closest;
use std::cmp::Ordering;

pub use crate::relate::PreparedGeometry;
pub use geo_types::{coord, line_string, point, polygon, wkt, CoordFloat, CoordNum};

pub mod geometry;
pub use geometry::*;

/// This module includes all the functions of geometric calculations
pub mod algorithm;
mod geometry_cow;
mod types;
mod utils;
use crate::kernels::{RobustKernel, SimpleKernel};
pub(crate) use geometry_cow::GeometryCow;

#[cfg(test)]
#[macro_use]
extern crate approx;

#[macro_use]
extern crate log;

/// Mean radius of Earth in meters
/// This is the value recommended by the IUGG:
/// Moritz, H. (2000). Geodetic Reference System 1980. Journal of Geodesy, 74(1), 128–133. doi:10.1007/s001900050278
/// "Derived Geometric Constants: mean radius" (p133)
/// https://link.springer.com/article/10.1007%2Fs001900050278
/// https://sci-hub.se/https://doi.org/10.1007/s001900050278
/// https://en.wikipedia.org/wiki/Earth_radius#Mean_radius
const MEAN_EARTH_RADIUS: f64 = 6371008.8;

// Radius of Earth at the equator in meters (derived from the WGS-84 ellipsoid)
const EQUATORIAL_EARTH_RADIUS: f64 = 6_378_137.0;

// Radius of Earth at the poles in meters (derived from the WGS-84 ellipsoid)
const POLAR_EARTH_RADIUS: f64 = 6_356_752.314_245;

// Flattening of the WGS-84 ellipsoid - https://en.wikipedia.org/wiki/Flattening
const EARTH_FLATTENING: f64 =
    (EQUATORIAL_EARTH_RADIUS - POLAR_EARTH_RADIUS) / EQUATORIAL_EARTH_RADIUS;

/// A prelude which re-exports the traits for manipulating objects in this
/// crate. Typically imported with `use geo::prelude::*`.
pub mod prelude {
    pub use crate::algorithm::*;
}

/// A common numeric trait used for geo algorithms
///
/// Different numeric types have different tradeoffs. `geo` strives to utilize generics to allow
/// users to choose their numeric types. If you are writing a function which you'd like to be
/// generic over all the numeric types supported by geo, you probably want to constrain
/// your function input to `GeoFloat`. For methods which work for integers, and not just floating
/// point, see [`GeoNum`].
///
/// # Examples
///
/// ```
/// use geo::{GeoFloat, MultiPolygon, Polygon, Point};
///
/// // An admittedly silly method implementation, but the signature shows how to use the GeoFloat trait
/// fn farthest_from<'a, T: GeoFloat>(point: &Point<T>, polygons: &'a MultiPolygon<T>) -> Option<&'a Polygon<T>> {
///     polygons.iter().fold(None, |accum, next| {
///         match accum {
///             None => Some(next),
///             Some(farthest) => {
///                 use geo::{euclidean_distance::EuclideanDistance};
///                 if next.euclidean_distance(point) > farthest.euclidean_distance(point) {
///                     Some(next)
///                 } else {
///                     Some(farthest)
///                 }
///             }
///         }
///     })
/// }
/// ```
pub trait GeoFloat:
    GeoNum + num_traits::Float + num_traits::Signed + num_traits::Bounded + float_next_after::NextAfter
{
}
impl<T> GeoFloat for T where
    T: GeoNum
        + num_traits::Float
        + num_traits::Signed
        + num_traits::Bounded
        + float_next_after::NextAfter
{
}

/// A trait for methods which work for both integers **and** floating point
pub trait GeoNum: CoordNum {
    type Ker: Kernel<Self>;

    /// Return the ordering between self and other.
    ///
    /// For integers, this should behave just like [`Ord`].
    ///
    /// For floating point numbers, unlike the standard partial comparison between floating point numbers, this comparison
    /// always produces an ordering.
    ///
    /// See [f64::total_cmp](https://doc.rust-lang.org/src/core/num/f64.rs.html#1432) for details.
    fn total_cmp(&self, other: &Self) -> Ordering;
}

macro_rules! impl_geo_num_for_float {
    ($t: ident) => {
        impl GeoNum for $t {
            type Ker = RobustKernel;
            fn total_cmp(&self, other: &Self) -> Ordering {
                self.total_cmp(other)
            }
        }
    };
}
macro_rules! impl_geo_num_for_int {
    ($t: ident) => {
        impl GeoNum for $t {
            type Ker = SimpleKernel;
            fn total_cmp(&self, other: &Self) -> Ordering {
                self.cmp(other)
            }
        }
    };
}

// This is the list of primitives that we support.
impl_geo_num_for_float!(f32);
impl_geo_num_for_float!(f64);
impl_geo_num_for_int!(i16);
impl_geo_num_for_int!(i32);
impl_geo_num_for_int!(i64);
impl_geo_num_for_int!(i128);
impl_geo_num_for_int!(isize);

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

    #[test]
    fn total_ord_float() {
        assert_eq!(GeoNum::total_cmp(&3.0f64, &2.0f64), Ordering::Greater);
        assert_eq!(GeoNum::total_cmp(&2.0f64, &2.0f64), Ordering::Equal);
        assert_eq!(GeoNum::total_cmp(&1.0f64, &2.0f64), Ordering::Less);
        assert_eq!(GeoNum::total_cmp(&1.0f64, &f64::NAN), Ordering::Less);
        assert_eq!(GeoNum::total_cmp(&f64::NAN, &f64::NAN), Ordering::Equal);
        assert_eq!(GeoNum::total_cmp(&f64::INFINITY, &f64::NAN), Ordering::Less);
    }

    #[test]
    fn total_ord_int() {
        assert_eq!(GeoNum::total_cmp(&3i32, &2i32), Ordering::Greater);
        assert_eq!(GeoNum::total_cmp(&2i32, &2i32), Ordering::Equal);
        assert_eq!(GeoNum::total_cmp(&1i32, &2i32), Ordering::Less);
    }

    #[test]
    fn numeric_types() {
        let _n_i16 = Point::new(1i16, 2i16);
        let _n_i32 = Point::new(1i32, 2i32);
        let _n_i64 = Point::new(1i64, 2i64);
        let _n_i128 = Point::new(1i128, 2i128);
        let _n_isize = Point::new(1isize, 2isize);
        let _n_f32 = Point::new(1.0f32, 2.0f32);
        let _n_f64 = Point::new(1.0f64, 2.0f64);
    }
}