Make the watershed benchmarks use GeometryOps methods

[asinghvi17]

Compact code for the watershed benchmark is:

import GeometryOps as GO, 
    GeoInterface as GI, 
    GeometryBasics as GB, 
    LibGEOS as LG
import GeoJSON
# 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
# We use CairoMakie to visualize our results!
using CairoMakie, MakieThemes, GeoInterfaceMakie
# Finally, we import some general utility packages:
using Statistics, CoordinateTransformations


# 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"`
# TODO: this doesn't work with WellKnownGeometry.  Why?
wkt_geoms = LG.readgeom.(eachline("watersheds.wkt.txt"), (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)
@be GO.union($(go_vp[1]), $(go_vp[2]); target = $GI.PolygonTrait()) seconds=3
@be LG.union($(vancouver_polygons[1]), $(vancouver_polygons[2])) seconds=3

go_u = GO.union((go_vp[1]), (go_vp[2]); target = GI.PolygonTrait())[1] 
lg_u = LG.union((vancouver_polygons[1]), (vancouver_polygons[2]))

GO.equals(go_u, lg_u)


all_intersected = falses(length(go_vp))
accumulator = deepcopy(go_vp[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(go_vp; color = rand(RGBf, length(go_vp))))
display(poly(accumulator))
intersect_states = [deepcopy(all_intersected)]
@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] || !(GO.intersects(go_vp[idx], accumulator)))
            println("Assimilating $idx")
            result = GO.union(go_vp[idx], accumulator; target = GI.PolygonTrait())
            lr=length(result)
            if lr > 1
                @warn "Result length > 1"
            end
            accumulator = result[1]
            all_intersected[idx] = true
        end
    end
    push!(intersect_states, deepcopy(all_intersected))#display(poly(go_vp[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 

Dpwnload the below file
watersheds.wkt.txt

[asinghvi17]

Currently failing with:

Began iteration 1
Assimilating 6
Assimilating 7
Assimilating 118
Assimilating 148
Assimilating 158
Assimilating 159
Assimilating 165
Assimilating 167
Assimilating 169
Assimilating 170
Assimilating 172
Assimilating 173
Assimilating 174
ERROR: BoundsError: attempt to access 71-element Vector{GeometryOps.PolyNode{Float64}} at index [0]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] _classify_crossing!(::Type{Float64}, a_list::Vector{GeometryOps.PolyNode{…}}, b_list::Vector{GeometryOps.PolyNode{…}})
   @ GeometryOps ~/.julia/dev/GeometryOps/src/methods/clipping/clipping_processor.jl:327
 [3] _build_ab_list(::Type{…}, poly_a::GeoInterface.Wrappers.LinearRing{…}, poly_b::GeoInterface.Wrappers.LinearRing{…}, delay_cross_f::typeof(GeometryOps._union_delay_cross_f), delay_bounce_f::typeof(GeometryOps._union_delay_bounce_f))
   @ GeometryOps ~/.julia/dev/GeometryOps/src/methods/clipping/clipping_processor.jl:53
 [4] _union(::GeometryOps.TraitTarget{…}, ::Type{…}, ::GeoInterface.PolygonTrait, poly_a::GeoInterface.Wrappers.Polygon{…}, ::GeoInterface.PolygonTrait, poly_b::GeoInterface.Wrappers.Polygon{…})
   @ GeometryOps ~/.julia/dev/GeometryOps/src/methods/clipping/union.jl:47
 [5] #union#141
   @ ~/.julia/dev/GeometryOps/src/methods/clipping/union.jl:32 [inlined]
 [6] union
   @ ~/.julia/dev/GeometryOps/src/methods/clipping/union.jl:29 [inlined]
 [7] macro expansion
   @ ~/.julia/dev/GeometryOps/benchmarks/benchmarks.jl:117 [inlined]
Some type information was truncated. Use `show(err)` to see complete types.

[skygering]

@asinghvi17 The Polygon formed isn't valid according to LibGEOS. Basically, the blue polygon is connected at that "neck" portion by one point.

import LibGEOS as LG
wkt_geoms = LG.readgeom.(eachline("/Users/skylargering/Downloads/merging_error.wkt"), (LG.WKTReader(LG.get_global_context()),))
LG.isValid(wkt_geoms[1])   # returns false

So technically that interior space filled by the yellow polygon should be a hole, but it isn't. My code doesn't know how to deal with that. Once you use LibGEOS.makeValid, my code works fine (with another small bug fix 😅 ).

Previous:

After:

I guess this actually points to an earlier issue within the code that created an invalid polygon. However, to me, it isn't even clear that a hole should be able to share a point with the edge. This seems equally wrong to me...

[asinghvi17]

Huh! I didn't think to check that. I guess that another merging order might have averted this?

Yeah I wouldn't think the hole should be able to intersect the edge, but I don't see what else could be done here (it is a bit of an edge case, but also not, since this is probably a primary usecase to merge polygons!)

[skygering]

This might also have to do with what are called "split" candidate edges that make sure that self-intersecting polygons aren't returned by the Foster algorithm. Perhaps this will be resolved with those. Regardless, I would like that PR and this fix to be tied together. Perhaps we can still merge the "glues edges" PR without resolving this for the time being.