Update benchmark code

- allow arbitrary legend kws
- import GeoFormatTypes, WellKnownGeometry, CoordinateTransformations, ProgressMeter
- make USA benchmarks print better
- add benchmarks for coverage union
  - It turns out JTS (and also GEOS) use a special purpose noder and algorithm for coverage unions.  We should also implement that at some stage - but it's a bit tough without the OverlayNG methodology.  Maybe we can get an AI to implement that here :D

Author: asinghvi17

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)

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})

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)