Repository: JuliaGeo/LibGEOS.jl
Pull request #115 (2022-06-11 16:07:58)

implement STRtree
Happy querying! Unfortunately I can't make it work without the $ in cfunction, but maybe we can live like that.


Patch:
diff --git a/Project.toml b/Project.toml
index b56d96a..7511f36 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,7 +14,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 [compat]
 julia = "1.6"
 CEnum = "0.2, 0.3, 0.4"
-GEOS_jll = "3.9"
+GEOS_jll = "3.10"
 GeoInterface = "0.4, 0.5"
 
 [targets]
diff --git a/src/LibGEOS.jl b/src/LibGEOS.jl
index 761aef4..6c11b5f 100644
--- a/src/LibGEOS.jl
+++ b/src/LibGEOS.jl
@@ -43,6 +43,7 @@ export Point,
     snap,
     delaunayTriangulation,
     delaunayTriangulationEdges,
+    constrainedDelaunayTriangulation,
     disjoint,
     touches,
     intersects,
@@ -89,7 +90,9 @@ export Point,
     getYMax,
     minimumRotatedRectangle,
     getGeometry,
-    getGeometries
+    getGeometries,
+    STRtree,
+    query
 
 include("libgeos_api.jl")
 
@@ -189,5 +192,6 @@ include("geos_functions.jl")
 include("geos_types.jl")
 include("geos_operations.jl")
 include("geo_interface.jl")
+include("strtree.jl")
 include("deprecated.jl")
 end
diff --git a/src/geos_functions.jl b/src/geos_functions.jl
index 1964246..7a04694 100644
--- a/src/geos_functions.jl
+++ b/src/geos_functions.jl
@@ -642,6 +642,17 @@ function delaunayTriangulation(
     result
 end
 
+function constrainedDelaunayTriangulation(
+    ptr::GEOSGeom,
+    context::GEOSContext = _context,
+)
+    result = GEOSConstrainedDelaunayTriangulation_r(context.ptr, ptr)
+    if result == C_NULL
+        error("LibGEOS: Error in GEOSConstrainedDelaunayTriangulation")
+    end
+    result
+end
+
 # -----
 # Binary predicates - return 2 on exception, 1 on true, 0 on false
 # -----
diff --git a/src/geos_operations.jl b/src/geos_operations.jl
index 54150ef..10b4d70 100644
--- a/src/geos_operations.jl
+++ b/src/geos_operations.jl
@@ -182,6 +182,8 @@ for geom in (
         MultiLineString(delaunayTriangulation(obj.ptr, tol, true))
     @eval delaunayTriangulation(obj::$geom, tol::Real = 0.0) =
         GeometryCollection(delaunayTriangulation(obj.ptr, tol, false))
+    @eval constrainedDelaunayTriangulation(obj::$geom) =
+        GeometryCollection(constrainedDelaunayTriangulation(obj.ptr))
 end
 
 sharedPaths(obj1::LineString, obj2::LineString) =
diff --git a/src/strtree.jl b/src/strtree.jl
new file mode 100644
index 0000000..48928b5
--- /dev/null
+++ b/src/strtree.jl
@@ -0,0 +1,61 @@
+struct STRtree{T}
+    ptr::Ptr{GEOSSTRtree}  # void pointer to the tree
+    items::T  # any geometry for which an envelope can be derived
+end
+
+"""
+    STRtree(items; nodecapacity=10, context::GEOSContext=_context)
+
+Create a Sort Tile Recursive tree for fast intersection queries of objects for which an
+envelope is defined. The tree cannot be changed after its creation.
+
+## Arguments
+- `items`: The items to store in the tree.
+- `nodecapacity`: The maximum number of items that can be stored per tree node (default 10).
+- `context`: The GEOS context in which the tree should be created in. Defaults to the global
+  LibGEOS context.
+"""
+function STRtree(items; nodecapacity = 10, context::GEOSContext = _context)
+    tree = LibGEOS.GEOSSTRtree_create_r(context.ptr, nodecapacity)
+    for item in items
+        envptr = envelope(item).ptr
+        GEOSSTRtree_insert_r(context.ptr, tree, envptr, pointer_from_objref(item))
+    end
+    return STRtree(tree, items)
+end
+
+"""
+    query(tree::STRtree, geometry; context::GEOSContext=_context)
+
+Returns the objects within `tree`, whose envelope intersects the envelope of `geometry`.
+
+## Arguments
+- `tree`: The STRtree to query
+- `geometry`: The LibGEOS geometry (e.g. Polygon) to run the query for
+- `context`: The GEOS context. Defaults to the global LibGEOS context.
+"""
+function query(
+    tree::STRtree{T},
+    geometry::Geometry;
+    context::GEOSContext = _context,
+) where {T}
+    matches = eltype(T)[]
+
+    # called for each overlapping item in the tree, used to push the item to matches
+    function callback(item, userdata)::Ptr{Cvoid}
+        item_ = unsafe_pointer_to_objref(item)
+        userdata_ = unsafe_pointer_to_objref(userdata)
+        push!(userdata_, item_)
+        return C_NULL
+    end
+
+    cfun_callback = @cfunction($callback, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}))
+    GEOSSTRtree_query_r(
+        context.ptr,
+        tree.ptr,
+        geometry.ptr,
+        cfun_callback,
+        pointer_from_objref(matches),
+    )
+    return matches
+end


Test patch:
diff --git a/test/runtests.jl b/test/runtests.jl
index 74628a9..4042820 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -18,4 +18,5 @@ end
     include("test_geo_interface.jl")
     include("test_regressions.jl")
     include("test_invalid_geometry.jl")
+    include("test_strtree.jl")
 end
diff --git a/test/test_geos_operations.jl b/test/test_geos_operations.jl
index 4b1a711..89848dd 100644
--- a/test/test_geos_operations.jl
+++ b/test/test_geos_operations.jl
@@ -70,6 +70,22 @@ end
     g2 = delaunayTriangulationEdges(g1, 2.0)
     equivalent_to_wkt(g2, "MULTILINESTRING ((0 0, 10 10), (0 0, 10 0), (10 0, 10 10))")
 
+    # LibGEOS.constrainedDelaunayTriangulationTest
+    g1 = readgeom("POLYGON EMPTY")
+    g2 = constrainedDelaunayTriangulation(g1)
+    @test isEmpty(g1)
+    @test isEmpty(g2)
+    @test GeoInterface.geotype(g2) == :GeometryCollection
+
+    g1 = readgeom("POINT(0 0)")
+    g2 = constrainedDelaunayTriangulation(g1)
+    @test isEmpty(g2)
+    @test GeoInterface.geotype(g2) == :GeometryCollection
+
+    g1 = readgeom("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))")
+    g2 = constrainedDelaunayTriangulation(g1)
+    equivalent_to_wkt(g2, "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))")
+
     # GEOSDistanceTest
     g1 = readgeom("POINT(10 10)")
     g2 = readgeom("POINT(3 6)")
diff --git a/test/test_strtree.jl b/test/test_strtree.jl
new file mode 100644
index 0000000..d8fc4e9
--- /dev/null
+++ b/test/test_strtree.jl
@@ -0,0 +1,28 @@
+@testset "STRtree" begin
+    # front page examples
+    p1 = readgeom("POLYGON((0 0,1 0,1 1,0 0))")
+    p2 = readgeom("POLYGON((0 0,1 0,1 1,0 1,0 0))")
+    p3 = readgeom("POLYGON((2 0,3 0,3 1,2 1,2 0))")
+
+    # query polygon
+    p4 = readgeom("POLYGON((0.5 0.5, 1.5 0.5, 1.5 1.5, 0.5 1.5, 0.5 0.5))")
+
+    tree = STRtree((p1, p2, p3))
+    res = query(tree, p4)
+    @test res isa Vector{Polygon}
+    @test res == [p1, p2]
+
+    # heterogeneous tree with linestring and polygon
+    ls = readgeom("LINESTRING(0 0,1 0,1 1,0 0)")
+    tree = STRtree((p1, p2, p3, ls))
+    res = query(tree, p4)
+    @test res isa Vector{GeoInterface.AbstractGeometry}
+    @test res == [p1, p2, ls]
+
+    # Let's accidentally "drop" p2. Since it is in the tree it won't be finalized.
+    p2 = nothing
+    p5 = readgeom("POLYGON((0 0,1 0,1 1,0 1,0 0))") # equal to p2
+    res = query(tree, p4)
+    @test equals(res[1], p1)
+    @test equals(res[2], p5)
+end