mostly passing

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

@@ -306,6 +306,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,17 +1,13 @@
 import GeometryOps as GO, GeoInterface as GI
 using FlexiJoins, DataFrames
 
-@testset "Tuple-point DataDrame" begin
-    points = tuple.(rand(100), rand(100))
-    points_df = DataFrame(geometry = points)
-
-    @test_nowarn joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
-end
+points = tuple.(rand(100), rand(100))
+points_df = DataFrame(geometry = points)
 
 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_all_integrations "Polygon DataDrame" (pl, pu) begin
+@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/helpers.jl

@@ -42,16 +42,23 @@ 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)
-    _test_all_implementations_inner("", args, code)
+macro test_all_implementations(args, code::Expr)
+    _test_all_implementations_inner("", args, TEST_MODULES, code)
 end
-macro test_all_implementations(title::String, args, code)
-    _test_all_implementations_inner(string(title, " "), args, code)
+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, code)
+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
@@ -66,11 +73,21 @@ function _test_all_implementations_inner(title, args, code)
             end
         end
     end
-    quote
-        for mod in TEST_MODULES 
-            @testset "$($(isempty(title) ? "" : "$title : " ))$mod" begin
+    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

@@ -6,30 +6,32 @@ 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])
 
-@test_all_implementations (pt1, mpt1, l1, l2, l3, r1, r2, r3, p1, p2, p3, p4, p5, mp1, c1) begin
+# 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))
-    @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]

test/methods/area.jl

@@ -30,25 +30,44 @@ 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")
 
-@test_all_implementations (pt, empty_pt, mpt, empty_mpt, l1, empty_l, ml1, empty_ml, empty_l, r1, empty_r, p1, p2, p3, p4, empty_p, mp1, empty_mp, c, empty_c) begin 
+@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
+
+@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  
+
+@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(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
 
     # Polygons have non-zero area
     # CCW polygons have positive signed area
@@ -67,21 +86,8 @@ empty_c = LG.readgeom("GEOMETRYCOLLECTION EMPTY")
     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
 
     # 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
-
-
-    # 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
-
 end

test/methods/clipping/coverage.jl

@@ -6,7 +6,7 @@ cell_area = 400.0
 # # Points, lines, curves
 pt1 = (1.0, 0.0)
 mpt1 = GI.MultiPoint([pt1, (25.0, 0.0), (100.0, 50.0)])
-l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22, 5.0)])
+l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22.0, 5.0)])
 
 @test_all_implementations "Point" pt1 begin
     @test GO.coverage(pt1, cell_extremes...) == 0.0
@@ -21,7 +21,7 @@ end
 # # Polygons
 
 p1 = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)]])
-p2 = GI.Polygon([[(-10, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]])
+p2 = GI.Polygon([[(-10.0, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]])
 p3 = GI.Polygon([[(5.0, 5.0), (5.0, 15.0), (15.0, 15.0), (15.0, 5.0), (5.0, 5.0)]])
 p4 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (15.0, 25.0), (15.0, 5.0), (5.0, 5.0)]])
 p5 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (25.0, 25.0), (25.0, 5.0), (5.0, 5.0)]])

test/methods/clipping/cut.jl

@@ -9,7 +9,7 @@ p1 = GI.Polygon([r1])
 p2 = GI.Polygon([r1, r2])
 p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (10.0, 5.0), (10.0, 0.0), (0.0, 0.0)]])
 
-@test_all_implementations "Cut Polygon" (l1, l2, l3, l4, l5, r1, r2, p1, p2, p3), begin
+@test_all_implementations "Cut Polygon with Line" (l1, l2, l3, l4, l5, r1, r2, p1, p2, p3) [GeoInterface] begin
     # Cut convex polygon
     cut_polys = GO.cut(p1, l1)
     @test all(GO.equals.(

test/methods/distance.jl

@@ -27,7 +27,22 @@ mp1 = LG.MultiPolygon([p1, p2])
 
 c1 = LG.GeometryCollection([pt1, r1, p1])
 
-@test_all_implementations (pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8, pt9, pt10, pt11, mpt1, l1, r1, r2, r3, r4, r5, p1, p2, mp1#=, c1 =#) begin
+@test_all_implementations "Where LinearRing exists" (pt1, pt2, pt3, pt4, pt5, r1, r2, r3, r4, r5) [LibGEOS, GeoInterface] begin
+    # Point on linear ring
+    @test GO.distance(pt1, r1) == LG.distance(pt1, r1)
+    @test GO.distance(pt3, r1) == LG.distance(pt3, r1)
+    # Point outside of linear ring
+    @test GO.distance(pt5, r1) ≈ LG.distance(pt5, r1)
+    # Point inside of hole created by linear ring
+    @test GO.distance(pt3, r1) ≈ LG.distance(pt3, r1)
+    @test GO.distance(pt4, r1) ≈ LG.distance(pt4, r1)
+end
+
+@test_all_implementations "Where GeometryCollection exists" (pt1, c1) [LibGEOS, ArchGDAL, GeoInterface] begin
+    @test GO.distance(pt1, c1) == LG.distance(pt1, c1)
+end
+
+@test_all_implementations (pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8, pt9, pt10, pt11, mpt1, l1, p1, p2, mp1) begin
     # Point and Point
 
     # Distance from point to same point
@@ -54,15 +69,6 @@ c1 = LG.GeometryCollection([pt1, r1, p1])
 
     # Point and Ring
 
-    # Point on linear ring
-    @test GO.distance(pt1, r1) == LG.distance(pt1, r1)
-    @test GO.distance(pt3, r1) == LG.distance(pt3, r1)
-    # Point outside of linear ring
-    @test GO.distance(pt5, r1) ≈ LG.distance(pt5, r1)
-    # Point inside of hole created by linear ring
-    @test GO.distance(pt3, r1) ≈ LG.distance(pt3, r1)
-    @test GO.distance(pt4, r1) ≈ LG.distance(pt4, r1)
-
     # Point and Polygon
     # Point on polygon exterior edge
     @test GO.distance(pt1, p1) == LG.distance(pt1, p1)
@@ -107,6 +113,3 @@ c1 = LG.GeometryCollection([pt1, r1, p1])
     @test GO.signed_distance(pt11, mp1) ≈
         -(min(LG.distance(pt11, r2), LG.distance(pt11, r3), LG.distance(pt11, r4), LG.distance(pt11, r5)))
 end
-
-# Point and Geometry Collection: convert doesn't work yet
-@test GO.distance(pt1, c1) == LG.distance(pt1, c1)

test/methods/equals.jl

@@ -31,7 +31,7 @@ end
     r2 = LG.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]])
     r3 = GI.LinearRing([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0]])
 
-    @test_all_implementations (l1, l2, l3, r1, r2, r3) begin
+    @test_all_implementations (l1, l2, l3, r1, r2) begin
         # Equal lines
         @test GO.equals(l1, l1) == LG.equals(l1, l1)
         @test GO.equals(l2, l2) == LG.equals(l2, l2)
@@ -41,16 +41,20 @@ end
         # Equal rings
         @test GO.equals(r1, r1) == LG.equals(r1, r1)
         @test GO.equals(r2, r2) == LG.equals(r2, r2)
-        # Test equal rings without closing point
-        @test GO.equals(r2, r3)
-        @test GO.equals(r3, l3)
         # Different rings
         @test GO.equals(r1, r2) == GO.equals(r2, r1) == LG.equals(r1, r2)
         # Equal linear ring and line string
         @test GO.equals(r2, l3) == LG.equals(r2, l3)
         # Equal line string and line
         @test GO.equals(l1, GI.Line([(0.0, 0.0), (0.0, 10.0)]))
     end
+
+    # LibGEOS rejects rings that are not closed, and they are not eaual in GeometryBasics or ArchGDAL?
+    @test_all_implementations (r2, r3, l3) [GeoInterface] begin
+        # Test equal rings without closing point
+        @test GO.equals(r2, r3)
+        @test GO.equals(r3, l3)
+    end
 end
 
 @testset "Polygons/MultiPolygons" begin
@@ -104,7 +108,7 @@ end
     ])
     m3 = LG.MultiPolygon([p3])
 
-    @test_all_implementations "Polygons" (pt1, r1, p1, p2, p3, p4, p5, p6, p7, p8) begin
+    @test_all_implementations "Polygons" (pt1, r1, p1, p2, p3, p4, p5, p6, p7, p9) begin
         # Point and polygon aren't equal
         GO.equals(pt1, p1) == LG.equals(pt1, p1)
         # Linear ring and polygon aren't equal
@@ -116,8 +120,6 @@ end
         @test GO.equals(p2, p6) == LG.equals(p2, p6)
         # Equal but opposite winding orders
         @test GO.equals(p2, p7) == LG.equals(p2, p7)
-        # Equal but without closing point (implied)
-        @test GO.equals(p7, p8) 
         # Different polygons
         @test GO.equals(p1, p2) == LG.equals(p1, p2)
         # Equal polygons with holes
@@ -130,6 +132,11 @@ end
         @test GO.equals(p9, p9) == LG.equals(p9, p9)
     end
 
+    @test_all_implementations "Unclosed Polygons" (p7, p8) [GeometryBasics, GeoInterface] begin
+        # Equal but without closing point (implied)
+        @test GO.equals(p7, p8) 
+    end
+
     @test_all_implementations "MultiPolygons" (p1, m1, m2, m3) begin
         # Equal multipolygon
         @test GO.equals(m1, m1) == LG.equals(m1, m1)