add macro to test all implementations

src/methods/clipping/clipping_processor.jl

@@ -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]

src/methods/clipping/union.jl

@@ -62,9 +62,9 @@ function _union(
     if n_pieces == 0 # no crossing points, determine if either poly is inside the other
         a_in_b, b_in_a = _find_non_cross_orientation(a_list, b_list, ext_a, ext_b; exact)
         if a_in_b
-            push!(polys, GI.Polygon([tuples(ext_b)]))
+            push!(polys, GI.Polygon([_linearring(tuples(ext_b))]))
         elseif b_in_a
-            push!(polys,  GI.Polygon([tuples(ext_a)]))
+            push!(polys,  GI.Polygon([_linearring(tuples(ext_a))]))
         else
             push!(polys, tuples(poly_a))
             push!(polys, tuples(poly_b))
@@ -124,6 +124,7 @@ function _add_union_holes!(polys, a_in_b, b_in_a, poly_a, poly_b; exact)
         current_poly = n_a_holes > 0 ? ext_poly_b : poly_a
         # Loop over all holes in both original polygons
         for (i, ih) in enumerate(Iterators.flatten((GI.gethole(poly_a), GI.gethole(poly_b))))
+            ih = _linearring(ih)
             in_ext, _, _ = _line_polygon_interactions(ih, curr_exterior_poly; exact, closed_line = true)
             if !in_ext
                 #= if the hole isn't in the overlapping region between the two polygons, add

src/primitives.jl

@@ -260,6 +260,19 @@ end
     geoms = _maptasks(apply_to_geom, 1:GI.ngeom(geom), threaded)
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
+@inline function _apply(f::F, target::TraitTarget{<:PointTrait}, trait::GI.PolygonTrait, geom;
+    crs=GI.crs(geom), calc_extent=_False(), threaded
+)::(GI.geointerface_geomtype(trait)) where F
+    # We need to force rebuilding a LinearRing not a LineString
+    geoms = _maptasks(1:GI.ngeom(geom), threaded) do i
+        lr = GI.getgeom(geom, i)
+        points = map(GI.getgeom(lr)) do p
+            _apply(f, target, p; crs, calc_extent, threaded=_False())
+        end
+        _linearring(_apply_inner(lr, points, crs, calc_extent))
+    end
+    return _apply_inner(geom, geoms, crs, calc_extent)
+end
 function _apply_inner(geom, geoms, crs, calc_extent::_True)
     # Calculate the extent of the sub geometries
     extent = mapreduce(GI.extent, Extents.union, geoms)

src/transformations/correction/intersecting_polygons.jl

@@ -60,6 +60,7 @@ application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
 
 function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
     union_multipoly = tuples(multipoly)
+    @show union_multipoly
     n_polys = GI.npolygon(multipoly)
     if n_polys > 1
         keep_idx = trues(n_polys)  # keep track of sub-polygons to remove
@@ -141,4 +142,4 @@ function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
         keepat!(diff_multipoly.geom, keep_idx)
     end
     return diff_multipoly
-end
+end

src/utils.jl

@@ -126,3 +126,6 @@ function _point_in_extent(p, extent::Extents.Extent)
     (x1, x2), (y1, y2) = extent.X, extent.Y
     return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2
 end
+
+_linearring(geom::GI.LineString) = GI.LinearRing(parent(geom); extent=geom.extent, crs=geom.crs)
+_linearring(geom::GI.LinearRing) = geom

test/extensions/flexijoins.jl

@@ -1,22 +1,19 @@
 import GeometryOps as GO, GeoInterface as GI
 using FlexiJoins, DataFrames
 
-pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])])
-pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])])
-poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue])
-
 points = tuple.(rand(100), rand(100))
 points_df = DataFrame(geometry = points)
 
-@testset "Basic integration" begin
+pl = GI.Polygon([GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 0)])])
+pu = GI.Polygon([GI.LinearRing([(0, 0), (0, 1), (1, 1), (0, 0)])])
 
-    @test_nowarn joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
+@test_all_implementations "Polygon DataDrame" (pl, pu) begin
+    poly_df = DataFrame(geometry = [pl, pu], color = [:red, :blue])
     # Test that the join happened correctly
     joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
     @test all(GO.contains.((pl,), joined_df.geometry_1[joined_df.color .== :red]))
     @test all(GO.contains.((pu,), joined_df.geometry_1[joined_df.color .== :blue]))
     # Test that within also works
     @test_nowarn joined_df = FlexiJoins.innerjoin((points_df, poly_df), by_pred(:geometry, GO.within, :geometry))
-
 end
 

test/helpers.jl

@@ -0,0 +1,93 @@
+using Test, GeoInterface, ArchGDAL, GeometryBasics, LibGEOS
+
+const TEST_MODULES = [GeoInterface, ArchGDAL, GeometryBasics, LibGEOS]
+
+# Monkey-patch GeometryBasics to have correct methods.
+# TODO: push this up to GB!
+
+@eval GeometryBasics begin
+    # MultiGeometry ncoord implementations
+    GeoInterface.ncoord(::GeoInterface.MultiPolygonTrait, ::GeometryBasics.MultiPolygon{N}) where N = N
+    GeoInterface.ncoord(::GeoInterface.MultiLineStringTrait, ::GeometryBasics.MultiLineString{N}) where N = N
+    GeoInterface.ncoord(::GeoInterface.MultiPointTrait, ::GeometryBasics.MultiPoint{N}) where N = N
+    # LinearRing and LineString confusion
+    GeometryBasics.geointerface_geomtype(::GeoInterface.LinearRingTrait) = LineString
+    function GeoInterface.convert(::Type{LineString}, ::GeoInterface.LinearRingTrait, geom)
+        return GeoInterface.convert(LineString, GeoInterface.LineStringTrait(), geom) # forward to the linestring conversion method
+    end
+    # Line interface
+    GeometryBasics.geointerface_geomtype(::GeoInterface.LineTrait) = Line
+    function GeoInterface.convert(::Type{Line}, ::GeoInterface.LineTrait, geom)
+        p1, p2 = GeoInterface.getpoint(geom)
+        return Line(GeoInterface.convert(Point, GeoInterface.PointTrait(), p1), GeoInterface.convert(Point, GeoInterface.PointTrait(), p2))
+    end
+    # GeometryCollection interface - currently just a large Union
+    const _ALL_GB_GEOM_TYPES = Union{Point, Line, LineString, Polygon, MultiPolygon, MultiLineString, MultiPoint}
+    GeometryBasics.geointerface_geomtype(::GeoInterface.GeometryCollectionTrait) = Vector{_ALL_GB_GEOM_TYPES}
+    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
+end
+
+@eval LibGEOS begin
+    function GI.convert(
+        ::Type{GeometryCollection},
+        ::GeometryCollectionTrait,
+        geom;
+        context = get_global_context(),
+    )
+        return GeometryCollection(GI.convert.((LibGEOS,), GI.getgeom(geom)))
+    end
+end
+
+# Macro to run a block of `code` for multiple modules, 
+# using GeoInterface.convert for each var in `args`
+macro test_all_implementations(args, code::Expr)
+    _test_all_implementations_inner("", args, TEST_MODULES, code)
+end
+macro test_all_implementations(title::String, args, code::Expr)
+    _test_all_implementations_inner(title::String, args, TEST_MODULES, code)
+end
+macro test_all_implementations(args, modules, code::Expr)
+    _test_all_implementations_inner("", args, modules, code)
+end
+macro test_all_implementations(title::String, args, modules, code::Expr)
+    _test_all_implementations_inner(title, args, modules, code)
+end
+
+function _test_all_implementations_inner(title, args, modules, code)
+    args1 = esc(args)
+    code1 = esc(code)
+    modules1 = modules isa Expr ? modules.args : modules
+
+    let_expr = if args isa Symbol # Handle a single variable name
+        quote
+            let $args1 = GeoInterface.convert(mod, $args1)
+                $code1
+            end
+        end
+    else # Handle a tuple of variable names
+        quote
+            let ($args1 = map(g -> GeoInterface.convert(mod, g), $args1))
+                $code1
+            end
+        end
+    end
+    testsets = Expr(:block)
+
+    for mod in modules1
+        expr = quote
+            mod = $mod
+            @testset "$mod" begin
+                $let_expr
+            end
+        end
+        push!(testsets.args, expr)
+    end
+
+    quote 
+        @testset "$($title)" begin
+            $testsets
+        end
+    end
+end

test/methods/angles.jl

@@ -9,39 +9,43 @@ concave_coords = [(0.0, 0.0), (0.0, 1.0), (-1.0, 1.0), (-1.0, 2.0), (2.0, 2.0),
 l2 = GI.LineString(concave_coords)
 l3 = GI.LineString(concave_coords[1:(end - 1)])
 r1 = GI.LinearRing(concave_coords)
-r2 = GI.LinearRing(concave_coords[1:(end - 1)])
-r3 = GI.LinearRing([(1.0, 1.0), (1.0, 1.5), (1.5, 1.5), (1.5, 1.0), (1.0, 1.0)])
+r2 = GI.LinearRing([(1.0, 1.0), (1.0, 1.5), (1.5, 1.5), (1.5, 1.0), (1.0, 1.0)])
 concave_angles = [90.0, 270.0, 90.0, 90.0, 90.0, 90.0]
 
-p1 = GI.Polygon([r3])
+p1 = GI.Polygon([r2])
 p2 = GI.Polygon([[(0.0, 0.0), (0.0, 4.0), (3.0, 0.0), (0.0, 0.0)]])
 p3 = GI.Polygon([[(-3.0, -2.0), (0.0,0.0), (5.0, 0.0), (-3.0, -2.0)]])
 p4 = GI.Polygon([r1])
-p5 = GI.Polygon([r1, r3])
+p5 = GI.Polygon([r1, r2])
 
 mp1 = GI.MultiPolygon([p2, p3])
 c1 = GI.GeometryCollection([pt1, l2, p2])
 
-# Points and lines
-@test isempty(GO.angles(pt1))
-@test isempty(GO.angles(mpt1))
-@test isempty(GO.angles(l1))
-
-# LineStrings and Linear Rings
-@test all(isapprox.(GO.angles(l2), concave_angles, atol = 1e-3))
-@test all(isapprox.(GO.angles(l3), concave_angles[2:(end - 1)], atol = 1e-3))
-@test all(isapprox.(GO.angles(r1), concave_angles, atol = 1e-3))
-@test all(isapprox.(GO.angles(r2), concave_angles, atol = 1e-3))
-
-# Polygons
-p2_angles = [90.0, 36.8699, 53.1301]
-p3_angles = [19.6538, 146.3099, 14.0362]
-@test all(isapprox.(GO.angles(p1), [90.0 for _ in 1:4], atol = 1e-3))
-@test all(isapprox.(GO.angles(p2), p2_angles, atol = 1e-3))
-@test all(isapprox.(GO.angles(p3), p3_angles, atol = 1e-3))
-@test all(isapprox.(GO.angles(p4), concave_angles, atol = 1e-3))
-@test all(isapprox.(GO.angles(p5), vcat(concave_angles, [270.0 for _ in 1:4]), atol = 1e-3))
-
-# Multi-geometries
-@test all(isapprox.(GO.angles(mp1), [p2_angles; p3_angles], atol = 1e-3))
-@test all(isapprox.(GO.angles(c1), [concave_angles; p2_angles], atol = 1e-3))
+# Line is not a widely available geometry type
+@test_all_implementations "line angles" l1 [GeometryBasics, GeoInterface] begin
+    @test isempty(GO.angles(l1))
+end
+
+@test_all_implementations "angles" (pt1, mpt1, l2, l3, r1, p1, p2, p3, p4, p5, mp1, c1) begin
+    # Points and lines
+    @test isempty(GO.angles(pt1))
+    @test isempty(GO.angles(mpt1))
+
+    # LineStrings and Linear Rings
+    @test all(isapprox.(GO.angles(l2), concave_angles, atol = 1e-3))
+    @test all(isapprox.(GO.angles(l3), concave_angles[2:(end - 1)], atol = 1e-3))
+    @test all(isapprox.(GO.angles(r1), concave_angles, atol = 1e-3))
+
+    # Polygons
+    p2_angles = [90.0, 36.8699, 53.1301]
+    p3_angles = [19.6538, 146.3099, 14.0362]
+    @test all(isapprox.(GO.angles(p1), [90.0 for _ in 1:4], atol = 1e-3))
+    @test all(isapprox.(GO.angles(p2), p2_angles, atol = 1e-3))
+    @test all(isapprox.(GO.angles(p3), p3_angles, atol = 1e-3))
+    @test all(isapprox.(GO.angles(p4), concave_angles, atol = 1e-3))
+    @test all(isapprox.(GO.angles(p5), vcat(concave_angles, [270.0 for _ in 1:4]), atol = 1e-3))
+
+    # Multi-geometries
+    @test all(isapprox.(GO.angles(mp1), [p2_angles; p3_angles], atol = 1e-3))
+    @test all(isapprox.(GO.angles(c1), [concave_angles; p2_angles], atol = 1e-3))
+end

test/methods/area.jl

@@ -33,55 +33,64 @@ p4 = LG.Polygon([
 empty_p = LG.readgeom("POLYGON EMPTY")
 mp1 = LG.MultiPolygon([p2, p4])
 empty_mp = LG.readgeom("MULTIPOLYGON EMPTY")
-c = LG.GeometryCollection([p1, p2, r1, empty_l])
+c = LG.GeometryCollection([p1, p2, r1])
+c_with_epty_l = LG.GeometryCollection([p1, p2, r1, empty_l])
 empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY")
 
-# Points, lines, and rings have zero area
-@test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0
-@test GO.area(empty_pt) == LG.area(empty_pt) == 0
-@test GO.area(pt) isa Float64
-@test GO.signed_area(pt, Float32) isa Float32
-@test GO.signed_area(pt) isa Float64
-@test GO.area(pt, Float32) isa Float32
-@test GO.area(mpt) == GO.signed_area(mpt) == LG.area(mpt) == 0
-@test GO.area(empty_mpt) == LG.area(empty_mpt) == 0
-@test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0
-@test GO.area(empty_l) == LG.area(empty_l) == 0
-@test GO.area(ml1) == GO.signed_area(ml1) == LG.area(ml1) == 0
-@test GO.area(empty_ml) == LG.area(empty_ml) == 0
-@test GO.area(r1) == GO.signed_area(r1) == LG.area(r1) == 0
-@test GO.area(empty_r) == LG.area(empty_r) == 0
+@test_all_implementations "That handle empty geoms" (empty_pt, empty_mpt, empty_l, empty_ml, empty_l, empty_r, empty_p, empty_mp, empty_c) [LibGEOS, ArchGDAL] begin 
+    @test GO.area(empty_pt) == LG.area(empty_pt) == 0
+    @test GO.area(empty_mpt) == LG.area(empty_mpt) == 0
+    @test GO.area(empty_l) == LG.area(empty_l) == 0
+    @test GO.area(empty_ml) == LG.area(empty_ml) == 0
+    @test GO.area(empty_r) == LG.area(empty_r) == 0
+    # Empty polygon
+    @test GO.signed_area(empty_p) == 0
+    @test GO.area(empty_p) == LG.area(empty_p) == 0
+    # Empty multipolygon
+    @test GO.area(empty_mp) == LG.area(empty_mp) == 0
+    # Empty collection
+    @test GO.area(c_with_epty_l) == LG.area(c_with_epty_l)
+    @test GO.area(c_with_epty_l, Float32) isa Float32
+    @test GO.area(empty_c) == LG.area(empty_c) == 0
+end
 
-# Polygons have non-zero area
-# CCW polygons have positive signed area
-@test GO.area(p1) == GO.signed_area(p1) == LG.area(p1)
-@test GO.signed_area(p1) > 0
-# Float32 calculations
-@test GO.area(p1) isa Float64
-@test GO.area(p1, Float32) isa Float32
-# CW polygons have negative signed area
-a2 = LG.area(p2)
-@test GO.area(p2) == a2
-@test GO.signed_area(p2) == -a2
-# Winding order of holes doesn't affect sign of signed area
-@test GO.signed_area(p3) == -a2
-# Concave polygon correctly calculates area
-a4 = LG.area(p4)
-@test GO.area(p4) == a4
-@test GO.signed_area(p4) == -a4
-# Empty polygon
-@test GO.area(empty_p) == LG.area(empty_p) == 0
-@test GO.signed_area(empty_p) == 0
+@test_all_implementations "With GeometryCollection" c [LibGEOS, ArchGDAL, GeoInterface] begin 
+    # Geometry collection summed area
+    @test GO.area(c) == LG.area(c)
+    @test GO.area(c, Float32) isa Float32
+end  
 
-# Multipolygon calculations work
-@test GO.area(mp1) == a2 + a4
-@test GO.area(mp1, Float32) isa Float32
-# Empty multipolygon
-@test GO.area(empty_mp) == LG.area(empty_mp) == 0
+@test_all_implementations "all" (pt, mpt, l1, ml1, r1, p1, p2, p3, p4, mp1) begin 
+    # Points, lines, and rings have zero area
+    @test GO.area(pt) == GO.signed_area(pt) == LG.area(pt) == 0
+    @test GO.area(pt) isa Float64
+    @test GO.signed_area(pt, Float32) isa Float32
+    @test GO.signed_area(pt) isa Float64
+    @test GO.area(pt, Float32) isa Float32
+    @test GO.area(mpt) == GO.signed_area(mpt) == LG.area(mpt) == 0
+    @test GO.area(l1) == GO.signed_area(l1) == LG.area(l1) == 0
+    @test GO.area(ml1) == GO.signed_area(ml1) == LG.area(ml1) == 0
+    @test GO.area(r1) == GO.signed_area(r1) == LG.area(r1) == 0
 
+    # Polygons have non-zero area
+    # CCW polygons have positive signed area
+    @test GO.area(p1) == GO.signed_area(p1) == LG.area(p1)
+    @test GO.signed_area(p1) > 0
+    # Float32 calculations
+    @test GO.area(p1) isa Float64
+    @test GO.area(p1, Float32) isa Float32
+    # CW polygons have negative signed area
+    a2 = LG.area(p2)
+    @test GO.area(p2) == a2
+    @test GO.signed_area(p2) == -a2
+    # Winding order of holes doesn't affect sign of signed area
+    @test GO.signed_area(p3) == -a2
+    # Concave polygon correctly calculates area
+    a4 = LG.area(p4)
+    @test GO.area(p4) == a4
+    @test GO.signed_area(p4) == -a4
 
-# Geometry collection summed area
-@test GO.area(c) == LG.area(c)
-@test GO.area(c, Float32) isa Float32
-# Empty collection
-@test GO.area(empty_c) == LG.area(empty_c) == 0
+    # Multipolygon calculations work
+    @test GO.area(mp1) == a2 + a4
+    @test GO.area(mp1, Float32) isa Float32
+end