diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index e0ea5b7c2..97484404c 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/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]
diff --git a/src/methods/clipping/union.jl b/src/methods/clipping/union.jl
index cdde8c286..e4b4bda8e 100644
--- a/src/methods/clipping/union.jl
+++ b/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
diff --git a/src/primitives.jl b/src/primitives.jl
index aa3cfdfea..3a6af61a9 100644
--- a/src/primitives.jl
+++ b/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)
diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
index e326b6a16..2223f6b89 100644
--- a/src/transformations/correction/intersecting_polygons.jl
+++ b/src/transformations/correction/intersecting_polygons.jl
@@ -141,4 +141,4 @@ function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
         keepat!(diff_multipoly.geom, keep_idx)
     end
     return diff_multipoly
-end
\ No newline at end of file
+end
diff --git a/src/utils.jl b/src/utils.jl
index 82a1e664d..bf815db2e 100644
--- a/src/utils.jl
+++ b/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
diff --git a/test/extensions/flexijoins.jl b/test/extensions/flexijoins.jl
index 80e0432ac..0d6f4b287 100644
--- a/test/extensions/flexijoins.jl
+++ b/test/extensions/flexijoins.jl
@@ -1,22 +1,23 @@
-import GeometryOps as GO, GeoInterface as GI
-using FlexiJoins, DataFrames
+using Test
+using FlexiJoins
+using DataFrames
+import GeometryOps as GO
+import GeoInterface as GI
+using ..TestHelpers
+
+points = GI.MultiPoint(tuple.(rand(100), rand(100)))
 
 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
-
-    @test_nowarn joined_df = FlexiJoins.innerjoin((poly_df, points_df), by_pred(:geometry, GO.contains, :geometry))
+@testset_implementations "Polygon DataDrame" begin
+    points_df = DataFrame(geometry=collect(GI.getpoint($points)))
+    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 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
 
diff --git a/test/helpers.jl b/test/helpers.jl
new file mode 100644
index 000000000..ce5e11e5b
--- /dev/null
+++ b/test/helpers.jl
@@ -0,0 +1,119 @@
+module TestHelpers
+
+using Test, GeoInterface, ArchGDAL, GeometryBasics, LibGEOS
+
+export @test_implementations, @testset_implementations
+
+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
+
+
+# Macro to run a block of `code` for multiple modules,
+# using GeoInterface.convert for each var in `args`
+macro test_implementations(code::Expr)
+    _test_implementations_inner(TEST_MODULES, code)
+end
+macro test_implementations(modules::Union{Expr,Vector}, code::Expr)
+    _test_implementations_inner(modules, code)
+end
+
+function _test_implementations_inner(modules::Union{Expr,Vector}, code::Expr)
+    vars = Dict{Symbol,Symbol}()
+    code1 = _quasiquote!(code, vars)
+    modules1 = modules isa Expr ? modules.args : modules
+    tests = Expr(:block)
+
+    for mod in modules1
+        expr = Expr(:block)
+        for (var, genkey) in pairs(vars)
+            push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var)))
+        end
+        push!(expr.args, :(@test $code1))
+        push!(tests.args, expr)
+    end
+
+    return esc(tests)
+end
+
+# Macro to run a block of `code` for multiple modules,
+# using GeoInterface.convert for each var in `args`
+macro testset_implementations(code::Expr)
+    _testset_implementations_inner("", TEST_MODULES, code)
+end
+macro testset_implementations(arg, code::Expr)
+    if arg isa String || arg isa Expr && arg.head == :string
+        _testset_implementations_inner(arg, TEST_MODULES, code)
+    else
+        _testset_implementations_inner("", arg, code)
+    end
+end
+macro testset_implementations(title, modules::Union{Expr,Vector}, code::Expr)
+    _testset_implementations_inner(title, modules, code)
+end
+
+function _testset_implementations_inner(title, modules::Union{Expr,Vector}, code::Expr)
+    vars = Dict{Symbol,Symbol}()
+    code1 = _quasiquote!(code, vars)
+    modules1 = modules isa Expr ? modules.args : modules
+    testsets = Expr(:block)
+
+    for mod in modules1
+        expr = Expr(:block)
+        for (var, genkey) in pairs(vars)
+            push!(expr.args, :($genkey = $GeoInterface.convert($mod, $var)))
+        end
+        label = title == "" ? "$mod" : "$title: $mod"
+        push!(expr.args, :(@testset $label $code1))
+        push!(testsets.args, expr)
+    end
+
+    return esc(testsets)
+end
+
+# Taken from BenchmarkTools.jl
+_quasiquote!(ex, vars) = ex
+function _quasiquote!(ex::Expr, vars::Dict)
+    if ex.head === :($)
+        v = ex.args[1]
+        gen = if v isa Symbol 
+            haskey(vars, v) ? vars[v] : gensym(v)
+        else
+            gensym()
+        end
+        vars[v] = gen
+        return v 
+    elseif ex.head !== :quote
+        for i in 1:length(ex.args)
+            ex.args[i] = _quasiquote!(ex.args[i], vars)
+        end
+    end
+    return ex
+end
+
+end
diff --git a/test/methods/angles.jl b/test/methods/angles.jl
index 0f6c55160..4004642a5 100644
--- a/test/methods/angles.jl
+++ b/test/methods/angles.jl
@@ -1,5 +1,9 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import GeoInterface as GI
+import GeometryOps as GO 
+import GeometryBasics as GB 
+import LibGEOS as LG
+using ..TestHelpers
 
 pt1 = GI.Point((0.0, 0.0))
 mpt1 = GI.MultiPoint([pt1, pt1])
@@ -9,39 +13,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
+@testset_implementations "line angles" [GB, GI] begin
+    @test isempty(GO.angles($l1))
+end
+
+@testset_implementations "angles" 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
diff --git a/test/methods/area.jl b/test/methods/area.jl
index c1b4a495e..3b2c7ea01 100644
--- a/test/methods/area.jl
+++ b/test/methods/area.jl
@@ -1,5 +1,8 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import GeoInterface as GI
+import GeometryOps as GO 
+import LibGEOS as LG
+using ..TestHelpers
 
 pt = LG.Point([0.0, 0.0])
 empty_pt = LG.readgeom("POINT EMPTY")
@@ -33,55 +36,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
+@testset_implementations "That handle empty geoms" [LG] 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
+@testset "With GeometryCollection" 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
+@testset_implementations "all" 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
\ No newline at end of file
+    # Multipolygon calculations work
+    @test GO.area($mp1) == a2 + a4
+    @test GO.area($mp1, Float32) isa Float32
+end
diff --git a/test/methods/barycentric.jl b/test/methods/barycentric.jl
index 75319a88d..a6faa9bd9 100644
--- a/test/methods/barycentric.jl
+++ b/test/methods/barycentric.jl
@@ -1,30 +1,34 @@
 using Test
 import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
-using GeometryOps, GeoInterface, GeometryBasics
+using GeometryOps, GeoInterface, GeometryBasics, barycentric_coordinates, MeanValue
 
-@testset "Mean value coordinates" begin
-    @testset "Preserving return type" begin
-        @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2{BigFloat}(1,1)) isa Vector{BigFloat}
-        @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2f(1,1)) isa Vector{BigFloat} # keep the most precise type
-        @test barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) isa Vector{Float64}
-        @test barycentric_coordinates(MeanValue(), Point2{Float32}[(0,0), (1,0), (0,1)], Point2{Float32}(1,1)) isa Vector{Float32}
-    end
-    @testset "Triangle coordinates" begin
-        # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (0.5,0.5) are (0.5,0.25,0.25)
-        @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(0.5,0.5)) .== (0.5,0.25,0.25))
-        # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (1,1) are (-1,1,1)
-        @test all(barycentric_coordinates(MeanValue(), Point2f[(0,0), (1,0), (0,1)], Point2f(1,1)) .≈ (-1,1,1))
-        # Test that calculations with different number types (in this case Float64) are also accurate
-        @test all(barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) .≈ (-1,1,1))
+t1, p1 = GI.LineString([(0.,0.), (1.,0.), (0.,1.)]), (0.5, 0.5)
+t2, p2 = GI.LineString([(0.,0.), (1.,0.), (0.,1.)]), (1., 1.)
+t3, p3 = GI.LineString(GeometryBasics.Point2d[(0,0), (1,0), (0,1)]), (1., 1.)
+
+@testset "Triangle coordinates" begin
+    # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (0.5,0.5) are (0.5,0.25,0.25)
+    @test all(barycentric_coordinates(MeanValue(), t1, p1) .== (0.5,0.25,0.25))
+    # Test that the barycentric coordinates for (0,0), (1,0), (0,1), and (1,1) are (-1,1,1)
+    @test all(barycentric_coordinates(MeanValue(), t2, p2) .≈ (-1,1,1))
+    # Test that calculations with different number types (in this case Float64) are also accurate
+    @test all(barycentric_coordinates(MeanValue(), Point2{Float64}, Point2{Float64}(1,1)) .≈ (-1,1,1))
+end
+
+@testset "Preserving return type" begin
+    @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2{BigFloat}(1,1)) isa Vector{BigFloat}
+    @test barycentric_coordinates(MeanValue(), Point2{BigFloat}[(0,0), (1,0), (0,1)], Point2f(1,1)) isa Vector{BigFloat} # keep the most precise type
+    @test barycentric_coordinates(MeanValue(), Point2{Float64}[(0,0), (1,0), (0,1)], Point2{Float64}(1,1)) isa Vector{Float64}
+    @test barycentric_coordinates(MeanValue(), Point2{Float32}[(0,0), (1,0), (0,1)], Point2{Float32}(1,1)) isa Vector{Float32}
+end
+
+@testset "Tests for helper methods" begin
+    @testset "`t_value`" begin
+        @test GO.t_value(Point2f(0,0), Point2f(1,0), 1, 1) == 0
+        @test GO.t_value(Point2f(0, 1), Point2f(1, 0), 1, 2) == -0.5f0
     end
-    @testset "Tests for helper methods" begin
-        @testset "`t_value`" begin
-            @test GeometryOps.t_value(Point2f(0,0), Point2f(1,0), 1, 1) == 0
-            @test GeometryOps.t_value(Point2f(0, 1), Point2f(1, 0), 1, 2) == -0.5f0
-        end
-        @testset "`_det`" begin
-            @test GeometryOps._det((1,0), (0,1)) == 1f0
-            @test GeometryOps._det(Point2f(1, 2), Point2f(3, 4)) == -2.0f0
-        end
+    @testset "`_det`" begin
+        @test GO._det((1,0), (0,1)) == 1f0
+        @test GO._det(Point2f(1, 2), Point2f(3, 4)) == -2.0f0
     end
 end
diff --git a/test/methods/centroid.jl b/test/methods/centroid.jl
index 1a7b81a8c..224a96da5 100644
--- a/test/methods/centroid.jl
+++ b/test/methods/centroid.jl
@@ -1,116 +1,125 @@
 using Test
-import GeoInterface as GI, 
-    GeometryOps as GO, 
-    LibGEOS as LG,
-    ArchGDAL as AG
+import GeoInterface as GI
+import GeometryOps as GO 
+import LibGEOS as LG
+import ArchGDAL as AG
+using ..TestHelpers
 
 @testset "Lines/Rings" begin
-    # Basic line string
-    l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
-    c1, len1 = GO.centroid_and_length(l1)
-    c1_from_LG = LG.centroid(l1)
-    @test c1[1] ≈ GI.x(c1_from_LG)
-    @test c1[2] ≈ GI.y(c1_from_LG)
-    @test len1 ≈ 20.0
 
-    # Spiral line string
+    l1 = LG.LineString([[0.0, 0.0], [10.0, 0.0], [10.0, 10.0]])
     l2 = LG.LineString([[0.0, 0.0], [2.5, -2.5], [-5.0, -3.0], [-4.0, 6.0], [10.0, 10.0], [12.0, -14.56]])
-    c2, len2 = GO.centroid_and_length(l2)
-    c2_from_LG = LG.centroid(l2)
-    @test c2[1] ≈ GI.x(c2_from_LG) 
-    @test c2[2] ≈ GI.y(c2_from_LG)
-    @test len2 ≈ 59.3090856788928
-
-    # Test that non-closed line strings throw an error for centroid_and_area
-    @test_throws AssertionError GO.centroid_and_area(l2)
-
     # Basic linear ring - note that this still uses weighting by length
     r1 = LG.LinearRing([[0.0, 0.0], [3456.0, 7894.0], [6291.0, 1954.0], [0.0, 0.0]])
-    c3 = GO.centroid(r1)
-    c3_from_LG = LG.centroid(r1)
-    @test c3[1] ≈ GI.x(c3_from_LG) 
-    @test c3[2] ≈ GI.y(c3_from_LG)
+    
+    @testset_implementations "Basic line string" begin
+        c1_from_LG = LG.centroid($l1)
+        c1, len1 = GO.centroid_and_length($l1)
+        @test c1[1] ≈ GI.x(c1_from_LG)
+        @test c1[2] ≈ GI.y(c1_from_LG)
+        @test len1 ≈ 20.0
+    end
+
+    @testset_implementations "Spiral line string" begin
+        c2_from_LG = LG.centroid($l2)
+        c2, len2 = GO.centroid_and_length($l2)
+        @test c2[1] ≈ GI.x(c2_from_LG) 
+        @test c2[2] ≈ GI.y(c2_from_LG)
+        @test len2 ≈ 59.3090856788928
+
+        # Test that non-closed line strings throw an error for centroid_and_area
+        @test_throws AssertionError GO.centroid_and_area(l2)
+    end
+
+    @testset_implementations "Basic linear ring" begin
+        c3_from_LG = LG.centroid($r1)
+        c3 = GO.centroid($r1)
+        @test c3[1] ≈ GI.x(c3_from_LG) 
+        @test c3[2] ≈ GI.y(c3_from_LG)
+    end
 end
 
-@testset "Polygons" begin
-    # Basic rectangle
-    p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
-    c1 = GO.centroid(p1)
+# Basic rectangle
+p1 = AG.fromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))")
+# Concave c-shaped polygon
+p2 = LG.Polygon([[
+    [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
+    [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0],
+]])
+# Randomly generated polygon with lots of sides
+p3 = LG.Polygon([[
+    [14.567, 8.974], [13.579, 8.849], [12.076, 8.769], [11.725, 8.567],
+    [11.424, 6.451], [10.187, 7.712], [8.187, 6.795], [8.065, 8.096],
+    [6.827, 8.287], [7.1628, 8.9221], [5.8428, 10.468], [7.987, 11.734],
+    [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340],
+    [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974],
+]])
+# Polygon with one hole
+p4 = LG.Polygon([
+    [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+    [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
+])
+# Polygon with two holes
+p5 = LG.Polygon([
+    [[-10.0, -10.0], [-2.0, 0.0], [6.0, -10.0], [-10.0, -10.0]],
+    [[-8.0, -8.0], [-8.0, -7.0], [-4.0, -7.0], [-4.0, -8.0], [-8.0, -8.0]],
+    [[-3.0,-9.0], [3.0, -9.0], [3.0, -8.5], [-3.0, -8.5], [-3.0, -9.0]],
+])
+# Same polygon as P5 but using a GeoInterface polygon
+p6 = GI.Polygon([
+    [(-10.0, -10.0), (-2.0, 0.0), (6.0, -10.0), (-10.0, -10.0)],
+    [(-8.0, -8.0), (-8.0, -7.0), (-4.0, -7.0), (-4.0, -8.0), (-8.0, -8.0)],
+    [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)],
+])
+
+@testset_implementations "Polygons" begin
+    c1 = GO.centroid($p1)
     c1 .≈ (5, 5)
     @test GI.x(c1) ≈ 5
     @test GI.y(c1) ≈ 5
 
-    # Concave c-shaped polygon
-    p2 = LG.Polygon([[
-        [11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
-        [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0],
-    ]])
-    c2, area2 = GO.centroid_and_area(p2)
-    c2_from_LG = LG.centroid(p2)
+    c2, area2 = GO.centroid_and_area($p2)
+    c2_from_LG = LG.centroid($p2)
     @test c2[1] ≈ GI.x(c2_from_LG)
     @test c2[2] ≈ GI.y(c2_from_LG)
     @test area2 ≈ LG.area(p2)
 
-    # Randomly generated polygon with lots of sides
-    p3 = LG.Polygon([[
-        [14.567, 8.974], [13.579, 8.849], [12.076, 8.769], [11.725, 8.567],
-        [11.424, 6.451], [10.187, 7.712], [8.187, 6.795], [8.065, 8.096],
-        [6.827, 8.287], [7.1628, 8.9221], [5.8428, 10.468], [7.987, 11.734],
-        [7.989, 12.081], [8.787, 11.930], [7.568, 13.926], [9.330, 13.340],
-        [9.6817, 13.913], [10.391, 12.222], [12.150, 12.032], [14.567, 8.974],
-    ]])
-    c3, area3 = GO.centroid_and_area(p3)
-    c3_from_LG = LG.centroid(p3)
+    c3, area3 = GO.centroid_and_area($p3)
+    c3_from_LG = LG.centroid($p3)
     @test c3[1] ≈ GI.x(c3_from_LG)
     @test c3[2] ≈ GI.y(c3_from_LG)
-    @test area3 ≈ LG.area(p3)
+    @test area3 ≈ LG.area($p3)
 
-    # Polygon with one hole
-    p4 = LG.Polygon([
-        [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
-        [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
-    ])
-    c4, area4 = GO.centroid_and_area(p4)
-    c4_from_LG = LG.centroid(p4)
+    c4, area4 = GO.centroid_and_area($p4)
+    c4_from_LG = LG.centroid($p4)
     @test c4[1] ≈ GI.x(c4_from_LG)
     @test c4[2] ≈ GI.y(c4_from_LG)
-    @test area4 ≈ LG.area(p4)
+    @test area4 ≈ LG.area($p4)
 
-    # Polygon with two holes
-    p5 = LG.Polygon([
-        [[-10.0, -10.0], [-2.0, 0.0], [6.0, -10.0], [-10.0, -10.0]],
-        [[-8.0, -8.0], [-8.0, -7.0], [-4.0, -7.0], [-4.0, -8.0], [-8.0, -8.0]],
-        [[-3.0,-9.0], [3.0, -9.0], [3.0, -8.5], [-3.0, -8.5], [-3.0, -9.0]],
-    ])
-    c5 = GO.centroid(p5)
-    c5_from_LG = LG.centroid(p5)
+    c5 = GO.centroid($p5)
+    c5_from_LG = LG.centroid($p5)
     @test c5[1] ≈ GI.x(c5_from_LG)
     @test c5[2] ≈ GI.y(c5_from_LG)
-
-    # Same polygon as P5 but using a GeoInterface polygon
-    p6 = GI.Polygon([
-        [(-10.0, -10.0), (-2.0, 0.0), (6.0, -10.0), (-10.0, -10.0)],
-        [(-8.0, -8.0), (-8.0, -7.0), (-4.0, -7.0), (-4.0, -8.0), (-8.0, -8.0)],
-        [(-3.0, -9.0), (3.0, -9.0), (3.0, -8.5), (-3.0, -8.5), (-3.0, -9.0)],
-    ])
-    c6 = GO.centroid(p6)
+    c6 = GO.centroid($p6)
     @test all(c5 .≈ c6)
 end
-@testset "MultiPolygons" begin
-    # Combine polygons made above
-    m1 = LG.MultiPolygon([
-        [
-            [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
-            [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]],
-        ],
-        [
-            [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
-            [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
-        ]
-    ])
-    c1, area1 = GO.centroid_and_area(m1)
-    c1_from_LG = LG.centroid(m1)
+
+# Combine polygons made above
+m1 = LG.MultiPolygon([
+    [
+        [[11.0, 0.0], [11.0, 3.0], [14.0, 3.0], [14.0, 2.0], [12.0, 2.0],
+        [12.0, 1.0], [14.0, 1.0], [14.0, 0.0], [11.0, 0.0]],
+    ],
+    [
+        [[0.0, 0.0], [10.0, 0.0], [10.0, 10.0], [0.0, 10.0], [0.0, 0.0]],
+        [[2.3, 2.7], [2.5, 4.9], [4.1, 5.2], [4.2, 1.9], [2.3, 2.7]],
+    ]
+])
+
+@testset_implementations "MultiPolygons" begin
+    c1, area1 = GO.centroid_and_area($m1)
+    c1_from_LG = LG.centroid($m1)
     @test c1[1] ≈ GI.x(c1_from_LG)
     @test c1[2] ≈ GI.y(c1_from_LG)
-    @test area1 ≈ LG.area(m1)
+    @test area1 ≈ LG.area($m1)
 end
diff --git a/test/methods/clipping/coverage.jl b/test/methods/clipping/coverage.jl
index aa4b388df..7127bb4c6 100644
--- a/test/methods/clipping/coverage.jl
+++ b/test/methods/clipping/coverage.jl
@@ -1,64 +1,80 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
-using GeoInterface
-using GeoInterface.Extents
+import GeoInterface as GI
+import GeometryBasics as GB
+import GeometryOps as GO
+import LibGEOS as LG
+using ..TestHelpers
 
 cell_extremes = (0.0, 20.0, 0.0, 20.0)
-cell_extent = Extents.Extent(X = (0.0, 20.0), Y = (0.0, 20.0))
+cell_extent = GI.Extents.Extent(X = (0.0, 20.0), Y = (0.0, 20.0))
 cell_poly = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)]])
 cell_area = 400.0
 
-# # Points, lines, curves
 pt1 = (1.0, 0.0)
-@test GO.coverage(pt1, cell_extremes...) == 0.0
 mpt1 = GI.MultiPoint([pt1, (25.0, 0.0), (100.0, 50.0)])
-@test GO.coverage(mpt1, cell_extremes...) == 0.0
-l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22, 5.0)])
-@test GO.coverage(l1, cell_extremes...) == 0.0
+l1 = GI.LineString([(-1.0, 5.0), (5.0, 10.0), (10.0, -2.0), (22.0, 5.0)])
+
+@testset_implementations "Points, lines, curves" begin
+    @test GO.coverage($pt1, cell_extremes...) == 0.0
+    @test GO.coverage($mpt1, cell_extremes...) == 0.0
+    @test GO.coverage($l1, cell_extremes...) == 0.0
+end
 
 # # Polygons
 
-# polygon is the same as the cell (input is extent)
 p1 = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)]])
-@test GO.coverage(p1, cell_extent) == cell_area
-# polygon is the same as the cell (input is min/max values)
-@test GO.coverage(p1, cell_extremes...) == cell_area
-# polygon is bigger than the cell
-p2 = GI.Polygon([[(-10, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]])
-@test GO.coverage(p2, cell_extremes...) == cell_area
+p2 = GI.Polygon([[(-10.0, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]])
+p2b = GI.Polygon([[(-10, -10.0), (-10.0, 30.0), (30.0, 30.0), (300.0, -10.0), (-10.0, -10.0)]])
 # polygon is completely inside of cell
 p3 = GI.Polygon([[(5.0, 5.0), (5.0, 15.0), (15.0, 15.0), (15.0, 5.0), (5.0, 5.0)]])
-@test GO.coverage(p3, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p3))
-# polygon exits cell through one edge
 p4 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (15.0, 25.0), (15.0, 5.0), (5.0, 5.0)]])
-@test GO.coverage(p4, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p4))
 p5 = GI.Polygon([[(5.0, 5.0), (5.0, 25.0), (25.0, 25.0), (25.0, 5.0), (5.0, 5.0)]])
-@test GO.coverage(p5, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, p5))
-# polygon exits cell through multiple edges (north and east)
 p6 = GI.Polygon([[(20.8826, 6.4239), (15.9663, 2.3014), (8.6078, 2.0995), (2.6849, 6.4088),
     (0.8449, 12.7452), (3.0813, 19.1654), (9.1906, 23.2520), (15.5835, 22.9101),
     (20.9143, 18.5933), (20.8826, 6.4239)]])
-@test GO.coverage(p6, cell_extremes...) ≈ LG.area(LG.intersection(p6, cell_poly))
-# polygon exits cell through multiple edges (west and east)
 p7 = GI.Polygon([[(-5.0, 10.0), (-5.0, 25.0), (25.0, 25.0), (25.0, 10.0), (-5.0, 10.0)]])
-@test GO.coverage(p7, cell_extremes...) ≈ LG.area(LG.intersection(p7, cell_poly))
-# non-convex polygon split into two pieces (same wall)
 p8 =  GI.Polygon([[(-10.0, 15.0), (10.0, 15.0), (10.0, 12.0), (-5.0, 12.0), (-5.0, 9.0), (10.0, 9.0), (10.0, 6.0), (-10.0, 6.0), (-10.0, 15.0)]])
-@test GO.coverage(p8, cell_extremes...) ≈ LG.area(LG.intersection(p8, cell_poly))
-# counter-clockwise polygon
 p9 =  GI.Polygon([[(-10.0, 15.0), (-10.0, 6.0), (10.0, 6.0), (10.0, 9.0), (-5.0, 9.0), (-5.0, 12.0), (10.0, 12.0), (10.0, 15.0), (-10.0, 15.0),]])
-@test GO.coverage(p9, cell_extremes...) ≈ LG.area(LG.intersection(p9, cell_poly))
-# polygon with a hole
 p10 = GI.Polygon([[(0.0, 0.0), (0.0, 20.0), (20.0, 20.0), (20.0, 0.0), (0.0, 0.0)],
     [(10.0, 10.0), (10.0, 15.0), (15.0, 15.0), (15.0, 10.0), (10.0, 10.0)],
     [(7.0, 7.0), (5.0, 7.0), (5.0, 5.0), (7.0, 7.0)],
 ])
-@test GO.coverage(p10, cell_extremes...) ≈ LG.area(LG.intersection(p10, cell_poly))
 # non-convex polygon split into two pieces (different walls)
 p11 = GI.Polygon([[(-10.0, 15.0), (-10.0, -10.0), (15.0, -10.0), (15.0, 5.0), (10.0, 5.0),
     (10.0, -5.0), (-5.0, -5.0), (-5.0, 10.0), (5.0, 10.0), (5.0, 15.0), (-10.0, 15.0)]])
-@test GO.coverage(p11, cell_extremes...) ≈ LG.area(LG.intersection(p11, cell_poly))
 
-# # MultiPolygons
+@testset_implementations "Polygons" [GI, GB, LG] begin
+    # polygon is the same as the cell (input is min/max values)
+    @test GO.coverage($p1, cell_extent) == cell_area
+    @test GO.coverage($p1, cell_extremes...) == cell_area
+    @test GO.coverage($p10, cell_extremes...) ≈ LG.area(LG.intersection($p10, cell_poly))
+    # polygon is the same as the cell (input is extent)
+    @test GO.coverage($p1, cell_extent) == cell_area
+    # polygon is the same as the cell (input is min/max values)
+    @test GO.coverage($p1, cell_extremes...) == cell_area
+    # polygon is bigger than the cell
+    @test GO.coverage($p2, cell_extremes...) == cell_area
+    # polygon is bigger than the cell
+    @test_implementations GO.coverage($p2b, cell_extremes...) == cell_area
+    # polygon is completly inside of cell
+    @test GO.coverage($p3, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p3))
+    # polygon exits cell through one edge
+    @test GO.coverage($p4, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p4))
+    @test GO.coverage($p5, cell_extremes...) ≈ LG.area(LG.intersection(cell_poly, $p5))
+    # polygon exits cell through multiple edges (north and east)
+    @test GO.coverage($p6, cell_extremes...) ≈ LG.area(LG.intersection($p6, cell_poly))
+    # polygon exits cell through multiple edges (west and east)
+    @test GO.coverage($p7, cell_extremes...) ≈ LG.area(LG.intersection($p7, cell_poly))
+    # non-convex polygon split into two pieces (same wall)
+    @test GO.coverage($p8, cell_extremes...) ≈ LG.area(LG.intersection($p8, cell_poly))
+    # counter-clockwise polygon
+    @test GO.coverage($p9, cell_extremes...) ≈ LG.area(LG.intersection($p9, cell_poly))
+    # polygon with a hole
+    @test GO.coverage($p10, cell_extremes...) ≈ LG.area(LG.intersection($p10, cell_poly))
+    # non-convex polygon split into two pieces (differnet walls)
+    @test GO.coverage($p11, cell_extremes...) ≈ LG.area(LG.intersection($p11, cell_poly))
+end
+
+# TODO: MultiPolygons and GeometryCollection
 mp1 = GI.MultiPolygon([p1, p2]) == 2cell_area
-gc1 = GI.GeometryCollection([pt1, l1, p1]) == cell_area
\ No newline at end of file
+gc1 = GI.GeometryCollection([pt1, l1, p1]) == cell_area
diff --git a/test/methods/clipping/cut.jl b/test/methods/clipping/cut.jl
index b2e0630e2..7530bde2a 100644
--- a/test/methods/clipping/cut.jl
+++ b/test/methods/clipping/cut.jl
@@ -1,7 +1,10 @@
 using Test
+import GeoInterface as GI
+import GeometryOps as GO
+using ..TestHelpers
 
-import GeoInterface as GI, GeometryOps as GO
-
+# TODO this should also work with a LineString
+# Some geom packages don't have `Line`
 l1 = GI.Line([(5.0, -5.0), (5.0, 15.0)])
 l2 = GI.Line([(-1.0, 8.0), (2.0, 11.0)])
 l3 = GI.Line([(-1.0, 6.0), (11.0, 6.0)])
@@ -13,9 +16,8 @@ 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)]])
 
-@testset "Cut Polygon" begin
-    # Cut convex polygon
-    cut_polys = GO.cut(p1, l1)
+@testset_implementations "Cut convex polygon" begin
+    cut_polys = GO.cut($p1, l1)
     @test all(GO.equals.(
         cut_polys,
         [
@@ -24,7 +26,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1
         ],
     ))
 
-    cut_polys = GO.cut(p1, l2)
+    cut_polys = GO.cut($p1, l2)
     @test all(GO.equals.(
         cut_polys,
         [
@@ -34,7 +36,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1
     ))
 
     # Cut convex polygon with hole
-    cut_polys = GO.cut(p2, l1)
+    cut_polys = GO.cut($p2, l1)
     @test all(GO.equals.(
         cut_polys,
         [
@@ -42,7 +44,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1
             GI.Polygon([[(5.0, 0.0), (10.0, 0.0), (10.0, 10.0), (5.0, 10.0), (5.0, 4.0), (8.0, 4.0), (8.0, 2.0), (5.0, 2.0), (5.0, 0.0)]]),
         ],
     ))
-    cut_polys = GO.cut(p2, l3)
+    cut_polys = GO.cut($p2, l3)
     @test all(GO.equals.(
         cut_polys,
         [
@@ -52,7 +54,7 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1
     ))
 
     # Cut concave polygon into three pieces
-    cut_polys = GO.cut(p3, l3)
+    cut_polys = GO.cut($p3, l3)
     @test all(GO.equals.(
         cut_polys,
         [
@@ -63,8 +65,8 @@ p3 = GI.Polygon([[(0.0, 0.0), (0.0, 5.0), (2.5, 7.5), (5.0, 5.0), (7.5, 7.5), (1
     ))
 
     # Line doesn't cut through polygon
-    cut_polys = GO.cut(p1, l4)
-    @test all(GO.equals.(cut_polys, [p1]))
-    cut_polys = GO.cut(p1, l5)
-    @test all(GO.equals.(cut_polys, [p1]))
-end
\ No newline at end of file
+    cut_polys = GO.cut($p1, l4)
+    @test all(GO.equals.(cut_polys, [$p1]))
+    cut_polys = GO.cut($p1, l5)
+    @test all(GO.equals.(cut_polys, [$p1]))
+end
diff --git a/test/methods/clipping/intersection_points.jl b/test/methods/clipping/intersection_points.jl
index 28c5a00f4..4cfc7074a 100644
--- a/test/methods/clipping/intersection_points.jl
+++ b/test/methods/clipping/intersection_points.jl
@@ -1,5 +1,8 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import GeoInterface as GI
+import GeometryOps as GO
+import LibGEOS as LG
+using ..TestHelpers
 
 l1 = GI.LineString([(90000.0, 1000.0), (90000.0, 22500.0), (95000.0, 22500.0), (95000.0, 1000.0), (90000.0, 1000.0)])
 l2 = GI.LineString([(90000.0, 7500.0), (107500.0, 27500.0), (112500.0, 27500.0), (95000.0, 7500.0), (90000.0, 7500.0)])
@@ -10,25 +13,22 @@ l6 = GI.LineString([(0.0, 25000.0), (0.0, 29000.0), (20000.0, 29000.0), (20000.0
 
 p1, p2 = GI.Polygon([l1]), GI.Polygon([l2])
 
-# Three intersection points
-GO_l1_l2_mp = GI.MultiPoint(GO.intersection_points(l1, l2))
-LG_l1_l2_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l1, l2))))
-@test GO.equals(GO_l1_l2_mp, LG_l1_l2_mp)
+@testset_implementations begin
+    # Three intersection points
+    LG_l1_l2_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l1, $l2))))
+    @test GO.equals(GI.MultiPoint(GO.intersection_points($l1, $l2)), LG_l1_l2_mp)
 
-# Four intersection points with large intersection
-GO_l3_l4_mp = GI.MultiPoint(GO.intersection_points(l3, l4))
-LG_l3_l4_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l3, l4))))
-@test GO.equals(GO_l3_l4_mp, LG_l3_l4_mp)
+    # Four intersection points with large intersection
+    LG_l3_l4_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l3, $l4))))
+    @test GO.equals(GI.MultiPoint(GO.intersection_points($l3, $l4)), LG_l3_l4_mp)
 
-# Four intersection points with very small intersection
-GO_l5_l6_mp = GI.MultiPoint(GO.intersection_points(l5, l6))
-LG_l5_l6_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection(l5, l6))))
-@test GO.equals(GO_l5_l6_mp, LG_l5_l6_mp)
+    # Four intersection points with very small intersection
+    LG_l5_l6_mp = GI.MultiPoint(collect(GI.getpoint(LG.intersection($l5, $l6))))
+    @test GO.equals(GI.MultiPoint(GO.intersection_points($l5, $l6)), LG_l5_l6_mp)
 
-# Test that intersection points between lines and polygons is equivalent
-GO_p1_p2_mp = GI.MultiPoint(GO.intersection_points(p1, p2))
-@test GO.equals(GO_p1_p2_mp, GO_l1_l2_mp)
+    # Test that intersection points between lines and polygons is equivalent
+    @test GO.equals(GI.MultiPoint(GO.intersection_points($p1, $p2)), GI.MultiPoint(GO.intersection_points($l1, $l2)))
 
-# No intersection points between polygon and line
-GO_p1_l6 = GO.intersection_points(p1, l6)
-@test isempty(GO_p1_l6)
+    # No intersection points between polygon and line
+    @test isempty(GO.intersection_points($p1, $l6))
+end
diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index ccb538218..432650e94 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -1,5 +1,8 @@
 using Test
-import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG
+import GeometryOps as GO
+import GeoInterface as GI
+import LibGEOS as LG
+using ..TestHelpers
 
 # Test of polygon clipping
 p1 = GI.Polygon([[(0.0, 0.0), (5.0, 5.0), (10.0, 0.0), (5.0, -5.0), (0.0, 0.0)]])
@@ -202,9 +205,9 @@ end
 # Test clipping functions and print error message if tests fail
 function test_clipping(GO_f, LG_f, f_name)
     for (p1, p2, sg1, sg2, sdesc) in test_pairs
-        pass_test = compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
-        @test pass_test
-        !pass_test && println("\n↑ TEST INFO: $sg1 $f_name $sg2 - $sdesc \n\n")
+        @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin
+            @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) 
+        end
     end
     return
 end
diff --git a/test/methods/convex_hull.jl b/test/methods/convex_hull.jl
index 62c5d6a58..7caeb0c3d 100644
--- a/test/methods/convex_hull.jl
+++ b/test/methods/convex_hull.jl
@@ -1,5 +1,8 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import GeoInterface as GI
+import GeometryOps as GO 
+import LibGEOS as LG
+using ..TestHelpers
 
 @testset "Basic example" begin
     points = tuple.(rand(100), rand(100))
@@ -12,7 +15,7 @@ import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
     @test isempty(
         setdiff(
             collect(GO.flatten(GO.tuples, GI.PointTrait, hull)), 
-            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GEOS(), points)))
+            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GO.GEOS(), points)))
         )
     )
 end
@@ -33,7 +36,7 @@ end
     @test isempty(
         setdiff(
             collect(GO.flatten(GO.tuples, GI.PointTrait, double_hull)), 
-            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GEOS(), points)))
+            collect(GO.flatten(GO.tuples, GI.PointTrait, GO.convex_hull(GO.GEOS(), points)))
         )
     )
 end
diff --git a/test/methods/distance.jl b/test/methods/distance.jl
index 207feeb4c..a2a48c865 100644
--- a/test/methods/distance.jl
+++ b/test/methods/distance.jl
@@ -1,5 +1,9 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import ArchGDAL as AG
+import GeoInterface as GI
+import GeometryOps as GO 
+import LibGEOS as LG
+using ..TestHelpers
 
 pt1 = LG.Point([0.0, 0.0])
 pt2 = LG.Point([0.0, 1.0])
@@ -30,85 +34,88 @@ mp1 = LG.MultiPolygon([p1, p2])
 
 c1 = LG.GeometryCollection([pt1, r1, p1])
 
-# Point and Point
-
-# Distance from point to same point
-@test GO.distance(pt1, pt1) == LG.distance(pt1, pt1)
-# Distance from point to different point
-@test GO.distance(pt1, pt2) ≈ GO.distance(pt2, pt1) ≈ LG.distance(pt1, pt2)
-# Return types
-@test GO.distance(pt1, pt1) isa Float64
-@test GO.distance(pt1, pt1, Float32) isa Float32
-
-# Point and Line
-
-#Point on line vertex
-@test GO.distance(pt1, l1) ==  GO.distance(l1, pt1) == LG.distance(pt1, l1)
-# Point on line edge
-@test GO.distance(pt2, l1) ==  GO.distance(l1, pt2) == LG.distance(pt2, l1)
-# Point equidistant from both segments
-@test GO.distance(pt3, l1) ≈  GO.distance(l1, pt3) ≈ LG.distance(pt3, l1)
-# Point closer to one segment than another
-@test GO.distance(pt4, l1) ≈  GO.distance(l1, pt4) ≈ LG.distance(pt4, l1)
-# Return types
-@test GO.distance(pt1, l1) isa Float64
-@test GO.distance(pt1, l1, Float32) isa Float32
-
-# 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)
-@test GO.signed_distance(pt1, p1) == 0
-@test GO.distance(pt2, p1) == LG.distance(pt2, p1)
-# Point on polygon hole edge
-@test GO.distance(pt4, p1) == LG.distance(pt4, p1)
-@test GO.signed_distance(pt4, p1) == 0
-@test GO.distance(pt6, p1) == LG.distance(pt6, p1)
-# Point inside of polygon
-@test GO.distance(pt3, p1) == LG.distance(pt3, p1)
-@test GO.signed_distance(pt3, p1) ≈
-    -(min(LG.distance(pt3, r2), LG.distance(pt3, r3), LG.distance(pt3, r4)))
-@test GO.distance(pt7, p1) == LG.distance(pt7, p1)
-@test GO.signed_distance(pt7, p1) ≈
-    -(min(LG.distance(pt7, r2), LG.distance(pt7, r3), LG.distance(pt7, r4)))
-# Point outside of polygon exterior
-@test GO.distance(pt5, p1) ≈ LG.distance(pt5, p1)
-@test GO.signed_distance(pt5, p1) ≈ LG.distance(pt5, p1)
-# Point inside of polygon hole
-@test GO.distance(pt8, p1) ≈ LG.distance(pt8, p1)
-@test GO.signed_distance(pt8, p1) ≈ LG.distance(pt8, p1)
-@test GO.distance(pt9, p1) ≈ LG.distance(pt9, p1)
-# Return types
-@test GO.distance(pt1, p1) isa Float64
-@test GO.distance(pt1, p1, Float32) isa Float32
-
-# Point and MultiPoint
-@test GO.distance(pt4, mpt1) == LG.distance(pt4, mpt1)
-@test GO.distance(pt4, mpt1) isa Float64
-@test GO.distance(pt4, mpt1, Float32) isa Float32
-
-# Point and MultiPolygon
-# Point outside of either polygon
-@test GO.distance(pt5, mp1) ≈ LG.distance(pt5, mp1)
-@test GO.distance(pt10, mp1) ≈ LG.distance(pt10, mp1)
-# Point within one polygon
-@test GO.distance(pt3, mp1) == LG.distance(pt3, mp1)
-@test GO.signed_distance(pt3, mp1) ≈
-    -(min(LG.distance(pt3, r2), LG.distance(pt3, r3), LG.distance(pt3, r4), LG.distance(pt3, r5)))
-@test GO.distance(pt11, mp1) == LG.distance(pt11, mp1)
-@test GO.signed_distance(pt11, mp1) ≈
-    -(min(LG.distance(pt11, r2), LG.distance(pt11, r3), LG.distance(pt11, r4), LG.distance(pt11, r5)))
-
-# Point and Geometry Collection
-@test GO.distance(pt1, c1) == LG.distance(pt1, c1)
-
+@testset_implementations "Where LinearRing exists" [LG, GI] 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
+
+@testset_implementations "Where GeometryCollection exists" [LG, AG, GI] begin
+    @test GO.distance($pt1, c1) == LG.distance($pt1, c1)
+end
+
+@testset_implementations "Point and Point" begin
+    # Distance from point to same point
+    @test GO.distance($pt1, $pt1) == LG.distance($pt1, $pt1)
+    # Distance from point to different point
+    @test GO.distance($pt1, $pt2) ≈ GO.distance($pt2, $pt1) ≈ LG.distance($pt1, $pt2)
+    # Return types
+    @test GO.distance($pt1, $pt1) isa Float64
+    @test GO.distance($pt1, $pt1, Float32) isa Float32
+end
+
+@testset_implementations "Point and Line" begin
+    #Point on line vertex
+    @test GO.distance($pt1, $l1) == GO.distance($l1, $pt1) == LG.distance($pt1, $l1)
+    # Point on line edge
+    @test GO.distance($pt2, $l1) == GO.distance($l1, $pt2) == LG.distance($pt2, $l1)
+    # Point equidistant from both segments
+    @test GO.distance($pt3, $l1) ≈  GO.distance($l1, $pt3) ≈  LG.distance($pt3, $l1)
+    # Point closer to one segment than another
+    @test GO.distance($pt4, $l1) ≈  GO.distance($l1, $pt4) ≈  LG.distance($pt4, $l1)
+    # Return types
+    @test GO.distance($pt1, $l1) isa Float64
+    @test GO.distance($pt1, $l1, Float32) isa Float32
+end
+
+@testset_implementations "Point and Polygon" begin
+    # Point on polygon exterior edge
+    @test GO.distance($pt1, $p1) == LG.distance($pt1, $p1)
+    @test GO.signed_distance($pt1, $p1) == 0
+    @test GO.distance($pt2, $p1) == LG.distance($pt2, $p1)
+    # Point on polygon hole edge
+    @test GO.distance($pt4, $p1) == LG.distance($pt4, $p1)
+    @test GO.signed_distance($pt4, $p1) == 0
+    @test GO.distance($pt6, $p1) == LG.distance($pt6, $p1)
+    # Point inside of polygon
+    @test GO.distance($pt3, $p1) == LG.distance($pt3, $p1)
+    @test GO.signed_distance($pt3, $p1) ≈
+        -(min(LG.distance($pt3, r2), LG.distance($pt3, r3), LG.distance($pt3, r4)))
+    @test GO.distance($pt7, $p1) == LG.distance($pt7, $p1)
+    @test GO.signed_distance($pt7, $p1) ≈
+        -(min(LG.distance($pt7, r2), LG.distance($pt7, r3), LG.distance($pt7, r4)))
+    # Point outside of polygon exterior
+    @test GO.distance($pt5, $p1) ≈ LG.distance($pt5, $p1)
+    @test GO.signed_distance($pt5, $p1) ≈ LG.distance($pt5, $p1)
+    # Point inside of polygon hole
+    @test GO.distance($pt8, $p1) ≈ LG.distance($pt8, $p1)
+    @test GO.signed_distance($pt8, $p1) ≈ LG.distance($pt8, $p1)
+    @test GO.distance($pt9, $p1) ≈ LG.distance($pt9, $p1)
+    # Return types
+    @test GO.distance($pt1, $p1) isa Float64
+    @test GO.distance($pt1, $p1, Float32) isa Float32
+end
+
+@testset_implementations "Point and MultiPoint" begin
+    @test GO.distance($pt4, $mpt1) == LG.distance(pt4, mpt1)
+    @test GO.distance($pt4, $mpt1) isa Float64
+    @test GO.distance($pt4, $mpt1, Float32) isa Float32
+end
+
+@testset_implementations "Point and MultiPolygon" begin
+    # Point outside of either polygon
+    @test GO.distance($pt5, mp1) ≈ LG.distance($pt5, mp1)
+    @test GO.distance($pt10, mp1) ≈ LG.distance($pt10, mp1)
+    # Point within one polygon
+    @test GO.distance($pt3, mp1) == LG.distance($pt3, mp1)
+    @test GO.signed_distance($pt3, mp1) ≈
+        -(min(LG.distance($pt3, r2), LG.distance($pt3, r3), LG.distance($pt3, r4), LG.distance($pt3, r5)))
+    @test GO.distance($pt11, mp1) == LG.distance($pt11, mp1)
+    @test GO.signed_distance($pt11, mp1) ≈
+        -(min(LG.distance($pt11, r2), LG.distance($pt11, r3), LG.distance($pt11, r4), LG.distance($pt11, r5)))
+end
diff --git a/test/methods/equals.jl b/test/methods/equals.jl
index 1fd725f36..38fa45c04 100644
--- a/test/methods/equals.jl
+++ b/test/methods/equals.jl
@@ -1,134 +1,149 @@
 using Test
-import GeoInterface as GI, GeometryOps as GO, LibGEOS as LG
+import GeoInterface as GI
+import GeometryBasics as GB 
+import GeometryOps as GO 
+import LibGEOS as LG
+using ..TestHelpers
 
-@testset "Points/MultiPoints" begin
-    p1 = LG.Point([0.0, 0.0])
-    p2 = LG.Point([0.0, 1.0])
+p1 = LG.Point([0.0, 0.0])
+p2 = LG.Point([0.0, 1.0])
+mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]])
+mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
+mp3 = LG.MultiPoint([p2])
+
+@testset_implementations "Points/MultiPoints" begin
     # Same points
-    @test GO.equals(p1, p1) == LG.equals(p1, p1)
-    @test GO.equals(p2, p2) == LG.equals(p2, p2)
+    @test GO.equals($p1, $p1) == LG.equals($p1, $p1)
+    @test GO.equals($p2, $p2) == LG.equals($p2, $p2)
     # Different points
-    @test GO.equals(p1, p2) == LG.equals(p1, p2)
+    @test GO.equals($p1, $p2) == LG.equals($p1, $p2)
 
-    mp1 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]])
-    mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
-    mp3 = LG.MultiPoint([p2])
     # Same points
-    @test GO.equals(mp1, mp1) == LG.equals(mp1, mp1)
-    @test GO.equals(mp2, mp2) == LG.equals(mp2, mp2)
+    @test GO.equals($mp1, $mp1) == LG.equals($mp1, $mp1)
+    @test GO.equals($mp2, $mp2) == LG.equals($mp2, $mp2)
     # Different points
-    @test GO.equals(mp1, mp2) == LG.equals(mp1, mp2)
-    @test GO.equals(mp1, p1) == LG.equals(mp1, p1)
+    @test GO.equals($mp1, $mp2) == LG.equals($mp1, $mp2)
+    @test GO.equals($mp1, $p1) == LG.equals($mp1, $p1)
     # Point and multipoint
-    @test GO.equals(p2, mp3) == LG.equals(p2, mp3)
+    @test GO.equals($p2, $mp3) == LG.equals($p2, $mp3)
 end
 
-@testset "Lines/Rings" begin
-    l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]])
-    l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]])
+l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]])
+l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]])
+l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]])
+r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]])
+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]])
+
+@testset_implementations "LineStrings" begin
     # Equal lines
-    @test GO.equals(l1, l1) == LG.equals(l1, l1)
-    @test GO.equals(l2, l2) == LG.equals(l2, l2)
+    @test GO.equals($l1, $l1) == LG.equals($l1, $l1)
+    @test GO.equals($l2, $l2) == LG.equals($l2, $l2)
     # Different lines
-    @test GO.equals(l1, l2) == GO.equals(l2, l1) == LG.equals(l1, l2)
+    @test GO.equals($l1, $l2) == GO.equals($l2, $l1) == LG.equals($l1, $l2)
 
-    r1 = LG.LinearRing([[0.0, 0.0], [5.0, 5.0], [10.0, 0.0], [5.0, -5.0], [0.0, 0.0]])
-    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]])
-    l3 = LG.LineString([[3.0, 0.0], [8.0, 5.0], [13.0, 0.0], [8.0, -5.0], [3.0, 0.0]])
     # 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)
+    @test GO.equals($r1, $r1) == LG.equals($r1, $r1)
+    @test GO.equals($r2, $r2) == LG.equals($r2, $r2)
     # Different rings
-    @test GO.equals(r1, r2) == GO.equals(r2, r1) == LG.equals(r1, r2)
+    @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)
+    @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)]))
+    @test GO.equals($l1, GI.Line([(0.0, 0.0), (0.0, 10.0)]))
 end
 
-@testset "Polygons/MultiPolygons" begin
-    pt1 = LG.Point([0.0, 0.0])
-    r1 = GI.LinearRing([(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)])
-    p1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]])
-    p2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]])
-    p3 = LG.Polygon(
-        [
-            [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-            [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
-        ]
-    )
-    p4 = LG.Polygon(
-        [
-            [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-            [[16.0, 1.0], [16.0, 11.0], [25.0, 11.0], [25.0, 1.0], [16.0, 1.0]]
-        ]
-    )
-    p5 = LG.Polygon(
-        [
-            [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-            [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]],
-            [[11.0, 1.0], [11.0, 2.0], [12.0, 2.0], [12.0, 1.0], [11.0, 1.0]]
-        ]
-    )
-    p6 = GI.Polygon([[(6, 6), (6, 1), (1, 1), (1, 6), (6, 6)]])
-    p7 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1), (6, 6)]])
-    p8 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1)]])
+# LibGEOS rejects rings that are not closed, and they are not eaual in GeometryBasics or ArchGDAL?
+@testset_implementations "Rings" [GI] begin
+    # Test equal rings without closing point
+    @test GO.equals($r2, $r3)
+    @test GO.equals($r3, $l3)
+end
+
+pt1 = LG.Point([0.0, 0.0])
+r1 = GI.LinearRing([(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)])
+p1 = GI.Polygon([[(0, 0), (0, 5), (5, 5), (5, 0), (0, 0)]])
+p2 = GI.Polygon([[(1, 1), (1, 6), (6, 6), (6, 1), (1, 1)]])
+p3 = LG.Polygon(
+    [
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
+    ]
+)
+p4 = LG.Polygon(
+    [
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[16.0, 1.0], [16.0, 11.0], [25.0, 11.0], [25.0, 1.0], [16.0, 1.0]]
+    ]
+)
+p5 = LG.Polygon(
+    [
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]],
+        [[11.0, 1.0], [11.0, 2.0], [12.0, 2.0], [12.0, 1.0], [11.0, 1.0]]
+    ]
+)
+p6 = GI.Polygon([[(6, 6), (6, 1), (1, 1), (1, 6), (6, 6)]])
+p7 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1), (6, 6)]])
+p8 = GI.Polygon([[(6, 6), (1, 6), (1, 1), (6, 1)]])
+p9 = LG.Polygon(
+    [[
+        [-53.57208251953125, 28.287451910503744],
+        [-53.33038330078125, 28.29228897739706],
+        [-53.34136962890625, 28.430052892335723],
+        [-53.57208251953125, 28.287451910503744],
+    ]]
+)
+m1 = LG.MultiPolygon([
+    [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]],
+    [
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
+    ]
+])
+m2 = LG.MultiPolygon([
+    [
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
+    ],
+    [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]]
+])
+m3 = LG.MultiPolygon([p3])
+
+@testset_implementations "Polygons" begin
     # Point and polygon aren't equal
-    GO.equals(pt1, p1) == LG.equals(pt1, p1)
+    GO.equals($pt1, $p1) == LG.equals($pt1, $p1)
     # Linear ring and polygon aren't equal
-    @test GO.equals(r1, p1) == LG.equals(r1, p1)
+    @test GO.equals($r1, $p1) == LG.equals($r1, $p1)
     # Equal polygon
-    @test GO.equals(p1, p1) == LG.equals(p1, p1)
-    @test GO.equals(p2, p2) == LG.equals(p2, p2)
+    @test GO.equals($p1, $p1) == LG.equals($p1, $p1)
+    @test GO.equals($p2, $p2) == LG.equals($p2, $p2)
     # Equal but offset polygons
-    @test GO.equals(p2, p6) == LG.equals(p2, p6)
+    @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) 
+    @test GO.equals($p2, $p7) == LG.equals($p2, $p7)
     # Different polygons
-    @test GO.equals(p1, p2) == LG.equals(p1, p2)
+    @test GO.equals($p1, $p2) == LG.equals($p1, $p2)
     # Equal polygons with holes
-    @test GO.equals(p3, p3) == LG.equals(p3, p3)
+    @test GO.equals($p3, $p3) == LG.equals($p3, $p3)
     # Same exterior, different hole
-    @test GO.equals(p3, p4) == LG.equals(p3, p4)
+    @test GO.equals($p3, $p4) == LG.equals($p3, $p4)
     # Same exterior and first hole, has an extra hole
-    @test GO.equals(p3, p5) == LG.equals(p3, p5)
-
-    p9 = LG.Polygon(
-        [[
-            [-53.57208251953125, 28.287451910503744],
-            [-53.33038330078125, 28.29228897739706],
-            [-53.34136962890625, 28.430052892335723],
-            [-53.57208251953125, 28.287451910503744],
-        ]]
-    )
+    @test GO.equals($p3, $p5) == LG.equals($p3, $p5)
     # Complex polygon
-    @test GO.equals(p9, p9) == LG.equals(p9, p9)
+    @test GO.equals($p9, $p9) == LG.equals($p9, $p9)
+end
 
-    m1 = LG.MultiPolygon([
-        [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]],
-        [
-            [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-            [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
-        ]
-    ])
-    m2 = LG.MultiPolygon([
-        [
-            [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-            [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
-        ],
-        [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]]
-    ])
+@testset_implementations "Unclosed Polygons" [GB, GI] begin
+    # Equal but without closing point (implied)
+    @test GO.equals($p7, $p8) 
+end
+
+@testset_implementations "MultiPolygons" begin
     # Equal multipolygon
-    @test GO.equals(m1, m1) == LG.equals(m1, m1)
+    @test GO.equals($m1, $m1) == LG.equals($m1, $m1)
     # Equal multipolygon with different order
-    @test GO.equals(m2, m2) == LG.equals(m2, m2)
+    @test GO.equals($m2, $m2) == LG.equals($m2, $m2)
     # Equal polygon to multipolygon
-    m3 = LG.MultiPolygon([p3])
-    @test GO.equals(p1, m3) == LG.equals(p1, m3)
-end
\ No newline at end of file
+    @test GO.equals($p1, $m3) == LG.equals($p1, $m3)
+end
diff --git a/test/methods/geom_relations.jl b/test/methods/geom_relations.jl
index b2f3f05df..0435c10fe 100644
--- a/test/methods/geom_relations.jl
+++ b/test/methods/geom_relations.jl
@@ -1,6 +1,7 @@
 import GeometryOps as GO
 import GeoInterface as GI
 import LibGEOS as LG
+using ..TestHelpers
 
 # Tests of DE-9IM Methods
 pt1 = LG.Point([0.0, 0.0])
@@ -154,19 +155,18 @@ 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"),
+    # (gc1, ml1, "gc1", "ml1", "Make sure collection works with multi-geom"),
 ]
 
-function test_geom_relation(GO_f, LG_f, f_name; swap_points = false)
+function test_geom_relation(GO_f, LG_f, f_name; swap_points=false)
     for (g1, g2, sg1, sg2, sdesc) in test_pairs
-        if swap_points
-            g1, g2 = g2, g1
-            sg1, sg2 = sg2, sg1
+        @testset "$sg1 $sg2 $sdesc" begin
+            if swap_points
+                @test_implementations (GO_f($g2, $g1) == LG_f($g2, $g1))
+            else
+                @test_implementations GO_f($g1, $g2) == LG_f($g1, $g2)
+            end
         end
-        go_val = GO_f(g1, g2)
-        lg_val = LG_f(g1, g2)
-        @test go_val == lg_val
-        go_val != lg_val && println("\n↑ TEST INFO: $sg1 $f_name $sg2 - $sdesc \n\n")
     end
 end
 
@@ -180,118 +180,124 @@ end
 
 
 @testset "Overlaps" begin
-    @testset "Points/MultiPoints" begin
-        p1 = LG.Point([0.0, 0.0])
-        p2 = LG.Point([0.0, 1.0])
+    p1 = LG.Point([0.0, 0.0])
+    p2 = LG.Point([0.0, 1.0])
+
+    mp1 = LG.MultiPoint([[0.0, 1.0], [4.0, 4.0]])
+    mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]])
+    mp3 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
+    mp4 = LG.MultiPoint([
+        [-36.05712890625, 26.480407161007275],
+        [-35.7220458984375, 27.137368359795584],
+        [-35.13427734375, 26.83387451505858],
+        [-35.4638671875, 27.254629577800063],
+        [-35.5462646484375, 26.86328062676624],
+        [-35.3924560546875, 26.504988828743404],
+    ])
+    mp5 = GI.MultiPoint([
+        [-35.4638671875, 27.254629577800063],
+        [-35.5462646484375, 26.86328062676624],
+        [-35.3924560546875, 26.504988828743404],
+        [-35.2001953125, 26.12091815959972],
+        [-34.9969482421875, 26.455820238459893],
+    ])
+
+    @testset_implementations "Points/MultiPoints"  begin
         # Two points can't overlap
-        @test GO.overlaps(p1, p1) == LG.overlaps(p1, p2)
-    
-        mp1 = LG.MultiPoint([[0.0, 1.0], [4.0, 4.0]])
-        mp2 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0]])
-        mp3 = LG.MultiPoint([[0.0, 1.0], [2.0, 2.0], [3.0, 3.0]])
+        @test GO.overlaps($p1, $p1) == LG.overlaps($p1, $p2)
         # No shared points, doesn't overlap
-        @test GO.overlaps(p1, mp1) == LG.overlaps(p1, mp1)
+        @test GO.overlaps($p1, $mp1) == LG.overlaps($p1, $mp1)
         # One shared point, does overlap
-        @test GO.overlaps(p2, mp1) == LG.overlaps(p2, mp1)
+        @test GO.overlaps($p2, $mp1) == LG.overlaps($p2, $mp1)
         # All shared points, doesn't overlap
-        @test GO.overlaps(mp1, mp1) == LG.overlaps(mp1, mp1)
+        @test GO.overlaps($mp1, $mp1) == LG.overlaps($mp1, $mp1)
         # Not all shared points, overlaps
-        @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2)
+        @test GO.overlaps($mp1, $mp2) == LG.overlaps($mp1, $mp2)
         # One set of points entirely inside other set, doesn't overlap
-        @test GO.overlaps(mp2, mp3) == LG.overlaps(mp2, mp3)
+        @test GO.overlaps($mp2, $mp3) == LG.overlaps($mp2, $mp3)
         # Not all points shared, overlaps
-        @test GO.overlaps(mp1, mp3) == LG.overlaps(mp1, mp3)
-    
-        mp1 = LG.MultiPoint([
-            [-36.05712890625, 26.480407161007275],
-            [-35.7220458984375, 27.137368359795584],
-            [-35.13427734375, 26.83387451505858],
-            [-35.4638671875, 27.254629577800063],
-            [-35.5462646484375, 26.86328062676624],
-            [-35.3924560546875, 26.504988828743404],
-        ])
-        mp2 = GI.MultiPoint([
-            [-35.4638671875, 27.254629577800063],
-            [-35.5462646484375, 26.86328062676624],
-            [-35.3924560546875, 26.504988828743404],
-            [-35.2001953125, 26.12091815959972],
-            [-34.9969482421875, 26.455820238459893],
-        ])
+        @test GO.overlaps($mp1, $mp3) == LG.overlaps($mp1, $mp3)
         # Some shared points, overlaps
-        @test GO.overlaps(mp1, mp2) == LG.overlaps(mp1, mp2)
-        @test GO.overlaps(mp1, mp2) == GO.overlaps(mp2, mp1)
+        @test GO.overlaps($mp1, $mp2) == LG.overlaps($mp1, $mp2)
+        @test GO.overlaps($mp1, $mp2) == GO.overlaps($mp2, $mp1)
     end
     
-    @testset "Lines/Rings" begin
-        l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]])
-        l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]])
-        l3 = LG.LineString([[0.0, -10.0], [0.0, 3.0]])
-        l4 = LG.LineString([[5.0, -5.0], [5.0, 5.0]])
+    l1 = LG.LineString([[0.0, 0.0], [0.0, 10.0]])
+    l2 = LG.LineString([[0.0, -10.0], [0.0, 20.0]])
+    l3 = LG.LineString([[0.0, -10.0], [0.0, 3.0]])
+    l4 = LG.LineString([[5.0, -5.0], [5.0, 5.0]])
+    # Linear rings that intersect but don't overlap
+    r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]])
+    r2 = LG.LinearRing([[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]])
+
+    @testset_implementations "Lines/Rings" begin
         # Line can't overlap with itself
-        @test GO.overlaps(l1, l1) == LG.overlaps(l1, l1)
+        @test GO.overlaps($l1, $l1) == LG.overlaps($l1, $l1)
         # Line completely within other line doesn't overlap
-        @test GO.overlaps(l1, l2) == GO.overlaps(l2, l1) == LG.overlaps(l1, l2)
+        @test GO.overlaps($l1, $l2) == GO.overlaps($l2, $l1) == LG.overlaps($l1, $l2)
         # Overlapping lines
-        @test GO.overlaps(l1, l3) == GO.overlaps(l3, l1) == LG.overlaps(l1, l3)
+        @test GO.overlaps($l1, $l3) == GO.overlaps($l3, $l1) == LG.overlaps($l1, $l3)
         # Lines that don't touch
-        @test GO.overlaps(l1, l4) == LG.overlaps(l1, l4)
-        # Linear rings that intersect but don't overlap
-        r1 = LG.LinearRing([[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]])
-        r2 = LG.LinearRing([[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]])
-        @test GO.overlaps(r1, r2) == LG.overlaps(r1, r2)
+        @test GO.overlaps($l1, $l4) == LG.overlaps($l1, $l4)
+        @test LG.overlaps($r1, $r2) == LG.overlaps($r1, $r2)
     end
     
-    @testset "Polygons/MultiPolygons" begin
-        p1 = LG.Polygon([[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]])
-        p2 = LG.Polygon([
+    p1 = LG.Polygon([[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]])
+    p2 = LG.Polygon([
+        [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
+        [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
+    ])
+    p3 = LG.Polygon([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]])
+    p4 = LG.Polygon([[[20.0, 5.0], [20.0, 10.0], [18.0, 10.0], [18.0, 5.0], [20.0, 5.0]]])
+    p5 = LG.Polygon(
+        [[
+            [-53.57208251953125, 28.287451910503744],
+            [-53.33038330078125, 28.29228897739706],
+            [-53.34136352890625, 28.430052892335723],
+            [-53.57208251953125, 28.287451910503744],
+        ]]
+    )
+    m1 = LG.MultiPolygon([
+        [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]],
+        [
             [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
             [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
-        ])
+        ]
+    ])
+
+    @testset_implementations "Polygons/MultiPolygons" begin
         # Test basic polygons that don't overlap
-        @test GO.overlaps(p1, p2) == LG.overlaps(p1, p2)
-        @test !GO.overlaps(p1, (1, 1))
-        @test !GO.overlaps((1, 1), p2)
+        @test GO.overlaps($p1, $p2) == LG.overlaps($p1, $p2)
+        @test !GO.overlaps($p1, (1, 1))
+        @test !GO.overlaps((1, 1), $p2)
     
-        p3 = LG.Polygon([[[1.0, 1.0], [1.0, 6.0], [6.0, 6.0], [6.0, 1.0], [1.0, 1.0]]])
         # Test basic polygons that overlap
-        @test GO.overlaps(p1, p3) == LG.overlaps(p1, p3)
+        @test GO.overlaps($p1, $p3) == LG.overlaps($p1, $p3)
     
-        p4 = LG.Polygon([[[20.0, 5.0], [20.0, 10.0], [18.0, 10.0], [18.0, 5.0], [20.0, 5.0]]])
         # Test one polygon within the other
-        @test GO.overlaps(p2, p4) == GO.overlaps(p4, p2) == LG.overlaps(p2, p4)
+        @test GO.overlaps($p2, $p4) == GO.overlaps($p4, $p2) == LG.overlaps($p2, $p4)
     
-        p5 = LG.Polygon(
-            [[
-                [-53.57208251953125, 28.287451910503744],
-                [-53.33038330078125, 28.29228897739706],
-                [-53.34136352890625, 28.430052892335723],
-                [-53.57208251953125, 28.287451910503744],
-            ]]
-        )
         # Test equal polygons
-        @test GO.overlaps(p5, p5) == LG.overlaps(p5, p5)
+        @test GO.overlaps($p5, $p5) == LG.overlaps($p5, $p5)
     
-        # Test multipolygons
-        m1 = LG.MultiPolygon([
-            [[[0.0, 0.0], [0.0, 5.0], [5.0, 5.0], [5.0, 0.0], [0.0, 0.0]]],
-            [
-                [[10.0, 0.0], [10.0, 20.0], [30.0, 20.0], [30.0, 0.0], [10.0, 0.0]],
-                [[15.0, 1.0], [15.0, 11.0], [25.0, 11.0], [25.0, 1.0], [15.0, 1.0]]
-            ]
-        ])
         # Test polygon that overlaps with multipolygon
-        @test GO.overlaps(m1, p3) == LG.overlaps(m1, p3)
+        @test GO.overlaps($m1, $p3) == LG.overlaps($m1, $p3)
         # Test polygon in hole of multipolygon, doesn't overlap
-        @test GO.overlaps(m1, p4) == LG.overlaps(m1, p4)
+        @test GO.overlaps($m1, $p4) == LG.overlaps($m1, $p4)
     end
 end
+
 @testset "Crosses" begin
 	line6 = GI.LineString([(1.0, 1.0), (1.0, 2.0), (1.0, 3.0), (1.0, 4.0)])
+    line7 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
+    line8 = GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])
 	poly7 = GI.Polygon([[(-1.0, 2.0), (3.0, 2.0), (3.0, 3.0), (-1.0, 3.0), (-1.0, 2.0)]])
 
-	@test GO.crosses(GI.LineString([(-2.0, 2.0), (4.0, 2.0)]), line6) == true
-	@test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), poly7) == true
-	@test GO.crosses(GI.MultiPoint([(1.0, 2.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == true
-	@test GO.crosses(GI.MultiPoint([(1.0, 0.0), (12.0, 12.0)]), GI.LineString([(1, 1), (1, 2), (1, 3), (1, 4)])) == false
-	@test GO.crosses(GI.LineString([(-2.0, 2.0), (-4.0, 2.0)]), poly7) == false
-end
\ No newline at end of file
+    @testset_implementations begin
+        @test GO.crosses(GI.LineString([(-2.0, 2.0), (4.0, 2.0)]), $line6) == true
+        @test GO.crosses(GI.LineString([(0.5, 2.5), (1.0, 1.0)]), $poly7) == true
+        @test GO.crosses(GI.MultiPoint([(1.0, 2.0), (12.0, 12.0)]), $line7) == true
+        @test GO.crosses(GI.MultiPoint([(1.0, 0.0), (12.0, 12.0)]), $line8) == false
+        @test GO.crosses(GI.LineString([(-2.0, 2.0), (-4.0, 2.0)]), $poly7) == false
+    end
+end
diff --git a/test/methods/orientation.jl b/test/methods/orientation.jl
index cc5c7b5c6..fef53ebd4 100644
--- a/test/methods/orientation.jl
+++ b/test/methods/orientation.jl
@@ -1,34 +1,36 @@
 using Test, GeometryOps
-import GeoInterface as GI, GeometryOps as GO
+import GeoInterface as GI
+import GeometryOps as GO
+using ..TestHelpers
 
-@testset "Orientation" begin
-    line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]])
-    line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]])
-    line3 = GI.LineString([
-        (-111.544189453125, 24.186847428521244),
-        (-110.687255859375, 24.966140159912975),
-        (-110.4510498046875, 24.467150664739002),
-        (-109.9951171875, 25.180087808990645)
-    ])
-	line4 = GI.LineString([
-        (-111.4617919921875, 24.05148034322011),
-        (-110.8795166015625, 24.681961205014595),
-        (-110.841064453125, 24.14174098050432),
-        (-109.97863769531249, 24.617057340809524)
-    ])
+line1 = GI.LineString([[9.170356, 45.477985], [9.164434, 45.482551], [9.166644, 45.484003]])
+line2 = GI.LineString([[9.169356, 45.477985], [9.163434, 45.482551], [9.165644, 45.484003]])
+line3 = GI.LineString([
+    (-111.544189453125, 24.186847428521244),
+    (-110.687255859375, 24.966140159912975),
+    (-110.4510498046875, 24.467150664739002),
+    (-109.9951171875, 25.180087808990645)
+])
+line4 = GI.LineString([
+    (-111.4617919921875, 24.05148034322011),
+    (-110.8795166015625, 24.681961205014595),
+    (-110.841064453125, 24.14174098050432),
+    (-109.97863769531249, 24.617057340809524)
+])
 
-    # @test isparallel(line1, line2) == true
-	# @test isparallel(line3, line4) == false
-
-	poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]])
-	poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]])
+poly1 = GI.Polygon([[[0, 0], [1, 0], [1, 1], [0.5, 0.5], [0, 1], [0, 0]]])
+poly2 = GI.Polygon([[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]])
 
-	@test isconcave(poly1) == true
-	@test isconcave(poly2) == false
+l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]])
+l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]])
 
-	l1 = GI.LineString([[0, 0], [1, 1], [1, 0], [0, 0]])
-	l2 = GI.LineString([[0, 0], [1, 0], [1, 1], [0, 0]])
+@testset_implementations "Orientation" begin
+    # @test isparallel(line1, line2) == true
+    # @test isparallel(line3, line4) == false
+    
+    @test isconcave($poly1) == true
+    @test isconcave($poly2) == false
 
-	@test isclockwise(l1) == true
-	@test isclockwise(l2) == false
+    @test isclockwise($l1) == true
+    @test isclockwise($l2) == false
 end
diff --git a/test/methods/polygonize.jl b/test/methods/polygonize.jl
index 6927f01d2..d16a55c0a 100644
--- a/test/methods/polygonize.jl
+++ b/test/methods/polygonize.jl
@@ -1,5 +1,6 @@
 using GeometryOps, GeoInterface, Test
 import OffsetArrays, DimensionalData, Rasters
+using ..TestHelpers
 
 # Missing holes throw a warning, so testing there are
 # no warnings in a range of randomisation is one way to test 
@@ -66,4 +67,4 @@ end
         @test GO.equals(data_mp, evil_mp)
         @test GI.crs(evil_mp) == GI.crs(evil)
     end
-end
\ No newline at end of file
+end
diff --git a/test/primitives.jl b/test/primitives.jl
index 6661a4eea..837916683 100644
--- a/test/primitives.jl
+++ b/test/primitives.jl
@@ -1,9 +1,15 @@
 using Test
 
-import GeoInterface as GI, GeometryOps as GO, GeometryBasics as GB
+import ArchGDAL as AG
+import GeometryBasics as GB 
+import GeometryOps as GO 
+import GeoInterface as GI
+import LibGEOS as LG
 import Proj
-import Shapefile, DataFrames
+import Shapefile
+import DataFrames
 using Downloads: download
+using ..TestHelpers
 
 pv1 = [(1, 2), (3, 4), (5, 6), (1, 2)]
 pv2 = [(3, 4), (5, 6), (6, 7), (3, 4)]
@@ -13,13 +19,14 @@ poly = GI.Polygon([lr1, lr2])
 
 @testset "apply" begin
 
-    flipped_poly = GO.apply(GI.PointTrait, poly) do p
-        (GI.y(p), GI.x(p))
+    @testset_implementations "apply flip" begin
+        flipped_poly = GO.apply(GI.PointTrait, $poly) do p
+            (GI.y(p), GI.x(p))
+        end
+        @test flipped_poly == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
+                                          GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
     end
 
-    @test flipped_poly == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
-                                      GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
-
     @testset "Tables.jl support" begin
         # check to account for missing data
         missing_or_equal(x, y) = (ismissing(x) && ismissing(y)) || (x == y)
@@ -73,9 +80,8 @@ end
 
 @testset "flatten" begin
     very_wrapped = [[GI.FeatureCollection([GI.Feature(poly; properties=(;))])]]
-    @test collect(GO.flatten(GI.PointTrait, very_wrapped)) == vcat(pv1, pv2)
-    @test collect(GO.flatten(GI.LinearRingTrait, [poly])) == [lr1, lr2]
-    @test collect(GO.flatten(GI.LinearRingTrait, [poly])) == [lr1, lr2]
+    @test GO._tuple_point.(GO.flatten(GI.PointTrait, very_wrapped)) == vcat(pv1, pv2)
+    @test collect(GO.flatten(GI.AbstractCurveTrait, [poly])) == [lr1, lr2]
     @test collect(GO.flatten(GI.x, GI.PointTrait, very_wrapped)) == first.(vcat(pv1, pv2))
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index c671a87b7..d4d25637a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,49 +1,35 @@
 using GeometryOps
 using Test
+using SafeTestsets
 
-using GeometryOps.GeoInterface
-using GeometryOps.GeometryBasics
-using GeoInterface.Extents: Extents
-using ArchGDAL
-using LibGEOS
-using Random, Distributions
-using Proj
+include("helpers.jl")
 
-const GI = GeoInterface
-const AG = ArchGDAL
-const LG = LibGEOS
-const GO = GeometryOps
-
-@testset "GeometryOps.jl" begin
-    @testset "Primitives" begin include("primitives.jl") end
-    # Methods
-    @testset "Angles" begin include("methods/angles.jl") end
-    @testset "Area" begin include("methods/area.jl") end
-    @testset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
-    @testset "Orientation" begin include("methods/orientation.jl") end
-    @testset "Centroid" begin include("methods/centroid.jl") end
-    @testset "Convex Hull" begin include("methods/convex_hull.jl") end
-    @testset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
-    @testset "Distance" begin include("methods/distance.jl") end
-    @testset "Equals" begin include("methods/equals.jl") end
-    # Clipping
-    @testset "Coverage" begin include("methods/clipping/coverage.jl") end
-    @testset "Cut" begin include("methods/clipping/cut.jl") end
-    @testset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end
-    @testset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end
-    # Transformations
-    @testset "Embed Extent" begin include("transformations/extent.jl") end
-    @testset "Reproject" begin include("transformations/reproject.jl") end
-    @testset "Flip" begin include("transformations/flip.jl") end
-    @testset "Simplify" begin include("transformations/simplify.jl") end
-    @testset "Segmentize" begin include("transformations/segmentize.jl") end
-    @testset "Transform" begin include("transformations/transform.jl") end
-    @testset "Geometry correction" begin 
-        include("transformations/correction/geometry_correction.jl")
-        include("transformations/correction/closed_ring.jl") 
-        include("transformations/correction/intersecting_polygons.jl")
-    end
-    # Extensions
-    @testset "FlexiJoins" begin include("extensions/flexijoins.jl") end
-    @testset "LibGEOS" begin include("extensions/libgeos.jl") end
-end
+@safetestset "Primitives" begin include("primitives.jl") end
+# Methods
+@safetestset "Angles" begin include("methods/angles.jl") end
+@safetestset "Area" begin include("methods/area.jl") end
+# @safetestset "Barycentric coordinate operations" begin include("methods/barycentric.jl") end
+@safetestset "Orientation" begin include("methods/orientation.jl") end
+@safetestset "Centroid" begin include("methods/centroid.jl") end
+@safetestset "Convex Hull" begin include("methods/convex_hull.jl") end
+@safetestset "DE-9IM Geom Relations" begin include("methods/geom_relations.jl") end
+@safetestset "Distance" begin include("methods/distance.jl") end
+@safetestset "Equals" begin include("methods/equals.jl") end
+# Clipping
+@safetestset "Coverage" begin include("methods/clipping/coverage.jl") end
+@safetestset "Cut" begin include("methods/clipping/cut.jl") end
+@safetestset "Intersection Point" begin include("methods/clipping/intersection_points.jl") end
+@safetestset "Polygon Clipping" begin include("methods/clipping/polygon_clipping.jl") end
+# Transformations
+@safetestset "Embed Extent" begin include("transformations/extent.jl") end
+@safetestset "Reproject" begin include("transformations/reproject.jl") end
+@safetestset "Flip" begin include("transformations/flip.jl") end
+@safetestset "Simplify" begin include("transformations/simplify.jl") end
+@safetestset "Segmentize" begin include("transformations/segmentize.jl") end
+@safetestset "Transform" begin include("transformations/transform.jl") end
+@safetestset "Geometry Correction" begin include("transformations/correction/geometry_correction.jl") end
+@safetestset "Closed Rings" begin include("transformations/correction/closed_ring.jl")  end
+@safetestset "Intersecting Polygons" begin include("transformations/correction/intersecting_polygons.jl") end
+# Extensions
+@safetestset "FlexiJoins" begin include("extensions/flexijoins.jl") end
+@safetestset "LibGEOS" begin include("extensions/libgeos.jl") end
diff --git a/test/transformations/correction/closed_ring.jl b/test/transformations/correction/closed_ring.jl
index 939d88c0f..90819f40f 100644
--- a/test/transformations/correction/closed_ring.jl
+++ b/test/transformations/correction/closed_ring.jl
@@ -1,12 +1,16 @@
 using Test
+import ArchGDAL as AG
+import GeoInterface as GI
+import GeometryBasics as GB
+import GeometryOps as GO
+using ..TestHelpers
 
-import GeoInterface as GI, GeometryOps as GO
+open_rectangle = GI.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])])
 
-open_rectangle = GI.Wrappers.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])])
-
-@testset "Closed Ring correction" begin
-    closed_rectangle = GO.ClosedRing()(open_rectangle)
-    @test GI.npoint(closed_rectangle) == GI.npoint(open_rectangle) + 1 # test that the rectangle is closed
+# LibGEOS fails in GI.convert if the ring is not closed
+@testset_implementations "Closed Ring correction" [GI, AG, GB] begin
+    closed_rectangle = GO.ClosedRing()($open_rectangle)
+    @test GI.npoint(closed_rectangle) == GI.npoint($open_rectangle) + 1 # test that the rectangle is closed
     @test GI.getpoint(closed_rectangle.geom[1], 1) == GI.getpoint(closed_rectangle.geom[1], GI.npoint(closed_rectangle))
     @test all(GO.flatten(GI.PointTrait, closed_rectangle) .== GO.flatten(GI.PointTrait, GO.ClosedRing()(closed_rectangle)))
-end
\ No newline at end of file
+end
diff --git a/test/transformations/correction/geometry_correction.jl b/test/transformations/correction/geometry_correction.jl
index e2a1a48ac..de8d51a09 100644
--- a/test/transformations/correction/geometry_correction.jl
+++ b/test/transformations/correction/geometry_correction.jl
@@ -1,12 +1,16 @@
 using Test
+import ArchGDAL as AG
+import GeoInterface as GI
+import GeometryBasics as GB
+import GeometryOps as GO
+using ..TestHelpers
 
-import GeoInterface as GI, GeometryOps as GO
+open_rectangle = GI.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])])
 
-open_rectangle = GI.Wrappers.Polygon([collect.([(0, 0), (10, 0), (10, 10), (0, 10)])])
-
-@testset "Fix method" begin
-    closed_rectangle = GO.fix(open_rectangle)
-    @test GI.npoint(closed_rectangle) == GI.npoint(open_rectangle) + 1 # test that the rectangle is closed
+# LibGEOS fails in GI.convert if the ring is not closed
+@testset_implementations "fix open rectangle" [GI, AG, GB] begin
+    closed_rectangle = GO.fix($open_rectangle)
+    @test GI.npoint(closed_rectangle) == GI.npoint($open_rectangle) + 1 # test that the rectangle is closed
     @test GI.getpoint(closed_rectangle.geom[1], 1) == GI.getpoint(closed_rectangle.geom[1], GI.npoint(closed_rectangle))
     @test all(GO.flatten(GI.PointTrait, closed_rectangle) .== GO.flatten(GI.PointTrait, GO.fix(closed_rectangle)))
-end
\ No newline at end of file
+end
diff --git a/test/transformations/correction/intersecting_polygons.jl b/test/transformations/correction/intersecting_polygons.jl
index 94cb31261..dfac723bf 100644
--- a/test/transformations/correction/intersecting_polygons.jl
+++ b/test/transformations/correction/intersecting_polygons.jl
@@ -1,7 +1,7 @@
 using Test
-
 import GeoInterface as GI
 import GeometryOps as GO
+using ..TestHelpers
 
 p1 = GI.Polygon([[(0.0, 0.0), (3.0, 0.0), (3.0, 3.0), (0.0, 3.0), (0.0, 0.0)]])
 # (p1, p2) -> one polygon inside of the other, sharing an edge
@@ -25,55 +25,61 @@ mp6 = GI.MultiPolygon([  # four interlocking polygons forming a hole
     [[(-5.0, -5.0), (-5.0, 15.0), (0.0, 15.0), (0.0, -5.0), (-5.0, -5.0)]],
 ])
 
-union_fixed_mp1 = GO.UnionIntersectingPolygons()(mp1)
-@test GI.npolygon(union_fixed_mp1) == 1
-@test GO.equals(GI.getpolygon(union_fixed_mp1, 1), p1)
+@testset_implementations begin
+    union_fixed_mp1 = GO.UnionIntersectingPolygons()($mp1)
+    @test GI.npolygon(union_fixed_mp1) == 1
+    @test GO.equals(GI.getpolygon(union_fixed_mp1, 1), $p1)
 
-diff_fixed_mp1 = GO.DiffIntersectingPolygons()(mp1)
-@test GO.equals(diff_fixed_mp1, union_fixed_mp1)
+    diff_fixed_mp1 = GO.DiffIntersectingPolygons()($mp1)
+    @test GO.equals(diff_fixed_mp1, union_fixed_mp1)
 
-union_fixed_mp2 = GO.UnionIntersectingPolygons()(mp2)
-@test GI.npolygon(union_fixed_mp2) == 1
-@test GO.equals(GI.getpolygon(union_fixed_mp2, 1), p1)
+    union_fixed_mp2 = GO.UnionIntersectingPolygons()($mp2)
+    @test GI.npolygon(union_fixed_mp2) == 1
+    @test GO.equals(GI.getpolygon(union_fixed_mp2, 1), $p1)
 
-diff_fixed_mp2 = GO.DiffIntersectingPolygons()(mp2)
-@test GO.equals(diff_fixed_mp2, union_fixed_mp2)
+    diff_fixed_mp2 = GO.DiffIntersectingPolygons()($mp2)
+    @test GO.equals(diff_fixed_mp2, union_fixed_mp2)
 
-union_fixed_mp3 = GO.UnionIntersectingPolygons()(mp3)
-@test GI.npolygon(union_fixed_mp3) == 1
-@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3)))
-diff_polys = GO.difference(union_fixed_mp3, mp3; target = GI.PolygonTrait(), fix_multipoly = nothing)
-@test sum(GO.area, diff_polys; init = 0.0) == 0
+    union_fixed_mp3 = GO.UnionIntersectingPolygons()($mp3)
+    @test GI.npolygon(union_fixed_mp3) == 1
+    @test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(mp3)))
+    diff_polys = GO.difference(union_fixed_mp3, $mp3; 
+        target = GI.PolygonTrait(), fix_multipoly = nothing
+    )
+    @test sum(GO.area, diff_polys; init = 0.0) == 0
 
-diff_fixed_mp3 = GO.DiffIntersectingPolygons()(mp3)
-@test GI.npolygon(diff_fixed_mp3) == 3
-@test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3)))
+    diff_fixed_mp3 = GO.DiffIntersectingPolygons()($mp3)
+    @test GI.npolygon(diff_fixed_mp3) == 3
+    @test all((GO.coveredby(p, union_fixed_mp3) for p in GI.getpolygon(diff_fixed_mp3)))
 
-union_fixed_mp4 = GO.UnionIntersectingPolygons()(mp4)
-@test GI.npolygon(union_fixed_mp4) == 2
-@test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), p4)) ||
-    (GO.equals(GI.getpolygon(union_fixed_mp4, 2), p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), p4))
+    union_fixed_mp4 = GO.UnionIntersectingPolygons()($mp4)
+    @test GI.npolygon(union_fixed_mp4) == 2
+    @test (GO.equals(GI.getpolygon(union_fixed_mp4, 1), $p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 2), $p4)) ||
+        (GO.equals(GI.getpolygon(union_fixed_mp4, 2), $p1) && GO.equals(GI.getpolygon(union_fixed_mp4, 1), $p4))
 
-diff_fixed_mp4 = GO.DiffIntersectingPolygons()(mp4)
-@test GO.equals(diff_fixed_mp4, union_fixed_mp4)
+    diff_fixed_mp4 = GO.DiffIntersectingPolygons()($mp4)
+    @test GO.equals(diff_fixed_mp4, union_fixed_mp4)
 
-union_fixed_mp5 = GO.UnionIntersectingPolygons()(mp5)
-@test GI.npolygon(union_fixed_mp5) == 1
-@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(mp5)))
-diff_polys = GO.difference(union_fixed_mp5, mp5; target = GI.PolygonTrait(), fix_multipoly = nothing)
-@test sum(GO.area, diff_polys; init = 0.0) == 0
+    union_fixed_mp5 = GO.UnionIntersectingPolygons()($mp5)
+    @test GI.npolygon(union_fixed_mp5) == 1
+    @test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon($mp5)))
+    diff_polys = GO.difference(union_fixed_mp5, $mp5; target = GI.PolygonTrait(), fix_multipoly = nothing)
+    @test sum(GO.area, diff_polys; init = 0.0) == 0
 
-diff_fixed_mp5 = GO.DiffIntersectingPolygons()(mp5)
-@test GI.npolygon(diff_fixed_mp5) == 3
-@test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5)))
+    diff_fixed_mp5 = GO.DiffIntersectingPolygons()($mp5)
+    @test GI.npolygon(diff_fixed_mp5) == 3
+    @test all((GO.coveredby(p, union_fixed_mp5) for p in GI.getpolygon(diff_fixed_mp5)))
 
-union_fixed_mp6 = GO.UnionIntersectingPolygons()(mp6)
-@test GI.npolygon(union_fixed_mp6) == 1
-@test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1
-@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6)))
-diff_polys = GO.difference(union_fixed_mp6, mp6; target = GI.PolygonTrait(), fix_multipoly = nothing)
-@test sum(GO.area, diff_polys; init = 0.0) == 0
+    union_fixed_mp6 = GO.UnionIntersectingPolygons()($mp6)
+    @test GI.npolygon(union_fixed_mp6) == 1
+    @test GI.nhole(GI.getpolygon(union_fixed_mp6, 1)) == 1
+    @test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(mp6)))
+    diff_polys = GO.difference(union_fixed_mp6, $mp6; 
+        target = GI.PolygonTrait(), fix_multipoly = nothing
+    )
+    @test sum(GO.area, diff_polys; init = 0.0) == 0
 
-diff_fixed_mp6 = GO.DiffIntersectingPolygons()(mp6)
-@test GI.npolygon(diff_fixed_mp6) == 4
-@test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6)))
+    diff_fixed_mp6 = GO.DiffIntersectingPolygons()($mp6)
+    @test GI.npolygon(diff_fixed_mp6) == 4
+    @test all((GO.coveredby(p, union_fixed_mp6) for p in GI.getpolygon(diff_fixed_mp6)))
+end
diff --git a/test/transformations/extent.jl b/test/transformations/extent.jl
index e93390347..87026e9eb 100644
--- a/test/transformations/extent.jl
+++ b/test/transformations/extent.jl
@@ -1,15 +1,17 @@
 using Test
- 
-import GeoInterface as GI, GeometryOps as GO
+import GeoInterface as GI
+import GeometryOps as GO
 using GeoInterface: Extents
+using ..TestHelpers
 
-@testset "embed_extent" begin
-    poly = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
-                       GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
+poly = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
+                   GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
 
-    ext_poly = GO.embed_extent(poly)
+@testset_implementations "embed_extent" begin
+    ext_poly = GO.embed_extent($poly)
     lr1, lr2 = GI.getgeom(ext_poly)
     @test ext_poly.extent == Extents.Extent(X=(1, 6), Y=(2, 7))
     @test lr1.extent == Extents.Extent(X=(1, 5), Y=(2, 6))
     @test lr2.extent == Extents.Extent(X=(3, 6), Y=(4, 7))
 end
+
diff --git a/test/transformations/flip.jl b/test/transformations/flip.jl
index e977c7504..f02083ff5 100644
--- a/test/transformations/flip.jl
+++ b/test/transformations/flip.jl
@@ -1,12 +1,11 @@
 using Test
- 
-import GeoInterface as GI, GeometryOps as GO
+import GeoInterface as GI
+import GeometryOps as GO
+using ..TestHelpers
 
-@testset "flip" begin
-    geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
-                       GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
-
-
-    @test GO.flip(geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
-                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
+geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
+                   GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
+@testset_implementations "flip" begin
+    @test GO.flip($geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
+                                        GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
 end
diff --git a/test/transformations/reproject.jl b/test/transformations/reproject.jl
index ac7d22faf..52044a6da 100644
--- a/test/transformations/reproject.jl
+++ b/test/transformations/reproject.jl
@@ -1,38 +1,43 @@
 using Test
- 
-import GeoInterface as GI, GeometryOps as GO
 using GeoFormatTypes
+import GeoInterface as GI
+import GeometryOps as GO
 import Proj
+using ..TestHelpers
 
-@testset "reproject" begin
-    ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
-    ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
-    hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
+ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
+ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
+hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
 
-    # Set up a regular tranformation of the points for reference
-    source_crs = convert(Proj.CRS, EPSG(4326))
-    target_crs = convert(Proj.CRS, EPSG(3857))
-    trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
+# Set up a regular tranformation of the points for reference
+source_crs = convert(Proj.CRS, EPSG(4326))
+target_crs = convert(Proj.CRS, EPSG(3857))
+trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
 
-    polygon1 = GI.Polygon([ring1])
-    polygon2 = GI.Polygon([ring2, hole2])
-    multipolygon = GI.MultiPolygon([polygon1, polygon2])
+polygon1 = GI.Polygon([ring1])
+polygon2 = GI.Polygon([ring2, hole2])
+multipolygon = GI.MultiPolygon([polygon1, polygon2])
 
-    ref_points3857 = map(GI.getpoint(multipolygon)) do p
-        trans([GI.x(p), GI.y(p)])
-    end
+ref_points3857 = map(GI.getpoint(multipolygon)) do p
+    trans([GI.x(p), GI.y(p)])
+end
+
+_xy(p) = GI.x(p), GI.y(p)
 
-    multipolygon3857 = GO.reproject(multipolygon, EPSG(4326), EPSG(3857))
-    multipolygon4326 = GO.reproject(multipolygon3857; target_crs=EPSG(4326))
-    points4326_1 = collect(GI.getpoint(multipolygon))
-    points4326_2 = GI.getcoord.(GI.getpoint(multipolygon4326))
-    points3857 = GI.getcoord.(GI.getpoint(multipolygon3857))
+@testset_implementations "reproject" begin
+    multipolygon3857 = GO.reproject($multipolygon, EPSG(4326), EPSG(3857))
+    multipolygon4326 = GO.reproject($multipolygon; source_crs=EPSG(4326), target_crs=EPSG(4326))
+    points4326_1 = collect(GI.getpoint($multipolygon))
+    points4326_2 = collect.(GI.getcoord.(GI.getpoint(multipolygon4326)))
+    points3857 = collect.(GI.getcoord.(GI.getpoint(multipolygon3857)))
 
     # Comparison to regular `trans` on points
-    @test all(map((p1, p2) -> all(map(isapprox, p1, p2)), ref_points3857, points3857))
+    @test map(ref_points3857, points3857) do p1, p2
+        all(map(isapprox, _xy(p1), _xy(p2))) 
+    end |> all
 
     # Round trip comparison
-    @test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points4326_1, points4326_2))
+    @test all(map((p1, p2) -> all(map(isapprox, _xy(p1), _xy(p2))), points4326_1, points4326_2))
 
     # Embedded crs check
     @test GI.crs(multipolygon3857) == EPSG(3857)
@@ -73,6 +78,5 @@ import Proj
         GO.reproject(multipolygon4326; target_crs=utm32_wkt)
 
     GO.reproject(multipolygon4326; target_crs=ProjString("+proj=moll"))
-
 end
 
diff --git a/test/transformations/segmentize.jl b/test/transformations/segmentize.jl
index 3a5d18333..5095c4c84 100644
--- a/test/transformations/segmentize.jl
+++ b/test/transformations/segmentize.jl
@@ -1,6 +1,8 @@
 using Test
-
-import GeometryOps as GO, GeoInterface as GI
+using Proj
+import GeometryOps as GO
+import GeoInterface as GI
+using ..TestHelpers
 
 @testset "Segmentation on multiple geometry levels" begin
     ls = GI.LineString([(0, 0), (1, 1), (2, 2), (3, 3)])
@@ -8,60 +10,60 @@ import GeometryOps as GO, GeoInterface as GI
     p = GI.Polygon([lr])
     mp = GI.MultiPolygon([p, p, p])
     mls = GI.MultiLineString([ls, ls, ls])
-    
-    @testset "LinearSegments" begin
 
-        @test GO.segmentize(ls; max_distance = 0.5) isa GI.LineString
-        @test GO.segmentize(lr; max_distance = 0.5) isa GI.LinearRing
+    @testset_implementations "LinearSegments" begin
+        @test GO.segmentize($ls; max_distance = 0.5) isa GI.LineString
+        if GI.trait($lr) isa GI.LinearRingTrait
+            @test GO.segmentize($lr; max_distance = 0.5) isa GI.LinearRing
+        end
         # Test that linear rings are closed after segmentization
-        segmentized_ring = GO.segmentize(lr; max_distance = 0.5)
-        @test GI.getpoint(segmentized_ring, 1) == GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring))
+        segmentized_ring = GO.segmentize($lr; max_distance = 0.5)
+        @test GI.getpoint(segmentized_ring, 1) == 
+            GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring))
 
-        @test GO.segmentize(p; max_distance = 0.5) isa GI.Polygon
-        @test GO.segmentize(mp; max_distance = 0.5) isa GI.MultiPolygon
-        @test GI.ngeom(GO.segmentize(mp; max_distance = 0.5)) == 3
+        @test GO.segmentize($p; max_distance = 0.5) isa GI.Polygon
+        @test GO.segmentize($mp; max_distance = 0.5) isa GI.MultiPolygon
+        @test GI.ngeom(GO.segmentize($mp; max_distance = 0.5)) == 3
 
         # Now test multilinestrings
-        @test GO.segmentize(mls; max_distance = 0.5) isa GI.MultiLineString
-        @test GI.ngeom(GO.segmentize(mls; max_distance = 0.5)) == 3
+        @test GO.segmentize($mls; max_distance = 0.5) isa GI.MultiLineString
+        @test GI.ngeom(GO.segmentize($mls; max_distance = 0.5)) == 3
     end
 
-    @testset "GeodesicSegments" begin
-
-        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), ls) isa GI.LineString
-        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), lr) isa GI.LinearRing
+    @testset_implementations "GeodesicSegments" begin
+        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $ls) isa GI.LineString
+        if GI.trait($lr) isa GI.LinearRingTrait
+            @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr) isa GI.LinearRing
+        end
         # Test that linear rings are closed after segmentization
-        segmentized_ring = GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), lr)
+        segmentized_ring = GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr)
         @test GI.getpoint(segmentized_ring, 1) == GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring))
 
-        p = GI.Polygon([lr])
-        mp = GI.MultiPolygon([p, p, p])
-        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), p) isa GI.Polygon
-        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mp) isa GI.MultiPolygon
-        @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mp)) == 3
+        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $p) isa GI.Polygon
+        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp) isa GI.MultiPolygon
+        @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp)) == 3
 
         # Now test multilinestrings
-        mls = GI.MultiLineString([ls, ls, ls])
-        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mls) isa GI.MultiLineString
-        @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), mls)) == 3
+        @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls) isa GI.MultiLineString
+        @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls)) == 3
     end
 
 end
 
-@testset "LinearSegments" begin
-    lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
-    ct = GO.centroid(lr)
-    ar = GO.area(lr)
+lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+@testset_implementations "LinearSegments" begin
+    ct = GO.centroid($lr)
+    ar = GO.area($lr)
     for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10))
-        segmentized = GO.segmentize(GO.LinearSegments(; max_distance), lr)
+        segmentized = GO.segmentize(GO.LinearSegments(; max_distance), $lr)
         @test all(GO.centroid(segmentized) .≈ ct)
         @test GO.area(segmentized) ≈ ar
     end
 end
 
-@testset "GeodesicSegments" begin
-    lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)])
+@testset_implementations "GeodesicSegments" begin
     for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) .* 900
-        @test_nowarn segmentized = GO.segmentize(GO.GeodesicSegments(; max_distance), lr)
+        @test_nowarn segmentized = GO.segmentize(GO.GeodesicSegments(; max_distance), $lr)
     end
-end
\ No newline at end of file
+end
diff --git a/test/transformations/simplify.jl b/test/transformations/simplify.jl
index dd59ef56c..1661c566a 100644
--- a/test/transformations/simplify.jl
+++ b/test/transformations/simplify.jl
@@ -1,6 +1,9 @@
 using Test
 import GeoJSON, JLD2
-import GeometryOps as GO, GeoInterface as GI, LibGEOS as LG
+import GeometryOps as GO
+import GeoInterface as GI
+import LibGEOS as LG
+using ..TestHelpers
 
 datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data"))
 @testset "RadialDistance and VisvalingamWhyatt" begin
@@ -8,6 +11,7 @@ datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data"))
     fc2 = GeoJSON.read(joinpath(datadir, "simplify2.geojson"))
     fcs = [fc for i in 1:100]
 
+    # TODO: @testset_implementations doesn't handle feature collections yet
     for T in (GO.RadialDistance, GO.VisvalingamWhyatt)
         @test length(collect(GO.flatten(GI.PointTrait, GO.simplify(T(number=10), fc)))) == 10
         @test length(collect(GO.flatten(GI.PointTrait, GO.simplify(T(ratio=0.5), fc)))) == 39 # Half of 78
@@ -15,20 +19,25 @@ datadir = realpath(joinpath(dirname(pathof(GO)), "../test/data"))
         GO.simplify(T(tol=0.001), fcs; threaded=true, calc_extent=true)
     end
 end
+
 @testset "DouglasPeucker" begin
     poly_coords = JLD2.jldopen(joinpath(datadir, "complex_polygons.jld2"))["verts"][1:4]
     for c in poly_coords
         npoints = length(c[1])
         poly = LG.Polygon(c)
-        lg_vals = GI.coordinates(LG.simplify(poly, 100.0))[1]
-        reduced_npoints = length(lg_vals)
-        @test all(GI.coordinates(GO.simplify(poly; tol = 100.0))[1] .== lg_vals)
-        @test all(GI.coordinates(GO.simplify(poly; number = reduced_npoints))[1] .== lg_vals)
-        @test all(GI.coordinates(GO.simplify(poly; ratio = (reduced_npoints/npoints)))[1] .== lg_vals)
+        @testset_implementations "Polygon coords match LibGEOS simplify" begin
+            lg_vals = GI.coordinates(LG.simplify($poly, 100.0))[1]
+            reduced_npoints = length(lg_vals)
+            @test all(GI.coordinates(GO.simplify($poly; tol = 100.0))[1] .== lg_vals)
+            @test all(GI.coordinates(GO.simplify($poly; number = reduced_npoints))[1] .== lg_vals)
+            @test all(GI.coordinates(GO.simplify($poly; ratio = (reduced_npoints/npoints)))[1] .== lg_vals)
+        end
     end
     # Ensure last point isn't removed with curve
     c = poly_coords[1]
-    line = LG.LineString(c[1])
-    lg_vals = GI.coordinates(LG.simplify(line, 100.0))
-    @test all(GI.coordinates(GO.simplify(line; tol = 100.0)) .== lg_vals)
+    linestring = LG.LineString(c[1])
+    @testset_implementations "LineString coords match LibGEOS simplify" begin
+        lg_vals = GI.coordinates(LG.simplify($linestring, 100.0))
+        @test all(GI.coordinates(GO.simplify($linestring; tol = 100.0)) .== lg_vals)
+    end
 end
diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl
index f3d5b50ef..40a700e65 100644
--- a/test/transformations/transform.jl
+++ b/test/transformations/transform.jl
@@ -1,13 +1,15 @@
 using Test
- 
-import GeoInterface as GI, GeometryOps as GO
 using CoordinateTransformations
+import GeoInterface as GI
+import GeometryOps as GO
+using ..TestHelpers
 
-@testset "transform" begin
-    geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
-                       GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
+geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
+                   GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
+
+@testset_implementations "transform" begin
     translated = GI.Polygon([GI.LinearRing([[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]]), 
                              GI.LinearRing([[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]])])
     f = CoordinateTransformations.Translation(3.5, 1.5)
-    @test GO.transform(f, geom) == translated
+    @test GO.transform(f, $geom) == translated
 end