juliageo / geometryops.jl Goto Github PK

View Code? Open in Web Editor NEW

18.0 4.0 3.0 235.77 MB

GeoInterface-based geometry operations

Home Page: https://juliageo.org/GeometryOps.jl/

License: MIT License

Julia 100.00%

geometry geospatial julia

geometryops.jl's Introduction

Warning

This package is still under heavy development! Use with care.

GeometryOps.jl is a package for geometric calculations on (primarily 2D) geometries.

The driving idea behind this package is to unify all the disparate packages for geometric calculations in Julia, and make them GeoInterface.jl-compatible. We are focusing primarily on 2/2.5D geometries for now. All methods in this package will consume any geometry which is compatible with GeoInterface - see its integrations page for more info on that!

Most of the use cases are driven by GIS and similar Earth data workflows, so this might be a bit specialized towards that, but methods should always be general to any coordinate space.

We welcome contributions, either as pull requests or discussion on issues!

Methods

GeometryOps tries to offer most of the basic geometry operations you'd need, implemented in pure Julia and accepting any GeoInterface.jl compatible type.

General geometry methods (OGC methods): equals, extent, distance, crosses, contains, intersects, etc
Targeted function application over large nested geometries (apply) and reduction over geometries (applyreduce)
- Both apply and applyreduce consume arbitrary tables as well, like DataFrames!
signed_area, centroid, distance, etc for valid geometries
Line and polygon simplification (simplify)
Polygon clipping, intersection, difference and union
Generalized barycentric coordinates in polygons (barycentric_coordinates)
Projection of geometries between coordinate reference systems using Proj.jl
Polygonization of raster images by contour detection (polygonize)
Segmentization/densification of geometry, both linearly and by geodesic paths (segmentize)

See the "API" page in the docs for a more complete list!

How to navigate the docs

GeometryOps' docs are divided into three main sections: tutorials, explanations and source code.
Documentation and examples for many functions can be found in the source code section, since we use literate programming in GeometryOps.

Tutorials are meant to teach the fundamental concepts behind GeometryOps, and how to perform certain operations.
Explanations usually contain little code, and explain in more detail how GeometryOps works.
Source code usually contains explanations and examples at the top of the page, followed by annotated source code from that file.

Performance comparison to other packages

From the wonderful vector-benchmark,

More benchmarks coming soon!

Planned additions

Buffering, hulls (convex and otherwise)
Checks for valid geometries (empty linestrings, null points, etc) (#14)
Operations on spherical (non-Euclidean) geometry (#17)

geometryops.jl's People

Contributors

Stargazers

Watchers

Forkers

caltech-octo luntan-maker pablosanjose

geometryops.jl's Issues

`centroid` fails on degenerate polygon

MWE:

import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG

julia> p1 = GI.Polygon([[[0, 0], [1, 0], [2, 0], [0, 0]]])
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}([[0, 0], [1, 0], [2, 0], [0, 0]], nothing, nothing)], nothing, nothing)

julia> p2 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [1, 0], [0, 0]]])
GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Vector{Int64}}, Nothing, Nothing}([[0, 0], [1, 0], [1, 1], [1, 0], [0, 0]], nothing, nothing)], nothing, nothing)

julia> GO.centroid(p1)
(NaN, NaN)

julia> GO.centroid(p2)
(NaN, NaN)

julia> LG.centroid(p1)
POINT (1 0)

julia> LG.centroid(p2)
POINT (0.75 0.25)

Comparison with other libraries

https://github.com/kadyb/vector-benchmark has a bunch of comparisons in the vein of the Julia performance plot, for basic vector geometry ops (no polygon ops, just distance, contains, intersect, reproject/transform, buffer, and IO).

I opened a PR kadyb/vector-benchmark#12 to add GeometryOps benchmarks here as well, but even if that doesn't get merged there we can still run those benchmarks, potentially every release or commit to master (or just when manually requested) over here:

If anyone here knows of other good benchmark sources for geometries, please let me know!

to_points

In the utils.jl file, the function to_points is not working properly. In addition to fixing it, it would be helpful to add an optional argument "repeat_last" which is a boolean. This boolean would specify whether or not to_points returns an array where the first and the last point are the same or not. If "repeat_last" = true, then if the first and last point are not already the same, to_point adds the first point again to the end of the array.

Error-reduced online accumulators

https://github.com/JuliaMath/AccurateArithmetic.jl has this, it would be pretty useful to have support for these in applyreduce for example.

We can even use the accumulators in https://github.com/joshday/OnlineStats.jl as an example.

This would need explicit support and some form of manipulation using areas.

Move Proj.jl to an extension

Would be good to cut the dependencies Proj introduces.

Inexact calculations effects on clippings

Here is an example of inexact calculations causing issues:

scale_fac = 1e5
c1 = [[-28083.868447876892, -58059.13401805979], [-9833.052704767595, -48001.726711609794], [-16111.439295815226, -2.856614689791036e-11], [-76085.95770326033, -2.856614689791036e-11], [-28083.868447876892, -58059.13401805979]]
c1 = [c ./ scale_fac for c in c1]
c2 = [[-53333.333333333336, 0.0], [0.0, 0.0], [0.0, -80000.0], [-60000.0, -80000.0], [-53333.333333333336, 0.0]]
c2 = [c ./ scale_fac for c in c2]
p1 = GI.Polygon([c1])
p2 = GI.Polygon([c2])
GO.intersection(p1, p2; target = GI.PolygonTrait(), fix_multipoly = nothing)

When scale_fac = 1 this doesn't work, but when scale_fac = 1e5 it does work...

Basically, there are two edges over one another on top of the x-axis, but one of the lines is slightly below the x-axis. I am working on a "patch" for this type of issue, but we might want a stronger fix at some point.

WARNING: both GeometryOps and Base export "contains"

I'm not sure how others typically deal with this but does this mean that GeometryOps should not export contains or GeometryOps needs to add a method to Base.contains?

`unsafe_applyreduce` using `GI.getpoint` or `GI.gethole`

@rafaqz mentioned that GI.getpoint or GI.gethole on any geometry can be substantially faster in some cases (Shapefile.jl comes to mind). One can implement an unsafe_applyreduce, which does not guarantee any ordering, but still gives you the points or rings to operate on.

This would operate only on LinearRing, LineString, Point traits.

Simplification of multi-geometries

Currently, when GeometryOps simplifies multi-geometries, it just simplifies each sub-geometry and then reconstructs. However, LibGEOS suggests that entire sub-geometries can be removed. See the example below:

import GeoInterface as GI
import GeometryOps as GO
import LibGEOS as LG

mp_coords = [
    [[[24332.448443033114, 78156.5044920507], [24335.48046004802, 78155.02337099149], [24332.100142171057, 78154.15319217058], [24332.448443033114, 78156.5044920507]]],
    [[[24721.947085257176, 80785.92121346148], [32813.95278360643, 88654.68891218181], [34573.68704928234, 87616.90783669954], [33538.70270722534, 80524.1641959386], [24414.252883053858, 78175.30137343572], [24333.793108370985, 78165.5820218308], [24721.947085257176, 80785.92121346148]]],
]

lg_vals = GI.coordinates(LG.simplify(mp, 100.0))
go_vals= GI.coordinates(GO.simplify(mp; tol = 100.0))

lg_vals is the coordiantes of a polygon, while go_vals is a multipolygon.

Tutorial: adjusting geometry with MakieDraw

Adding this kind of tutorial would be cool, so we can accumulate things people can do with the Julia geo ecosystem here

Specifying `target` all the time in polygon methods is painful

Can we have it default to the greater of the two types given?

signed_distance returns different answer if polygon is explicitly closed or not

import GeometryOps as GO
import GeoInterface as GI

pt = (1.1, 1.0)
rect_open = GI.Polygon([[(1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0)]])
rect_closed = GI.Polygon([[(1.0, 1.0), (1.0, 2.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)]])

GO.signed_distance(pt, rect_open) # return 0.1
GO.signed_distance(pt, rect_closed) # return 0.0

`rebuild` mega-issue

Rebuild doesn't work with GI.Point (because it tries to pass extent as a kwarg). We could special case this behaviour either on our end or on the GeoInterface end.

julia> GO.rebuild((GI.Wrappers.Point(0, 0)), (1, 2, ))
ERROR: MethodError: no method matching (GeoInterface.Wrappers.Point{false, false})(::Tuple{Int64, Int64}; crs::Nothing, extent::Nothing)

Closest candidates are:
  (GeoInterface.Wrappers.Point{Z, M})(::T; crs) where {Z, M, T} got unsupported keyword argument "extent"
   @ GeoInterface ~/.julia/dev/GeoInterface/src/wrappers.jl:272
  (GeoInterface.Wrappers.Point{Z, M})(::Real, ::Real, ::Real...; crs) where {Z, M} got unsupported keyword argument "extent"
   @ GeoInterface ~/.julia/dev/GeoInterface/src/wrappers.jl:278

I'm not sure if this was intended behaviour but rebuild(GB.Point2(0), (1, 2, 3, 4, 5, 6, 7)) creates a Point7 instead of erroring.

Add *hull methods

Like convexhull, it would be nice to do a bunch of these algorithms for completeness https://en.wikipedia.org/wiki/Convex_hull_algorithms

Mostly they have code examples already in the wiki pages

Find arclength to point/find point at arclength

Kind of similar thought to segmentize. Should also have linear and geodesic backends.

TagBot trigger issue

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Spherical geometry

It could be useful to have some spherical geometry (read: geometry on the Earth surface).

In this pull request I implemented a formula to calculate the area of any spherical polygon, which of course sits awkwardly in Rasters.jl.

Some things we might want to implement are (please suggest more):

Shortest distance between two points on a sphere
Bearing
Polygon area

This page has some relevant math

Handle GeoInterface.jl compatible tables

Probably we should allow GeoInterface.jl compatible tables in apply ? (Tables.jl compatible with a :geometry column or something else specified by GeoInterface.geometrycolumns).

But probably it should also return a Table? to keep the "returns a similar object" pattern.

Could we default to using a DataFrame and put DataFrames.jl part in an extension?

This means applyreduce could work with just Tables.jl but apply would need DataFrames.jl

Geometry validation methods

We need some functions like Rs sf_make_valid to validate geometry correctness. R is just using GEOS for this.

Currently we are in a wild-west situation where anything reads, writes and plots and will just randomly fail when it hits methods where validity matters.

force2d, forcez, forcem, forcezm

Shapely has these and it seems pretty useful when dealing with nightmare datasets.

forcez(geoms; default = 0)
forcezm(geoms; defaultz = 1, defaultm = 0)

or something

buffering

These may be useful algorithms:

https://www.tandfonline.com/doi/full/10.1080/10095020.2012.747643

https://ieeexplore.ieee.org/abstract/document/9350239

Features needed for GeoMakie

Cutting arbitrary geometries by a line/linestring (antimeridian cutting and potentially also interpolation)
Arclength interpolation (#71)
Line simplification (given by simplify)
Point of intersection of linestring with rectangle
Direction of linestring / direction in projection / something which can tell me which direction I should push my ticks in

For more general Makie usage I look forward to integrating barycentric rendering into CairoMakie.

Polygon repeated last points

Are we assuming polygons and their rings have a repeated last point? Or are we always checking?

If they don't, should we add an extra last point or just try to make the algorithm work regardless of the repeated last point?

Let all polygon ops start with a bounding box query

If the bounding boxes of the two polygons don't intersect, we can immediately short-circuit the evaluation of any polygon set operation. Would give up to 1000x speedup in some cases, and of course more for really complex polygons.

Point in polygon check failing

For a polygon with multiple holes, the point in polygon check seems to be failing.

import GeometryOps as GO, GeoInterface as GI
using GeometryBasics

rings = [[Point2f.(r*cos(θ), r*sin(θ)) for θ in LinRange(0, 2π, 90)] for r in 5:-1:1]
p = GI.Polygon(rings)
poly(p)

xmin, ymin = minimum(ax.finallimits[])
xmax, ymax = maximum(ax.finallimits[])
points = [Point2f(x, y) for x in LinRange(xmin, xmax, 500), y in LinRange(ymin, ymax, 500)]
f, a, p = heatmap(xmin..xmax, ymin..ymax, GeometryOps.within.(points, (p,)))
poly!(a, p; alpha = 0.6)
f

"Prepared geometry" and other such wrapper types

A lot of packages use "prepared geometries", where points and edges are preprocessed in some way and stored, to reduce processing times for various algorithms.

GeoInterface wrappers kind of already have this with embedded extents. It could be more useful if we had a wrapper class like:

AbstractPreparedGeometry{NumType, CacheTypes, Z, M}

where CacheTypes is a Union of all the cache types the geometry holds in a named tuple, like STRTree or MonotoneChain, etc. Algorithms would also probably have to opt in or out of handling these, and we would need a prepare interface to generate and load such preparations.

maxfreu/STRTrees.jl already exists and works well, so that would be a good starting point!

The algorithm support might be more of an effort, but creating such geometries with cache is low hanging fruit, and ought to be fairly easy to implement.

Clipping on `GO.touches(p1, p2)`

Consider 2 polygons for which GO.touches(p1, p2) == true. In this case, the intersection is generally understood to be a linestring.

Two questions arise:
(a) GeometryOps currently returns an empty polygon for intersect(p1, p2). Should we allow it to return a LineString if the target is set to GeometryCollection? If so, how would that work with the current Foster algorithm?
(b) Do we need to check intersection dimensionality first before constructing the intersection?

Prepared geometry

"Prepared geometry" is basically a regular geometry with an attached cache, which can be a sorted edge list, an Rtree, a monotone chain decomposition, etc.

We should create prepared geometry wrapper structs which hold a geometry + its cache. It may also contain a field verifying whether the geometry has been validated or not. (This will require special handling in e.g. apply, probably).

Preparations would also require an interface, both to create them, and a trait like interface where we can use multiple libraries' implementations (LibSpatialIndex.jl, SortTileRecursiveTree.jl, NearestNeighbors.jl, etc.)

The wrapper would also have to wrap feature collections, so we can use rtrees defined in things like .gpkg files or elsewhere for spatial joins etc.

Lazy `apply`

I have a few uses for lazy apply, so having that would be cool. For example: one could pass a lazy-apply'ed geometry to another apply function.

This is pretty low priority - just a thought!

Clipping Type Instabilities

I don't know if #69 takes care of this but, I am still seeing type instability, both from tuples, and also from the target.

To reproduce use:

import GeoInterface as GI
import GeometryOps as GO
using ProfileView
using Cthulhu

p25 = GI.Polygon([[(0.0, 0.0), (5.0, 0.0), (5.0, 8.0), (0.0, 8.0), (0.0, 0.0)],
    [(4.0, 0.5), (4.5, 0.5), (4.5, 3.5), (4.0, 3.5), (4.0, 0.5)],
    [(2.0, 4.0), (4.0, 4.0), (4.0, 6.0), (2.0, 6.0), (2.0, 4.0)]])
p26 = GI.Polygon([
    [(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (6.0, 5.0), (6.0, 3.0), (3.0, 3.0), (3.0, 1.0)],
    [(3.5, 5.5), (6.0, 5.5), (6.0, 6.5), (3.5, 6.5), (3.5, 5.5)],
    [(5.5, 1.5), (5.5, 2.5), (3.5, 2.5), (3.5, 1.5), (5.5, 1.5)]])

GO.intersection(p25, p26; target = GI.PolygonTrait)

f(n, p1, p2) = for _ in 1:n GO.intersection(p1, p2; target = GI.PolygonTrait) end
ProfileView.@profview f(10000, p25, p26)

The profile won't look just like the one above until #72 is merged, but the type unstable bar will still be there.

Add a `applyreduce` primative

For things like area it would be good to just reduce accross whole collections. Like apply but it doesnt reconstruct anything and has a second reducing op

Documentor

We have a lot of warnings from Documenter.jl when generating the documentation. I fixed some of them in #50, but not all of them.

Union between more than two polygons

How would this be implemented? This is a pretty big blocker to seamlessly merging many polygons into, for example, a multipolygon.

This would also enable unions between multipolygons and polygons.

LibGEOS does this pretty seamlessly, so we probably should as well :D

Predicates

Following on here #3

Using exact predicates will be something important to consider, making use of ExactPredicates.jl being the best choice. There is an issue of having to convert to Tuples, but this isn't an issue unless points are represented as (non-static) vectors - you already access individual coordinates anyway, so there is no difference in allocations. Also importantly, the performance loss of considering exact predicates is nothing major - for coordinates far from degeneracy, the fast path in ExactPredicates.jl makes it the same as using standard determinant definitions. For point sets where you end up needing to use the more complicated paths in ExactPredicates.jl, the performance does slow down - but that's exactly where you need it!

As an example, my original implementations in DelaunayTriangulation.jl just used determinant definitions for the predicates. You quickly run into infinite loops or e.g. triangles with the incorrect orientation. Some good examples demonstrating this importance here https://www.sciencedirect.com/science/article/pii/S0925772107000697 if you are curious. See also the first chapter here.

For the function is_point_on_segment I mentioned in the PR, where you need an exact predicate is in the computation of cross - this is just the definition of the orient predicate. It's really difficult to detect points on line segments without exact predicates, and when it matters it really matters that we detect collinear points in certain geometric algorithms (at least the ones I use anyway). You can see how I determine whether a point is on a line segment with ExactPredicates.jl here (which assumes already that the points are collinear, detected via point_position_relative_to_line (same as orient).

Lastly: You might be interested in considering three-valued logic rather than Bools, as argued e.g. here. I implement this with a Certificate system in DelaunayTriangulation.jl, but that's rather specific; ExactPredicates.jl uses (-1, 0, 1). It just turns out to be nice to be able to distinguish between cases like e.g. on, left, right. (Certificates were needed since some predicates even needed > 3 outcomes! like classifying the type of an intersection or specifying how a line intersects a triangle).

Anyway, don't have to take all of this on but thought you might appreciate some input. I had to think a lot about this at the time of initially implementing DelaunayTriangulation.jl.

Try to back most things by staticvectors instead of tuples

They have better dispatches on math functions and coordinate transformations.

What should GeometryOps do?

The driving idea behind this package is to unify all the disparate packages for geometric calculations in Julia, and make them GeoInterface.jl-compatible. We seem to be focusing primarily on 2/2.5D geometries for now.

Most of the usecases are driven by GIS and similar Earth data workflows, so this might be a bit specialized towards that, but methods should always be general to any coordinate space.

Methods we want to implement

For my part, I definitely want:

Signed area
Signed distance to a polygon/mesh
Centroid
Polylabel (shouldn't be too hard to do)

@rafaqz has contributed some really awesome code which allows people to generically go into the tree and act on GeoInterface compatible geometries.

@skygering - you seemed to be interested primarily in clipping and intersection, so those could be some methods we can add as well.

@DanielVandH has the excellent package DanielVandH/DelaunayTriangulation.jl which implements some similar stuff, we could use/yank things from there or just depend on it.

The package is built on JuliaGeometry/GeometryBasics.jl so we can use that as well.

Tip

Take a look at the meeting notes for more details!

Planned feature: barycentric coordinates on arbitrary polygons (without holes)

I'm in the process of writing a mesh rasterizer for CairoMakie which can rasterize quads on the CPU, and uses true quads instead of splitting into triangles. A part of this is obtaining barycentric coordinates for interpolation.

Barycentric coordinates are basically just a decomposition from Cartesian space into a (possibly degenerate) space of dimension N (where N is the number of vertices of the polygon), and the values are the relative nearness of the point to each of the vertices. Barycentric coordinates always sum to 1.

Creating this software rasterizer would allow us to really fiddle with meshes and do some pretty interesting stuff, like custom color mixing, working in LAB instead of RGB space, etc. It also allows us to smoothly (and correctly) interpolate values across an arbitrary polygon. See this example from CGAL:

There are some pretty useful resources here:

https://doc.cgal.org/latest/Barycentric_coordinates_2/index.html

Wachspress + mean value coords.pdf
How to compute mean value coordinates.pdf
Bloomenthal2022Barycentric.pdf

I'll probably have a decently working implementation of this up in a week or so.

`reproject` will not work with Proj transformations and threaded=true

A new Proj context must be created for each thread.

https://github.com/alex-s-gardner/FastGeoProjections.jl/blob/d19c373c89cc7ef2a7e0b14861e231d32e26cebf/src/epsg2epsg.jl#L16-L28 has code to handle this which we should adapt.

Implement what we can, wrap what we can't

We should add a method GEOS() to nearly every method here, which literally just calls LibGEOS.(...) so that we can expose a correct and working example of the feature. This would be huge for adoption, since the GO interface is a lot more Julian and easy to use (plus accepts more types) than LG!

`applyfoldl` and `applyfoldr` like `mapfoldl` and `mapfoldr`

I have to calculate transformed extents, and having these would be great so I don't materialize the transformed geometry!

The use would be something like:

applyfoldl(transformation, _update_extent, GI.PointTrait(), geometry; init = Extents.Extent(...))
```.

Wrap `Contour.jl` as a backend for `polygonize`

This would be neat to have as an extension package, and could even be used to tie into Rasters etc for handling if GeoInterface ever gets a raster API.

Contour.jl is used in the Makie contour recipe.

Geometry return types

I wanted to start a conversation about what return types of functions that create geometries should be. For example, simplify, if given a LibGEOS.Polygon as an input, still returns a GeoInterface.Polygon. This can be kinda irritating for users. This problem becomes even more pronounced with things like union, difference, and intersection, which my intern is working on a PR for (probably in the next week). Should they all return GI shapes? Should they return coordinates (and a trait type maybe)? Should they rebuild the same type as the input? What do we think is best?

Would love your opinions @rafaqz @asinghvi17.

import LibGEOS as LG
import GeoInterface as GI
import GeometryOps as GO
using BenchmarkTools
pLG = LG.Polygon([
    [
        [0.0, 5.0], [2.0, 2.0], [5.0, 2.0], [2.0, -2.0], [5.0, -5.0],
        [0.0, -2.0], [-5.0, -5.0], [-2.0, -2.0], [-5.0, 2.0], [-2.0, 2.0],
        [0.0, 5.0],
    ],
])

pGI = GI.Polygon([
    [
        (0.0, 5.0), (2.0, 2.0), (5.0, 2.0), (2.0, -2.0), (5.0, -5.0),
        (0.0, -2.0), (-5.0, -5.0), (-2.0, -2.0), (-5.0, 2.0), (-2.0, 2.0),
        (0.0, 5.0),
    ],
])

GO.equals(pLG, pGI)  # true

@benchmark LG.area(pLG)
@benchmark GO.area(pLG)
@benchmark GO.area(pGI)

I got the following:

So it seems like the method itself isn't bad, but it does NOT play nice with LibGEOS.

import everything from GeoInterface.jl and Extents.jl so users can just use `GO.foo`

Union of touching geometries is incorrect

MWE:

using GeometryBasics
import GeometryOps as GO, GeoInterface as GI
p1 = GI.Polygon([decompose(Point2f, Rect2f(0,0,1,1))]) |> GO.fix
p2 = GI.Polygon([decompose(Point2f, Rect2f(1,1,1,1))]) |> GO.fix

julia> GO.intersects(p1, p2)
true

julia> GO.touches(p1, p2)
true

julia> GO.union(p1, p2; target = GI.PolygonTrait)
2-element Vector{GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}}:
 GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(0.0, 0.0), (1.0, 0.0), (0.0, 1.0), (1.0, 1.0), (0.0, 0.0)], nothing, nothing)], nothing, nothing)
 GeoInterface.Wrappers.Polygon{false, false, Vector{GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}}, Nothing, Nothing}(GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}[GeoInterface.Wrappers.LinearRing{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(1.0, 1.0), (2.0, 1.0), (1.0, 2.0), (2.0, 2.0), (1.0, 1.0)], nothing, nothing)], nothing, nothing)

This should be a single polygon, that represents a rectangle of height 2 and width 1.

Optimizations

Functions that could be parallelized (please add as we go so that we can optimize later):

Signed area could conceivably be multithreaded

Use Aqua.jl

I just noticed we have ambiguities in crosses and centroid.

I will be a little annoying to add Aqua because GeometryBasics.jl also has ambiguties and Aqua catches them.

(We also have unbound args and stae deps so its worth adding anyway)