Add a TGGeometry extension for fast inexact geometric predicates

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  

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"
@@ -67,4 +71,4 @@ Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 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", "Test"]

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

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)