Implement a GEOS algorithm in an extension

Project.toml

@@ -1,23 +1,26 @@
 name = "GeometryOps"
 uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 authors = ["Anshul Singhvi <anshulsinghvi@gmail.com> and contributors"]
-version = "0.1.6"
+version = "0.1.7"
 
 [deps]
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
+LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 
 [extensions]
+GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsFlexiJoinsExt = "FlexiJoins"
 GeometryOpsProjExt = "Proj"
 
@@ -27,6 +30,7 @@ ExactPredicates = "2.2.8"
 FlexiJoins = "0.1.30"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7"
+LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
 SortTileRecursiveTree = "0.1"
@@ -40,17 +44,20 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
+NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "FlexiJoins", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "OffsetArrays", "Shapefile", "Test"]
+test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "SafeTestsets", "Shapefile", "Test"]

ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl

@@ -0,0 +1,24 @@
+module GeometryOpsLibGEOSExt
+
+import GeometryOps as GO, LibGEOS as LG
+import GeometryOps: GI
+
+import GeometryOps: GEOS, enforce
+
+using GeometryOps
+# The filter statement is required because in Julia, each module has its own versions of these
+# functions, which serve to evaluate or include code inside the scope of the module.
+# However, if you import those from another module (which you would with `all=true`),
+# that creates an ambiguity which causes a warning during precompile/load time.
+# In order to avoid this, we filter out these special functions.
+for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(GeometryOps; all = true))
+    @eval using GeometryOps: $name
+end
+
+include("buffer.jl")
+include("segmentize.jl")
+include("simplify.jl")
+
+include("simple_overrides.jl")
+
+end

ext/GeometryOpsLibGEOSExt/buffer.jl

@@ -0,0 +1,33 @@
+const _GEOS_CAPSTYLE_LOOKUP = Dict{Symbol, LG.GEOSBufCapStyles}(
+    :round => LG.GEOSBUF_CAP_ROUND,
+    :flat => LG.GEOSBUF_CAP_FLAT,
+    :square => LG.GEOSBUF_CAP_SQUARE,
+)
+
+const _GEOS_JOINSTYLE_LOOKUP = Dict{Symbol, LG.GEOSBufJoinStyles}(
+    :round => LG.GEOSBUF_JOIN_ROUND,
+    :mitre => LG.GEOSBUF_JOIN_MITRE,
+    :bevel => LG.GEOSBUF_JOIN_BEVEL,
+)
+
+to_cap_style(style::Symbol) = _GEOS_CAPSTYLE_LOOKUP[style]
+to_cap_style(style::LG.GEOSBufCapStyles) = style
+to_cap_style(num::Integer) = num
+
+to_join_style(style::Symbol) = _GEOS_JOINSTYLE_LOOKUP[style]
+to_join_style(style::LG.GEOSBufJoinStyles) = style
+to_join_style(num::Integer) = num
+
+function GO.buffer(alg::GEOS, geometry, distance)
+    # The reason we use apply here is so that this also works with featurecollections,
+    # tables, vectors of geometries, etc!
+    return apply(TraitTarget{GI.AbstractGeometryTrait}(), geometry) do geom
+        LG.bufferWithStyle(
+            GI.convert(LG, geom), distance; 
+            quadsegs = get(alg, :quadsegs, 8),
+            endCapStyle = to_cap_style(get(alg, :endCapStyle, :round)),
+            joinStyle = to_join_style(get(alg, :joinStyle, :round)),
+            mitreLimit = get(alg, :mitreLimit, 5.0),
+        )
+    end
+end

ext/GeometryOpsLibGEOSExt/segmentize.jl

@@ -0,0 +1,30 @@
+# # Segmentize
+import GeometryOps: segmentize, apply
+
+#=
+This file implements the LibGEOS segmentization method for GeometryOps.
+=#
+
+function _segmentize_geos(geom::LG.AbstractGeometry, max_distance)
+    context = LG.get_context(geom)
+    result = LG.GEOSDensify_r(context, geom, max_distance)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSDensify")
+    end
+    return LG.geomFromGEOS(result, context)
+end
+
+_segmentize_geos(geom, max_distance) = _segmentize_geos(GI.convert(LG, geom), max_distance)
+
+# 2 behaviours:
+# - enforce: enforce the presence of a kwargs
+# - fetch: fetch the value of a kwargs, or return a default value
+@inline function GO.segmentize(alg::GEOS, geom; threaded::Union{Bool, GO.BoolsAsTypes} = _False())
+    max_distance = enforce(alg, :max_distance, GO.segmentize)
+    return GO.apply(
+        Base.Fix2(_segmentize_geos, max_distance), 
+        GO.TraitTarget(GI.GeometryCollectionTrait(), GI.MultiPolygonTrait(), GI.PolygonTrait(), GI.MultiLineStringTrait(), GI.LineStringTrait(), GI.LinearRingTrait(), GI.MultiPointTrait(), GI.PointTrait()),
+        geom; 
+        threaded
+    )
+end

ext/GeometryOpsLibGEOSExt/simple_overrides.jl

@@ -0,0 +1,66 @@
+#=
+# Simple overrides
+This file contains simple overrides for GEOS, essentially only those
+functions which have direct counterparts in LG and only
+require conversion before calling.
+=#
+# ## Polygon set operations
+# ### Difference
+function GO.difference(::GEOS, geom_a, geom_b; target=nothing)
+    return LG.difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Union
+function GO.union(::GEOS, geom_a, geom_b; target=nothing)
+    return LG.union(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Intersection
+function GO.intersection(::GEOS, geom_a, geom_b; target=nothing)
+    return LG.intersection(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Symmetric difference
+function GO.symdifference(::GEOS, geom_a, geom_b; target=nothing)
+    return LG.symmetric_difference(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+
+# ## DE-9IM boolean methods
+# ### Equals
+function GO.equals(::GEOS, geom_a, geom_b)
+    return LG.equals(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Disjoint
+function GO.disjoint(::GEOS, geom_a, geom_b)
+    return LG.disjoint(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Touches
+function GO.touches(::GEOS, geom_a, geom_b)
+    return LG.touches(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Crosses
+function GO.crosses(::GEOS, geom_a, geom_b)
+    return LG.crosses(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Within
+function GO.within(::GEOS, geom_a, geom_b)
+    return LG.within(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Contains
+function GO.contains(::GEOS, geom_a, geom_b)
+    return LG.contains(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Overlaps
+function GO.overlaps(::GEOS, geom_a, geom_b)
+    return LG.overlaps(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Covers
+function GO.covers(::GEOS, geom_a, geom_b)
+    return LG.covers(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### CoveredBy
+function GO.coveredby(::GEOS, geom_a, geom_b)
+    return LG.coveredby(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+# ### Intersects
+function GO.intersects(::GEOS, geom_a, geom_b)
+    return LG.intersects(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+end
+

ext/GeometryOpsLibGEOSExt/simplify.jl

@@ -0,0 +1,31 @@
+# Address potential ambiguities
+GO._simplify(::GI.PointTrait, ::GO.GEOS, geom; kw...) = geom
+GO._simplify(::GI.MultiPointTrait, ::GO.GEOS, geom; kw...) = geom
+
+function GO._simplify(::GI.AbstractGeometryTrait, alg::GO.GEOS, geom; kwargs...)
+    method = get(alg, :method, :TopologyPreserve)
+    @assert haskey(alg.params, :tol) """
+        The `:tol` parameter is required for the GEOS algorithm in `simplify`, 
+        but it was not provided.  
+
+        Provide it by passing `GEOS(; tol = ...,) as the algorithm.
+        """
+    tol = alg.params.tol
+    if method == :TopologyPreserve
+        return LG.topologyPreserveSimplify(GI.convert(LG, geom), tol)
+    elseif method == :DouglasPeucker
+        return LG.simplify(GI.convert(LG, geom), tol)
+    else
+        error("Invalid method passed to `GO.simplify(GEOS(...), ...)`: $method. Please use :TopologyPreserve or :DouglasPeucker")
+    end
+end
+
+function GO._simplify(trait::GI.AbstractCurveTrait, alg::GO.GEOS, geom; kw...)
+    Base.invoke(
+        GO._simplify,
+        Tuple{GI.AbstractGeometryTrait, GO.GEOS, typeof(geom)},
+        trait, alg, geom; 
+        kw...
+    )
+end
+

src/GeometryOps.jl

@@ -18,12 +18,15 @@ const GB = GeometryBasics
 const TuplePoint{T} = Tuple{T, T} where T <: AbstractFloat
 const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T
 
+include("types.jl")
 include("primitives.jl")
 include("utils.jl")
+include("not_implemented_yet.jl")
 
 include("methods/angles.jl")
 include("methods/area.jl")
 include("methods/barycentric.jl")
+include("methods/buffer.jl")
 include("methods/centroid.jl")
 include("methods/distance.jl")
 include("methods/equals.jl")
@@ -70,6 +73,7 @@ function __init__()
     # Handle all available errors!
     Base.Experimental.register_error_hint(_reproject_error_hinter, MethodError)
     Base.Experimental.register_error_hint(_geodesic_segments_error_hinter, MethodError)
+    Base.Experimental.register_error_hint(_buffer_error_hinter, MethodError)
 end
 
 end

src/methods/buffer.jl

@@ -0,0 +1,27 @@
+#=
+# Buffer
+
+Buffering a geometry means computing the region `distance` away from it, and returning
+that region as the new geometry.
+
+As of now, we only support `GEOS` as the backend, meaning that LibGEOS must be loaded.
+=#
+
+function buffer(geometry, distance; kwargs...)
+    buffered = buffer(GEOS(; kwargs...), geometry, distance)
+    return tuples(buffered)
+end
+
+# Below is an error handler similar to the others we have for e.g. segmentize,
+# which checks if there is a method error for the geos backend.
+
+
+# Add an error hint for GeodesicSegments if Proj is not loaded!
+function _buffer_error_hinter(io, exc, argtypes, kwargs)
+    if isnothing(Base.get_extension(GeometryOps, :GeometryOpsLibGEOSExt)) && exc.f == buffer && first(exc.argtypes) == GEOS
+        print(io, "\n\nThe `buffer` method requires the LibGEOS.jl package to be explicitly loaded.\n")
+        print(io, "You can do this by simply typing ")
+        printstyled(io, "using LibGEOS"; color = :cyan, bold = true)
+        println(io, " in your REPL, \nor otherwise loading LibGEOS.jl via using or import.")
+    end
+end

src/not_implemented_yet.jl

@@ -0,0 +1,9 @@
+# # Not implemented yet
+# All of the functions in this file are not implemented in Julia yet.
+# Some of them may have implementations in LibGEOS which we can use
+# via an extension, but there is no native-Julia implementation for them.
+
+function symdifference end
+function buffer end
+function convexhull end
+function concavehull end

src/primitives.jl

@@ -53,52 +53,6 @@ TraitTarget
 
 =#
 
-
-#=
-We pass `threading` and `calc_extent` as types, not simple boolean values.  
-
-This is to help compilation - with a type to hold on to, it's easier for 
-the compiler to separate threaded and non-threaded code paths.
-
-Note that if we didn't include the parent abstract type, this would have been really 
-type unstable, since the compiler couldn't tell what would be returned!
-
-We had to add the type annotation on the `_booltype(::Bool)` method for this reason as well.
-=#
-abstract type BoolsAsTypes end
-struct _True <: BoolsAsTypes end
-struct _False <: BoolsAsTypes end
-
-@inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False()
-@inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
-
-"""
-    TraitTarget{T}
-
-This struct holds a trait parameter or a union of trait parameters.
-
-It is primarily used for dispatch into methods which select trait levels, 
-like `apply`, or as a parameter to `target`.
-
-## Constructors
-```julia
-TraitTarget(GI.PointTrait())
-TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like
-TraitTarget(TraitTarget(...))
-# There are also type based constructors available, but that's not advised.
-TraitTarget(GI.PointTrait)
-TraitTarget(Union{GI.LineStringTrait, GI.LinearRingTrait})
-# etc.
-```
-
-"""
-struct TraitTarget{T} end
-TraitTarget(::Type{T}) where T = TraitTarget{T}()
-TraitTarget(::T) where T<:GI.AbstractTrait = TraitTarget{T}()
-TraitTarget(::TraitTarget{T}) where T = TraitTarget{T}()
-TraitTarget(::Type{<:TraitTarget{T}}) where T = TraitTarget{T}()
-TraitTarget(traits::GI.AbstractTrait...) = TraitTarget{Union{map(typeof, traits)...}}()
-
 const THREADED_KEYWORD = "- `threaded`: `true` or `false`. Whether to use multithreading. Defaults to `false`."
 const CRS_KEYWORD = "- `crs`: The CRS to attach to geometries. Defaults to `nothing`."
 const CALC_EXTENT_KEYWORD = "- `calc_extent`: `true` or `false`. Whether to calculate the extent. Defaults to `false`."

src/transformations/segmentize.jl

@@ -174,7 +174,7 @@ end
 function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False())
     @assert method.max_distance > 0 "`max_distance` should be positive and nonzero!  Found $(method.max_distance)."
     segmentize_function = Base.Fix1(_segmentize, method)
-    return apply(segmentize_function, Union{GI.LinearRingTrait, GI.LineStringTrait}, geom; threaded)
+    return apply(segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded)
 end
 
 _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom))