diff --git a/GeometryOpsCore/src/types/algorithm.jl b/GeometryOpsCore/src/types/algorithm.jl
index 5fb21c4a3..4e2f171e2 100644
--- a/GeometryOpsCore/src/types/algorithm.jl
+++ b/GeometryOpsCore/src/types/algorithm.jl
@@ -42,8 +42,80 @@ MyIndependentAlgorithm(m::Manifold; kw1 = 1, kw2 = "hello") = MyIndependentAlgor
 
 export Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm
 
+"""
+    abstract type Algorithm{M <: Manifold}
+
+The abstract supertype for all GeometryOps algorithms.  
+These define how to perform a particular [`Operation`](@ref).
+
+An algorithm may be associated with one or many [`Manifold`](@ref)s.  
+It may either have the manifold as a field, or have it as a static parameter 
+(e.g. `struct GEOS <: Algorithm{Planar}`).
+
+## Interface
+
+All `Algorithm`s must implement the following methods:
+
+- `rebuild(alg, manifold::Manifold)` Rebuild algorithm `alg` with a new manifold 
+  as passed in the second argument.  This may error and throw a [`WrongManifoldException`](@ref)
+  if the manifold is not compatible with that algorithm.
+- `manifold(alg::Algorithm)` Return the manifold associated with the algorithm.
+- `best_manifold(alg::Algorithm, input)`: Return the best manifold for that algorithm, in the absence of
+  any other context.  WARNING: this may change in future and is not stable!
+
+The actual implementation is left to the implementation of that particular [`Operation`](@ref).
+
+## Notable subtypes
+
+- [`AutoAlgorithm`](@ref): Tells the [`Operation`](@ref) receiving 
+  it to automatically select the best algorithm for its input data.
+- [`ManifoldIndependentAlgorithm`](@ref): An abstract supertype for an algorithm that works on any manifold.
+  The manifold must be stored in the algorithm for a `ManifoldIndependentAlgorithm`, and accessed via `manifold(alg)`.
+- [`SingleManifoldAlgorithm`](@ref): An abstract supertype for an algorithm that only works on a 
+  single manifold, specified in its type parameter.  `SingleManifoldAlgorithm{Planar}` is a special case
+  that does not have to store its manifold, since that doesn't contain any information.  All other 
+  `SingleManifoldAlgorithm`s must store their manifold, since they do contain information.
+- [`NoAlgorithm`](@ref): A type that indicates no algorithm is to be used, essentially the equivalent
+  of `nothing`.
+"""
 abstract type Algorithm{M <: Manifold} end
 
+"""
+    manifold(alg::Algorithm)::Manifold
+
+Return the manifold associated with the algorithm.  
+
+May be any subtype of [`Manifold`](@ref).
+"""
+function manifold end
+
+# The below definition is a special case, since [`Planar`](@ref) has no contents, being a 
+# singleton struct.
+# If that changes in the future, then this method must be deleted.
+manifold(::Algorithm{<: Planar}) = Planar()
+
+"""
+    best_manifold(alg::Algorithm, input)::Manifold
+    
+Return the best [`Manifold`](@ref) for the algorithm `alg` based on the given `input`.
+
+May be any subtype of [`Manifold`](@ref).
+"""
+function best_manifold end
+
+# ## Implementation of basic algorithm types
+
+# ### `AutoAlgorithm`
+
+"""
+    AutoAlgorithm{T, M <: Manifold}(manifold::M, x::T)
+
+Indicates that the [`Operation`](@ref) should automatically select the best algorithm for
+its input data, based on the passed in manifold (may be an [`AutoManifold`](@ref)) and data 
+`x`.
+
+The actual implementation is left to the implementation of that particular [`Operation`](@ref).
+"""
 struct AutoAlgorithm{T, M <: Manifold} <: Algorithm{M} 
     manifold::M
     x::T
@@ -52,18 +124,33 @@ end
 AutoAlgorithm(m::Manifold; kwargs...) = AutoAlgorithm(m, kwargs)
 AutoAlgorithm(; kwargs...) = AutoAlgorithm(AutoManifold(), kwargs)
 
+manifold(a::AutoAlgorithm) = a.manifold
+rebuild(a::AutoAlgorithm, m::Manifold) = AutoAlgorithm(m, a.x)
+
+
+# ### `ManifoldIndependentAlgorithm`
+
+"""
+    abstract type ManifoldIndependentAlgorithm{M <: Manifold} <: Algorithm{M}
+
+The abstract supertype for a manifold-independent algorithm, i.e., one which may work on any manifold.
 
+The manifold is stored in the algorithm for a `ManifoldIndependentAlgorithm`, and accessed via `manifold(alg)`.
+"""
 abstract type ManifoldIndependentAlgorithm{M <: Manifold} <: Algorithm{M} end
 
-abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M} end
 
-struct NoAlgorithm{M <: Manifold} <: Algorithm{M} 
-    m::M
-end
+# ### `SingleManifoldAlgorithm`
 
-NoAlgorithm() = NoAlgorithm(Planar()) # TODO: add a NoManifold or AutoManifold type?
-# Maybe AutoManifold
-# and then we have DD.format like materialization  
+"""
+    abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M}
+
+The abstract supertype for a single-manifold algorithm, i.e., one which is known to only work 
+on a single manifold.
+
+The manifold may be accessed via `manifold(alg)`.
+"""
+abstract type SingleManifoldAlgorithm{M <: Manifold} <: Algorithm{M} end
 
 function (Alg::Type{<: SingleManifoldAlgorithm{M}})(m::M; kwargs...) where {M}
     # successful - the algorithm is designed for this manifold
@@ -76,3 +163,24 @@ function (Alg::Type{<: SingleManifoldAlgorithm{M}})(m::Manifold; kwargs...) wher
     # throw a WrongManifoldException and be done with it
     throw(WrongManifoldException{typeof(m), M, Alg}())
 end
+
+
+# ### `NoAlgorithm`
+
+"""
+    NoAlgorithm(manifold)
+
+A type that indicates no algorithm is to be used, essentially the equivalent
+of `nothing`.
+
+Stores a manifold within itself.
+"""
+struct NoAlgorithm{M <: Manifold} <: Algorithm{M} 
+    m::M
+end
+
+NoAlgorithm() = NoAlgorithm(Planar()) # TODO: add a NoManifold or AutoManifold type?
+
+manifold(a::NoAlgorithm) = a.m
+# Maybe AutoManifold
+# and then we have DD.format like materialization  
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index fcdd6bb59..de398abd3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -8,6 +8,7 @@ CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df"
 ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
+Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -21,17 +22,20 @@ Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
+TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 
 [extensions]
 GeometryOpsFlexiJoinsExt = "FlexiJoins"
 GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsProjExt = "Proj"
+GeometryOpsTGGeometryExt = "TGGeometry"
 
 [compat]
 CoordinateTransformations = "0.5, 0.6"
 DataAPI = "1"
 DelaunayTriangulation = "1.0.4"
 ExactPredicates = "2.2.8"
+Extents = "0.1.5"
 FlexiJoins = "0.1.30"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.2"
@@ -42,8 +46,9 @@ LinearAlgebra = "1"
 Proj = "1"
 SortTileRecursiveTree = "0.1"
 Statistics = "1"
+TGGeometry = "0.1"
 Tables = "1"
-julia = "1.9"
+julia = "1.10"
 
 [extras]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
@@ -53,7 +58,6 @@ DimensionalData = "0703355e-b756-11e9-17c0-8b28908087d0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
-GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
@@ -64,7 +68,8 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
+TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "SafeTestsets", "Shapefile", "Test"]
+test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "SafeTestsets", "Shapefile", "TGGeometry", "Test"]
diff --git a/benchmarks/benchmark_plots.jl b/benchmarks/benchmark_plots.jl
index b0f5432a4..a3ac53d03 100644
--- a/benchmarks/benchmark_plots.jl
+++ b/benchmarks/benchmark_plots.jl
@@ -73,6 +73,7 @@ function plot_trials(
         legend_orientation = :horizontal,
         legend_halign = 1.0,
         legend_valign = -0.25,
+        legend_kws = (;),
     )
 
     xs, ys, labels = [], [], []
@@ -118,7 +119,8 @@ function plot_trials(
             tellheight = legend_position isa Union{Tuple, Makie.Automatic} && (legend_position isa Makie.Automatic || length(legend_position) != 3) && legend_orientation == :horizontal, 
             halign = legend_halign, 
             valign = legend_valign, 
-            orientation = legend_orientation
+            orientation = legend_orientation,
+            legend_kws...,
         )
         ax.xtickcolor[] = ax.xgridcolor[]
         ax
diff --git a/benchmarks/benchmarks.jl b/benchmarks/benchmarks.jl
index af5f15ed9..438f60f23 100644
--- a/benchmarks/benchmarks.jl
+++ b/benchmarks/benchmarks.jl
@@ -10,18 +10,23 @@
 import GeometryOps as GO, 
     GeoInterface as GI, 
     GeometryBasics as GB, 
-    LibGEOS as LG
-import GeoJSON, NaturalEarth
+    LibGEOS as LG,
+    GeoFormatTypes as GFT
+import GeoJSON, NaturalEarth, WellKnownGeometry
+using CoordinateTransformations: Translation, LinearMap
 # In order to benchmark, we'll actually use the new [Chairmarks.jl](https://github.com/lilithhafner/Chairmarks.jl), 
 # since it's significantly faster than BenchmarkTools.  To keep benchmarks organized, though, we'll still use BenchmarkTools' 
 # `BenchmarkGroup` structure.
 using Chairmarks
 import BenchmarkTools: BenchmarkGroup
+using ProgressMeter
 # We use CairoMakie to visualize our results!
 using CairoMakie, MakieThemes, GeoInterfaceMakie
 # Finally, we import some general utility packages:
 using Statistics, CoordinateTransformations
 
+include("benchmark_plots.jl")
+
 # We also set up some utility functions for later on.
 """
 Returns LibGEOS and GeometryOps' native geometries, 
@@ -155,7 +160,7 @@ circle_difference_suite = circle_suite["difference"] = BenchmarkGroup(["title:Ci
 circle_intersection_suite = circle_suite["intersection"] = BenchmarkGroup(["title:Circle intersection", "subtitle:Tested on a regular circle"])
 circle_union_suite = circle_suite["union"] = BenchmarkGroup(["title:Circle union", "subtitle:Tested on a regular circle"])
 
-n_points_values = round.(Int, exp10.(LinRange(1, 4, 10)))
+n_points_values = round.(Int, exp10.(LinRange(0.7, 6, 15)))
 @time for n_points in n_points_values
     circle = GI.Wrappers.Polygon([tuple.((cos(θ) for θ in LinRange(0, 2π, n_points)), (sin(θ) for θ in LinRange(0, 2π, n_points)))])
     closed_circle = GO.ClosedRing()(circle)
@@ -197,29 +202,117 @@ f
 # Now, we get to benchmarking:
 
 
-usa_o_lg, usa_o_go = lg_and_go(usa_poly)
-usa_r_lg, usa_r_go = lg_and_go(usa_reflected)
+usa_o_lg, usa_o_go = lg_and_go(usa_poly);
+usa_r_lg, usa_r_go = lg_and_go(usa_reflected);
 
 # First, we'll test union:
-printstyled("LibGEOS"; color = :red, bold = true)
-println()
-@be LG.union($usa_o_lg, $usa_r_lg) seconds=5
-printstyled("GeometryOps"; color = :blue, bold = true)
-println()
-@be GO.union($usa_o_go, $usa_r_go; target = GI.PolygonTrait) seconds=5
-
-# Next, intersection:
-printstyled("LibGEOS"; color = :red, bold = true)
-println()
-@be LG.intersection($usa_o_lg, $usa_r_lg) seconds=5
-printstyled("GeometryOps"; color = :blue, bold = true)
-println()
-@be GO.intersection($usa_o_go, $usa_r_go; target = GI.PolygonTrait) seconds=5
-
-# and finally the difference:
-printstyled("LibGEOS"; color = :red, bold = true)
-println()
-@be LG.difference(usa_o_lg, usa_r_lg) seconds=5
-printstyled("GeometryOps"; color = :blue, bold = true)
-println()
-@be go_diff = GO.difference(usa_o_go, usa_r_go; target = GI.PolygonTrait) seconds=5
+begin
+    printstyled("Union"; color = :green, bold = true)
+    println()
+    printstyled("LibGEOS"; color = :red, bold = true)
+    println()
+    display(@be LG.union($usa_o_lg, $usa_r_lg) seconds=5)
+    printstyled("GeometryOps"; color = :blue, bold = true)
+    println()
+    display(@be GO.union($usa_o_go, $usa_r_go; target = GI.PolygonTrait) seconds=5)
+    println()
+    # Next, intersection:
+    printstyled("Intersection"; color = :green, bold = true)
+    println()
+    printstyled("LibGEOS"; color = :red, bold = true)
+    println()
+    display(@be LG.intersection($usa_o_lg, $usa_r_lg) seconds=5)
+    printstyled("GeometryOps"; color = :blue, bold = true)
+    println()
+    display(@be GO.intersection($usa_o_go, $usa_r_go; target = GI.PolygonTrait) seconds=5)
+    # and finally the difference:
+    printstyled("Difference"; color = :green, bold = true)
+    println()
+    printstyled("LibGEOS"; color = :red, bold = true)
+    println()
+    display(@be LG.difference(usa_o_lg, usa_r_lg) seconds=5)
+    printstyled("GeometryOps"; color = :blue, bold = true)
+    println()
+    display(@be GO.difference(usa_o_go, usa_r_go; target = GI.PolygonTrait) seconds=5)
+end
+
+
+
+
+# ## Vancouver watershed benchmarks
+#=
+
+Vancouver Island has ~1,300 watersheds.  LibGEOS uses this exact data
+in their tests, so we'll use it in ours as well!
+
+We'll start by loading the data, and then we'll use it to benchmark various operations.
+
+=#
+
+# The CRS for this file is EPSG:3005, or as a PROJ string,
+# `"+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"`
+# TODO: this doesn't work with WellKnownGeometry.  Why?
+
+watersheds = mktempdir() do dir
+    cd(dir) do
+        wkt_gz = download("https://github.com/pramsey/geos-performance/raw/refs/heads/master/data/watersheds.wkt.gz", "watersheds.wkt.gz")
+        run(`gunzip watersheds.wkt.gz`)
+        return [
+            GO.tuples(GFT.WellKnownText(GFT.Geom(), line))
+            for line in eachline("watersheds.wkt")
+        ]
+    end
+end
+ 
+watershed_polygons = only.(GI.getgeom.(watersheds))
+
+import SortTileRecursiveTree as STR
+tree = STR.STRtree(watershed_polygons)
+query_result = STR.query(tree, GI.extent(watershed_polygons[1]))
+
+GO.intersects.((watershed_polygons[1],), watershed_polygons[query_result])
+
+@be GO.union($(watershed_polygons[1]), $(watershed_polygons[2]); target = $GI.PolygonTrait())
+@be LG.union($(watershed_polygons[1] |> GI.convert(LG)), $(watershed_polygons[2] |> GI.convert(LG)))
+
+function union_coverage(intersection_f::IF, union_f::UF, polygons::Vector{T}; showplot = true, showprogress = true) where {T, IF, UF}
+    tree = STR.STRtree(polygons)
+    all_intersected = falses(length(polygons))
+    accumulator = polygons[1]
+    all_intersected[1] = true
+    i = 1
+
+    (showprogress && (prog = Progress(length(all_intersected))))
+
+    for i in 1:length(polygons)
+        query_result = STR.query(tree, GI.extent(accumulator))
+        for idx in query_result
+            if !(all_intersected[idx] || !(intersection_f(polygons[idx], accumulator)))
+                result = union_f(polygons[idx], accumulator)
+                accumulator = result
+                all_intersected[idx] = true
+                (showprogress && next!(prog))
+            end
+        end
+        showplot && display(poly(view(polygons, all_intersected); color = rand(RGBf, sum(all_intersected))), axis = (; title = "$(GI.trait(accumulator) isa GI.PolygonTrait ? "Polygon" : "MultiPolygon with $(GI.ngeom(accumulator)) individual polys")"))
+        all(all_intersected) && break # if done then finish
+    end
+
+    return accumulator
+end
+
+@time union_coverage(LG.intersects, LG.union, watershed_polygons .|> GI.convert(LG); showplot = false, showprogress = true)
+
+@time union_coverage(GO.intersects, (x, y) -> (GO.union(x, y; target = GI.MultiPolygonTrait())), watershed_polygons; showplot = false, showprogress = true)
+
+
+using GADM
+
+# austria is landlocked and will form a coverage
+# something like India will not -- because it has islands
+ind_fc = GADM.get("AUT"; depth = 1)
+ind_polys = GI.geometry.(GI.getfeature(ind_fc)) |> x -> GO.tuples(x; calc_extent = true)
+
+
+
+@time union_coverage(GO.intersects, (x, y) -> (GO.union(x, y; target = GI.MultiPolygonTrait())), ind_polys; showplot = true, showprogress = true)
diff --git a/benchmarks/geometry_providers.jl b/benchmarks/geometry_providers.jl
index 9f917eeb8..80715160f 100644
--- a/benchmarks/geometry_providers.jl
+++ b/benchmarks/geometry_providers.jl
@@ -139,8 +139,8 @@ end
 
 function Makie.convert_arguments(::Makie.SampleBased, labels::AbstractVector{<: AbstractString}, bs::AbstractVector{<: Chairmarks.Benchmark})
     ts = map(b -> getproperty.(b.samples, :time), bs)
-    labels = 
-    return flatten
+    labels = reduce(vcat, (fill(l, length(t)) for (l, t) in zip(labels, ts)))
+    return flatten(labels), flatten(ts)
 end
 
 function Makie.convert_arguments(::Type{Makie.Errorbars}, xs, bs::AbstractVector{<: Chairmarks.Benchmark})
diff --git a/benchmarks/polygon_scaling.jl b/benchmarks/polygon_scaling.jl
index f38cfcc23..f1c74a5e3 100644
--- a/benchmarks/polygon_scaling.jl
+++ b/benchmarks/polygon_scaling.jl
@@ -7,7 +7,7 @@ p26 = GI.Polygon([[(3.0, 1.0), (8.0, 1.0), (8.0, 7.0), (3.0, 7.0), (3.0, 5.0), (
 
 suite = BenchmarkGroup(["title:Polygon intersection timing","subtitle:Single polygon, densified"])
 
-for max_distance in exp10.(LinRange(-1, 1.5, 10))
+for max_distance in exp10.(LinRange(-2, 1.5, 10))
     p25s = GO.segmentize(p25; max_distance)
     p26s = GO.segmentize(p26; max_distance)
     n_verts = GI.npoint(p25s)
diff --git a/ext/GeometryOpsTGGeometryExt/GeometryOpsTGGeometryExt.jl b/ext/GeometryOpsTGGeometryExt/GeometryOpsTGGeometryExt.jl
new file mode 100644
index 000000000..2ea0f0ac0
--- /dev/null
+++ b/ext/GeometryOpsTGGeometryExt/GeometryOpsTGGeometryExt.jl
@@ -0,0 +1,12 @@
+module GeometryOpsTGGeometryExt
+
+using GeometryOps: TG
+import GeometryOps as GO
+
+using TGGeometry
+
+for jl_fname in TGGeometry.TG_PREDICATES
+    @eval GO.$jl_fname(::TG, geom1, geom2) = TGGeometry.$jl_fname(geom1, geom2)
+end
+
+end
\ No newline at end of file
diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl
index 030597023..c0d169f6e 100644
--- a/src/GeometryOps.jl
+++ b/src/GeometryOps.jl
@@ -6,6 +6,7 @@ import GeometryOpsCore
 import GeometryOpsCore: 
                 TraitTarget,
                 Manifold, Planar, Spherical, Geodesic, AutoManifold, WrongManifoldException,
+                manifold, best_manifold,
                 Algorithm, AutoAlgorithm, ManifoldIndependentAlgorithm, SingleManifoldAlgorithm, NoAlgorithm,
                 BoolsAsTypes, True, False, booltype, istrue,
                 TaskFunctors, 
diff --git a/src/types.jl b/src/types.jl
index 3a89f0f9f..48ddbebb8 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,14 +1,21 @@
 #=
 # Types
 
-This file defines some fundamental types used in GeometryOps.
+This file defines some types used in GeometryOps.  
 
 !!! warning
-    Unlike in other Julia packages, only some types are defined in this file, not all. 
-    This is because we define types in the files where they are used, to make it easier to understand the code.
+    Many type definitions are in `GeometryOpsCore`, not here.  Look there for the definitions of the basic types like `Manifold`, `Algorithm`, etc.
+
+
+## Naming
+
+We force all our external algorithm types to be uppercase, to make them different
+from the package names.  This is really relevant for the `PROJ` algorithm, since 
+the Julia package is called `Proj.jl`.  If we called our type `Proj`, it would
+conflict with the package's name - as we saw with GeoFormatTypes' `GeoJSON` and `GeoJSON.jl`.
 
 =#
-export GEOS
+export GEOS, TG, PROJ
 #=
 
 ## `GEOS`
@@ -23,6 +30,77 @@ useful for two reasons:
 
 =#
 
+# ## C-library planar algorithms
+"""
+    abstract type CLibraryPlanarAlgorithm <: GeometryOpsCore.SingleManifoldAlgorithm{Planar} end
+
+This is a type which extends `GeometryOpsCore.SingleManifoldAlgorithm{Planar}`,
+and is used as an abstract supertype for some C library based algorithms.
+
+The type requires that algorithm structs be arranged as:
+```
+struct MyAlgorithm <: CLibraryPlanarAlgorithm
+    manifold::Planar
+    params::NamedTuple
+end
+```
+
+Then you get a nice constructor for free, as well as the 
+`get(alg, key, value)` and `get(alg, key) do ...` syntax.
+Plus the [`enforce`](@ref) method, which will check that given keyword arguments
+are present.
+"""
+abstract type CLibraryPlanarAlgorithm <: GeometryOpsCore.SingleManifoldAlgorithm{Planar} end
+
+function (::Type{T})(; params...) where {T <: CLibraryPlanarAlgorithm}
+    nt = NamedTuple(params)
+    return T(nt)
+end
+(T::Type{<: CLibraryPlanarAlgorithm})(::Planar, params::NamedTuple) = T(params)
+
+manifold(alg::CLibraryPlanarAlgorithm) = Planar()
+best_manifold(alg::CLibraryPlanarAlgorithm, input) = Planar()
+
+# Rebuild methods with manifolds are here.
+function rebuild(alg::T, m::Planar) where {T <: CLibraryPlanarAlgorithm}
+    return T(alg.params) # TODO: should this not rebuild at all, then, since nothing will change?
+end
+function rebuild(alg::T, m::AutoManifold) where {T <: CLibraryPlanarAlgorithm}
+    return T(alg.params) # "rebuild" as a planar algorithm.
+end
+function rebuild(alg::T, m::M) where {T <: CLibraryPlanarAlgorithm, M <: Manifold}
+    throw(GeometryOpsCore.WrongManifoldException{M, Planar, T}("The algorithm `$(typeof(alg))` is only compatible with planar manifolds."))
+end
+# Rebuild methods for parameters are here.  This ends up being quite useful really.
+rebuild(alg::T, params::NamedTuple) where {T <: CLibraryPlanarAlgorithm} = T(params)
+rebuild(alg::T; params...) where {T <: CLibraryPlanarAlgorithm} = T(NamedTuple(params))
+
+# These are definitions for convenience, so we don't have to type out 
+# `alg.params` every time.
+
+Base.get(alg::CLibraryPlanarAlgorithm, key, value) = Base.get(alg.params, key, value)
+Base.get(f::Function, alg::CLibraryPlanarAlgorithm, key) = Base.get(f, alg.params, key)
+
+"""
+    enforce(alg::CLibraryPlanarAlgorithm, kw::Symbol, f)
+
+Enforce the presence of a keyword argument in a `GEOS` algorithm, and return `alg.params[kw]`.
+
+Throws an error if the key is not present, and mentions `f` in the error message (since there isn't 
+a good way to get the name of the function that called this method).
+
+This applies to all `CLibraryPlanarAlgorithm` types, like [`GEOS`](@ref) and [`TG`](@ref).
+"""
+function enforce(alg::CLibraryPlanarAlgorithm, kw::Symbol, f)
+    if haskey(alg.params, kw)
+        return alg.params[kw]
+    else
+        error("$(f) requires a `$(kw)` keyword argument to the `GEOS` algorithm, which was not provided.")
+    end
+end
+
+# ## GEOS - call into LibGEOS.jl
+
 """
     GEOS(; params...)
 
@@ -33,31 +111,72 @@ Dispatch is generally carried out using the names of the keyword arguments.
 For example, `segmentize` will only accept a `GEOS` struct with only a 
 `max_distance` keyword, and no other.
 
-It's generally a lot slower than the native Julia implementations, since
+It's generally somewhat slower than the native Julia implementations, since
 it must convert to the LibGEOS implementation and back - so be warned!
+
+## Extended help
+
+This uses the [LibGEOS.jl](https://github.com/JuliaGeometry/LibGEOS.jl) package,
+which is a Julia wrapper around the C library GEOS (https://trac.osgeo.org/geos).
 """
-struct GEOS
+struct GEOS <: CLibraryPlanarAlgorithm # SingleManifoldAlgorithm{Planar}
     params::NamedTuple
 end
 
-function GEOS(; params...)
-    nt = NamedTuple(params)
-    return GEOS(nt)
+# ## TG - call into TGGeometry.jl
+
+"""
+    TG(; params...)
+
+A struct which instructs the method it's passed to as an algorithm
+to use the appropriate TG function via `TGGeometry.jl` for the operation.
+
+It's generally a lot faster than the native Julia implementations, but only
+supports planar manifolds / operations.  Also, it only supports geometric predicates,
+specifically the ones which the underlying `tg` library supports.  These are:
+
+[`equals`](@ref), [`intersects`](@ref), [`disjoint`](@ref), [`contains`](@ref), 
+[`within`](@ref), [`covers`](@ref), [`coveredby`](@ref), and [`touches`](@ref).
+
+## Extended help
+
+This uses the [TGGeometry.jl](https://github.com/JuliaGeo/TGGeometry.jl) package,
+which is a Julia wrapper around the `tg` C library (https://github.com/tidwall/tg).
+"""
+struct TG <: CLibraryPlanarAlgorithm
+    params::NamedTuple
 end
-# These are definitions for convenience, so we don't have to type out 
-# `alg.params` every time.
-Base.get(alg::GEOS, key, value) = Base.get(alg.params, key, value)
-Base.get(f::Function, alg::GEOS, key) = Base.get(f, alg.params, key)
+
+# ## PROJ - call into Proj.jl
 
 """
-    enforce(alg::GO.GEOS, kw::Symbol, f)
+    PROJ(; params...)
 
-Enforce the presence of a keyword argument in a `GEOS` algorithm, and return `alg.params[kw]`.
+A struct which instructs the method it's passed to as an algorithm
+to use the appropriate PROJ function via `Proj.jl` for the operation.
 
-Throws an error if the key is not present, and mentions `f` in the error message (since there isn't 
-a good way to get the name of the function that called this method).
+## Extended help
+
+This is the default algorithm for [`reproject`](@ref), and will also be the default algorithm for 
+operations on geodesics like [`area`](@ref) and [`arclength`](@ref).
 """
-function enforce(alg::GEOS, kw::Symbol, f)
+struct PROJ{M <: Manifold} <: Algorithm{M}
+    manifold::M
+    params::NamedTuple
+end
+
+PROJ(; params...) = PROJ(Planar(), NamedTuple(params))
+PROJ(m::Manifold) = PROJ(m, NamedTuple())
+
+manifold(alg::PROJ) = alg.manifold
+rebuild(alg::PROJ, m::Manifold) = PROJ(m, alg.params)
+rebuild(alg::PROJ, params::NamedTuple) = PROJ(alg.manifold, params)
+
+# We repeat these functions here because PROJ does not subtype `CLibraryPlanarAlgorithm`.
+
+Base.get(alg::PROJ, key, value) = Base.get(alg.params, key, value)
+Base.get(f::Function, alg::PROJ, key) = Base.get(f, alg.params, key)
+function enforce(alg::PROJ, kw::Symbol, f)
     if haskey(alg.params, kw)
         return alg.params[kw]
     else
diff --git a/test/extensions/tggeometry.jl b/test/extensions/tggeometry.jl
new file mode 100644
index 000000000..2f5b70396
--- /dev/null
+++ b/test/extensions/tggeometry.jl
@@ -0,0 +1,60 @@
+using Test
+
+import GeometryOps as GO, GeoInterface as GI
+import TGGeometry
+
+# This test file is only really certifying that the extension is working...
+# we rely on TGGeometry.jl's tests to verify correctness, since it's not an exact
+# library and a lot of our tests do check exactness.
+
+# That's why it isn't included in the main polygon test suite.
+
+point = (0.0, 0.0)
+multipoint = GI.MultiPoint([(0.0, 0.0), (1.0, 1.0)])
+linestring = GI.LineString([(0.0, 0.0), (1.0, 1.0)])
+multilinestring = GI.MultiLineString([[(0.0, 0.0), (1.0, 1.0)], [(2.0, 2.0), (3.0, 3.0)]])
+polygon = GI.Polygon([[(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0), (0.0, 0.0)]])
+multipolygon = GI.MultiPolygon([polygon, GO.transform(p -> (GI.x(p), GI.y(p) + 1.0), polygon)])
+
+_xplus5(g) = GO.transform(p -> (GI.x(p) + 5.0, GI.y(p)), g)
+
+disjoint_point = _xplus5(point)
+disjoint_multipoint = _xplus5(multipoint)
+disjoint_linestring = _xplus5(linestring)
+disjoint_multilinestring = _xplus5(multilinestring)
+disjoint_polygon = _xplus5(polygon)
+disjoint_multipolygon = _xplus5(multipolygon)
+
+@testset "Internal consistency with TGGeometry.jl" begin
+    for funsym in TGGeometry.TG_PREDICATES
+        for geom1 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+            for geom2 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+                @testset let predicate = funsym, geom1 = geom1, geom2 = geom2
+                    @test Base.getproperty(TGGeometry, predicate)(geom1, geom2) == Base.getproperty(GO, predicate)(GO.TG(), geom1, geom2)
+                end
+            end
+            for geom2 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+                @testset let predicate = funsym, geom1 = geom1, geom2 = geom2
+                    @test Base.getproperty(TGGeometry, predicate)(geom1, geom2) == Base.getproperty(GO, predicate)(GO.TG(), geom1, geom2)
+                end
+            end
+        end
+    end
+end
+
+@testset "Consistency with GeometryOps algorithms for simple cases" begin
+    for funsym in TGGeometry.TG_PREDICATES
+        for geom1 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+            for geom2 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+                @testset let predicate = funsym, geom1 = geom1, geom2 = geom2
+                    @test Base.getproperty(TGGeometry, predicate)(geom1, geom2) == Base.getproperty(GO, predicate)(GO.TG(), geom1, geom2)
+                end
+            end
+            for geom2 in (point, multipoint, linestring, multilinestring, polygon, multipolygon)
+                @testset let predicate = funsym, geom1 = geom1, geom2 = geom2
+                    @test Base.getproperty(TGGeometry, predicate)(geom1, geom2) == Base.getproperty(GO, predicate)(GO.TG(), geom1, geom2)
+                end
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/external/jts/jts_testset_reader.jl b/test/external/jts/jts_testset_reader.jl
new file mode 100644
index 000000000..57e71c218
--- /dev/null
+++ b/test/external/jts/jts_testset_reader.jl
@@ -0,0 +1,123 @@
+using XML
+
+import WellKnownGeometry
+import GeoFormatTypes as GFT
+import GeometryOps as GO
+import GeoInterface as GI
+
+"""
+    jts_wkt_to_geom(wkt::String)
+
+Convert a JTS WKT string to a GeometryOps geometry, via WellKnownGeometry.jl and GO.tuples.
+
+The reason this exists is because WellKnownGeometry doesn't work well with newlines in subsidiary geometries,
+so this sanitizes the input before parsing and converting.
+"""
+function jts_wkt_to_geom(wkt::String)
+    sanitized_wkt = join(strip.(split(wkt, "\n")), "")
+    geom = GFT.WellKnownText(GFT.Geom(), sanitized_wkt)
+    return GO.tuples(geom)
+end
+
+struct TestItem{T}
+    operation::String
+    arg1::GO.GI.Wrappers.WrapperGeometry
+    arg2::GO.GI.Wrappers.WrapperGeometry
+    expected_result::T
+end
+
+Base.show(io::IO, ::MIME"text/plain", item::TestItem) = print(io, "TestItem(operation = $(item.operation), expects $(GI.trait(item.expected_result)))")
+Base.show(io::IO, item::TestItem) = show(io, MIME"text/plain"(), item)
+
+struct Case
+    description::String
+    geom_a::GO.GI.Wrappers.WrapperGeometry
+    geom_b::GO.GI.Wrappers.WrapperGeometry
+    items::Vector{TestItem}
+end
+
+function load_test_cases(filepath::String)
+    doc = read(filepath, XML.Node) # lazy parsing
+    run = only(children(doc))
+    test_cases = Case[]
+    for case in children(run)
+        if tag(case) != "case"
+            continue
+        end
+        push!(test_cases, parse_case(case))
+    end
+    return test_cases
+end
+
+function parse_case(case::XML.Node)
+    description = value(only(children(case.children[1])))
+    a = jts_wkt_to_geom(value(only(children(case.children[2]))))
+    b = jts_wkt_to_geom(value(only(children(case.children[3]))))
+
+    items = TestItem[]
+    for item in children(case)[4:end]
+        ops = children(item)
+        for op in ops
+            op_attrs = XML.attributes(op)
+            operation = op_attrs["name"]
+            arg1 = op_attrs["arg1"]
+            arg2 = op_attrs["arg2"]
+            expected_result = jts_wkt_to_geom(value(only(op.children)))
+            push!(items, TestItem(operation, lowercase(arg1) == "a" ? a : b, lowercase(arg2) == "a" ? a : b, expected_result))
+        end
+    end
+    return Case(description, a, b, items)
+end
+
+
+testfile = "/Users/anshul/.julia/dev/geo/jts/modules/tests/src/test/resources/testxml/general/TestOverlayAA.xml"
+
+cases = load_test_cases(testfile)
+
+using Test
+
+for case in cases
+    @testset "$(case.description)" begin
+        for item in case.items
+            @testset "$(item.operation)" begin
+                result = if item.operation == "union"
+                     GO.union(item.arg1, item.arg2; target = GO.PolygonTrait())
+                elseif item.operation == "difference"
+                    GO.difference(item.arg1, item.arg2; target = GO.PolygonTrait())
+                elseif item.operation == "intersection"
+                    GO.intersection(item.arg1, item.arg2; target = GO.PolygonTrait())
+                elseif item.operation == "symdifference"
+                    continue
+                else
+                    continue
+                end
+
+                finalresult = if length(result) == 0
+                    nothing
+                elseif length(result) == 1
+                    only(result)
+                else
+                    GI.MultiPolygon(result)
+                end
+
+                if isnothing(finalresult)
+                    @warn("No result")
+                    continue
+                end
+
+                if GI.geomtrait(item.expected_result) isa Union{GI.MultiPolygonTrait, GI.PolygonTrait}
+                    difference_in_areas = GO.area(GO.difference(finalresult, item.expected_result; target = GO.PolygonTrait()))
+                    if difference_in_areas > 1
+                        @warn("Difference in areas: $(difference_in_areas)")
+                        f, a, p = poly(finalresult; label = "Final result (GO)", axis = (; title = "$(case.description) - $(item.operation)"))
+                        poly!(a, item.expected_result; label = "Expected result (JTS)")
+                        display(f)
+                    end
+                    # @test difference_in_areas < 1
+                else
+                    @test_broken GO.equals(finalresult, item.expected_result)
+                end
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/external/tg/tg_testset_reader.jl b/test/external/tg/tg_testset_reader.jl
new file mode 100644
index 000000000..52f9cc4f3
--- /dev/null
+++ b/test/external/tg/tg_testset_reader.jl
@@ -0,0 +1,65 @@
+using JSON3
+import GeoFormatTypes as GFT, GeoInterface as GI, GeometryOps as GO
+import TGGeometry
+import LibGEOS
+import WellKnownGeometry
+
+
+using Test
+
+TG_PRED_SYMBOL_TO_FUNCTION = Dict(
+    [sym => getproperty(GO, sym) for sym in TGGeometry.TG_PREDICATES]
+)
+
+TG_IGNORE_LIST = Set([:crosses, :overlaps, :equals])
+
+testsets = JSON3.read(read(joinpath(@__DIR__, "data", "mpoly.jsonc")))
+
+# # assume GEOS is the ultimate source of truth
+
+testset = testsets[4]
+geoms = GO.tuples.(GFT.WellKnownText.((GFT.Geom(),), testset.geoms))
+
+
+
+function run_testsets(json_file_path)
+    testsets = JSON3.read(read(json_file_path))
+    for i in eachindex(testsets)
+        testset = testsets[i]
+        geoms = GO.tuples.(GFT.WellKnownText.((GFT.Geom(),), testset.geoms))
+        @testset let testset_index = i
+            for (predname, results) in testset.predicates
+                !haskey(TG_PRED_SYMBOL_TO_FUNCTION, predname) && continue
+                predname in TG_IGNORE_LIST && continue
+                predicate_f = TG_PRED_SYMBOL_TO_FUNCTION[predname]
+
+                @testset let predicate = predname
+                    go_result = try
+                        predicate_f(geoms[1], geoms[2])
+                    catch e
+                        println("Error: $e")
+                        @test_broken false
+                        continue
+                    end
+                    lg_result = predicate_f(GO.GEOS(), geoms[1], geoms[2])
+                    # tg_result = predicate_f(GO.TG(), geoms[1], geoms[2])
+
+                    expected = first(results) == "T"
+
+                    @test go_result == expected
+                    @test lg_result == expected
+                    # @test tg_result == expected
+                end
+            end
+        end
+    end
+end
+
+@testset "All TG testsets" begin
+    for file in filter(endswith(".jsonc"), readdir(joinpath(@__DIR__, "data"); join = true))
+        filename = splitext(basename(file))[1]
+        @testset "$filename" begin
+            run_testsets(file)
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 07b08fea2..3d5465c96 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,7 @@ include("helpers.jl")
     @safetestset "Manifold" begin include("core/manifold.jl") end
 end
 
+@safetestset "Types" begin include("types.jl") end
 @safetestset "Primitives" begin include("primitives.jl") end
 # Methods
 @safetestset "Angles" begin include("methods/angles.jl") end
@@ -39,3 +40,4 @@ end
 # Extensions
 @safetestset "FlexiJoins" begin include("extensions/flexijoins.jl") end
 @safetestset "LibGEOS" begin include("extensions/libgeos.jl") end
+@safetestset "TGGeometry" begin include("extensions/tggeometry.jl") end
\ No newline at end of file
diff --git a/test/types.jl b/test/types.jl
new file mode 100644
index 000000000..9e18e6259
--- /dev/null
+++ b/test/types.jl
@@ -0,0 +1,148 @@
+using Test
+using GeometryOps
+using GeometryOpsCore
+
+using GeometryOps: Planar, Spherical, AutoManifold, rebuild, manifold, best_manifold, enforce, get
+using GeometryOps: CLibraryPlanarAlgorithm, GEOS, TG, PROJ
+
+@testset "CLibraryPlanarAlgorithm" begin
+    # Test that it's a subtype of SingleManifoldAlgorithm{Planar}
+    @test CLibraryPlanarAlgorithm <: GeometryOpsCore.SingleManifoldAlgorithm{Planar}
+
+    # Test that it requires params::NamedTuple
+    struct TestAlg <: CLibraryPlanarAlgorithm
+        params::NamedTuple
+    end
+
+    # Test constructor with keyword arguments
+    alg = TestAlg(; a=1, b=2)
+    @test alg.params == (; a=1, b=2)
+
+    # Test constructor with NamedTuple
+    alg = TestAlg((; a=1, b=2))
+    @test alg.params == (; a=1, b=2)
+
+    # Test null constructor
+    alg = TestAlg()
+    @test alg.params == NamedTuple()
+
+    # Test manifold methods
+    @test manifold(alg) == Planar()
+    @test best_manifold(alg, nothing) == Planar()
+
+    # Test rebuild methods
+    @test rebuild(alg, Planar()) == TestAlg(alg.params)
+    @test rebuild(alg, AutoManifold()) == TestAlg(alg.params)
+    @test_throws GeometryOpsCore.WrongManifoldException rebuild(alg, Spherical())
+    @test rebuild(alg, (; c=3, d=4)) == TestAlg((; c=3, d=4))
+    @test rebuild(alg; c=3, d=4) == TestAlg((; c=3, d=4))
+
+    # Test get methods
+    @test Base.get(alg, :a, 0) == 0
+    @test Base.get(alg, :c, 0) == 0
+
+    # Test enforce method
+    @test_throws ErrorException enforce(alg, :a, "test")
+    @test_throws ErrorException enforce(alg, :c, "test")
+end
+
+@testset "GEOS" begin
+    # Test that it's a subtype of CLibraryPlanarAlgorithm
+    @test GEOS() isa CLibraryPlanarAlgorithm
+
+    # Test null constructor
+    alg = GEOS()
+    @test alg.params == NamedTuple()
+    @test manifold(alg) == Planar()
+
+    # Test constructor
+    alg = GEOS(; a=1, b=2)
+    @test alg.params == (; a=1, b=2)
+
+    # Test manifold methods
+    @test manifold(alg) == Planar()
+    @test best_manifold(alg, nothing) == Planar()
+
+    # Test rebuild methods
+    @test rebuild(alg, Planar()) == GEOS(alg.params)
+    @test rebuild(alg, AutoManifold()) == GEOS(alg.params)
+    @test_throws GeometryOpsCore.WrongManifoldException rebuild(alg, Spherical())
+    @test rebuild(alg, (; c=3, d=4)) == GEOS((; c=3, d=4))
+    @test rebuild(alg; c=3, d=4) == GEOS((; c=3, d=4))
+
+    # Test get methods
+    @test Base.get(alg, :a, 0) == 1
+    @test Base.get(alg, :c, 0) == 0
+
+    # Test enforce method
+    @test enforce(alg, :a, "test") == 1
+    @test_throws ErrorException enforce(alg, :c, "test")
+end
+
+@testset "TG" begin
+    # Test that it's a subtype of CLibraryPlanarAlgorithm
+    @test TG <: CLibraryPlanarAlgorithm
+
+    # Test null constructor
+    alg = TG()
+    @test alg.params == NamedTuple()
+    @test manifold(alg) == Planar()
+
+    # Test constructor
+    alg = TG(; a=1, b=2)
+    @test alg.params == (; a=1, b=2)
+
+    # Test manifold methods
+    @test manifold(alg) == Planar()
+    @test best_manifold(alg, nothing) == Planar()
+
+    # Test rebuild methods
+    @test rebuild(alg, Planar()) == TG(alg.params)
+    @test rebuild(alg, AutoManifold()) == TG(alg.params)
+    @test_throws GeometryOpsCore.WrongManifoldException rebuild(alg, Spherical())
+    @test rebuild(alg, (; c=3, d=4)) == TG((; c=3, d=4))
+    @test rebuild(alg; c=3, d=4) == TG((; c=3, d=4))
+
+    # Test get methods
+    @test Base.get(alg, :a, 0) == 1
+    @test Base.get(alg, :c, 0) == 0
+
+    # Test enforce method
+    @test enforce(alg, :a, "test") == 1
+    @test_throws ErrorException enforce(alg, :c, "test")
+end
+
+@testset "PROJ" begin
+    # Test that it's a subtype of Algorithm
+    @test PROJ <: GeometryOpsCore.Algorithm
+
+    # Test null constructor
+    alg = PROJ()
+    @test alg.manifold == Planar()
+    @test alg.params == NamedTuple()
+
+    # Test constructors
+    alg1 = PROJ(; a=1, b=2)
+    @test alg1.manifold == Planar()
+    @test alg1.params == (; a=1, b=2)
+
+    alg2 = PROJ(Spherical())
+    @test alg2.manifold == Spherical()
+    @test alg2.params == NamedTuple()
+
+    # Test manifold method
+    @test manifold(alg1) == Planar()
+    @test manifold(alg2) == Spherical()
+
+    # Test rebuild methods
+    @test rebuild(alg1, Spherical()) == PROJ(Spherical(), alg1.params)
+    @test rebuild(alg1, (; c=3, d=4)) == PROJ(Planar(), (; c=3, d=4))
+
+    # Test get methods
+    @test Base.get(alg1, :a, 0) == 1
+    @test Base.get(alg1, :c, 0) == 0
+
+    # Test enforce method
+    @test enforce(alg1, :a, "test") == 1
+    @test_throws ErrorException enforce(alg1, :c, "test")
+end 
\ No newline at end of file