Wrap concavehull

Project.toml

@@ -26,14 +26,16 @@ GeoInterface = "1.6"
 Makie = "0.23, 0.24"
 Plots = "1"
 RecipesBase = "1"
+TaskLocalValues = "0.1"
 Test = "<0.0.1,1"
 julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+TaskLocalValues = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Makie", "Plots", "RecipesBase", "Test"]
+test = ["Aqua", "Makie", "Plots", "RecipesBase", "TaskLocalValues", "Test"]

src/geo_interface.jl

@@ -330,6 +330,8 @@ for f in (
     :envelope,
     :minimumRotatedRectangle,
     :convexhull,
+    :concavehull,
+    :concavehull_of_polygons,
     :boundary,
     :unaryUnion,
     :pointOnSurface,

src/geos_functions.jl

@@ -677,6 +677,53 @@ function convexhull(obj::Geometry, context::GEOSContext = get_context(obj))
     geomFromGEOS(result, context)
 end
 
+"""
+    concavehull(geom, ratio; bylength = false, allow_holes = true)
+
+Compute the concave hull of a geometry, according to the chi-shape generated from 
+`geom` with the given `ratio`.
+
+If `bylength` is true, the chi-shape is generated from the length of the geometry, 
+calling `GEOSConcaveHullByLength`; otherwise it is generated from the area, calling 
+`GEOSConcaveHull`.
+
+If `allow_holes` is true, the concave hull may have holes.  If false, then holes are disallowed.
+
+See also [`convexhull`](@ref), [`concavehull_of_polygons`](@ref).
+"""
+function concavehull(obj::Geometry; ratio::Union{Real, Nothing} = nothing, length_ratio::Union{Real, Nothing} = nothing, allow_holes::Bool = true, context::GEOSContext = get_context(obj))
+    result = if isnothing(ratio) && isnothing(length_ratio)
+        throw(ArgumentError("Either `ratio` or `length_ratio` must be provided to `LibGEOS.concavehull`, none were provided."))
+    elseif !isnothing(ratio) && !isnothing(length_ratio)
+        throw(ArgumentError("Only one of `ratio` or `length_ratio` may be provided to `LibGEOS.concavehull`; got both."))
+    elseif !isnothing(ratio) && isnothing(length_ratio)
+        GEOSConcaveHull_r(context, obj, ratio, allow_holes)
+    elseif isnothing(ratio) && !isnothing(length_ratio)
+        GEOSConcaveHullByLength_r(context, obj, ratio, allow_holes)
+    end
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSConcaveHull")
+    end
+    geomFromGEOS(result, context)
+end
+
+"""
+    concavehull_of_polygons(obj::Geometry, ratio; tight = false, allow_holes = true)
+
+Compute the concave hull of a geometry, according to the chi-shape generated from 
+`obj` with the given `ratio`.
+
+If `tight` is true, the chi-shape is generated from the tightest fit of the geometry, 
+otherwise it is generated from the area.
+"""
+function concavehull_of_polygons(obj::Geometry, ratio::Real; tight::Bool = false, allow_holes::Bool = true, context::GEOSContext = get_context(obj))
+    result = GEOSConcaveHullOfPolygons_r(context, obj, ratio, tight, allow_holes)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSConcaveHullOfPolygons")
+    end
+    geomFromGEOS(result, context)
+end
+
 function difference(
     obj1::Geometry,
     obj2::Geometry,

test/test_geo_interface.jl

@@ -2,6 +2,7 @@ using Test, GeoInterface, LibGEOS, Extents
 import Makie, Plots
 const GI = GeoInterface
 const LG = LibGEOS
+import LibGEOS: Point
 
 @testset "Geo interface" begin
     pt = LibGEOS.Point(1.0, 2.0)
@@ -300,8 +301,8 @@ const LG = LibGEOS
     @test GeoInterface.geomtrait(geomcollection) == GeometryCollectionTrait()
     @test GeoInterface.is3d(geomcollection) == false
     @test GeoInterface.testgeometry(geomcollection)
-
     @testset "Conversion" begin
+        _concavehull(x) = LG.concavehull(x; ratio = 0.3)
         one_arg_functions = (
             LG.area,
             LG.geomLength,
@@ -317,6 +318,7 @@ const LG = LibGEOS
             LG.delaunayTriangulationEdges,
             LG.delaunayTriangulation,
             LG.constrainedDelaunayTriangulation,
+            _concavehull,
             # these have different signatures
             # LG.simplify, LG.topologyPreserveSimplify,
         )
@@ -357,39 +359,56 @@ const LG = LibGEOS
         @test geom isa MultiPoint
         @test GeoInterface.coordinates(geom) == coords
         for f in one_arg_functions
-            @test f(LibGEOS.MultiPoint(coords)) == f(GeoInterface.MultiPoint(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.MultiPoint(coords)) == f(GeoInterface.MultiPoint(coords))
+            end
         end
         coords2 = [[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]
         for f in two_arg_functions
-            @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) ==
-                  f(GeoInterface.LineString(coords), GeoInterface.LineString(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) ==
+                      f(GeoInterface.LineString(coords), GeoInterface.LineString(coords))
+            end
         end
 
         coords = [[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]]
         geom = GeoInterface.convert(LineString, GeoInterface.LineString(coords))
         @test geom isa LineString
         @test GeoInterface.coordinates(geom) == coords
         for f in one_arg_functions
-            @test f(LibGEOS.LineString(coords)) == f(GeoInterface.LineString(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.LineString(coords)) == f(GeoInterface.LineString(coords))
+            end
         end
         coords2 = [[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]
         for f in two_arg_functions
-            @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) ==
-                  f(GeoInterface.LineString(coords), GeoInterface.LineString(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.LineString(coords), LibGEOS.LineString(coords)) ==
+                      f(GeoInterface.LineString(coords), GeoInterface.LineString(coords))
+            end
         end
 
         coords = [[[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]]]
         geom = GeoInterface.convert(MultiLineString, GeoInterface.MultiLineString(coords))
         @test geom isa MultiLineString
         @test GeoInterface.coordinates(geom) == coords
         for f in one_arg_functions
-            @test f(LibGEOS.MultiLineString(coords)) ==
-                  f(GeoInterface.MultiLineString(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.MultiLineString(coords)) ==
+                      f(GeoInterface.MultiLineString(coords))
+            end
         end
         coords2 = [[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]]
         for f in two_arg_functions
-            @test f(LibGEOS.MultiLineString(coords), LibGEOS.MultiLineString(coords2)) ==
-                  f(GeoInterface.MultiLineString(coords), LibGEOS.MultiLineString(coords2))
+            @testset let current_function = f
+                @test f(
+                    LibGEOS.MultiLineString(coords),
+                    LibGEOS.MultiLineString(coords2),
+                ) == f(
+                    GeoInterface.MultiLineString(coords),
+                    LibGEOS.MultiLineString(coords2),
+                )
+            end
         end
 
         coords = [[[0.0, 0], [0.0, 10], [10.0, 10], [10.0, 0], [0.0, 0]]]
@@ -400,12 +419,16 @@ const LG = LibGEOS
         @test GeoInterface.nhole(geom) == 0
         @test GeoInterface.coordinates(geom) == coords
         for f in one_arg_functions
-            @test f(LibGEOS.Polygon(coords)) == f(GeoInterface.Polygon(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.Polygon(coords)) == f(GeoInterface.Polygon(coords))
+            end
         end
         coords2 = [[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]]
         for f in two_arg_functions
-            @test f(LibGEOS.Polygon(coords), LibGEOS.Polygon(coords2)) ==
-                  f(GeoInterface.Polygon(coords), LibGEOS.Polygon(coords2))
+            @testset let current_function = f
+                @test f(LibGEOS.Polygon(coords), LibGEOS.Polygon(coords2)) ==
+                      f(GeoInterface.Polygon(coords), LibGEOS.Polygon(coords2))
+            end
         end
 
         pgeom = LibGEOS.prepareGeom(geom)
@@ -417,12 +440,17 @@ const LG = LibGEOS
         @test geom isa MultiPolygon
         @test GeoInterface.coordinates(geom) == coords
         for f in one_arg_functions
-            @test f(LibGEOS.MultiPolygon(coords)) == f(GeoInterface.MultiPolygon(coords))
+            @testset let current_function = f
+                @test f(LibGEOS.MultiPolygon(coords)) ==
+                      f(GeoInterface.MultiPolygon(coords))
+            end
         end
         coords2 = [[[[0.0, 10], [0.5, 10], [20.0, 20], [10.0, 10], [0.0, 10]]]]
         for f in two_arg_functions
-            @test f(LibGEOS.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) ==
-                  f(GeoInterface.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2))
+            @testset let current_function = f
+                @test f(LibGEOS.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2)) ==
+                      f(GeoInterface.MultiPolygon(coords), LibGEOS.MultiPolygon(coords2))
+            end
         end
     end
 

test/test_geos_functions.jl

@@ -1,6 +1,7 @@
 using Test
 using LibGEOS
 import GeoInterface
+import GeoInterface as GI
 using Extents
 
 @testset "WKTWriter" begin
@@ -956,7 +957,21 @@ end
     @test LibGEOS.reverse(readgeom("LINESTRING(0 0, 1 1)")) ==
           readgeom("LINESTRING(1 1, 0 0)")
 
+    # NOTE: maximum inscribed circle is not an exact algorithm,
+    # and may be susceptible to both floating point error
+    # and changes to internal hashing order of the quadtree.
+    # So we need to test for what we really care about here:
+    # that the center is in the right place, and that the radius is correct.
+    # For example, between GEOS v3.13 and v3.14, the test broke
+    # because it had been checking that the returned linestring was precisely equal to a reference,
+    # but the side switched from the left to the right.  Still a completely correct
+    # and valid result - but not the same as the reference and it settled in a different cell.
+    # The algorithm LibGEOS and jts use here is the same one that Mapbox uses
+    # as well as Polylabel.jl in Julia (soon to be integrated into GeometryOps).
     geo = readgeom("POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))")
     mic = LibGEOS.maximumInscribedCircle(geo, 1e-4)
-    @test mic == readgeom("LINESTRING (0.5000000000000001 0.5, 1 0.5)")
+    @test GI.coordinates(GI.getpoint(mic, 1)) ≈ [0.5, 0.5] atol = 1e-5
+    x1, y1 = GI.coordinates(GI.getpoint(mic, 1))
+    x2, y2 = GI.coordinates(GI.getpoint(mic, 2))
+    @test hypot(x1 - x2, y1 - y2) ≈ 0.5 atol = 1e-5
 end

test/test_geos_types.jl

@@ -1,3 +1,4 @@
+import TaskLocalValues
 
 @testset "open_issue_if_conversion_makes_sense" begin
     polygon = readgeom("POLYGON EMPTY")
@@ -390,10 +391,10 @@ end
 @testset "Multi threading" begin
     function f91(n)
         # adapted from https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1267732709
-        contexts = [LibGEOS.GEOSContext() for i = 1:Threads.nthreads()]
+        contexts = TaskLocalValues.TaskLocalValue{LibGEOS.GEOSContext}(() -> LibGEOS.GEOSContext())
         p = [[[-1.0, -1], [+1, -1], [+1, +1], [-1, +1], [-1, -1]]]
-        Threads.@threads :static for i = 1:n
-            ctx = contexts[Threads.threadid()-Threads.nthreads(:interactive)]
+        Threads.@threads for i = 1:n
+            ctx = contexts[]
             g1 = LibGEOS.Polygon(p, ctx)
             g2 = LibGEOS.Polygon(p, ctx)
             for j = 1:n
@@ -407,10 +408,10 @@ end
     @testset "clone" begin
         function f(n)
             # adapted from https://github.com/JuliaGeo/LibGEOS.jl/issues/91#issuecomment-1267732709
-            contexts = [LibGEOS.GEOSContext() for i = 1:Threads.nthreads()]
+            contexts = TaskLocalValues.TaskLocalValue{LibGEOS.GEOSContext}(() -> LibGEOS.GEOSContext())
             p = LibGEOS.Polygon([[[-1.0, -1], [+1, -1], [+1, +1], [-1, +1], [-1, -1]]])
-            Threads.@threads :static for i = 1:n
-                ctx = contexts[Threads.threadid()-Threads.nthreads(:interactive)]
+            Threads.@threads for i = 1:n
+                ctx = contexts[]
                 g1 = LibGEOS.clone(p, ctx)
                 g2 = LibGEOS.clone(p, ctx)
                 for j = 1:n