Merge branch 'main' into dd

.github/workflows/CI.yml

@@ -35,6 +35,15 @@ jobs:
       - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
+      - uses: julia-actions/julia-processcoverage@v1
+        with:
+          directories: src,ext
+      - name: Upload coverage reports to Codecov
+        uses: codecov/[email protected]
+        with:
+          files: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          slug: JuliaGeo/GeometryOps.jl
   docs:
     name: Documentation
     runs-on: ubuntu-latest
@@ -55,6 +64,7 @@ jobs:
           install-package: false
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
       - run: |
           julia --project=docs -e '
             using Documenter: DocMeta, doctest

Project.toml

@@ -1,10 +1,11 @@
 name = "GeometryOps"
 uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 authors = ["Anshul Singhvi <[email protected]> and contributors"]
-version = "0.1.10"
+version = "0.1.11"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -28,6 +29,7 @@ GeometryOpsProjExt = "Proj"
 
 [compat]
 CoordinateTransformations = "0.5, 0.6"
+DelaunayTriangulation = "1.0.4"
 ExactPredicates = "2.2.8"
 DimensionalData = "0.27"
 FlexiJoins = "0.1.30"

README.md

@@ -1,7 +1,7 @@
 <img width="400" alt="GeometryOps.jl" src="https://github.com/JuliaGeo/GeometryOps.jl/assets/32143268/92c5526d-23a9-4e01-aee0-2fcea99c5001">
 
 ![Lifecycle:Experimental](https://img.shields.io/badge/Lifecycle-Experimental-339999)
-[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/v0.1.10/)
+[![Stable](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/stable/)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaGeo.github.io/GeometryOps.jl/dev/)
 [![Build Status](https://github.com/JuliaGeo/GeometryOps.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaGeo/GeometryOps.jl/actions/workflows/CI.yml?query=branch%3Amain)
 

benchmarks/benchmarks.jl

@@ -11,7 +11,7 @@ import GeometryOps as GO,
     GeoInterface as GI, 
     GeometryBasics as GB, 
     LibGEOS as LG
-import GeoJSON
+import GeoJSON, NaturalEarth
 # In order to benchmark, we'll actually use the new [Chairmarks.jl](https://github.com/lilithhafner/Chairmarks.jl), 
 # since it's significantly faster than BenchmarkTools.  To keep benchmarks organized, though, we'll still use BenchmarkTools' 
 # `BenchmarkGroup` structure.
@@ -184,7 +184,7 @@ usa_difference_suite = usa_poly_suite["difference"] = BenchmarkGroup(["title:USA
 usa_intersection_suite = usa_poly_suite["intersection"] = BenchmarkGroup(["title:USA intersection", "subtitle:Tested on CONUS"])
 usa_union_suite = usa_poly_suite["union"] = BenchmarkGroup(["title:USA union", "subtitle:Tested on CONUS"])
 
-fc = GeoJSON.read(read(download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/ca96624a56bd078437bca8184e78163e5039ad19/geojson/ne_10m_admin_0_countries.geojson")))
+fc = NaturalEarth.naturalearth("admin_0_countries", 10)
 usa_multipoly = fc.geometry[findfirst(==("United States of America"), fc.NAME)] |> x -> GI.convert(LG, x) |> LG.makeValid |> GO.tuples
 
 usa_poly = GI.getgeom(usa_multipoly, findmax(GO.area.(GI.getgeom(usa_multipoly)))[2]) # isolate the poly with the most area
@@ -194,7 +194,7 @@ f, a, p = plot(usa_poly; label = "Original", axis = (; aspect = DataAspect()));
 axislegend(a)
 f
 
-# Now, we  get to benchmarking:
+# Now, we get to benchmarking:
 
 
 usa_o_lg, usa_o_go = lg_and_go(usa_poly)

docs/src/tutorials/creating_geometry.md

@@ -186,7 +186,7 @@ Read the land `MultiPolygon`s as a `GeoJSON.FeatureCollection`.
 land_geo = GeoJSON.read(land_path)
 ````
 
-We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destinaton` CRS. For GeoMakie, `source` is the CRS of the input and `dest` is the CRS you want to visualize in.
+We then need to create a figure with a `GeoAxis` that can handle the projection between `source` and `destination` CRS. For GeoMakie, `source` is the CRS of the input and `dest` is the CRS you want to visualize in.
 
 ````@example creating_geometry
 fig = Figure(size=(1000, 500));
@@ -267,7 +267,7 @@ y = r .* (k + 1) .* sin.(ϴ) .- r .* sin.((k + 1) .* ϴ);
 ring4 = GI.LinearRing(Point.(zip(x, y)))
 ````
 
-But this time when we create the `Polygon` we beed to specify the `CRS` at the time of creation, making it a geospatial polygon
+But this time when we create the `Polygon` we need to specify the `CRS` at the time of creation, making it a geospatial polygon
 
 ````@example creating_geometry
 geopoly1 = GI.Polygon([ring4], crs = source_crs1)

docs/src/tutorials/geodesic_paths.md

@@ -0,0 +1,17 @@
+# Geodesic paths
+
+Geodesic paths are paths computed on an ellipsoid, as opposed to a plane.  
+
+```@example geodesic
+import GeometryOps as GO, GeoInterface as GI
+using CairoMakie, GeoMakie
+
+
+IAH = (-95.358421, 29.749907)
+AMS = (4.897070, 52.377956)
+
+
+fig, ga, _cp = lines(GeoMakie.coastlines(); axis = (; type = GeoAxis))
+lines!(ga, GO.segmentize(GO.GeodesicSegments(; max_distance = 100_000), GI.LineString([IAH, AMS])); color = Makie.wong_colors()[2])
+fig
+```

ext/GeometryOpsLibGEOSExt/simple_overrides.jl

@@ -64,3 +64,17 @@ function GO.intersects(::GEOS, geom_a, geom_b)
     return LG.intersects(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
 end
 
+# ## Convex hull
+function GO.convex_hull(::GEOS, geoms)
+    return LG.convexhull(
+        LG.MultiPoint(
+            collect(
+                GO.flatten(
+                    x -> GI.convert(LG, x), 
+                    GI.PointTrait, 
+                    geoms
+                )
+            )
+        )
+    )
+end

ext/GeometryOpsProjExt/reproject.jl

@@ -14,7 +14,7 @@ function reproject(geom;
         end
 
         # If its still nothing, error
-        isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword"))
+        isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword"))
 
         # Otherwise reproject
         reproject(geom, source_crs, target_crs; kw...)

src/GeometryOps.jl

@@ -4,13 +4,14 @@ module GeometryOps
 
 using GeoInterface
 using GeometryBasics
-import Tables
 using LinearAlgebra, Statistics
+
+import Tables
 import GeometryBasics.StaticArrays
+import DelaunayTriangulation # for convex hull and triangulation
 import ExactPredicates
 import Base.@kwdef
-
-using GeoInterface.Extents: Extents
+import GeoInterface.Extents: Extents
 
 const GI = GeoInterface
 const GB = GeometryBasics
@@ -28,6 +29,7 @@ include("methods/area.jl")
 include("methods/barycentric.jl")
 include("methods/buffer.jl")
 include("methods/centroid.jl")
+include("methods/convex_hull.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
 include("methods/clipping/predicates.jl")

src/methods/angles.jl

@@ -39,7 +39,7 @@ This is computed differently for different geometries:
     - The angles of a point is an empty vector.
     - The angles of a single line segment is an empty vector.
     - The angles of a linestring or linearring is a vector of angles formed by the curve.
-    - The angles of a polygin is a vector of vectors of angles formed by each ring.
+    - The angles of a polygon is a vector of vectors of angles formed by each ring.
     - The angles of a multi-geometry collection is a vector of the angles of each of the
         sub-geometries as defined above.
 
@@ -83,7 +83,7 @@ function _angles(::Type{T}, ::GI.LinearRingTrait, geom; interior = true) where T
 end
 
 #= The angles of a polygon is a vector of polygon angles. Note that if there are holes
-within the polyogn, the angles will be listed after the exterior ring angles in order of the
+within the polygon, the angles will be listed after the exterior ring angles in order of the
 holes. All angles, including the hole angles, are interior angles of the polygon.=#
 function _angles(::Type{T}, ::GI.PolygonTrait, geom) where T
     angles = _angles(T, GI.LinearRingTrait(), GI.getexterior(geom); interior = true)

src/methods/area.jl

@@ -29,7 +29,7 @@ lines!(
 f
 ```
 The points are ordered in a counterclockwise fashion, which means that the signed area
-is negative.  If we reverse the order of the points, we get a postive area.
+is negative.  If we reverse the order of the points, we get a positive area.
 ```@example rect
 GO.signed_area(rect)  # -1.0
 ```
@@ -76,7 +76,7 @@ end
     signed_area(geom, [T = Float64])::T
 
 Returns the signed area of a single geometry, based on winding order. 
-This is computed slighly differently for different geometries:
+This is computed slightly differently for different geometries:
 
     - The signed area of a point is always zero.
     - The signed area of a curve is always zero.

src/methods/centroid.jl

@@ -45,7 +45,7 @@ LineString or LinearRing when they are closed, for example as the interior hole
 of a polygon.
 
 The helper functions centroid_and_length and centroid_and_area are made
-availible just in case the user also needs the area or length to decrease
+available just in case the user also needs the area or length to decrease
 repeat computation.
 =#
 """

src/methods/clipping/clipping_processor.jl

@@ -74,10 +74,10 @@ poly_a, including its intersection points with poly_b. The information stored in
 PolyNode is needed for clipping using the Greiner-Hormann clipping algorithm.
     
 Note: After calling this function, a_list is not fully formed because the neighboring
-indicies of the intersection points in b_list still need to be updated. Also we still have
+indices of the intersection points in b_list still need to be updated. Also we still have
 not update the entry and exit flags for a_list.
     
-The a_idx_list is a list of the indicies of intersection points in a_list. The value at
+The a_idx_list is a list of the indices of intersection points in a_list. The value at
 index i of a_idx_list is the location in a_list where the ith intersection point lies.
 =#
 function _build_a_list(::Type{T}, poly_a, poly_b; exact) where T
@@ -179,7 +179,7 @@ creates a vector of PolyNodes to represent poly_b. The information stored in eac
 is needed for clipping using the Greiner-Hormann clipping algorithm.
     
 Note: after calling this function, b_list is not fully updated. The entry/exit flags still
-need to be updated. However, the neightbor value in a_list is now updated.
+need to be updated. However, the neighbor value in a_list is now updated.
 =#
 function _build_b_list(::Type{T}, a_idx_list, a_list, n_b_intrs, poly_b) where T
     # Sort intersection points by insertion order in b_list
@@ -609,8 +609,8 @@ _get_poly_type(::Type{T}) where T =
 #=
     _find_non_cross_orientation(a_list, b_list, a_poly, b_poly; exact)
 
-For polygns with no crossing intersection points, either one polygon is inside of another,
-or they are seperate polygons with no intersection (other than an edge or point).
+For polygons with no crossing intersection points, either one polygon is inside of another,
+or they are separate polygons with no intersection (other than an edge or point).
 
 Return two booleans that represent if a is inside b (potentially with shared edges / points)
 and visa versa if b is inside of a.
@@ -643,6 +643,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
     for i in 1:n_polys
         n_new_per_poly = 0
         for curr_hole in Iterators.map(tuples, hole_iterator) # loop through all holes
+            curr_hole = _linearring(curr_hole)
             # loop through all pieces of original polygon (new pieces added to end of list)
             for j in Iterators.flatten((i:i, (n_polys + 1):(n_polys + n_new_per_poly)))
                 curr_poly = return_polys[j]
@@ -655,7 +656,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                         append!(remove_poly_idx, falses(n_new_pieces))
                         n_new_per_poly += n_new_pieces
                     end
-                    if !on_ext && !out_ext  # hole is completly within exterior
+                    if !on_ext && !out_ext  # hole is completely within exterior
                         push!(curr_poly.geom, new_hole)
                     else  # hole is partially within and outside of polygon's exterior
                         new_polys = difference(curr_poly_ext, new_hole_poly, T; target=GI.PolygonTrait())
@@ -669,7 +670,7 @@ function _add_holes_to_polys!(::Type{T}, return_polys, hole_iterator, remove_pol
                             n_new_per_poly += n_new_polys
                         end
                     end
-                # polygon is completly within hole
+                # polygon is completely within hole
                 elseif coveredby(curr_poly_ext, GI.Polygon(StaticArrays.SVector(curr_hole)))
                     remove_poly_idx[j] = true
                 end
@@ -687,7 +688,7 @@ end
 
 The new hole is combined with any existing holes in curr_poly. The holes can be combined
 into a larger hole if they are intersecting. If this happens, then the new, combined hole is
-returned with the orignal holes making up the new hole removed from curr_poly. Additionally,
+returned with the original holes making up the new hole removed from curr_poly. Additionally,
 if the combined holes form a ring, the interior is added to the return_polys as a new
 polygon piece. Additionally, holes leftover after combination will be checked for it they
 are in the "main" polygon or in one of these new pieces and moved accordingly. 
@@ -750,7 +751,7 @@ function _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
                     continue
                 else
                     p3 = p
-                    # check if p2 is approximatly on the edge formed by p1 and p3 - remove if so
+                    # check if p2 is approximately on the edge formed by p1 and p3 - remove if so
                     if Predicates.orient(p1, p2, p3; exact = _False()) == 0
                         remove_idx[i - 1] = true
                     end

src/methods/clipping/coverage.jl

@@ -4,7 +4,7 @@ export coverage
 ## What is coverage?
 
 Coverage is the amount of geometry area within a bounding box defined by the minimum and
-maximum x and y-coordiantes of that bounding box, or an Extent containing that information.
+maximum x and y-coordinates of that bounding box, or an Extent containing that information.
 
 To provide an example, consider this rectangle:
 ```@example rect

src/methods/clipping/cut.jl

@@ -98,7 +98,7 @@ of cut geometry in Vector{Vector{Tuple}} format.
 
 Note: degenerate cases where intersection points are vertices do not work right now. =#
 function _cut(::Type{T}, geom, line, geom_list, intr_list, n_intr_pts; exact) where T
-    # Sort and catagorize the intersection points
+    # Sort and categorize the intersection points
     sort!(intr_list, by = x -> geom_list[x].fracs[2])
     _flag_ent_exit!(GI.LineTrait(), line, geom_list; exact)
     # Add first point to output list

src/methods/clipping/difference.jl

@@ -83,7 +83,7 @@ function _difference(
             end
         end
     end
-    # Remove uneeded collinear points on same edge
+    # Remove unneeded collinear points on same edge
     _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
     return polys
 end
@@ -147,7 +147,7 @@ end
 #= Multipolygon with multipolygon difference - note that all intersection regions between
 sub-polygons of `multipoly_a` and sub-polygons of `multipoly_b` will be removed from the
 corresponding sub-polygon of `multipoly_a`. Unless specified with `fix_multipoly = nothing`,
-`multipolygon_a` will be validated using the given (defauly is `UnionIntersectingPolygons()`)
+`multipolygon_a` will be validated using the given (default is `UnionIntersectingPolygons()`)
 correction. =#
 function _difference(
     target::TraitTarget{GI.PolygonTrait}, ::Type{T},
@@ -169,7 +169,7 @@ function _difference(
         else
             difference(GI.MultiPolygon(polys), poly_b; target, fix_multipoly)
         end
-        #= One multipoly_a has been completly covered (and thus removed) there is no need to
+        #= One multipoly_a has been completely covered (and thus removed) there is no need to
         continue taking the difference =#
         isempty(polys) && break
     end

src/methods/clipping/intersection.jl

@@ -5,7 +5,7 @@ export intersection, intersection_points
     Enum LineOrientation
 Enum for the orientation of a line with respect to a curve. A line can be
 `line_cross` (crossing over the curve), `line_hinge` (crossing the endpoint of the curve),
-`line_over` (colinear with the curve), or `line_out` (not interacting with the curve).
+`line_over` (collinear with the curve), or `line_out` (not interacting with the curve).
 """
 @enum LineOrientation line_cross=1 line_hinge=2 line_over=3 line_out=4
 
@@ -86,7 +86,7 @@ function _intersection(
         hole_iterator = Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b)))
         _add_holes_to_polys!(T, polys, hole_iterator, remove_idx; exact)
     end
-    # Remove uneeded collinear points on same edge
+    # Remove unneeded collinear points on same edge
     _remove_collinear_points!(polys, remove_idx, poly_a, poly_b)
     return polys
 end
@@ -128,7 +128,7 @@ function _intersection(
 end
 
 #= Multipolygon with polygon intersection is equivalent to taking the intersection of the
-poylgon with the multipolygon and thus simply switches the order of operations and calls the
+polygon with the multipolygon and thus simply switches the order of operations and calls the
 above method. =#
 _intersection(
     target::TraitTarget{GI.PolygonTrait}, ::Type{T},
@@ -197,7 +197,7 @@ intersection_points(geom_a, geom_b, ::Type{T} = Float64) where T <: AbstractFloa
     _intersection_points(T, GI.trait(geom_a), geom_a, GI.trait(geom_b), geom_b)
 
 
-#= Calculates the list of intersection points between two geometries, inlcuding line
+#= Calculates the list of intersection points between two geometries, including line
 segments, line strings, linear rings, polygons, and multipolygons. =#
 function _intersection_points(::Type{T}, ::GI.AbstractTrait, a, ::GI.AbstractTrait, b; exact = _True()) where T
     # Initialize an empty list of points
@@ -384,7 +384,7 @@ end
 
 #= If lines defined by (a1, a2) and (b1, b2) meet at one point that is not an endpoint of
 either segment, they form a crossing intersection with a singular intersection point. That 
-point is caculated by finding the fractional distance along each segment the point occurs
+point is calculated by finding the fractional distance along each segment the point occurs
 at (α, β). If the point is too close to an endpoint to be distinct, the point shares a value
 with the endpoint, but with a non-zero and non-one fractional value. If the intersection
 point calculated is outside of the envelope of the two segments due to floating point error,
@@ -407,9 +407,9 @@ function _find_cross_intersection(::Type{T}, a1, a2, b1, b2, a_ext, b_ext) where
     β = _clamped_frac(Δbax * Δay - Δbay * Δax, a_cross_b, eps(T))
 
     #= Intersection will be where a1 + α * Δa = b1 + β * Δb. However, due to floating point
-    innacurracies, α and β calculations may yeild different intersection points. Average
+    inaccuracies, α and β calculations may yield different intersection points. Average
     both points together to minimize difference from real value, as long as segment isn't 
-    vertical or horizontal as this will almost certianly lead to the point being outside the
+    vertical or horizontal as this will almost certainly lead to the point being outside the
     envelope due to floating point error. Also note that floating point limitations could
     make intersection be endpoint if α≈0 or α≈1.=#
     x = if Δax == 0