georust / geo

Geospatial primitives and algorithms for Rust

geo

Geospatial Primitives, Algorithms, and Utilities

Chat or ask questions on Discord

The geo crate provides geospatial primitive types such as Point, LineString, and Polygon, and provides algorithms and operations such as:

• Area and centroid calculation
• Simplification and convex hull operations
• Euclidean and Haversine distance measurement
• Intersection checks
• Affine transforms such as rotation and translation
• All DE-9IM spatial predicates such as contains, crosses, and touches.

Please refer to the documentation for a complete list.

The primitive types also provide the basis for other functionality in the Geo ecosystem, including:

• Coordinate transformation and projection
• Serialization to and from GeoJSON and WKT
• Geocoding
• Working with GPS data

Example

// primitives
use geo::{line_string, polygon};

// algorithms
use geo::ConvexHull;

// An L shape
let poly = polygon![
    (x: 0.0, y: 0.0),
    (x: 4.0, y: 0.0),
    (x: 4.0, y: 1.0),
    (x: 1.0, y: 1.0),
    (x: 1.0, y: 4.0),
    (x: 0.0, y: 4.0),
    (x: 0.0, y: 0.0),
];

// Calculate the polygon's convex hull
let hull = poly.convex_hull();

assert_eq!(
    hull.exterior(),
    &line_string![
        (x: 4.0, y: 0.0),
        (x: 4.0, y: 1.0),
        (x: 1.0, y: 4.0),
        (x: 0.0, y: 4.0),
        (x: 0.0, y: 0.0),
        (x: 4.0, y: 0.0),
    ]
);

Contributing

Contributions are welcome! Have a look at the issues, and open a pull request if you'd like to add an algorithm or some functionality.

License

Licensed under either of

• Apache License, Version 2.0 (LICENSE-APACHE or http://www.apache.org/licenses/LICENSE-2.0)
• MIT license (LICENSE-MIT or http://opensource.org/licenses/MIT)

at your option.

Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted for inclusion in the work by you, as defined in the Apache-2.0 license, shall be dual licensed as above, without any additional terms or conditions.

geo's Issues

Floating point comparison test failure

On current master:

---- algorithm::distance::test::line_segment_distance_test stdout ----
        thread 'algorithm::distance::test::line_segment_distance_test' panicked at 'assertion failed: `(left == right)` (left: `2.048590078926335`, right: `2.0485900789263356`)', src/algorithm/distance.rs:194
note: Run with `RUST_BACKTRACE=1` for a backtrace.

The difference is that the last digit of left is not included in right. Exact comparison of floating points is always somewhat iffy - what do you think about moving to a "close-enough" comparison?

I can implement the sample code from here: http://floating-point-gui.de/errors/comparison/ as a function, and replace the assert_eq! with assert!(float_eq(actual, expected))

point types

how do we want to handle the difference between x,y, x,y,z, x,y,z,w (someone could conceivably want a 4th digit for something like time). Possibilities

• WKT style: use templates or something to have distinct PointXY, PoingXYZ, types a draw back is that if things assume 2d points having a 3d point could break stuff
• GeoJSON style: make it a vector functions would have to double check that that enough dimensions where present and the type system wouldn't be able to help

Stop reexporting struct types in 'centroid.rs'

We don't need the pub here.

This results in the structs getting added within the module which is not what we want: https://georust.github.io/rust-geo/geo/algorithm/centroid/index.html

Figure out if it's possible to make the algorithm traits generically commutative

e.g. if we implement impl Intersects<Polygon> for LineString, is it possible to write a generic trait implementation such that we automatically get impl Intersects<LineString> for Polygon implemented?

don't mix x/y and lat/lon

The Point struct has both x/y and lat/lon getter methods. This can be confusing when both are used in the same method.

Additionally i think using lat/lon is very dangerous in general. My experience is that many people think x = lat and y = lon ...

How best to calculate Point to Polygon distance?

I'm working on an algorithm that requires calculating the minimum signed distance from a Point to a Polygon exterior (negative implies that the point lies inside the Polygon / polygon.contains(point) is true).
I can't see a good way to do this, currently – I've implemented point-line-distance in simplify, but that's for segments, and there's currently no way to get Polygon exterior ring segments. Or is there? If anyone can think of an approach here, or point to an existing implementation I'd be very grateful (and I'll implement it ASAP)

Spatial Reference System of geometries

Spatial data is always represented in a spatial / coordinate reference system. While most data, e.g. from GPS sensors, is represented in WGS84 (LatLon) this is not always the case.

Currently geometries in rust-geo have no information about a coordinate reference system. This allows to compare a LatLon- Point with an UTM- Point.

In JTS and GEOS each Geometry has a CRS-id field. Algorithms have to check if the CRS of two geometries are the same. This could also be done for geometries in rust-geo.

Maybe it would also be possible to use the advanced type-system of Rust to add a generic CRS-type to geometries....

Create `Union` algorithm trait

Similar to #80

Geometry type system design propersal

I'm not native English speaker, if you have any question, just point it out and I'll try my best to make it clear.

After we reached agreement to design the type system against some open standard, it's time to go further.

My idea is to make an enum Geometry as the entry point of our type system. It will be like the Json enumeration in rustc-serialize, because in Rust, enum is a good way for polymorphism. When we are parsing geometry objects from Geojson or WKT, we might have no idea about the concrete type of them so we can use Geometry. Concretely, we will have a FromWkt trait and a ToWkt trait for I/O on Geometry without knowing the exact type. We can use pattern matching then to deal with different types. So other operations can also be impled on Geometry, like DE-9IM. It computes relationship between two objects (Point-LineString, Point-Point, LineString-Polygon or whatever).

The value of enum will be the structs for different types. We still need concrete types because some operations or attributes are unique to some type.

enum Geometry {
  Point (Point),
  LineString (LineString),
  Polygon (Polygon),

  ...
}

We many have some functions like as_point(&self) -> Option<Point> to convert it from Geometry to Point.

I will be working on a simple implementation for this. Any question is welcomed.

Seemingly unnecessary bounds on `Translate`

Is the FromPrimitive bound still necessary for the Translate trait's blanket impl? It doesn't look like MapCoords requires FromPrimitive, and I can't see anything which explicitly requires it.

The code in question:

impl<T, G> Translate<T> for G
    where T: Float + FromPrimitive,
        G: MapCoords<T, T, Output=G>
{
    fn translate(&self, xoff: T, yoff: T) -> Self {
        self.map_coords(&|&(x, y)| (x + xoff, y + yoff))
    }
}

This isn't really an issue in practice because all types implementing num_traits::Float are also FromPrimitive, but it feels a bit odd to add the extra constraint.

---

Also, I don't think the & in front of the closure is necessary. MapCoords::map_coords() should be generic over anything which is Fn(&(T, T)) -> (NT, NT) so you can skip the extra syntax and dynamic dispatch, although fixing that now would be a breaking change...

Idea: Implement constraints on geometry types

Examples:

• LineString: panics if fewer than two points
• Polygon: panics if any rings has fewer than three points
• Polygon should check if first+last points are the same?
• Polygon should ensure interiors are inside exterior? Maybe not this is one as it might be expensive.

`centroid` returns incorrect result when multiplogon consists of polygon in different clockwise

Consider:
A Multi polygon consists of two polygon [p(0., 0.), p(2., 0.), p(2., 2.), p(0., 2.), p(0., 0.)] (clockwise) and [p(0., 0.), p(-2., 0.), p(-2., 2.), p(0., 2.), p(0., 0.)](anti-clockwise), whose centroid is (0,1).

Current computation gives an incorrect result (Nan, Nan).

Perpendicular distance and LineString simplification

Is there any interest in a simplify method for LineStrings, using Douglas–Peucker? I've just finished implementing it in a separate crate, and then realised it might be useful here. This would also supply a perpendicular distance (from a LineString to a Point) method…

Possible bug with polygon to polygon intersections.

I needed to divide a region (polygon) to 50x50 meter cells, what I did was looping through the width and height of the region's bounding box and create a cell for every step then test it via bigger_polygon.intersects(&cell_polygon). If I'm not doing stupid mistake it seems like a there is bug on calculating intersections.

And here is the visualization, as we can see red square shouldn't return true for intersection test.

Here is the failing test case. POINTS array are black polygon in the visual and small_polygon is the red square.

#[cfg(test)]
mod test {
    use geo::Polygon;
    use geo::intersects::Intersects;

    #[test]
    fn interstects() {
        let big_polygon = Polygon::new(POINTS.to_vec().into(), vec![]);

        // Just a 50x50 meters rect
        let small_polygon = Polygon::new(vec![
            ( 127.0072984994905, 37.577366522199405 ),
            ( 127.0078658729216, 37.577366520841686 ),
            ( 127.0078658729216, 37.57691686003872 ),
            ( 127.0072984994905, 37.57691686139645 )
        ].into(), vec![]);

        // It shouldn't intersect with bigger polygon
        assert_eq!(big_polygon.intersects(&small_polygon), false);
    }

    pub const POINTS: [(f64, f64); 26] = [
        ( 126.99294090270998, 37.5772698028868 ),
        ( 126.99240982532503, 37.57736758665555 ),
        ( 126.99149787425996, 37.5764237553521 ),
        ( 126.99141204357149, 37.57567548419531 ),
        ( 126.99364364147186, 37.57203605817081 ),
        ( 126.99377775192261, 37.57054367201618 ),
        ( 126.99527978897093, 37.57054792387107 ),
        ( 126.99573040008544, 37.56655107315775 ),
        ( 127.00274169445039, 37.56686572662952 ),
        ( 127.00982272624971, 37.56547103626546 ),
        ( 127.01042890548705, 37.5659387707002 ),
        ( 127.01180756092073, 37.56629594774544 ),
        ( 127.01261758804321, 37.56694226375953 ),
        ( 127.01248347759248, 37.56781818309174 ),
        ( 127.009693980217,   37.568872674254486 ),
        ( 127.00952231884003, 37.570369345756944 ),
        ( 127.010595202446,   37.57149182964705 ),
        ( 127.00900733470918, 37.573056475933804 ),
        ( 127.00767695903781, 37.57298844851792 ),
        ( 127.00883567333221, 37.571253728417325 ),
        ( 126.99647605419162, 37.57070949417681 ),
        ( 126.99636876583098, 37.57202755463073 ),
        ( 126.99677646160124, 37.572771610715066 ),
        ( 126.99616760015486, 37.574859179812385 ),
        ( 126.99551582336427, 37.57638974318991 ),
        ( 126.99294090270998, 37.5772698028868 )
    ];
}

Add a line/line-segment type?

Does it make sense to have a separate type for a Line, like https://github.com/paulmach/go.geo/blob/master/line.go (in that it is treated as a LineString for GeoJSON/WKT)?

In particular, I would like this feature for strictly typed conversions using Postgres's builtin Line Segment type (using https://github.com/sfackler/rust-postgres). I had heard that it might be the plan to represent geo::types as traits; in which case this could either be its own trait or an implementation of LineString. (Aside: if you've written up the current direction of georust somewhere, I'd love a link!).

Determine fate of `Geometry` enum

What is the point of this enum? Doesn't seem like it serves much of a purpose. What are qualifications for a struct/variant to be added to it?

Better README

The README could do with some care and attention. Without wanting to duplicate the docs excessively, would it be useful to list some of the present primitives and functionality?

Polygon centroid calculation should account for interior holes

It doesn't at the moment, meaning that rotation around a polygon centroid (see #107) is only possible for simple polygons. If anyone has suggestions for existing clean algorithm implementations, I'm happy to rewrite the existing implementation.

Convex Hull is oriented clockwise

The convex hull implementation outputs a clockwise-oriented exterior ring. This isn't strictly wrong, but many other algorithms assume a counter-clockwise exterior ring. The simplest fix here is to reverse hull with into_iter().rev().collect(), but I'm not sure what its overhead is – it'd be better to fix the implementation to order the edge segments correctly.

Proposal: Introduce `geo::GeoBuf` struct, alter `geo::Geo` struct

Idea

Currently, there exists just geo::Geo, which is a structure containing owned data. I propose we rename this to geo::GeoBuf and introduce a new structure containing borrowed data called geo::Geo. Table of similar analogies:

Owned	Borrowed
Vec	slice
String	str
PathBuf	Path
GeoBuf	Geo

Reasoning

Currently, if there exists a method centerpoint on Geo, and I want to apply this operation on a Wkt struct (from rust-wkt), I have to convert my Wkt struct to a Geo, and since Geo is owned, this results in memory allocations occurring.

With my proposal above, we would just need to interpret the Wkt struct as the new borrowed Geo, which results in a Geo struct that references the existing data in the Wkt struct, resulting in no new memory allocations.

Implementation

struct Geo will be renamed to struct GeoBuf

struct Geo will be a new struct with borrowed fields. Not too sure how this will look right now, but it's do-able.

Deref<Target=Geo> will be implemented on GeoBuf

Polylabel

I've implemented the very nice Mapbox Polylabel algorithm for optimised polygon label placement, here: https://github.com/urschrei/polylabel-rs

I think it should be part of georust. If you agree, should it be:

• an impl on Polygon within rust-geo
  • Should it have the same name, or a different name?
• a separate crate, like Polyline

geo structs

basic geo types need to be defined, some basics

• point, line and poygons
• how to deal with multi versions, separate type or points/line/polygon always store the type in a vector
• minimum number of distinct points in a line (geojson is 2, shapefiles is 1).
• minimum number of distinct points in a polygon (geojson is 3 plus first and last must be the same so 4 total, esri does not have first and last be the same).
• winding order, does it matter (geojson no, esrijson clockwise is a polygon and counterclockwise is a hole).
• if winding order doesn't matter how do you deal with polygons that cross into multiple hemispheres
• how do you deal with holes in polygons (geojson has a polygon object be an outer ring followed by zero or more inner rings, esrijson uses winding order).
• how do you deal with self intersecting polygons.
• how do you deal with polygons which intersect their own holes.
• do we deal with typologies
• how many dimensions
• geometry collections (esri does not touch this)
• projections and that whole can of shit ...

Implement `Intersects<Polygon<T>> for Polygon`

https://github.com/georust/rust-geo/blob/master/src/algorithm/intersects.rs

https://georust.github.io/rust-geo/geo/algorithm/intersects/trait.Intersects.html

Implement validity checks for each geometry type

See discussion in #123

Validity checks are described here: https://postgis.net/docs/using_postgis_dbmanagement.html#OGC_Validity

e.g. Point::is_simple, LineString::is_simple, etc.

Prelude

A common Rust pattern is to expose some top-level prelude module which re-exports all the commonly-used traits in a module.

Seeing as a lot of the functionality for various operations is separated into their various submodules in geo::algorithms, and you have to manually import each trait, would it be possible to expose a geo::prelude module which people can glob import (e.g. use geo::prelude::*)?

Change COORD_PRECISION from f64 to f32.

https://github.com/georust/rust-geo/blob/e23eaa69749c1279d7e41cf38e6335fcec324b23/src/types.rs#L7

f64 is probably overkill for a value of 1e-1.

Provide basic operations for geometry types

http://turfjs.org/static/docs/ has a nice list of operations. Might be good to have a discussion about whether to have these in a separate library is more appropriate

Introduce trait for accessing geometries

In the course of refactoring rust-postgis to use rust-geo geometries, I'm looking at needed conversions between different georust geometries. One goal would be to store a rust-gdal geometry with rust-postgis as efficient as possible.

One possible solution would be using iterators returning a Point trait instead of accessing the structs directly. If rust-gdal would implement this trait, storing a geometry with rust-postgis could be done without any conversion. And if the geo algorithms would use this trait, they could be applied directly on geometries implementing this trait. It should also solve the ownership problem in #21.

A simplified trait could be

pub trait PointType {
    fn x(&self) -> f64;
    fn y(&self) -> f64;
    //...
}

When this trait is implemented for geo::Point and geo::LineString implements IntoIterator with Item=&PointType, you can iterate like

let xsum = linestring.into_iter().fold(0., |sum, p| sum + p.x());

Relicense under dual MIT/Apache-2.0

Why?

The MIT license requires reproducing countless copies of the same copyright
header with different names in the copyright field, for every MIT library in
use. The Apache license does not have this drawback, and has protections from
patent trolls and an explicit contribution licensing clause. However, the
Apache license is incompatible with GPLv2. This is why Rust is dual-licensed as
MIT/Apache (the "primary" license being Apache, MIT only for GPLv2 compat), and
doing so would be wise for this project. This also makes this crate suitable
for inclusion in the Rust standard distribution and other project using dual
MIT/Apache.

How?

To do this, get explicit approval from each contributor of copyrightable work
(as not all contributions qualify for copyright) and then add the following to
your README:

## License

Licensed under either of
 * Apache License, Version 2.0 ([LICENSE-APACHE](LICENSE-APACHE) or http://www.apache.org/licenses/LICENSE-2.0)
 * MIT license ([LICENSE-MIT](LICENSE-MIT) or http://opensource.org/licenses/MIT)
at your option.

### Contribution

Unless you explicitly state otherwise, any contribution intentionally submitted
for inclusion in the work by you shall be dual licensed as above, without any
additional terms or conditions.

and in your license headers, use the following boilerplate (based on that used in Rust):

// Copyright (c) 2015 t developers
// Licensed under the Apache License, Version 2.0
// <LICENSE-APACHE or
// http://www.apache.org/licenses/LICENSE-2.0> or the MIT
// license <LICENSE-MIT or http://opensource.org/licenses/MIT>,
// at your option. All files in the project carrying such
// notice may not be copied, modified, or distributed except
// according to those terms.

And don't forget to update the license metadata in your Cargo.toml!

Contributor checkoff

• @frewsxcv
• @sunng87
• @Turbo87
• @mgax
• @groteworld

Using Simplification algorithm on a polygon can result in an invalid polygon

#135 (comment)

Maybe the algorithm should sometimes return an Option<T> since it can sometimes fail.

Change struct property names: Point {lat, lng} -> {x, y}

https://github.com/georust/rust-geo/blob/38fe34109bb7dd2001c0ea9153b9c1dca98a44e8/src/point.rs

It's more concise and more aligned with existing specifications. Don't worry about Y and Z for now

f32 vs f64

Hello,

I've noticed that when computing haversine distance between points yields wildly different results when using f32 vs f64 for coordinate types.

extern crate geo;

use geo::Point;
use geo::algorithm::haversine_distance::HaversineDistance;

#[warn(dead_code)]
fn distance_32(){
    let p = Point::new(-77.036585f32, 38.897448f32);
    let dist = p.haversine_distance(&Point::new(-77.009080f32, 38.889825f32));
    println!("Distance = {:?}", dist)
}

#[warn(dead_code)]
fn distance_64(){
    let p = Point::new(-77.036585f64, 38.897448f64);
    let dist = p.haversine_distance(&Point::new(-77.009080f64, 38.889825f64));
    println!("Distance = {:?}", dist)
}


fn main() {

    distance_32();
    distance_64();
}

Output

Distance (f32) = 3114.3247
Distance (f64) = 2529.6506398715096

Is this normal?

Thanks!

Generic numeric type for all structs

We should allow the user to use whatever numeric type they want (e.g. f64, f32, i64). Just need to parameterize it.

Polygon-LineString containment check seems to be wrong

Given an approximately circular Polygon, and an approximately circular concentric LineString:

The containment check should return true - it doesn't.

Polygon
LineString

Dealing with Validity

See previous discussions in #56 and #122
First, what do we mean when we talk about validity? The OGC OpenGIS spec discussed in the PostGIS docs gives a good overview.

• It's probably not desirable to have validating constructors by default, as validation can incur a significant performance overhead, but they should probably be available, either as separate constructors, or with a validate=True argument, although the latter could cause a lot of headaches in terms of revising calls to the existing constructors, and dealing with invalid geometries; e.g. both of the simplification algorithms can return invalid geometries, so those traits would always have to return a Result (they should arguably already be doing this)
• Should we check for validity before carrying out certain operations? This is what Shapely does, raising an error if invalid geometries are detected. I'm not sure what JTS / GEOS do.
• I think the two points above directly imply that every geometry operation must return a Result.

Possible bug in distance between point and polygon

Hi everybody, I'm starting to use this library, and while experimenting I found a possible bug.

I tried this both on the pinned version on crate.io and the master branch here.

Basically in some cases asking for the distance between a polygon and a point that lies on the polygon's exterior ring will not return 0.

Here is a basic main.rs to reproduce the problem.

extern crate geo;

use geo::{Point, Polygon, LineString};
use geo::distance::Distance;


fn main() {
    let exterior = LineString(vec![
                              Point::new(0., 0.),
                              Point::new(0., 0.0004),
                              Point::new(0.0004, 0.0004),
                              Point::new(0.0004, 0.),
                              Point::new(0., 0.)]);

    let poly = Polygon::new(exterior.clone(), vec![]);

    let good_point_1 = Point::new(0.0001, 0.0001);
    let good_point_2 = Point::new(0.0002, 0.);
    let bugged_point = Point::new(0.0001, 0.);

    assert_eq!(poly.distance(&good_point_1), 0.);
    assert_eq!(poly.distance(&good_point_2), 0.);
    // This fails!
    assert_eq!(poly.distance(&bugged_point), 0.);
}

As you see the point (0.0001, 0.) should be inside the polygon (which goes from x=0 to x=0.0004), but the distance result is instead 0.00000000000000000006776263578034403

Could you confirm this is not the expected behaviour?

Update MapCoords to lazily return an iterator instead of always collecting

See my comment here: #136 (comment)

Coordinate vs. Point

what is the purpose of having a Coordinate and a Point type? it seems to me that having Point should be sufficient.

Rebrand?

From #rust-geo on IRC:

6:11 PM <@dremonkey> also I have a good name for rust-geo if interested
6:11 PM <@dremonkey> in changing
6:12 PM <@dremonkey> noticed a lot of rust libs play off of elements or iron... so was thinking lodestone would be kinda fun

Not suggesting we need to, but when (if?) this library ever matures and we ever want to rebrand, ideas can be shared here.

Specifically with the 'lodestone' idea above, I think it's very creative, though we lose a bit of the implicit description of the library (unless there is some large overlap in software developers and geologists I'm not aware of)

semver breakage

geo-0.4.7 depends on serde-1, which is a breaking change. geo should be 0.5.

Document what the parameters to the Polygon struct are

pub struct Polygon<T>(pub LineString<T>, pub Vec<LineString<T>>) where T: Float;

The first item is the exterior ring. The second item is a vector of interior rings. We should either add a doc comment documenting this, or change the struct to something like:

pub struct Polygon<T> {
    pub exterior: pub LineString<T>,
    pub interior: Vec<LineString<T>>
} where T: Float;

Remove unused imports in 'src/algorithm/simplify.rs'

Signed and Unsigned Area

Area returned by current algorithm is signed, which may lead to confusing usage.

Here is what we'd like to propose for enhancement :

• implement UnsingedArea
• Complete the doc

Coordinates are stored as Floats. Could they also be integers?

Currently, the fundamental "unit" is num_traits::Float, i.e. all x/ys (or lat/longs) are represented as Float. And I wonder if this needs to be? I have a project which does things and eventually projects/rasterises geometry to a integer grid, and there's no 100% easy way to store this. Yes I can use Floats, and then truncate at the very end, but that seems subpar. I'm not suggesting removing Float, but expanding what can be stored in/as a Geometry/etc.

Obviously some algorithms wouldn't work on integers (e.g. Area calculation), and the Area trait would still use T: Float.

I have a branch locally which does this, and it's that much of a change to support PrimInt and Float.

What do yous think?

"keywords" typo in Cargo.toml

http://doc.crates.io/manifest.html#package-metadata

https://github.com/georust/rust-geo/blob/master/Cargo.toml#L10

Woops

Implement `Intersects<Bbox<T>> for Polygon`

https://github.com/georust/rust-geo/blob/master/src/algorithm/intersects.rs

https://georust.github.io/rust-geo/geo/algorithm/intersects/trait.Intersects.html

Intersection algorithm says a line intersects with a point not on it

I'm using the following code:

extern crate geo;
use geo::*;
use geo::algorithm::intersects::Intersects;
fn main() {
    let ln = Line::new(Point::new(51.3561738,0.1089333), Point::new(51.3555692,0.1097743));
    let pt = Point::new(51.3566007,0.1083048);
    assert!(ln.intersects(&pt));
}

which fails the assertion:

thread 'main' panicked at 'assertion failed: ln.intersects(&pt)', src/main.rs:7:4
note: Run with `RUST_BACKTRACE=1` for a backtrace.

As far as I can see, the point (51.3566007,0.1083048) is not on the line segment given (although it is on the line, i.e. what you get if you extend the line segment to infinity in both directions):

Implement spherical distance

https://github.com/georust/rust-geo/blob/master/src/algorithm/distance.rs

Currently it looks like rust-geo calculates distance with Euclidean distance, so it's in degrees. That can be useful, but is not typical for geo math: instead, Haversine distance or Vincenty distance is used - the latter is slightly more accurate but also somewhat slower.

Happy to contribute a patch if people agree with this direction.

Create `Intersection` algorithm trait

It might be the case (in the case of intersecting two polygons) that the operation will return a polygon or a linestring or a point or nothing, the return type will need to reflect this. Maybe something like this?

trait Intersection<...> {
    fn intersection(&self, other) -> Option<IntersectionResult> {
        ...
    }
}

// not a fan of "Result" here since it sounds similar to std::result::Result
enum IntersectionResult {   
    Polygon(Polygon),
    LineString(LineString),
    Point(Point),
}

Note that this is different from the Intersects trait which simply returns a bool indicating if two geometries intersect at any point

Extract cross product operation

The cross product of two points is implemented inline in several places. We should implement a cross product method directly on Point.

See also: #132 (comment)