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
#![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
//!
//! ## 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
//!
//! ## Distance
//!
//! - **[`EuclideanDistance`]**: Calculate the minimum euclidean distance between geometries
//! - **[`GeodesicDistance`]**: Calculate the minimum geodesic distance between geometries using the algorithm presented in _Algorithms for geodesics_ by Charles Karney (2013)
//! - **[`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)
//! - **[`HaversineDistance`]**: Calculate the minimum geodesic distance between geometries using the haversine formula
//! - **[`RhumbDistance`]**: Calculate the length of a rhumb line connecting the two geometries
//! - **[`VincentyDistance`]**: Calculate the minimum geodesic distance between geometries using Vincenty’s formula
//!
//! ## Length
//!
//! - **[`EuclideanLength`]**: Calculate the euclidean length of a geometry
//! - **[`GeodesicLength`]**: Calculate the geodesic length of a geometry using the algorithm presented in _Algorithms for geodesics_ by Charles Karney (2013)
//! - **[`HaversineLength`]**: Calculate the geodesic length of a geometry using the haversine formula
//! - **[`RhumbLength`]**: Calculate the length of a geometry assuming it's composed of rhumb lines
//! - **[`VincentyLength`]**: Calculate the geodesic length of a geometry using Vincenty’s formula
//!
//! ## 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
//!
//! - **[`HaversineBearing`]**: Calculate the bearing between points using great circle calculations.
//! - **[`GeodesicBearing`]**: Calculate the bearing between points on a [geodesic](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid)
//! - **[`RhumbBearing`]**: Calculate the angle from north of the rhumb line connecting two points.
//! - **[`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
//!
//! ## Similarity
//!
//! - **[`FrechetDistance`]**: Calculate the similarity between [`LineString`]s using the Fréchet distance
//!
//! ## 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).
//!
//! ## 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.
//! - **[`Densify`]**: Densify linear geometry components by interpolating points
//! - **[`DensifyHaversine`]**: Densify spherical geometry by interpolating points on a sphere
//! - **[`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
//! - **[`HaversineIntermediate`]**: Calculate intermediate points on a sphere along a great-circle line
//! - **[`RhumbDestination`]**: Given a start point, bearing, and distance, calculate the destination point on a sphere assuming travel along a rhumb line
//! - **[`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.
//!
//! # Features
//!
//! The following optional [Cargo features] are available:
//!
//! - `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
//! - `use-proj`: Enables coordinate conversion and transformation of `Point` geometries using the [`proj` crate]
//! - `use-serde`: Allows geometry types to be serialized and deserialized with [Serde]
//!
//! # 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]
//!
//! [`geo-types`]: https://crates.io/crates/geo-types
//! [`proj` crate]: https://github.com/georust/proj
//! [geojson crate]: https://crates.io/crates/geojson
//! [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
//! [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 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. It should match our set of implementations for HasKernel since GeoNum
// depends on HasKernel.
impl_geo_num_for_float!(f32);
impl_geo_num_for_float!(f64);
impl_geo_num_for_int!(i64);
impl_geo_num_for_int!(i32);
impl_geo_num_for_int!(i16);
impl_geo_num_for_int!(isize);
#[cfg(has_i128)]
impl_geo_num_for_int!(i128);

#[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);
    }
}