diff --git a/Project.toml b/Project.toml
index 64c916de6..ec62f505a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -48,7 +48,7 @@ GeometryOpsCore = "=0.1.5"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
-SortTileRecursiveTree = "0.1"
+SortTileRecursiveTree = "0.1.2"
 Statistics = "1"
 TGGeometry = "0.1"
 Tables = "1"
diff --git a/docs/src/experiments/regridding/regridding.jl b/docs/src/experiments/regridding/regridding.jl
new file mode 100644
index 000000000..2a4f0c8fc
--- /dev/null
+++ b/docs/src/experiments/regridding/regridding.jl
@@ -0,0 +1,118 @@
+import GeoInterface as GI, GeometryOps as GO
+using SortTileRecursiveTree: STRtree
+using SparseArrays: spzeros
+using Extents
+
+using CairoMakie, GeoInterfaceMakie
+
+include("sphericalpoints.jl")
+
+function area_of_intersection_operator(grid1, grid2; nodecapacity1 = 10, nodecapacity2 = 10)
+    area_of_intersection_operator(GO.Planar(), grid1, grid2; nodecapacity1 = nodecapacity1, nodecapacity2 = nodecapacity2)
+end
+
+function area_of_intersection_operator(m::GO.Manifold, grid1, grid2; nodecapacity1 = 10, nodecapacity2 = 10) # grid1 and grid2 are both vectors of polygons
+    A = spzeros(Float64, length(grid1), length(grid2))
+    # Prepare STRtrees for the two grids, to speed up intersection queries
+    # we may want to separately tune nodecapacity if one is much larger than the other.  
+    # specifically we may want to tune leaf node capacity via Hilbert packing while still 
+    # constraining inner node capacity.  But that can come later.
+    tree1 = STRtree(grid1; nodecapacity = nodecapacity1) 
+    tree2 = STRtree(grid2; nodecapacity = nodecapacity2)
+    # Do the dual query, which is the most efficient way to do this,
+    # by iterating down both trees simultaneously, rejecting pairs of nodes that do not intersect.
+    # when we find an intersection, we calculate the area of the intersection and add it to the result matrix.
+    GO.SpatialTreeInterface.do_dual_query(Extents.intersects, tree1, tree2) do i1, i2
+        p1, p2 = grid1[i1], grid2[i2]
+        # may want to check if the polygons intersect first, 
+        # to avoid antimeridian-crossing multipolygons viewing a scanline.
+        intersection_polys = try # can remove this now, got all the errors cleared up in the fix.
+            # At some future point, we may want to add the manifold here
+            # but for right now, GeometryOps only supports planar polygons anyway.
+            GO.intersection(p1, p2; target = GI.PolygonTrait())
+        catch e
+            @error "Intersection failed!" i1 i2
+            rethrow(e)
+        end
+
+        area_of_intersection = GO.area(m, intersection_polys)
+        if area_of_intersection > 0
+            A[i1, i2] += area_of_intersection
+        end
+    end
+
+    return A
+end
+
+grid1 = begin
+    gridpoints = [(i, j) for i in 0:2, j in 0:2]
+    [GI.Polygon([GI.LinearRing([gridpoints[i, j], gridpoints[i, j+1], gridpoints[i+1, j+1], gridpoints[i+1, j], gridpoints[i, j]])]) for i in 1:size(gridpoints, 1)-1, j in 1:size(gridpoints, 2)-1] |> vec
+end
+
+grid2 = begin
+    diamondpoly = GI.Polygon([GI.LinearRing([(0, 1), (1, 2), (2, 1), (1, 0), (0, 1)])])
+    trianglepolys = GI.Polygon.([
+        [GI.LinearRing([(0, 0), (1, 0), (0, 1), (0, 0)])],
+        [GI.LinearRing([(0, 1), (0, 2), (1, 2), (0, 1)])],
+        [GI.LinearRing([(1, 2), (2, 1), (2, 2), (1, 2)])],
+        [GI.LinearRing([(2, 1), (2, 0), (1, 0), (2, 1)])],
+    ])
+    [diamondpoly, trianglepolys...]
+end
+
+A = area_of_intersection_operator(grid1, grid2)
+
+# Now, let's perform some interpolation!
+area1 = vec(sum(A, dims=2))
+# test: @assert area1 == GO.area.(grid1)
+area2 = vec(sum(A, dims=1))
+# test: @assert area2 == GO.area.(grid2)
+
+values_on_grid2 = [0, 0, 5, 0, 0]
+poly(grid2; color = values_on_grid2, strokewidth = 2, strokecolor = :red)
+
+values_on_grid1 = A * values_on_grid2 ./ area1
+@assert sum(values_on_grid1 .* area1) == sum(values_on_grid2 .* area2)
+poly(grid1; color = values_on_grid1, strokewidth = 2, strokecolor = :blue)
+
+values_back_on_grid2 = A' * values_on_grid1 ./ area2
+@assert sum(values_back_on_grid2 .* area2) == sum(values_on_grid2 .* area2)
+poly(grid2; color = values_back_on_grid2, strokewidth = 2, strokecolor = :green)
+# We can see here that some data has diffused into the central diamond cell of grid2,
+# since it was overlapped by the top left cell of grid1.
+
+
+using SpeedyWeather
+using GeoMakie
+
+SpeedyWeatherGeoMakieExt = Base.get_extension(SpeedyWeather, :SpeedyWeatherGeoMakieExt)
+
+grid1 = rand(OctaHEALPixGrid, 5 + 100)
+grid2 = rand(FullGaussianGrid, 4 + 100)
+
+faces1 = SpeedyWeatherGeoMakieExt.get_faces(grid1)
+faces2 = SpeedyWeatherGeoMakieExt.get_faces(grid2)
+
+polys1 = GI.Polygon.(GI.LinearRing.(eachcol(faces1))) .|> GO.CutAtAntimeridianAndPoles() .|> GO.fix
+polys2 = GI.Polygon.(GI.LinearRing.(eachcol(faces2))) .|> GO.CutAtAntimeridianAndPoles() .|> GO.fix
+
+A = @time area_of_intersection_operator(polys1, polys2)
+
+p1 = polys1[93]
+p2 = polys2[105]
+
+f, a, p = poly(p1)
+poly!(a, p2)
+f
+
+# bug found in Foster Hormann tracing
+# but geos also does the same thing
+boxpoly = GI.Polygon([GI.LinearRing([(0, 0), (2, 0), (2, 2), (0, 2), (0, 0)])])
+diamondpoly = GI.Polygon([GI.LinearRing([(0, 1), (1, 2), (2, 1), (1, 0), (0, 1)])])
+
+diffpoly = GO.difference(boxpoly, diamondpoly; target = GI.PolygonTrait()) |> only
+cutpolys = GO.cut(diffpoly, GI.Line([(0, 0), (4, 0)])) # even cut misbehaves!
+
+
+
+
diff --git a/docs/src/experiments/regridding/sphericalpoints.jl b/docs/src/experiments/regridding/sphericalpoints.jl
new file mode 100644
index 000000000..4444a1ca9
--- /dev/null
+++ b/docs/src/experiments/regridding/sphericalpoints.jl
@@ -0,0 +1,116 @@
+import GeoInterface as GI
+import GeometryOps as GO
+import LinearAlgebra
+import LinearAlgebra: dot, cross
+
+## Get the area of a LinearRing with coordinates in radians
+struct SphericalPoint{T <: Real}
+	data::NTuple{3, T}
+end
+SphericalPoint(x, y, z) = SphericalPoint((x, y, z))
+
+# define the 4 basic mathematical operators elementwise on the data tuple
+Base.:+(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .+ q.data)
+Base.:-(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .- q.data)
+Base.:*(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data .* q.data)
+Base.:/(p::SphericalPoint, q::SphericalPoint) = SphericalPoint(p.data ./ q.data)
+# Define sum on a SphericalPoint to sum across its data
+Base.sum(p::SphericalPoint) = sum(p.data)
+
+# define dot and cross products
+LinearAlgebra.dot(p::SphericalPoint, q::SphericalPoint) = sum(p * q)
+function LinearAlgebra.cross(a::SphericalPoint, b::SphericalPoint)
+	a1, a2, a3 = a.data
+    b1, b2, b3 = b.data
+	SphericalPoint((a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1))
+end
+
+# Using Eriksson's formula for the area of spherical triangles: https://www.jstor.org/stable/2691141
+# This melts down only when two points are antipodal, which in our case will not happen.
+function _unit_spherical_triangle_area(a, b, c)
+    #t = abs(dot(a, cross(b, c)))
+    #t /= 1 + dot(b,c) + dot(c, a) + dot(a, b)
+    t = abs(dot(a, (cross(b - a, c - a))) / dot(b + a, c + a))
+    2*atan(t)
+end
+
+_lonlat_to_sphericalpoint(p) = _lonlat_to_sphericalpoint(GI.x(p), GI.y(p))
+function _lonlat_to_sphericalpoint(lon, lat)
+    lonsin, loncos = sincosd(lon)
+    latsin, latcos = sincosd(lat)
+    x = latcos * loncos
+    y = latcos * lonsin
+    z = latsin
+    return SphericalPoint(x,y,z)
+end
+
+
+
+# Extend area to spherical
+
+# TODO: make this the other way around, but that can wait.
+GO.area(m::GO.Planar, geoms) = GO.area(geoms)
+
+function GO.area(m::GO.Spherical, geoms)
+    return GO.applyreduce(+, GI.PolygonTrait(), geoms) do poly
+        GO.area(m, GI.PolygonTrait(), poly)
+    end * ((-m.radius^2)/ 2) # do this after the sum, to increase accuracy and minimize calculations.
+end
+
+function GO.area(m::GO.Spherical, ::GI.PolygonTrait, poly)
+    area = abs(_ring_area(m, GI.getexterior(poly)))
+    for interior in GI.gethole(poly)
+        area -= abs(_ring_area(m, interior))
+    end
+    return area
+end
+
+function _ring_area(m::GO.Spherical, ring)
+    # convert ring to a sequence of SphericalPoints
+    points = _lonlat_to_sphericalpoint.(GI.getpoint(ring))[1:end-1] # deliberately drop the closing point
+    p1, p2, p3 = points[1], points[2], points[3]
+    # For a spherical polygon, we can compute the area by splitting it into triangles
+    # and summing their areas. We use the first point as a common vertex for all triangles.
+    area = 0.0
+    # Sum areas of triangles formed by first point and consecutive pairs of points
+    np = length(points)
+    for i in 1:np
+        p1, p2, p3 = p2, p3, points[i]
+        area += _unit_spherical_triangle_area(p1, p2, p3)
+    end
+    return area
+end
+
+function _ring_area(m::GO.Spherical, ring)
+    # Convert ring points to spherical coordinates
+    points = GO.tuples(ring).geom
+    
+    # Remove last point if it's the same as first (closed ring)
+    if points[end] == points[1]
+        points = points[1:end-1]
+    end
+    
+    n = length(points)
+    if n < 3
+        return 0.0
+    end
+
+    area = 0.0
+
+    # Use L'Huilier's formula to sum the areas of spherical triangles
+    # formed by first point and consecutive pairs of points
+    for i in 1:n
+        p1, p2, p3 = points[mod1(i-1, n)], points[mod1(i, n)], points[mod1(i+1, n)]
+        area += sind(GI.y(p2)) * (GI.x(p3) - GI.x(p1))
+    end
+    
+    return area
+end
+
+
+
+
+# Test the area calculation
+p1 = GI.Polygon([GI.LinearRing(Point2f[(0, 0), (1, 0), (0, 1), (0, 0)] .- (Point2f(0.5, 0.5),))])
+
+GO.area(GO.Spherical(), p1)
\ No newline at end of file
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index db7efd2ee..b2faaf35e 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -21,12 +21,15 @@ using GeoInterface
 using GeometryBasics
 using LinearAlgebra, Statistics
 
+using GeometryBasics.StaticArrays
+
 import Tables, DataAPI
-import GeometryBasics.StaticArrays
 import DelaunayTriangulation # for convex hull and triangulation
 import ExactPredicates
 import Base.@kwdef
 import GeoInterface.Extents: Extents
+import SortTileRecursiveTree
+import SortTileRecursiveTree: STRtree
 
 const GI = GeoInterface
 const GB = GeometryBasics
@@ -44,6 +47,12 @@ include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
 
 using .LoopStateMachine, .SpatialTreeInterface
 
+include("utils/NaturalIndexing.jl")
+using .NaturalIndexing
+
+
+# Load utility modules in
+using .NaturalIndexing, .SpatialTreeInterface, .LoopStateMachine
 
 include("methods/angles.jl")
 include("methods/area.jl")
@@ -84,6 +93,7 @@ include("transformations/forcedims.jl")
 include("transformations/correction/geometry_correction.jl")
 include("transformations/correction/closed_ring.jl")
 include("transformations/correction/intersecting_polygons.jl")
+include("transformations/correction/cut_at_antimeridian.jl")
 
 # Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`.
 for name in names(GeoInterface)
diff --git a/src/methods/clipping/clipping_processor.jl b/src/methods/clipping/clipping_processor.jl
index 94bf4e859..41b30e1de 100644
--- a/src/methods/clipping/clipping_processor.jl
+++ b/src/methods/clipping/clipping_processor.jl
@@ -157,6 +157,264 @@ function _build_ab_list(alg::FosterHormannClipping, ::Type{T}, poly_a, poly_b, d
     return a_list, b_list, a_idx_list
 end
 
+"The number of vertices past which we should use a STRtree for edge intersection checking."
+const GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS = 32
+# Fallback convenience method so we can just pass the algorithm in
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    alg::FosterHormannClipping{M, A}, f_on_each_a::FA, f_after_each_a::FAAfter, f_on_each_maybe_intersect::FI, poly_a, poly_b, _t::Type{T} = Float64
+) where {FA, FAAfter, FI, T, M, A}
+    return foreach_pair_of_maybe_intersecting_edges_in_order(alg.manifold, alg.accelerator, f_on_each_a, f_after_each_a, f_on_each_maybe_intersect, poly_a, poly_b, T)
+end
+
+"""
+    foreach_pair_of_maybe_intersecting_edges_in_order(
+        manifold::M, accelerator::A,
+        f_on_each_a::FA,
+        f_after_each_a::FAAfter,
+        f_on_each_maybe_intersect::FI,
+        geom_a,
+        geom_b,
+        ::Type{T} = Float64
+    ) where {FA, FAAfter, FI, T, M <: Manifold, A <: IntersectionAccelerator}
+
+Decompose `geom_a` and `geom_b` into edge lists (unsorted), and then, logically, 
+perform the following iteration:
+
+```julia
+for (a_edge, i) in enumerate(eachedge(geom_a))
+    f_on_each_a(a_edge, i)
+    for (b_edge, j) in enumerate(eachedge(geom_b))
+        if may_intersect(a_edge, b_edge)
+            f_on_each_maybe_intersect(a_edge, b_edge)
+        end
+    end
+    f_after_each_a(a_edge, i)
+end
+```
+
+This may not be the exact acceleration that is performed - but it is 
+the logical sequence of events.  It also uses the `accelerator`, 
+and can automatically choose the best one based on an internal heuristic
+if you pass in an [`AutoAccelerator`](@ref).  
+
+For example, the `SingleSTRtree` accelerator is used along
+with extent thinning to avoid unnecessary edge intersection 
+checks in the inner loop.
+
+"""
+function foreach_pair_of_maybe_intersecting_edges_in_order(
+    manifold::M, accelerator::A, f_on_each_a::FA, f_after_each_a::FAAfter, f_on_each_maybe_intersect::FI, poly_a, poly_b, _t::Type{T} = Float64
+) where {FA, FAAfter, FI, T, M <: Manifold, A <: IntersectionAccelerator}
+    # TODO: dispatch on manifold
+    # this is suitable for planar
+    # but spherical / geodesic will need s2 support at some point,
+    # or -- even now -- just buffering
+    na = GI.npoint(poly_a)
+    nb = GI.npoint(poly_b)
+
+    accelerator = if accelerator isa AutoAccelerator
+        if na < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS && nb < GEOMETRYOPS_NO_OPTIMIZE_EDGEINTERSECT_NUMVERTS
+            NestedLoop()
+        else
+            SingleSTRtree()
+        end
+    else
+        accelerator
+    end
+
+    if accelerator isa NestedLoop
+        # if we don't have enough vertices in either of the polygons to merit a tree,
+        # then we can just do a simple nested loop
+        # this becomes extremely useful in e.g. regridding, 
+        # where we know the polygon will only ever have a few vertices.
+        # This is also applicable to any manifold, since the checking is done within
+        # the loop.
+        # First, loop over "each edge" in poly_a
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            isnothing(f_on_each_a) ||f_on_each_a(a1t, i)
+            for (j, (b1t, b2t)) in enumerate(eachedge(poly_b, T))
+                b1t == b2t && continue
+                LoopStateMachine.@controlflow f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # this should be aware of manifold by construction.
+            end
+            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+        end
+        # And we're done!
+    elseif accelerator isa SingleSTRtree
+        # This is the "middle ground" case - run only a strtree 
+        # on poly_b without doing so on poly_a.
+        # This is less complex than running a dual tree traversal,
+        # and reduces the overhead of constructing an edge list and tree on poly_a.
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        edges_b, indices_b = to_edgelist(ext_a, poly_b, T)
+        if isempty(edges_b) && !isnothing(f_on_each_a) && !isnothing(f_after_each_a)
+            # shortcut - nothing can possibly intersect
+            # so we just call f_on_each_a for each edge in poly_a
+            for i in 1:GI.npoint(poly_a)-1
+                pt = _tuple_point(GI.getpoint(poly_a, i), T)
+                f_on_each_a(pt, i)
+                f_after_each_a(pt, i)
+            end
+            return nothing
+        end
+
+        # This is the STRtree generated from the edges of poly_b
+        tree_b = STRtree(edges_b)
+
+        # this is a pre-allocation that will store the resuits of the query into tree_b
+        query_result = Int[] 
+        
+        # Loop over each vertex in poly_a
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            l1 = GI.Line(SVector{2}(a1t, a2t))
+            ext_l = GI.extent(l1)
+            # l = GI.Line(SVector{2}(a1t, a2t); extent=ext_l) # this seems to be unused - TODO remove
+            isnothing(f_on_each_a) || f_on_each_a(a1t, i)
+            # Query the STRtree for any edges in b that may intersect this edge
+            # This is sorted because we want to pretend we're doing the same thing
+            # as the nested loop above, and iterating through poly_b in order.
+            if Extents.intersects(ext_l, ext_b)
+                empty!(query_result)
+                SortTileRecursiveTree.query!(query_result, tree_b.rootnode, ext_l) # this is already sorted and uniqueified in STRtree.
+                # Loop over the edges in b that might intersect the edges in a
+                for j in query_result
+                    b1t, b2t = edges_b[j].geom
+                    b1t == b2t && continue
+                    # Manage control flow if the function returns a LoopStateMachine.Action
+                    # like Break(), Continue(), or Return()
+                    # This allows the function to break out of the loop early if it wants
+                    # without being syntactically inside the loop.
+                    LoopStateMachine.@controlflow f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), indices_b[j])) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
+                end
+            end
+            isnothing(f_after_each_a) || f_after_each_a(a1t, i)
+        end
+    elseif accelerator isa SingleNaturalTree
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        edges_b = to_edgelist(poly_b, T)
+
+        b_tree = NaturalIndexing.NaturalIndex(edges_b)
+
+        for (i, (a1t, a2t)) in enumerate(eachedge(poly_a, T))
+            a1t == a2t && continue
+            ext_l = Extents.Extent(X = minmax(a1t[1], a2t[1]), Y = minmax(a1t[2], a2t[2]))
+            isnothing(f_on_each_a) || f_on_each_a(a1t, i)
+            # Query the STRtree for any edges in b that may intersect this edge
+            # This is sorted because we want to pretend we're doing the same thing
+            # as the nested loop above, and iterating through poly_b in order.
+            if Extents.intersects(ext_l, ext_b)
+                # Loop over the edges in b that might intersect the edges in a
+                SpatialTreeInterface.depth_first_search(Base.Fix1(Extents.intersects, ext_l), b_tree) do j
+                    b1t, b2t = edges_b[j].geom
+                    b1t == b2t && return LoopStateMachine.Continue()
+                    # LoopStateMachine control is managed outside the loop, by the depth_first_search function.
+                    return f_on_each_maybe_intersect(((a1t, a2t), i), ((b1t, b2t), j)) # note the indices_b[j] here - we are using the index of the edge in the original edge list, not the index of the edge in the STRtree.
+                end
+            end
+        end
+
+    elseif accelerator isa DoubleNaturalTree
+        edges_a = to_edgelist(poly_a, T)
+        edges_b = to_edgelist(poly_b, T)
+
+        tree_a = NaturalIndexing.NaturalIndex(edges_a)
+        tree_b = NaturalIndexing.NaturalIndex(edges_b)
+
+        last_a_idx = 0
+
+        SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_edge_idx, b_edge_idx
+            a1t, a2t = edges_a[a_edge_idx].geom
+            b1t, b2t = edges_b[b_edge_idx].geom
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_on_each_a)
+                    for i in (last_a_idx+1):(a_edge_idx-1)
+                        f_on_each_a((edges_a[i].geom[1]), i)
+                        !isnothing(f_after_each_a) && f_after_each_a((edges_a[i].geom[1]), i)
+                    end
+                end
+                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
+            end
+
+            f_on_each_maybe_intersect(((a1t, a2t), a_edge_idx), ((b1t, b2t), b_edge_idx))
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_after_each_a)
+                    f_after_each_a(a1t, a_edge_idx)
+                end
+                last_a_idx = a_edge_idx
+            end
+        end
+
+        @show last_a_idx
+
+        if last_a_idx == 0 # the query did not find any intersections
+            if !isnothing(f_on_each_a) && isnothing(f_after_each_a)
+                return
+            else
+                for (i, edge) in enumerate(edges_a)
+                    !isnothing(f_on_each_a) && f_on_each_a(edge.geom[1], i)
+                    !isnothing(f_after_each_a) && f_after_each_a(edge.geom[1], i)
+                end
+            end
+        elseif last_a_idx < length(edges_a)
+            # the query terminated early - this will almost always be the case.
+            if !isnothing(f_on_each_a) && isnothing(f_after_each_a)
+                return
+            else
+                for (i, edge) in zip(last_a_idx+1:length(edges_a), view(edges_a, last_a_idx+1:length(edges_a)))
+                    !isnothing(f_on_each_a) && f_on_each_a(edge.geom[1], i)
+                    !isnothing(f_after_each_a) && f_after_each_a(edge.geom[1], i)
+                end
+            end
+        end
+    elseif accelerator isa ThinnedDoubleNaturalTree
+        ext_a, ext_b = GI.extent(poly_a), GI.extent(poly_b)
+        mutual_extent = Extents.intersection(ext_a, ext_b)
+
+        edges_a, indices_a = to_edgelist(mutual_extent, poly_a, T)
+        edges_b, indices_b = to_edgelist(mutual_extent, poly_b, T)
+
+        tree_a = NaturalIndexing.NaturalIndex(edges_a)
+        tree_b = NaturalIndexing.NaturalIndex(edges_b)
+
+        last_a_idx = 1
+
+        SpatialTreeInterface.dual_depth_first_search(Extents.intersects, tree_a, tree_b) do a_thinned_idx, b_thinned_idx
+            a_edge_idx = indices_a[a_thinned_idx]
+            b_edge_idx = indices_b[b_thinned_idx]
+
+            a1t, a2t = edges_a[a_thinned_idx].geom
+            b1t, b2t = edges_b[b_thinned_idx].geom
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_on_each_a)
+                    for i in last_a_idx:(a_edge_idx-1)
+                        f_on_each_a(a1t, a_edge_idx)
+                        !isnothing(f_after_each_a) && f_after_each_a(a1t, a_edge_idx)
+                    end
+                end
+                !isnothing(f_on_each_a) && f_on_each_a(a1t, a_edge_idx)
+            end
+
+            f_on_each_maybe_intersect(((a1t, a2t), a_edge_idx), ((b1t, b2t), b_edge_idx))
+
+            if last_a_idx < a_edge_idx
+                if !isnothing(f_after_each_a)
+                    f_after_each_a(a1t, a_edge_idx)
+                end
+                last_a_idx = a_edge_idx
+            end
+        end
+    else
+        error("Unsupported accelerator type: $accelerator.  FosterHormannClipping only supports NestedLoop() or SingleSTRtree().")
+    end
+
+    return nothing
+
+end
+
 #=
     _build_a_list(::Type{T}, poly_a, poly_b) -> (a_list, a_idx_list)
 
@@ -176,89 +434,101 @@ function _build_a_list(alg::FosterHormannClipping{M, A}, ::Type{T}, poly_a, poly
     a_list = PolyNode{T}[]  # list of points in poly_a
     sizehint!(a_list, n_a_edges)
     a_idx_list = Vector{Int}()  # finds indices of intersection points in a_list
-    a_count = 0  # number of points added to a_list
-    n_b_intrs = 0
-    # Loop through points of poly_a
-    local a_pt1
-    for (i, a_p2) in enumerate(GI.getpoint(poly_a))
-        a_pt2 = (T(GI.x(a_p2)), T(GI.y(a_p2)))
-        if i <= 1 || (a_pt1 == a_pt2)  # don't repeat points
-            a_pt1 = a_pt2
-            continue
-        end
-        # Add the first point of the edge to the list of points in a_list
-        new_point = PolyNode{T}(;point = a_pt1)
+    local a_count::Int = 0  # number of points added to a_list
+    local n_b_intrs::Int = 0
+    local prev_counter::Int = 0
+
+    function on_each_a(a_pt, i)
+        new_point = PolyNode{T}(;point = a_pt)
         a_count += 1
         push!(a_list, new_point)
-        # Find intersections with edges of poly_b
-        local b_pt1
         prev_counter = a_count
-        for (j, b_p2) in enumerate(GI.getpoint(poly_b))
-            b_pt2 = _tuple_point(b_p2, T)
-            if j <= 1 || (b_pt1 == b_pt2)  # don't repeat points
-                b_pt1 = b_pt2
-                continue
-            end
-            # Determine if edges intersect and how they intersect
-            line_orient, intr1, intr2 = _intersection_point(T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
-            if line_orient != line_out  # edges intersect
-                if line_orient == line_cross  # Intersection point that isn't a vertex
-                    int_pt, fracs = intr1
+        return nothing
+    end
+
+    function after_each_a(a_pt, i)
+        # Order intersection points by placement along edge using fracs value
+        if prev_counter < a_count
+            Δintrs = a_count - prev_counter
+            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
+            sort!(inter_points, by = x -> x.fracs[1])
+        end
+        return nothing
+    end
+
+    function on_each_maybe_intersect(((a_pt1, a_pt2), i), ((b_pt1, b_pt2), j))
+        if (b_pt1 == b_pt2)  # don't repeat points
+            b_pt1 = b_pt2
+            return
+        end
+        # Determine if edges intersect and how they intersect
+        line_orient, intr1, intr2 = _intersection_point(alg.manifold, T, (a_pt1, a_pt2), (b_pt1, b_pt2); exact)
+        if line_orient != line_out  # edges intersect
+            if line_orient == line_cross  # Intersection point that isn't a vertex
+                int_pt, fracs = intr1
+                new_intr = PolyNode{T}(;
+                    point = int_pt, inter = true, neighbor = j, # j is now equivalent to old j-1
+                    crossing = true, fracs = fracs,
+                )
+                a_count += 1
+                n_b_intrs += 1
+                push!(a_list, new_intr)
+                push!(a_idx_list, a_count)
+            else
+                (_, (α1, β1)) = intr1
+                # Determine if a1 or b1 should be added to a_list
+                add_a1 = α1 == 0 && 0 ≤ β1 < 1
+                a1_β = add_a1 ? β1 : zero(T)
+                add_b1 = β1 == 0 && 0 < α1 < 1
+                b1_α = add_b1 ? α1 : zero(T)
+                # If lines are collinear and overlapping, a second intersection exists
+                if line_orient == line_over
+                    (_, (α2, β2)) = intr2
+                    if α2 == 0 && 0 ≤ β2 < 1
+                        add_a1, a1_β = true, β2
+                    end
+                    if β2 == 0 && 0 < α2 < 1
+                        add_b1, b1_α = true, α2
+                    end
+                end
+                # Add intersection points determined above
+                if add_a1
+                    n_b_intrs += a1_β == 0 ? 0 : 1
+                    a_list[prev_counter] = PolyNode{T}(;
+                        point = a_pt1, inter = true, neighbor = j,
+                        fracs = (zero(T), a1_β),
+                    )
+                    push!(a_idx_list, prev_counter)
+                end
+                if add_b1
                     new_intr = PolyNode{T}(;
-                        point = int_pt, inter = true, neighbor = j - 1,
-                        crossing = true, fracs = fracs,
+                        point = b_pt1, inter = true, neighbor = j,
+                        fracs = (b1_α, zero(T)),
                     )
                     a_count += 1
-                    n_b_intrs += 1
                     push!(a_list, new_intr)
                     push!(a_idx_list, a_count)
-                else
-                    (_, (α1, β1)) = intr1
-                    # Determine if a1 or b1 should be added to a_list
-                    add_a1 = α1 == 0 && 0 ≤ β1 < 1
-                    a1_β = add_a1 ? β1 : zero(T)
-                    add_b1 = β1 == 0 && 0 < α1 < 1
-                    b1_α = add_b1 ? α1 : zero(T)
-                    # If lines are collinear and overlapping, a second intersection exists
-                    if line_orient == line_over
-                        (_, (α2, β2)) = intr2
-                        if α2 == 0 && 0 ≤ β2 < 1
-                            add_a1, a1_β = true, β2
-                        end
-                        if β2 == 0 && 0 < α2 < 1
-                            add_b1, b1_α = true, α2
-                        end
-                    end
-                    # Add intersection points determined above
-                    if add_a1
-                        n_b_intrs += a1_β == 0 ? 0 : 1
-                        a_list[prev_counter] = PolyNode{T}(;
-                            point = a_pt1, inter = true, neighbor = j - 1,
-                            fracs = (zero(T), a1_β),
-                        )
-                        push!(a_idx_list, prev_counter)
-                    end
-                    if add_b1
-                        new_intr = PolyNode{T}(;
-                            point = b_pt1, inter = true, neighbor = j - 1,
-                            fracs = (b1_α, zero(T)),
-                        )
-                        a_count += 1
-                        push!(a_list, new_intr)
-                        push!(a_idx_list, a_count)
-                    end
                 end
             end
-            b_pt1 = b_pt2
         end
-        # Order intersection points by placement along edge using fracs value
-        if prev_counter < a_count
-            Δintrs = a_count - prev_counter
-            inter_points = @view a_list[(a_count - Δintrs + 1):a_count]
-            sort!(inter_points, by = x -> x.fracs[1])
+        return nothing
+    end
+
+    # do the iteration but in an accelerated way
+    # this is equivalent to (but faster than)
+    #=
+    ```julia
+    for ((a1, a2), i) in eachedge(poly_a)
+        on_each_a(a1, i)
+        for ((b1, b2), j) in eachedge(poly_b)
+            on_each_maybe_intersect(((a1, a2), i), ((b1, b2), j))
         end
-        a_pt1 = a_pt2
+        after_each_a(a1, i)
     end
+    ```
+    =#
+    foreach_pair_of_maybe_intersecting_edges_in_order(alg, on_each_a, after_each_a, on_each_maybe_intersect, poly_a, poly_b, T)
+
     return a_list, a_idx_list, n_b_intrs
 end
 
diff --git a/src/methods/clipping/intersection.jl b/src/methods/clipping/intersection.jl
index 47d82077c..7af70d740 100644
--- a/src/methods/clipping/intersection.jl
+++ b/src/methods/clipping/intersection.jl
@@ -236,11 +236,13 @@ function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.Abstr
     # Check if the geometries extents even overlap
     Extents.intersects(GI.extent(a), GI.extent(b)) || return result
     # Create a list of edges from the two input geometries
-    edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
+    # edges_a, edges_b = map(sort! ∘ to_edges, (a, b))
     # Loop over pairs of edges and add any unique intersection points to results
-    for a_edge in edges_a, b_edge in edges_b
-        line_orient, intr1, intr2 = _intersection_point(T, a_edge, b_edge; exact)
-        line_orient == line_out && continue  # no intersection points
+    # TODO: add intersection acceleration here.
+
+    function f_on_each_maybe_intersect((a_edge, a_idx), (b_edge, b_idx))
+        line_orient, intr1, intr2 = _intersection_point(manifold, T, a_edge, b_edge; exact)
+        line_orient == line_out && return LoopStateMachine.Action(:continue) # use LoopStateMachine.Continue() to skip this edge - in this case it doesn't matter but you could use it to e.g. break once you found the first intersecting point.
         pt1, _ = intr1
         push!(result, pt1)  # if not line_out, there is at least one intersection point
         if line_orient == line_over # if line_over, there are two intersection points
@@ -248,6 +250,20 @@ function _intersection_points(manifold::M, accelerator::A, ::Type{T}, ::GI.Abstr
             push!(result, pt2)
         end
     end
+
+    # iterate over each pair of intersecting edges only,
+    # calling `f_on_each_maybe_intersect` for each pair 
+    # that may intersect.
+    foreach_pair_of_maybe_intersecting_edges_in_order(
+        manifold, accelerator, 
+        nothing, # f_on_each_a
+        nothing, # f_after_each_a
+        f_on_each_maybe_intersect, # f_on_each_maybe_intersect
+        a,
+        b,
+        T
+    )
+    
     #= TODO: We might be able to just add unique points with checks on the α and β values
     returned from `_intersection_point`, but this would be different for curves vs polygons
     vs multipolygons depending on if the shape is closed. This then wouldn't allow using the
diff --git a/src/transformations/correction/closed_ring.jl b/src/transformations/correction/closed_ring.jl
index b8076b821..7f4c7e905 100644
--- a/src/transformations/correction/closed_ring.jl
+++ b/src/transformations/correction/closed_ring.jl
@@ -46,7 +46,7 @@ See also [`GeometryCorrection`](@ref).
 """
 struct ClosedRing <: GeometryCorrection end
 
-application_level(::ClosedRing) = GI.PolygonTrait
+application_level(::ClosedRing) = TraitTarget(GI.PolygonTrait())
 
 function (::ClosedRing)(::GI.PolygonTrait, polygon)
     exterior = _close_linear_ring(GI.getexterior(polygon))
diff --git a/src/transformations/correction/cut_at_antimeridian.jl b/src/transformations/correction/cut_at_antimeridian.jl
new file mode 100644
index 000000000..6a88040ac
--- /dev/null
+++ b/src/transformations/correction/cut_at_antimeridian.jl
@@ -0,0 +1,358 @@
+#=
+# Antimeridian Cutting
+=#
+export CutAtAntimeridianAndPoles # TODO: too wordy?
+
+#=
+This correction cuts the geometry at the antimeridian and the poles, and returns a fixed geometry.
+
+The implementation is translated from the implementation in https://github.com/gadomski/antimeridian, 
+which is a Python package.  Several ports of that algorithm exist in e.g. Go, Rust, etc.
+
+At some point we will have to go in and clean up the implementation, remove all the hardcoded code,
+and make it more efficient by using raw geointerface, and not allocating so much (perhaps, by passing allocs around).
+
+But for right now it works, that's all I need.
+=#
+
+"""
+    CutAtAntimeridianAndPoles(; kwargs...) <: GeometryCorrection
+
+This correction cuts the geometry at the antimeridian and the poles, and returns a fixed geometry.
+"""
+Base.@kwdef struct CutAtAntimeridianAndPoles <: GeometryCorrection
+    "The value at the north pole, in your angular units"
+    northpole::Float64 = 90.0 # TODO not used!!!
+    "The value at the south pole, in your angular units"
+    southpole::Float64 = -90.0 # TODO not used!!!
+    "The value at the left edge of the antimeridian, in your angular units"
+    left::Float64 = -180.0
+    "The value at the right edge of the antimeridian, in your angular units"
+    right::Float64 = 180.0
+    "The period of the cyclic / cylindrical coordinate system's x value, usually computed automatically so you don't have to provide it."
+    period::Float64 = right - left
+    "If the polygon is known to enclose the north pole, set this to true"
+    force_north_pole::Bool=false # TODO not used!!!
+    "If the polygon is known to enclose the south pole, set this to true"
+    force_south_pole::Bool=false # TODO not used!!!
+    "If true, use the great circle method to find the antimeridian crossing, otherwise use the flat projection method.  There is no reason to set this to be off."
+    great_circle = true
+end
+
+function Base.show(io::IO, cut::CutAtAntimeridianAndPoles)
+    print(io, "CutAtAntimeridianAndPoles(left=$(cut.left), right=$(cut.right))")
+end
+Base.show(io::IO, ::MIME"text/plain", cut::CutAtAntimeridianAndPoles) = Base.show(io, cut)
+
+application_level(::CutAtAntimeridianAndPoles) = TraitTarget(GI.PolygonTrait(), GI.LineStringTrait(), GI.MultiLineStringTrait(), GI.MultiPolygonTrait())
+
+function (c::CutAtAntimeridianAndPoles)(trait::Union{GI.PolygonTrait, GI.MultiPolygonTrait, GI.LineStringTrait, GI.MultiLineStringTrait}, geom)
+    return _AntimeridianHelpers.cut_at_antimeridian(trait, geom)
+end
+
+module _AntimeridianHelpers
+
+import GeoInterface as GI
+import ..GeometryOps as GO
+using ..GeometryOps: CutAtAntimeridianAndPoles
+
+# Custom cross product for 3D tuples
+function _cross(a::Tuple{Float64,Float64,Float64}, b::Tuple{Float64,Float64,Float64})
+    return (
+        a[2]*b[3] - a[3]*b[2],
+        a[3]*b[1] - a[1]*b[3],
+        a[1]*b[2] - a[2]*b[1]
+    )
+end
+
+# Convert spherical degrees to cartesian coordinates
+function spherical_degrees_to_cartesian(c::CutAtAntimeridianAndPoles, point::Tuple{Float64,Float64})::Tuple{Float64,Float64,Float64}
+    # TODO: handle non-degree domains somehow
+    lon, lat = point
+    slon, clon = sincosd(lon)
+    slat, clat = sincosd(lat)
+    return (
+        clon * clat,
+        slon * clat,
+        slat
+    )
+end
+
+# Calculate crossing latitude using great circle method
+function crossing_latitude_great_circle(c::CutAtAntimeridianAndPoles, start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64})::Float64
+    # Convert points to 3D vectors
+    p1 = spherical_degrees_to_cartesian(c, start)
+    p2 = spherical_degrees_to_cartesian(c, endpoint)
+    
+    # Cross product defines plane through both points
+    n1 = _cross(p1, p2)
+    
+    # Unit vector that defines the meridian plane
+    n2 = spherical_degrees_to_cartesian(c, (c.left, 0.0))
+    
+    # Intersection of planes defined by cross product
+    intersection = _cross(n1, n2)
+    norm = sqrt(sum(intersection .^ 2))
+    intersection = intersection ./ norm
+    
+    # Convert back to spherical coordinates (just need latitude)
+    round(asind(intersection[3]), digits=7)
+end
+
+# Calculate crossing latitude using flat projection method
+function crossing_latitude_flat(c::CutAtAntimeridianAndPoles, start::Tuple{Float64,Float64}, stop::Tuple{Float64,Float64})::Float64
+    lat_delta = stop[2] - start[2]
+    if stop[1] > 0
+        round(
+            start[2] + (c.right - start[1]) * lat_delta / (stop[1] + c.period - start[1]),
+            digits=7
+        )
+    else
+        round(
+            start[2] + (start[1] - c.left) * lat_delta / (start[1] + c.period - stop[1]),
+            digits=7
+        )
+    end
+end
+
+# Main crossing latitude calculation function
+function crossing_latitude(c::CutAtAntimeridianAndPoles, start::Tuple{Float64,Float64}, endpoint::Tuple{Float64,Float64}, great_circle::Bool)::Float64
+    abs(start[1]) == 180 && return start[2]
+    abs(endpoint[1]) == 180 && return endpoint[2]
+    
+    return great_circle ? crossing_latitude_great_circle(c, start, endpoint) : crossing_latitude_flat(c, start, endpoint)
+end
+
+# Normalize coordinates to ensure longitudes are between -180 and 180
+function normalize_coords(c::CutAtAntimeridianAndPoles, coords::Vector{Tuple{Float64,Float64}})::Vector{Tuple{Float64,Float64}}
+    normalized = deepcopy(coords)
+    all_on_antimeridian = true
+    
+    for i in eachindex(normalized)
+        point = normalized[i]
+        prev_point = normalized[mod1(i-1, length(normalized))]
+        
+        if isapprox(point[1], c.right)
+            if abs(point[2]) != c.northpole && isapprox(prev_point[1], c.left)
+                normalized[i] = (c.left, point[2])
+            else
+                normalized[i] = (c.right, point[2])
+            end
+        elseif isapprox(point[1], c.left)
+            if abs(point[2]) != c.northpole && isapprox(prev_point[1], c.right)
+                normalized[i] = (c.right, point[2])
+            else
+                normalized[i] = (c.left, point[2])
+            end
+        else
+            normalized[i] = (((point[1] - c.left) % c.period) + c.left, point[2])
+            all_on_antimeridian = false
+        end
+    end
+    
+    return all_on_antimeridian ? coords : normalized
+end
+
+# Segment a ring of coordinates at antimeridian crossings
+function segment_ring(c::CutAtAntimeridianAndPoles, coords::Vector{Tuple{Float64,Float64}}, great_circle::Bool)::Vector{Vector{Tuple{Float64,Float64}}}
+    segments = Vector{Vector{Tuple{Float64,Float64}}}()
+    current_segment = Tuple{Float64,Float64}[]
+    
+    for i in 1:length(coords)-1
+        start = coords[i]
+        endpoint = coords[i+1]
+        push!(current_segment, start)
+        
+        # Check for antimeridian crossing
+        if (endpoint[1] - start[1] > 180) && (endpoint[1] - start[1] != 360)  # left crossing
+            lat = crossing_latitude(c, start, endpoint, great_circle)
+            push!(current_segment, (-180.0, lat))
+            push!(segments, current_segment)
+            current_segment = [(180.0, lat)]
+        elseif (start[1] - endpoint[1] > 180) && (start[1] - endpoint[1] != 360)  # right crossing
+            lat = crossing_latitude(c, endpoint, start, great_circle)
+            push!(current_segment, (180.0, lat))
+            push!(segments, current_segment)
+            current_segment = [(-180.0, lat)]
+        end
+    end
+    
+    # Handle last point and segment
+    if isempty(segments)
+        return Vector{Vector{Tuple{Float64,Float64}}}()  # No crossings
+    elseif coords[end] == segments[1][1]
+        # Join polygons
+        segments[1] = vcat(current_segment, segments[1])
+    else
+        push!(current_segment, coords[end])
+        push!(segments, current_segment)
+    end
+    
+    return segments
+end
+
+# Check if a segment is self-closing
+function is_self_closing(c::CutAtAntimeridianAndPoles, segment::Vector{Tuple{Float64,Float64}})::Bool
+    is_right = segment[end][1] == c.right
+    return segment[1][1] == segment[end][1] && (
+        (is_right && segment[1][2] > segment[end][2]) ||
+        (!is_right && segment[1][2] < segment[end][2])
+    )
+end
+
+# Build polygons from segments
+function build_polygons(c::CutAtAntimeridianAndPoles, segments::Vector{Vector{Tuple{Float64,Float64}}})::Vector{GI.Polygon}
+    isempty(segments) && return GI.Polygon[]
+    
+    segment = pop!(segments)
+    is_right = segment[end][1] == c.right
+    candidates = Tuple{Union{Nothing,Int},Float64}[]
+    
+    if is_self_closing(c, segment)
+        push!(candidates, (nothing, segment[1][2]))
+    end
+    
+    for (i, s) in enumerate(segments)
+        if s[1][1] == segment[end][1]
+            if (is_right && s[1][2] > segment[end][2] && 
+                (!is_self_closing(c, s) || s[end][2] < segment[1][2])) ||
+               (!is_right && s[1][2] < segment[end][2] && 
+                (!is_self_closing(c, s) || s[end][2] > segment[1][2]))
+                push!(candidates, (i, s[1][2]))
+            end
+        end
+    end
+    
+    # Sort candidates
+    sort!(candidates, by=x->x[2], rev=!is_right)
+    
+    index = isempty(candidates) ? nothing : candidates[1][1]
+    
+    if !isnothing(index)
+        # Join segments and recurse
+        segment = vcat(segment, segments[index])
+        segments[index] = segment
+        return build_polygons(c, segments)
+    else
+        # Handle self-joining segment
+        polygons = build_polygons(c, segments)
+        if !all(p == segment[1] for p in segment)
+            push!(polygons, GI.Polygon([segment]))
+        end
+        return filter(polygons) do p
+            x1 = GI.x(first(GI.getpoint(p)))
+            !all(p -> GI.x(p) == x1, GI.getpoint(p))
+        end
+    end
+end
+
+# Main function to cut a polygon at the antimeridian
+cut_at_antimeridian(x) = cut_at_antimeridian(CutAtAntimeridianAndPoles(), GI.trait(x), x)
+cut_at_antimeridian(c::CutAtAntimeridianAndPoles, x) = cut_at_antimeridian(c, GI.trait(x), x)
+cut_at_antimeridian(t::GI.AbstractTrait, x) = cut_at_antimeridian(CutAtAntimeridianAndPoles(), t, x)
+
+function cut_at_antimeridian(
+    c::CutAtAntimeridianAndPoles,
+    ::GI.PolygonTrait,
+    polygon::T;
+    force_north_pole::Bool=false,
+    force_south_pole::Bool=false,
+    fix_winding::Bool=true,
+    great_circle::Bool=true
+) where {T}
+    # Get exterior ring
+    exterior = GO.tuples(GI.getexterior(polygon)).geom
+    exterior = normalize_coords(c, exterior)
+    
+    # Segment the exterior ring
+    segments = segment_ring(c, exterior, great_circle)
+    
+    if isempty(segments)
+        # No antimeridian crossing
+        if fix_winding && GO.isclockwise(GI.LinearRing(exterior))
+            coord_vecs = reverse.(getproperty.(GO.tuples.(GI.getring(polygon)), :geom))
+            return GI.Polygon(normalize_coords.((c,), coord_vecs))
+        end
+        return polygon
+    end
+    
+    # Handle holes
+    holes = Vector{Vector{Tuple{Float64,Float64}}}()
+    for hole_idx in 1:GI.nhole(polygon)
+        hole = GO.tuples(GI.gethole(polygon, hole_idx)).geom
+        hole_segments = segment_ring(c, hole, great_circle)
+        
+        if !isempty(hole_segments)
+            if fix_winding
+                unwrapped = [(x % 360, y) for (x, y) in hole]
+                if !GO.isclockwise(GI.LineString(unwrapped))
+                    hole_segments = segment_ring(c, reverse(hole), great_circle)
+                end
+            end
+            append!(segments, hole_segments)
+        else
+            push!(holes, hole)
+        end
+    end
+    
+    # Build final polygons
+    result_polygons = build_polygons(c, segments)
+    
+    # Add holes to appropriate polygons
+    for hole in holes
+        for (i, poly) in enumerate(result_polygons)
+            if GO.contains(poly, GI.Point(hole[1]))
+                rings = GI.getring(poly)
+                push!(rings, hole)
+                result_polygons[i] = GI.Polygon(rings)
+                break
+            end
+        end
+    end
+
+    filter!(poly -> !iszero(GO.area(poly)), result_polygons)
+    
+    return length(result_polygons) == 1 ? result_polygons[1] : GI.MultiPolygon(result_polygons)
+end
+
+function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.AbstractCurveTrait, line::T; great_circle::Bool=true) where {T}
+    coords = GO.tuples(line).geom
+    segments = segment_ring(c, coords, great_circle)
+    
+    if isempty(segments)
+        return line
+    else
+        return GI.MultiLineString(segments)
+    end
+end
+
+
+function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.MultiPolygonTrait, x; kwargs...)
+    results = GI.Polygon[]
+    for poly in GI.getgeom(x)
+        result = cut_at_antimeridian(c, GI.PolygonTrait(), poly; kwargs...)
+        if result isa GI.Polygon
+            push!(results, result)
+        elseif result isa GI.MultiPolygon
+            append!(results, result.geom)
+        end
+    end
+    return GI.MultiPolygon(results)
+end
+
+function cut_at_antimeridian(c::CutAtAntimeridianAndPoles, t::GI.MultiLineStringTrait, multiline::T; great_circle::Bool=true) where {T}
+    linestrings = Vector{Vector{Tuple{Float64,Float64}}}()
+    
+    for line in GI.getgeom(multiline)
+        fixed = cut_at_antimeridian(c, GI.LineStringTrait(), line; great_circle)
+        if fixed isa GI.LineString
+            push!(linestrings, GO.tuples(fixed).geom)
+        else
+            append!(linestrings, GO.tuples.(GI.getgeom(fixed)).geom)
+        end
+    end
+    
+    return GI.MultiLineString(linestrings)
+end
+
+end
\ No newline at end of file
diff --git a/src/transformations/correction/geometry_correction.jl b/src/transformations/correction/geometry_correction.jl
index 5cccb5ea2..a9e227f13 100644
--- a/src/transformations/correction/geometry_correction.jl
+++ b/src/transformations/correction/geometry_correction.jl
@@ -33,28 +33,49 @@ This abstract type represents a geometry correction.
 ## Interface
 
 Any `GeometryCorrection` must implement two functions:
-    * `application_level(::GeometryCorrection)::AbstractGeometryTrait`: This function should return the `GeoInterface` trait that the correction is intended to be applied to, like `PointTrait` or `LineStringTrait` or `PolygonTrait`.
-    * `(::GeometryCorrection)(::AbstractGeometryTrait, geometry)::(some_geometry)`: This function should apply the correction to the given geometry, and return a new geometry.
+  - `application_level(::GeometryCorrection)::TraitTarget`: This function should 
+  return the `GeoInterface` trait that the correction is intended to be applied to, 
+  like `PointTrait` or `LineStringTrait` or `PolygonTrait`.  It can also return a 
+  union of traits via `TraitTarget`, but that behaviour is a bit tricky...
+  - `(::GeometryCorrection)(::AbstractGeometryTrait, geometry)::(some_geometry)`: 
+  This function should apply the correction to the given geometry, and return a new 
+  geometry.
 """
 abstract type GeometryCorrection end
 
+# Make sure that geometry corrections are treated as scalars when broadcasting.
+Base.Broadcast.broadcastable(c::GeometryCorrection) = (c,)
+
 application_level(gc::GeometryCorrection) = error("Not implemented yet for $(gc)")
 
 (gc::GeometryCorrection)(geometry) = gc(GI.trait(geometry), geometry)
 
 (gc::GeometryCorrection)(trait::GI.AbstractGeometryTrait, geometry) = error("Not implemented yet for $(gc) and $(trait).")
 
-function fix(geometry; corrections = GeometryCorrection[ClosedRing(),], kwargs...)
-    traits = application_level.(corrections)
+function fix(geometry; corrections = GeometryCorrection[ClosedRing()], kwargs...)
+    final_geoms = geometry
+    # Iterate through the corrections and apply them to the input.
+    # This allocates a _lot_, especially when reconstructing tables,
+    # but it's the only fully general way to do this that I can think of.
+    for correction in corrections
+        final_geoms = apply(correction, application_level(correction), final_geoms; kwargs...)
+    end
+    #=
+    # This was the old implementation
+    application_levels = application_level.(corrections)
     final_geometry = geometry
-    for Trait in (GI.PointTrait, GI.MultiPointTrait, GI.LineStringTrait, GI.LinearRingTrait, GI.MultiLineStringTrait, GI.PolygonTrait, GI.MultiPolygonTrait)
-        available_corrections = findall(x -> x == Trait, traits)
+    for trait in (GI.PointTrait(), GI.MultiPointTrait(), GI.LineStringTrait(), GI.LinearRingTrait(), GI.MultiLineStringTrait(), GI.PolygonTrait(), GI.MultiPolygonTrait())
+        available_corrections = findall(x -> trait in x, application_levels)
         isempty(available_corrections) && continue
-        @debug "Correcting for $(Trait)"
+        @debug "Correcting for $(trait), with corrections: " available_corrections
         net_function = reduce(∘, corrections[available_corrections])
-        final_geometry = apply(net_function, Trait, final_geometry; kwargs...)
+        # TODO: this allocates too much, because it keeps reconstructing higher level geoms.
+        # We might want some way to embed the fixes in reconstruct/rebuild, which would imply a modified apply pipeline...
+        final_geometry = apply(net_function, trait, final_geometry; kwargs...)
     end
     return final_geometry
+    =#
+    return final_geoms
 end
 
 # ## Available corrections
diff --git a/src/transformations/correction/intersecting_polygons.jl b/src/transformations/correction/intersecting_polygons.jl
index 2223f6b89..73d857525 100644
--- a/src/transformations/correction/intersecting_polygons.jl
+++ b/src/transformations/correction/intersecting_polygons.jl
@@ -56,7 +56,7 @@ See also [`GeometryCorrection`](@ref).
 """
 struct UnionIntersectingPolygons <: GeometryCorrection end
 
-application_level(::UnionIntersectingPolygons) = GI.MultiPolygonTrait
+application_level(::UnionIntersectingPolygons) = TraitTarget(GI.MultiPolygonTrait())
 
 function (::UnionIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
     union_multipoly = tuples(multipoly)
@@ -99,7 +99,7 @@ See also [`GeometryCorrection`](@ref), [`UnionIntersectingPolygons`](@ref).
 """
 struct DiffIntersectingPolygons <: GeometryCorrection end
 
-application_level(::DiffIntersectingPolygons) = GI.MultiPolygonTrait
+application_level(::DiffIntersectingPolygons) = TraitTarget(GI.MultiPolygonTrait())
 
 function (::DiffIntersectingPolygons)(::GI.MultiPolygonTrait, multipoly)
     diff_multipoly = tuples(multipoly)
diff --git a/src/utils/NaturalIndexing.jl b/src/utils/NaturalIndexing.jl
new file mode 100644
index 000000000..080da277e
--- /dev/null
+++ b/src/utils/NaturalIndexing.jl
@@ -0,0 +1,245 @@
+module NaturalIndexing
+
+import GeoInterface as GI
+import Extents
+
+using ..SpatialTreeInterface
+
+import ..GeometryOps as GO # TODO: only needed for NaturallyIndexedRing, remove when that is removed.
+
+export NaturalIndex, NaturallyIndexedRing, prepare_naturally
+
+"""
+    NaturalLevel{E <: Extents.Extent}
+
+A level in the natural tree.  Stored in a vector in [`NaturalIndex`](@ref).
+
+- `extents` is a vector of extents of the children of the node
+"""
+struct NaturalLevel{E <: Extents.Extent}
+    extents::Vector{E} # child extents
+end
+
+Base.show(io::IO, level::NaturalLevel) = print(io, "NaturalLevel($(length(level.extents)) extents)")
+Base.show(io::IO, ::MIME"text/plain", level::NaturalLevel) = Base.show(io, level)
+
+"""
+    NaturalIndex{E <: Extents.Extent}
+
+A natural tree index.  Stored in a vector in [`NaturalIndex`](@ref).
+
+- `nodecapacity` is the "spread", number of children per node
+- `extent` is the extent of the tree
+- `levels` is a vector of [`NaturalLevel`](@ref)s
+"""
+struct NaturalIndex{E <: Extents.Extent}
+    nodecapacity::Int # "spread", number of children per node
+    extent::E
+    levels::Vector{NaturalLevel{E}}
+end
+
+Extents.extent(idx::NaturalIndex) = idx.extent
+
+function Base.show(io::IO, ::MIME"text/plain", idx::NaturalIndex)
+    println(io, "NaturalIndex with $(length(idx.levels)) levels and $(idx.nodecapacity) children per node")
+    println(io, "extent: $(idx.extent)")
+end
+function Base.show(io::IO, idx::NaturalIndex)
+    println(io, "NaturalIndex($(length(idx.levels)) levels, $(idx.extent))")
+end
+
+function NaturalIndex(geoms; nodecapacity = 32)
+    # Get the extent type initially (coord order, coord type, etc.)
+    # so that the construction is type stable.
+    e1 = GI.extent(first(geoms))
+    E = typeof(e1)
+    return NaturalIndex{E}(geoms; nodecapacity = nodecapacity)
+end
+function NaturalIndex(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+    # If we are passed a vector of extents - inflate immediately!
+    return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
+end
+
+function NaturalIndex{E}(geoms; nodecapacity = 32) where E <: Extents.Extent
+    # If passed a vector of geometries, and we know the type of the extent,
+    # then simply retrieve the extents so they can serve as the "last-level" 
+    # extents.
+    # Finally, call the lowest level method that performs inflation.
+    last_level_extents = GI.extent.(geoms)
+    return NaturalIndex{E}(last_level_extents; nodecapacity = nodecapacity)
+end
+# This is the main constructor that performs inflation.
+function NaturalIndex{E}(last_level_extents::Vector{E}; nodecapacity = 32) where E <: Extents.Extent
+    ngeoms = length(last_level_extents)
+    last_level = NaturalLevel(last_level_extents)
+
+    nlevels = _number_of_levels(nodecapacity, ngeoms)
+
+    levels = Vector{NaturalLevel{E}}(undef, nlevels)
+    levels[end] = last_level
+    # Iterate backwards, from bottom to top level,
+    # and build up the level extent vectors.
+    for level_index in (nlevels-1):(-1):1
+        prev_level = levels[level_index+1] # this is always instantiated, since we are iterating backwards
+        nrects = _number_of_keys(nodecapacity, nlevels - (level_index), ngeoms)
+        extents = [
+            begin
+                start = (rect_index - 1) * nodecapacity + 1
+                stop = min(start + nodecapacity - 1, length(prev_level.extents))
+                reduce(Extents.union, view(prev_level.extents, start:stop))
+            end
+            for rect_index in 1:nrects
+        ]
+        levels[level_index] = NaturalLevel(extents)
+    end
+
+    return NaturalIndex(nodecapacity, reduce(Extents.union, levels[1].extents), levels)
+
+end
+
+function _number_of_keys(nodecapacity::Int, level::Int, ngeoms::Int)
+    return ceil(Int, ngeoms / (nodecapacity ^ (level)))
+end
+
+"""
+    _number_of_levels(nodecapacity::Int, ngeoms::Int)
+
+Calculate the number of levels in a natural tree for a given number of geometries and node capacity.
+
+## How this works
+
+The number of keys in a level is given by `ngeoms / nodecapacity ^ level`.
+
+The number of levels is the smallest integer such that the number of keys in the last level is 1.
+So it goes - if that makes sense.
+"""
+function _number_of_levels(nodecapacity::Int, ngeoms::Int)
+    level = 1
+    while _number_of_keys(nodecapacity, level, ngeoms) > 1
+        level += 1
+    end
+    return level
+end
+
+
+# This is like a pointer to a node in the tree.
+"""
+    NaturalIndexNode{E <: Extents.Extent}
+
+A reference to a node in the natural tree.  Kind of like a tree cursor.
+
+- `parent_index` is a pointer to the parent index
+- `level` is the level of the node in the tree
+- `index` is the index of the node in the level
+- `extent` is the extent of the node
+"""
+struct NaturalIndexNode{E <: Extents.Extent}
+    parent_index::NaturalIndex{E}
+    level::Int
+    index::Int
+    extent::E
+end
+
+Extents.extent(node::NaturalIndexNode) = node.extent
+
+# What does SpatialTreeInterface require of trees?
+# - Parents completely cover their children
+# - `GI.extent(node)` returns `Extent` 
+#   - can mean that `Extents.extent(node)` returns the extent of the node
+# - `nchild(node)` returns the number of children of the node
+# - `getchild(node)` returns an iterator over all children of the node
+# - `getchild(node, i)` returns the i-th child of the node
+# - `isleaf(node)` returns a boolean indicating whether the node is a leaf
+# - `child_indices_extents(node)` returns an iterator over the indices and extents of the children of the node
+
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndex}) = true
+SpatialTreeInterface.isspatialtree(::Type{<: NaturalIndexNode}) = true
+
+function SpatialTreeInterface.nchild(node::NaturalIndexNode)
+    start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
+    stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
+    return stop_idx - start_idx + 1
+end
+
+function SpatialTreeInterface.getchild(node::NaturalIndexNode, i::Int)
+    child_index = (node.index - 1) * node.parent_index.nodecapacity + i
+    return NaturalIndexNode(
+        node.parent_index, 
+        node.level + 1, # increment level by 1
+        child_index, # index of this particular child
+        node.parent_index.levels[node.level+1].extents[child_index] # the extent of this child
+    )
+end
+
+# Get all children of a node
+function SpatialTreeInterface.getchild(node::NaturalIndexNode)
+    return (SpatialTreeInterface.getchild(node, i) for i in 1:SpatialTreeInterface.nchild(node))
+end
+
+SpatialTreeInterface.isleaf(node::NaturalIndexNode) = node.level == length(node.parent_index.levels) - 1
+
+function SpatialTreeInterface.child_indices_extents(node::NaturalIndexNode)
+    start_idx = (node.index - 1) * node.parent_index.nodecapacity + 1
+    stop_idx = min(start_idx + node.parent_index.nodecapacity - 1, length(node.parent_index.levels[node.level+1].extents))
+    return ((i, node.parent_index.levels[node.level+1].extents[i]) for i in start_idx:stop_idx)
+end
+
+# implementation for "root node" / top level tree
+
+SpatialTreeInterface.isleaf(node::NaturalIndex) = length(node.levels) == 1
+
+SpatialTreeInterface.nchild(node::NaturalIndex) = length(node.levels[1].extents)
+
+SpatialTreeInterface.getchild(node::NaturalIndex) = SpatialTreeInterface.getchild(NaturalIndexNode(node, 0, 1, node.extent))
+SpatialTreeInterface.getchild(node::NaturalIndex, i) = SpatialTreeInterface.getchild(NaturalIndexNode(node, 0, 1, node.extent), i)
+
+SpatialTreeInterface.child_indices_extents(node::NaturalIndex) = (i_ext for i_ext in enumerate(node.levels[1].extents))
+
+"""
+    NaturallyIndexedRing(points; nodecapacity = 32)
+
+A linear ring that contains a natural index.
+
+!!! warning
+    This will be removed in favour of prepared geometry - the idea here
+    is just to test what interface works best to store things in.
+"""
+struct NaturallyIndexedRing
+    points::Vector{Tuple{Float64, Float64}}
+    index::NaturalIndex{Extents.Extent{(:X, :Y), NTuple{2, NTuple{2, Float64}}}}
+end
+
+function NaturallyIndexedRing(points::Vector{Tuple{Float64, Float64}}; nodecapacity = 32)
+    index = NaturalIndex(GO.edge_extents(GI.LinearRing(points)); nodecapacity)
+    return NaturallyIndexedRing(points, index)
+end
+NaturallyIndexedRing(ring::NaturallyIndexedRing) = ring
+
+function GI.convert(::Type{NaturallyIndexedRing}, ::GI.LinearRingTrait, geom)
+    points = GO.tuples(geom).geom
+    return NaturallyIndexedRing(points)
+end
+
+Base.show(io::IO, ::MIME"text/plain", ring::NaturallyIndexedRing) = Base.show(io, ring)
+Base.show(io::IO, ring::NaturallyIndexedRing) = print(io, "NaturallyIndexedRing($(length(ring.points)) points) with index $(sprint(show, ring.index))")
+
+GI.ncoord(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = 2
+GI.is3d(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
+GI.ismeasured(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = false
+
+GI.ngeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = length(ring.points)
+GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing) = ring.points
+GI.getgeom(::GI.LinearRingTrait, ring::NaturallyIndexedRing, i::Int) = ring.points[i]
+
+Extents.extent(ring::NaturallyIndexedRing) = ring.index.extent
+
+GI.isgeometry(::Type{<: NaturallyIndexedRing}) = true
+GI.geomtrait(::NaturallyIndexedRing) = GI.LinearRingTrait()
+
+function prepare_naturally(geom)
+    return GO.apply(GI.PolygonTrait(), geom) do poly
+        return GI.Polygon([GI.convert(NaturallyIndexedRing, GI.LinearRingTrait(), ring) for ring in GI.getring(poly)])
+    end
+end
+
+end # module NaturalIndexing
\ No newline at end of file
diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
index d306d75b4..0066142f6 100644
--- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
+++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl
@@ -7,7 +7,7 @@ import GeoInterface as GI
 import AbstractTrees
 
 # public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
-export query, do_query
+export query
 export FlatNoTree
 
 # The spatial tree interface and its implementations are defined here.
diff --git a/test/methods/clipping/polygon_clipping.jl b/test/methods/clipping/polygon_clipping.jl
index 2fc4575fa..638b627bd 100644
--- a/test/methods/clipping/polygon_clipping.jl
+++ b/test/methods/clipping/polygon_clipping.jl
@@ -166,7 +166,7 @@ test_pairs = [
 const ϵ = 1e-10
 # Compare clipping results from GeometryOps and LibGEOS
 function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
-    GO_result_list = GO_f(p1, p2; target = GI.PolygonTrait())
+    
     LG_result_geom = LG_f(p1, p2)
     if LG_result_geom isa LG.GeometryCollection
         poly_list = LG.Polygon[]
@@ -175,38 +175,43 @@ function compare_GO_LG_clipping(GO_f, LG_f, p1, p2)
         end
         LG_result_geom = LG.MultiPolygon(poly_list)
     end
-    # Check if nothing is returned
-    if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
-        return true
-    end
-    # Check for unnecessary points
-    if sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom)
-        return false
-    end
-    # Make sure last point is repeated
-    for poly in GO_result_list
-        for ring in GI.getring(poly)
-            GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)) && return false
+
+    for _accelerator in (GO.AutoAccelerator(), GO.NestedLoop(), GO.SingleSTRtree(), GO.SingleNaturalTree(), #=GO.DoubleNaturalTree(), =# #=GO.ThinnedDoubleNaturalTree(), =# #=GO.DoubleSTRtree()=#)
+    @testset let accelerator = _accelerator # this is a ContextTestSet that is otherwise invisible but adds context to the testset
+        GO_result_list = GO_f(GO.FosterHormannClipping(accelerator), p1, p2; target = GI.PolygonTrait())
+        # Check if nothing is returned
+        if isempty(GO_result_list) && (LG.isEmpty(LG_result_geom) || LG.area(LG_result_geom) == 0)
+            @test true
+            continue
+        end
+        # Check for unnecessary points
+        @test !(sum(GI.npoint, GO_result_list; init = 0.0) > GI.npoint(LG_result_geom))
+        # Make sure last point is repeated
+        for poly in GO_result_list
+            for ring in GI.getring(poly)
+                @test !(GI.getpoint(ring, 1) != GI.getpoint(ring, GI.npoint(ring)))
+            end
         end
-    end
 
-    # Check if polygons cover the same area
-    local GO_result_geom
-    if length(GO_result_list) == 1
-        GO_result_geom = GO_result_list[1]
-    else
-        GO_result_geom = GI.MultiPolygon(GO_result_list)
-    end
-    diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
-    diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
-    return diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+        # Check if polygons cover the same area
+        local GO_result_geom
+        if length(GO_result_list) == 1
+            GO_result_geom = GO_result_list[1]
+        else
+            GO_result_geom = GI.MultiPolygon(GO_result_list)
+        end
+        diff_1_area = LG.area(LG.difference(GO_result_geom, LG_result_geom))
+        diff_2_area = LG.area(LG.difference(LG_result_geom, GO_result_geom))
+        @test diff_1_area ≤ ϵ && diff_2_area ≤ ϵ
+    end # testset
+    end # loop
 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
         @testset_implementations "$sg1 $f_name $sg2 - $sdesc" begin
-            @test compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) 
+            compare_GO_LG_clipping(GO_f, LG_f, $p1, $p2) # this executes tests internally
         end
     end
     return
diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl
index eb61f553b..3f78b05c3 100644
--- a/test/utils/SpatialTreeInterface.jl
+++ b/test/utils/SpatialTreeInterface.jl
@@ -4,6 +4,7 @@ using GeometryOps.SpatialTreeInterface
 using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
 using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search
 using GeometryOps.SpatialTreeInterface: FlatNoTree
+using GeometryOps.NaturalIndexing: NaturalIndex
 using Extents
 using SortTileRecursiveTree: STRtree
 using NaturalEarth
@@ -218,7 +219,15 @@ end
     test_find_point_in_all_countries(STRtree)
 end
 
-
+# Test NaturalIndex implementation
+@testset "STRtree" begin
+    test_basic_interface(NaturalIndex)
+    test_child_indices_extents(NaturalIndex)
+    test_query_functionality(NaturalIndex)
+    test_dual_query_functionality(NaturalIndex)
+    test_geometry_support(NaturalIndex)
+    test_find_point_in_all_countries(NaturalIndex)
+end
 
 # This testset is not used because Polylabel.jl has some issues.