[Experimental] Extent forwarding for predicates (#301)

Author: asinghvi17

ext/GeometryOpsLibGEOSExt/simple_overrides.jl

@@ -23,47 +23,27 @@ function GO.symdifference(::GEOS, geom_a, geom_b; target=nothing, calc_extent =
 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))
+# These are all the same so we loop over all names and eval them in
+for fn in (:equals, :disjoint, :touches, :crosses, :within, :contains, :overlaps, :covers, :coveredby, :intersects)
+    @eval begin
+        # The basic method for geometries
+        function GO.$fn(::GEOS, geom_a, geom_b)
+            return LG.$fn(GI.convert(LG, geom_a), GI.convert(LG, geom_b))
+        end
+        # Extents and geometries
+        function GO.$fn(alg::GEOS, geom_a::GO.Extents.Extent, geom_b)
+            return GO.$fn(alg, GO.extent_to_polygon(geom_a), geom_b)
+        end
+        function GO.$fn(alg::GEOS, geom_a, geom_b::GO.Extents.Extent)
+            return GO.$fn(alg, geom_a, GO.extent_to_polygon(geom_b))
+        end
+        # Pure extents - this should probably be some GEOSRect or something,
+        # but for now this works
+        function GO.$fn(alg::GEOS, geom_a::GO.Extents.Extent, geom_b::GO.Extents.Extent)
+            return GO.$fn(alg, GO.extent_to_polygon(geom_a), GO.extent_to_polygon(geom_b))
+        end
+    end
 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
-
 # ## Convex hull
 function GO.convex_hull(::GEOS, geoms)
     chull = LG.convexhull(

src/methods/geom_relations/coveredby.jl

@@ -271,3 +271,16 @@ function _coveredby(
     end
     return true
 end
+
+# Extent forwarding
+function _coveredby(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return _coveredby(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function _coveredby(t1, e1::Extents.Extent, t2, g2)
+    return _coveredby(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function _coveredby(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return Extents.coveredby(e1, e2)
+end
+
+

src/methods/geom_relations/crosses.jl

@@ -26,6 +26,25 @@ crosses(::LineStringTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_line(
 crosses(::PolygonTrait, g1, ::MultiPointTrait, g2)::Bool = multipoint_crosses_poly(g2, g1)
 crosses(::PolygonTrait, g1, ::LineStringTrait, g2)::Bool = line_crosses_poly(g2, g1)
 
+
+
+function crosses(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return crosses(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function crosses(t1, e1::Extents.Extent, t2, g2)
+    return crosses(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function crosses(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return false # two extents can never cross!
+    # TODO that's not quite true, if one extent is a line...
+end
+
+
+
+
+
+
+
 function multipoint_crosses_line(geom1, geom2)
     int_point = false
     ext_point = false
@@ -128,3 +147,5 @@ function _point_on_segment(point, (start, stop); exclude_boundary::Symbol=:none)
     end
     return false
 end
+
+

src/methods/geom_relations/disjoint.jl

@@ -255,3 +255,17 @@ function _disjoint(
     end
     return true
 end
+
+
+# Extent forwarding
+function _disjoint(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return _disjoint(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function _disjoint(t1, e1::Extents.Extent, t2, g2)
+    return _disjoint(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function _disjoint(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return Extents.disjoint(e1, e2)
+end
+
+

src/methods/geom_relations/overlaps.jl

@@ -227,6 +227,18 @@ function _overlaps(
     return seg_val == line_over && (!a_fully_within && !b_fully_within)
 end
 
+# Extent forwarding
+
+function _overlaps(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return _overlaps(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function _overlaps(t1, e1::Extents.Extent, t2, g2)
+    return _overlaps(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function _overlaps(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return Extents.overlaps(e1, e2)
+end
+
 #= TODO: Once overlaps is swapped over to use the geom relations workflow, can
 delete these helpers. =#
 

src/methods/geom_relations/touches.jl

@@ -293,4 +293,19 @@ function _touches(
         end
     end
     return has_touched
-end
+end
+
+# Extent forwarding
+
+
+function _touches(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return _touches(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function _touches(t1, e1::Extents.Extent, t2::GI.AbstractGeometryTrait, g2)
+    return _touches(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function _touches(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return Extents.touches(e1, e2)
+end
+
+

src/methods/geom_relations/within.jl

@@ -278,3 +278,19 @@ function _within(
     end
     return true
 end
+
+
+# Extent forwarding
+
+
+function _within(t1::GI.AbstractGeometryTrait, g1, t2, e::Extents.Extent)
+    return _within(t1, g1, GI.PolygonTrait(), extent_to_polygon(e))
+end
+function _within(t1, e1::Extents.Extent, t2, g2)
+    return _within(GI.PolygonTrait(), extent_to_polygon(e1), t2, g2)
+end
+function _within(t1, e1::Extents.Extent, t2, e2::Extents.Extent)
+    return Extents.within(e1, e2)
+end
+
+

src/utils/utils.jl

@@ -233,3 +233,40 @@ function lazy_edge_extents(geom)
     end
     for edge in eachedge(geom, Float64))
 end
+
+
+# Extent to polygon
+
+"""
+    extent_to_polygon(ext::Extents.Extent)
+
+Convert an extent to a polygon.
+
+# Examples
+
+```jldoctest
+import GeometryOps as GO, Extents
+    
+ext = Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0))
+GO.extent_to_polygon(ext)
+# output
+GeoInterface.Wrappers.Polygon{false, false}([GeoInterface.Wrappers.LinearRing([(1.0, 1.0), … (3) … , (1.0, 1.0)])])
+```
+"""
+function extent_to_polygon(ext::Extents.Extent{(:X, :Y)})
+    x1, x2 = ext.X
+    y1, y2 = ext.Y
+    return GI.Polygon(StaticArrays.@SVector[GI.LinearRing(StaticArrays.@SVector[(x1, y1), (x2, y1), (x2, y2), (x1, y2), (x1, y1)])])
+end
+
+
+function extent_to_polygon(ext::Extents.Extent{(:Y, :X)})
+    x1, x2 = ext.X
+    y1, y2 = ext.Y
+    return GI.Polygon(StaticArrays.@SVector[GI.LinearRing(StaticArrays.@SVector[(x1, y1), (x2, y1), (x2, y2), (x1, y2), (x1, y1)])])
+end
+
+
+
+
+

test/helpers.jl

@@ -1,5 +1,7 @@
 module TestHelpers
 
+import GeometryOps as GO
+
 using Test, GeoInterface, ArchGDAL, GeometryBasics, LibGEOS
 
 export @test_implementations, @testset_implementations
@@ -31,6 +33,45 @@ const TEST_MODULES = [GeoInterface, ArchGDAL, GeometryBasics, LibGEOS]
     function GeoInterface.convert(::Type{Vector{_ALL_GB_GEOM_TYPES}}, ::GeoInterface.GeometryCollectionTrait, geoms)
         return _ALL_GB_GEOM_TYPES[GeoInterface.convert(GeometryBasics, g) for g in GeoInterface.getgeom(geoms)]
     end
+
+    function GeoInterface.convert(
+        ::Type{GeometryBasics.LineString}, 
+        type::GeoInterface.LineStringTrait, 
+        geom::GeoInterface.Wrappers.LinearRing{false, false, StaticArraysCore.SVector{N, Tuple{Float64, Float64}}, Nothing, Nothing} where N
+        )
+        return GeoInterface.convert(LineString, GeoInterface.LineStringTrait(), collect(geom.geom))
+    end
+
+    function GeoInterface.convert(
+        ::Type{GeometryBasics.LineString}, 
+        type::GeoInterface.LineStringTrait, 
+        geom::GeoInterface.Wrappers.LinearRing{false, false, StaticArraysCore.SVector{N, Tuple{Float64, Float64}}, Nothing, Nothing} where N
+        )
+        return LineString(Point2{Float64}.(collect(geom.geom)))
+    end
+end
+
+
+@eval ArchGDAL begin
+    function GeoInterface.convert(
+        ::Type{T},
+        type::GeoInterface.PolygonTrait,
+        geom,
+    ) where {T<:IGeometry}
+        f = get(lookup_method, typeof(type), nothing)
+        isnothing(f) && error(
+            "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
+        )
+        poly = createpolygon()
+        foreach(GeoInterface.getring(geom)) do ring
+            xs = GeoInterface.x.(GeoInterface.getpoint(ring)) |> collect
+            ys = GeoInterface.y.(GeoInterface.getpoint(ring)) |> collect
+            subgeom = unsafe_createlinearring(xs, ys)
+            result = GDAL.ogr_g_addgeometrydirectly(poly, subgeom)
+            @ogrerr result "Failed to add linearring."
+        end
+        return poly
+    end
 end
 
 
@@ -52,7 +93,16 @@ function _test_implementations_inner(modules::Union{Expr,Vector}, code::Expr)
     for mod in modules1
         expr = Expr(:block)
         for (var, genkey) in pairs(vars)
-            push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var)))
+            push!(expr.args, :($genkey = if $var isa $(GeoInterface.Extents.Extent)
+                    if $mod in ($GeoInterface, $LibGEOS)
+                        $var
+                    else
+                        $GeoInterface.convert($mod, $(GO.extent_to_polygon)($var))
+                    end
+                else
+                    $GeoInterface.convert($mod, $var)
+                end
+            ))
         end
         push!(expr.args, :(@test $code1))
         push!(tests.args, expr)

test/methods/geom_relations.jl

@@ -1,6 +1,7 @@
 import GeometryOps as GO
 import GeoInterface as GI
 import LibGEOS as LG
+using Extents
 using ..TestHelpers
 
 # Tests of DE-9IM Methods
@@ -65,6 +66,12 @@ mp1 = LG.MultiPolygon([p1])
 mp2 = LG.MultiPolygon([p6, p7])
 gc1 = LG.GeometryCollection([pt1, l5, p6])
 
+ext1 = Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0))
+ext2 = Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5))
+ext3 = Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0))
+ext4 = Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0))  # Touches ext1
+
+
 test_pairs = [
     # Points and geometries
     (pt1, pt1, "pt1", "pt1", "Same point"),
@@ -86,7 +93,7 @@ test_pairs = [
     (pt6, p1, "pt6", "p1", "Point on hole edge"),
     (pt7, p1, "pt7", "p1", "Point inside of polygon hole"),
     (p1, pt5, "p1", "pt5", "Point inside of polygon (order swapped)"),
-    # # Lines and geometries
+    # Lines and geometries
     (l1, l1, "l1", "l1", "Same line"),
     (l2, l1, "l2", "l1", "L2 is one segment of l1"),
     (l3, l1, "l3", "l1", "L3 shares one segment with l1 and has one segment outside"),
@@ -155,7 +162,17 @@ test_pairs = [
     (mpt1, mpt3, "mpt1", "mpt3", "No shared points"),
     (ml1, ml2, "ml1", "ml2", "Lines in ml1 cross and touch ml2"),
     (mp1, mp2, "mp1", "mp2", "Polygons in mp1 are inside hole and overlap"),
-    # (gc1, ml1, "gc1", "ml1", "Make sure collection works with multi-geom"),
+    # Extent tests
+    (ext1, ext1, "ext1", "ext1", "Same extent"),
+    (ext1, ext2, "ext1", "ext2", "Overlapping extents"),
+    (ext1, ext3, "ext1", "ext3", "Disjoint extents"),
+    (ext1, ext4, "ext1", "ext4", "Touching extents"),
+    (ext1, pt1, "ext1", "pt1", "Point inside extent"),
+    (ext1, pt2, "ext1", "pt2", "Point outside extent"),
+    (ext1, l1, "ext1", "l1", "Line inside extent"),
+    (ext1, l2, "ext1", "l2", "Line crossing extent"),
+    (ext1, p1, "ext1", "p1", "Polygon same as extent"),
+    (ext1, p2, "ext1", "p2", "Polygon overlapping extent"),
 ]
 
 function test_geom_relation(GO_f, LG_f, f_name; swap_points=false)
@@ -170,13 +187,26 @@ function test_geom_relation(GO_f, LG_f, f_name; swap_points=false)
     end
 end
 
-@testset "Contains" begin test_geom_relation(GO.contains, LG.contains, "contains"; swap_points = true) end
-@testset "Covered By" begin test_geom_relation(GO.coveredby, LG.coveredby, "coveredby") end
-@testset "Covers" begin test_geom_relation(GO.covers, LG.covers, "covers"; swap_points = true) end
-@testset "Disjoint" begin test_geom_relation(GO.disjoint, LG.disjoint, "disjoint")end
-@testset "Intersect" begin test_geom_relation(GO.intersects, LG.intersects, "intersects") end
-@testset "Touches" begin test_geom_relation(GO.touches, LG.touches, "touches") end
-@testset "Within" begin test_geom_relation(GO.within, LG.within, "within") end
+
+function test_geom_relation(GO_f, alg::GO.Algorithm, f_name; swap_points=false)
+    for (g1, g2, sg1, sg2, sdesc) in test_pairs
+        @testset "$sg1 $sg2 $sdesc" begin
+            if swap_points
+                @test_implementations (GO_f($g2, $g1) == GO_f(alg, $g2, $g1))
+            else
+                @test_implementations GO_f($g1, $g2) == GO_f(alg, $g1, $g2)
+            end
+        end
+    end
+end
+
+@testset "Contains" begin test_geom_relation(GO.contains, GO.GEOS(), "contains"; swap_points = true) end
+@testset "Covered By" begin test_geom_relation(GO.coveredby, GO.GEOS(), "coveredby") end
+@testset "Covers" begin test_geom_relation(GO.covers, GO.GEOS(), "covers"; swap_points = true) end
+@testset "Disjoint" begin test_geom_relation(GO.disjoint, GO.GEOS(), "disjoint")end
+@testset "Intersect" begin test_geom_relation(GO.intersects, GO.GEOS(), "intersects") end
+@testset "Touches" begin test_geom_relation(GO.touches, GO.GEOS(), "touches") end
+@testset "Within" begin test_geom_relation(GO.within, GO.GEOS(), "within") end
 
 
 @testset "Overlaps" begin