From debedf7a87fd0b492787db7fe40bd5fa7acc952c Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 3 Mar 2025 19:15:25 -0500 Subject: [PATCH] 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 --- benchmarks/benchmark_plots.jl | 4 +- benchmarks/benchmarks.jl | 234 ++++++++++++++++--------------- benchmarks/geometry_providers.jl | 4 +- benchmarks/polygon_scaling.jl | 2 +- 4 files changed, 129 insertions(+), 115 deletions(-) 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..5dd0b646e 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, @@ -57,87 +62,6 @@ end plot_trials(circle_area_suite) -# ## 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? -wkt_geoms = LG.readgeom.(eachline("/Users/anshul/Downloads/watersheds.wkt"), (LG.WKTReader(LG.get_global_context()),)) -vancouver_polygons = GI.getgeom.(wkt_geoms, 1); #.|> GO.tuples; - -import SortTileRecursiveTree as STR -tree = STR.STRtree(vancouver_polygons) -query_result = STR.query(tree, GI.extent(vancouver_polygons[1])) - -GO.intersects.((vancouver_polygons[1],), vancouver_polygons[query_result]) - -go_vp = GO.tuples.(vancouver_polygons[1:2]) -@be GO.union($(go_vp[1]), $(go_vp[2]); target = $GI.PolygonTrait()) -@be LG.union($(vancouver_polygons[1]), $(vancouver_polygons[2])) - -all_intersected = falses(length(vancouver_polygons)) -accumulator = deepcopy(vancouver_polygons[1]) -all_intersected[1] = true -i = 1 -# query_result = STR.query(tree, GI.extent(accumulator)) -# for idx in query_result -# println("Assessing $idx") -# if !all_intersected[idx] && LG.intersects(vancouver_polygons[idx], accumulator) -# println("Assimilating $idx") -# result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#) -# # @show length(result) -# accumulator = result#[1] -# all_intersected[idx] = true -# end -# end -display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected)))) -display(poly(accumulator)) -@time while !(all(all_intersected)) && i < length(vancouver_polygons) - println("Began iteration $i") - query_result = STR.query(tree, GI.extent(accumulator)) - for idx in query_result - if !(all_intersected[idx] || !(LG.intersects(vancouver_polygons[idx], accumulator))) - println("Assimilating $idx") - result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#) - # @show length(result) - accumulator = result#[1] - all_intersected[idx] = true - end - end - display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected)))) - println("Finished iteration $i") - # wait_for_key("Press any key to continue to the next iteration.") - i += 1 -end - -fig = Figure() -ax = Axis(fig[1, 1]; title = "STRTree query for polygons", aspect = DataAspect()) -for (i, result_index) in enumerate(query_result) - poly!(ax, vancouver_polygons[result_index]; color = Makie.wong_colors()[i], label = "$result_index") -end -Legend(fig[1, 2], ax) -fig - - - -# TODO: -# - Vancouver watersheds: -# - Intersection on pre-buffered geometry -# - Polygon union by reduction (perhaps pre-sort by border order, so we don't end up with useless polygons) -# - Queries using STRTree.jl -# - Potentially using a prepared geometry based approach to multithreaded reductive joining -# - Implement multipolygon joining. How? Query intersection or touching for each individual geometry, -# and implement a - - ## Segmentization @@ -155,7 +79,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) @@ -175,9 +99,9 @@ n_points_values = round.(Int, exp10.(LinRange(1, 4, 10))) circle_union_suite["LibGEOS"][n_points] = @be LG.union($lg_circle_left, $lg_circle_right) end -plot_trials(circle_difference_suite) -plot_trials(circle_intersection_suite) -plot_trials(circle_union_suite) +plot_trials(circle_difference_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) +plot_trials(circle_intersection_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) +plot_trials(circle_union_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) usa_poly_suite = BenchmarkGroup() usa_difference_suite = usa_poly_suite["difference"] = BenchmarkGroup(["title:USA difference", "subtitle:Tested on CONUS"]) @@ -197,29 +121,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)