From 926b915ff9ea7d4b9bf6a6a9fd42627511c911b6 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 18:51:03 -0400 Subject: [PATCH 01/30] Move utils to `utils/utils.jl` and add edge list functionality --- src/{ => utils}/utils.jl | 87 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) rename src/{ => utils}/utils.jl (60%) diff --git a/src/utils.jl b/src/utils/utils.jl similarity index 60% rename from src/utils.jl rename to src/utils/utils.jl index 540a7211e3..9ea887bf74 100644 --- a/src/utils.jl +++ b/src/utils/utils.jl @@ -126,3 +126,90 @@ function _point_in_extent(p, extent::Extents.Extent) (x1, x2), (y1, y2) = extent.X, extent.Y return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 end + +#= +# `eachedge`, `to_edgelist` + +These functions are used to decompose geometries into lists of edges. +Currently they only work on linear rings. +=# + +""" + eachedge(geom, [::Type{T}]) + +Decompose a geometry into a list of edges. +Currently only works for LineString and LinearRing. + +Returns some iterator, which yields tuples of points. Each tuple is an edge. + +It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. +""" +eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) + +function eachedge(geom, ::Type{T}) where T + eachedge(GI.trait(geom), geom, T) +end + +# implementation for LineString and LinearRing +function eachedge(trait::GI.AbstractCurveTrait, geom, ::Type{T}) where T + return (_tuple_point.((GI.getpoint(geom, i), GI.getpoint(geom, i+1)), T) for i in 1:GI.npoint(geom)-1) +end + +# implementation for Polygon, MultiPolygon, MultiLineString, GeometryCollection +function eachedge(trait::GI.AbstractGeometryTrait, geom, ::Type{T}) where T + return Iterators.flatten((eachedge(r, T) for r in flatten(GI.AbstractCurveTrait, geom))) +end + +function eachedge(trait::GI.PointTrait, geom, ::Type{T}) where T + return ArgumentError("Can't get edges from points, $geom was a PointTrait.") +end + +function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T + return ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.") +end + +""" + to_edgelist(geom, [::Type{T}]) + +Convert a geometry into a vector of `GI.Line` objects with attached extents. +""" +to_edgelist(geom, ::Type{T}) where T = + [_lineedge(ps, T) for ps in eachedge(geom, T)] +function to_edgelist(ext::E, geom, ::Type{T}) where {E<:Extents.Extent,T} + edges_in = eachedge(geom, T) + l1 = _lineedge(first(edges_in), T) + edges_out = typeof(l1)[] + indices = Int[] + for (i, ps) in enumerate(edges_in) + l = _lineedge(ps, T) + if Extents.intersects(ext, GI.extent(l)) + push!(edges_out, l) + push!(indices, i) + end + end + return edges_out, indices +end + +function _lineedge(ps::Tuple, ::Type{T}) where T + l = GI.Line(SVector{2,NTuple{2, T}}(ps)) # TODO: make this flexible in dimension + e = GI.extent(l) + return GI.Line(l.geom; extent=e) +end + +function lazy_edgelist(geom, ::Type{T}) where T + (_lineedge(ps, T) for ps in eachedge(geom, T)) +end + +function edge_extents(geom) + return [begin + Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) + end + for edge in eachedge(geom, Float64)] +end + +function lazy_edge_extents(geom) + return (begin + Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) + end + for edge in eachedge(geom, Float64)) +end From 838932a419aa453b360f6b5d9688dc81ad0c6f22 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:27:13 -0400 Subject: [PATCH 02/30] Add LoopStateMachine and tests --- src/GeometryOps.jl | 10 +- .../LoopStateMachine/LoopStateMachine.jl | 106 ++++++++++++++++++ test/utils/LoopStateMachine.jl | 43 +++++++ 3 files changed, 158 insertions(+), 1 deletion(-) create mode 100644 src/utils/LoopStateMachine/LoopStateMachine.jl create mode 100644 test/utils/LoopStateMachine.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c0d169f6eb..dd7d4fcba4 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -36,9 +36,17 @@ const Edge{T} = Tuple{TuplePoint{T},TuplePoint{T}} where T include("types.jl") include("primitives.jl") -include("utils.jl") include("not_implemented_yet.jl") +include("utils/utils.jl") + +include("utils/LoopStateMachine/LoopStateMachine.jl") +using .LoopStateMachine + +include("utils/SpatialTreeInterface/SpatialTreeInterface.jl") +using .SpatialTreeInterface + + include("methods/angles.jl") include("methods/area.jl") include("methods/barycentric.jl") diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl new file mode 100644 index 0000000000..3de8239a44 --- /dev/null +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -0,0 +1,106 @@ +""" + LoopStateMachine + +Utilities for returning state from functions that run inside a loop. + +This is used in e.g clipping, where we may need to break or transition states. + +The main entry point is to return an [`Action`](@ref) from a function that +is wrapped in a `@controlflow f(...)` macro in a loop. When a known `Action` +(currently, `:continue`, `:break`, `:return`, or `:full_return` actions) is returned, +it is processed by the `@controlflow` macro, which allows the function to break out of the loop +early, continue to the next iteration, or return a value, basically a way to provoke syntactic +behaviour from a function called from a inside a loop, where you do not have access to that loop. + +## Example + +```julia +``` +""" +module LoopStateMachine + +export Action, @controlflow + +import ..GeometryOps as GO + +const ALL_ACTION_DESCRIPTIONS = """ +- `:continue`: continue to the next iteration of the loop. + This is the `continue` keyword in Julia. The contents of the action are not used. +- `:break`: break out of the loop. + This is the `break` keyword in Julia. The contents of the action are not used. +- `:return`: cause the function executing the loop to return with the wrapped value. +- `:full_return`: cause the function executing the loop to return `Action(:full_return, x)`. + This is very useful to terminate recursive funtions, like tree queries terminating after you + have found a single intersecting segment. +""" + +""" + Action(name::Symbol, [x]) + +Create an `Action` with the name `name` and optional contents `x`. + +`Action`s are returned from functions wrapped in a `@controlflow` macro, which +does something based on the return value of that function if it is an `Action`. + +## Available actions + +$ALL_ACTION_DESCRIPTIONS +""" +struct Action{T} + name::Symbol + x::T +end + +Action() = Action{Nothing}(:unnamed, nothing) +Action(x::T) where T = Action{:unnamed, T}(:unnamed, x) +Action(x::Symbol) = Action(x, nothing) + +function Base.show(io::IO, action::Action{T}) where T + print(io, "Action ", action.name) + if isnothing(action.x) + print(io, "()") + else + print(io, "(", action.x, ")") + end +end + +""" + @controlflow f(...) + +Process the result of `f(...)` and return the result if it's not an `Action`(@ref LoopStateMachine.Action). + +If it is an `Action`, then process it according to the following rules, and throw an error if it's not recognized. +`:continue`, `:break`, `:return`, or `:full_return` are valid actions. + +$ALL_ACTION_DESCRIPTIONS + +!!! warning + Only use this inside a loop, otherwise you'll get a syntax error, especially if you use `:continue` or `:break`. + +## Examples +""" +macro controlflow(expr) + varname = gensym("loop-state-machine-returned-value") + return quote + $varname = $(esc(expr)) + if $varname isa Action + if $varname.name == :continue + continue + elseif $varname.name == :break + break + elseif $varname.name == :return + return $varname.x + elseif $varname.name == :full_return + return $varname + else + throw(ArgumentError("Unknown action: $($varname.name)")) + end + else + $varname + end + end +end + +# You can define more actions as you desire. + +end \ No newline at end of file diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl new file mode 100644 index 0000000000..2523597454 --- /dev/null +++ b/test/utils/LoopStateMachine.jl @@ -0,0 +1,43 @@ +using Test +using GeometryOps.LoopStateMachine: @controlflow, Continue, Break + +@testset "Continue action" begin + count = 0 + f(i) = begin + count += 1 + if i == 3 + return Continue() + end + count += 1 + end + for i in 1:5 + @controlflow f(i) + end + @test count == 9 # Adds 1 for each iteration, but skips second +1 on i=3 +end + +@testset "Break action" begin + count = 0 + function f(i) + count += 1 + if i == 3 + return Break() + end + count += 1 + end + for i in 1:5 + @controlflow f(i) + end + @test count == 5 # Counts up to i=3, adding 2 for i=1,2 and 1 for i=3 +end + +@testset "Return value" begin + results = Int[] + for i in 1:3 + val = @controlflow begin + i * 2 + end + push!(results, val) + end + @test results == [2, 4, 6] +end From 34191219280f7580394aabcf2aad423842fb2cce Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:31:13 -0400 Subject: [PATCH 03/30] test all actions --- test/utils/LoopStateMachine.jl | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 2523597454..4ea8dfe132 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -1,12 +1,12 @@ using Test -using GeometryOps.LoopStateMachine: @controlflow, Continue, Break +using GeometryOps.LoopStateMachine: @controlflow, Action @testset "Continue action" begin count = 0 f(i) = begin count += 1 if i == 3 - return Continue() + return Action(:continue) end count += 1 end @@ -21,7 +21,7 @@ end function f(i) count += 1 if i == 3 - return Break() + return Action(:break) end count += 1 end @@ -31,6 +31,24 @@ end @test count == 5 # Counts up to i=3, adding 2 for i=1,2 and 1 for i=3 end +@testset "Return action" begin + f(i) = for j in 1:3 + i == j && @controlflow Action(:return, i) + end + @test f(1) == 1 + @test f(2) == 2 + @test f(3) == 3 +end + +@testset "Full return action" begin + f(i) = for j in 1:3 + i == j && @controlflow Action(:full_return, i) + end + @test f(1) == Action(:full_return, 1) + @test f(2) == Action(:full_return, 2) + @test f(3) == Action(:full_return, 3) +end + @testset "Return value" begin results = Int[] for i in 1:3 From d6c1ffc2c3f77da5cbd33f8290baff9374a49478 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:42:32 -0400 Subject: [PATCH 04/30] add SpatialTreeInterface implementation (but no tests yet, need to add those!) --- Project.toml | 2 + src/utils/SpatialTreeInterface/README.md | 47 ++++ .../SpatialTreeInterface.jl | 244 ++++++++++++++++++ 3 files changed, 293 insertions(+) create mode 100644 src/utils/SpatialTreeInterface/README.md create mode 100644 src/utils/SpatialTreeInterface/SpatialTreeInterface.jl diff --git a/Project.toml b/Project.toml index 8ee10bf442..d0eaaa7ff5 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi ", "Rafael Schouten Date: Wed, 16 Apr 2025 22:55:55 -0400 Subject: [PATCH 05/30] Define `isspatialtree` for FlatNoTree --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 2af8cfcd63..d1a08c756f 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -6,7 +6,7 @@ import Extents import GeoInterface as GI import AbstractTrees -# public isspatialtree, getchild, nchild, child_indices_extents, node_extent +# public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent export query, do_query # ## Interface @@ -225,6 +225,7 @@ struct FlatNoTree{T} geometries::T end +isspatialtree(::Type{<: FlatNoTree}) = true isleaf(tree::FlatNoTree) = true # NOTE: use pairs instead of enumerate here, so that we can support From 7b1b99b9bd53114ed7ecc621cca38bdea4310f28 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:56:05 -0400 Subject: [PATCH 06/30] Add basic tests for FlatNoTree --- test/utils/SpatialTreeInterface.jl | 77 ++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 test/utils/SpatialTreeInterface.jl diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl new file mode 100644 index 0000000000..5961f39819 --- /dev/null +++ b/test/utils/SpatialTreeInterface.jl @@ -0,0 +1,77 @@ +using Test + +using GeometryOps.SpatialTreeInterface +using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent +using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using Extents +using GeoInterface: GeoInterface as GI + +using GeometryOps.SpatialTreeInterface: FlatNoTree + +# Test FlatNoTree implementation +@testset "FlatNoTree" begin + # Test with a simple vector of extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = FlatNoTree(extents) + + @testset "Basic interface" begin + @test isleaf(tree) + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + end + + @testset "child_indices_extents" begin + # Test that we get the correct indices and extents + indices_extents = collect(child_indices_extents(tree)) + @test length(indices_extents) == 3 + @test indices_extents[1] == (1, extents[1]) + @test indices_extents[2] == (2, extents[2]) + @test indices_extents[3] == (3, extents[3]) + end + + @testset "Query functionality" begin + # Test query with a predicate that matches all + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == [1, 2, 3] + + # Test query with a predicate that matches none + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + # Test query with a specific extent predicate + search_extent = Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == [1, 2] # Should match first two extents + end + + @testset "Dual query functionality" begin + # Create two trees for dual query testing + tree1 = FlatNoTree([ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)) + ]) + tree2 = FlatNoTree([ + Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)), + Extents.Extent(X=(1.5, 2.5), Y=(1.5, 2.5)) + ]) + + # Test dual query with a predicate that matches all + all_pred = (x, y) -> true + results = Tuple{Int, Int}[] + do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) + @test sort(results) == [(1,1), (1,2), (2,1), (2,2)] + + # Test dual query with a specific predicate + intersects_pred = (x, y) -> Extents.intersects(x, y) + results = Tuple{Int, Int}[] + do_dual_query((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) + @test sort(results) == [(1,1), (2,1), (2,2)] + end +end + From 7dd401409cfc719965997fd758c93d4581bb86b5 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 22:56:27 -0400 Subject: [PATCH 07/30] ignore call notes since we don't want to commit them --- docs/.gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/.gitignore b/docs/.gitignore index 30c80518cb..468c60bf22 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -2,4 +2,6 @@ node_modules/ build/ package-lock.json src/source/ -Manifest.toml \ No newline at end of file +Manifest.toml + +src/call_notes.md \ No newline at end of file From aa8ea4370a67fcc4f899cd7f6aa5b53842001363 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 16 Apr 2025 23:15:29 -0400 Subject: [PATCH 08/30] test all of LoopStateMachine --- src/utils/LoopStateMachine/LoopStateMachine.jl | 9 +++++---- test/utils/LoopStateMachine.jl | 14 ++++++++++++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl index 3de8239a44..e01eea06b0 100644 --- a/src/utils/LoopStateMachine/LoopStateMachine.jl +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -52,15 +52,16 @@ struct Action{T} end Action() = Action{Nothing}(:unnamed, nothing) -Action(x::T) where T = Action{:unnamed, T}(:unnamed, x) +Action(x::T) where T = Action{T}(:unnamed, x) Action(x::Symbol) = Action(x, nothing) function Base.show(io::IO, action::Action{T}) where T - print(io, "Action ", action.name) + print(io, "Action") + print(io, "(:$(action.name)") if isnothing(action.x) - print(io, "()") + print(io, ")") else - print(io, "(", action.x, ")") + print(io, ", ",action.x, ")") end end diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 4ea8dfe132..99cdd99a7e 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -59,3 +59,17 @@ end end @test results == [2, 4, 6] end + +@testset "Show" begin + @test sprint(print, Action(:continue)) == "Action(:continue)" + @test sprint(print, Action(:break)) == "Action(:break)" + @test sprint(print, Action(:return, 1)) == "Action(:return, 1)" + @test sprint(print, Action(:full_return, 1)) == "Action(:full_return, 1)" +end + +@testset "Unnamed action" begin + @test sprint(print, Action()) == "Action(:unnamed)" + @test sprint(print, Action(1)) == "Action(:unnamed, 1)" + @test sprint(print, Action(:x)) == "Action(:x)" +end + From e9fb51502aff28dfa705a17503490a6c81973fe3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:44:46 -0400 Subject: [PATCH 09/30] generic tests for STI with FlatNoTree and STR --- test/utils/SpatialTreeInterface.jl | 109 ++++++++++++++++++++++++----- 1 file changed, 93 insertions(+), 16 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 5961f39819..33bde5a49d 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -3,28 +3,34 @@ using Test using GeometryOps.SpatialTreeInterface using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI +using SortTileRecursiveTree: STRtree -using GeometryOps.SpatialTreeInterface: FlatNoTree - -# Test FlatNoTree implementation -@testset "FlatNoTree" begin - # Test with a simple vector of extents - extents = [ - Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), - Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), - Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) - ] - tree = FlatNoTree(extents) - +# Generic test functions for spatial trees +function test_basic_interface(TreeType) @testset "Basic interface" begin + # Create a simple tree with one extent + extents = [Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0))] + tree = TreeType(extents) + @test isleaf(tree) @test isspatialtree(tree) @test isspatialtree(typeof(tree)) end +end +function test_child_indices_extents(TreeType) @testset "child_indices_extents" begin + # Create a tree with three extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = TreeType(extents) + # Test that we get the correct indices and extents indices_extents = collect(child_indices_extents(tree)) @test length(indices_extents) == 3 @@ -32,8 +38,18 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree @test indices_extents[2] == (2, extents[2]) @test indices_extents[3] == (3, extents[3]) end +end +function test_query_functionality(TreeType) @testset "Query functionality" begin + # Create a tree with three extents + extents = [ + Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)), + Extents.Extent(X=(2.0, 3.0), Y=(2.0, 3.0)) + ] + tree = TreeType(extents) + # Test query with a predicate that matches all all_pred = x -> true results = query(tree, all_pred) @@ -49,14 +65,16 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree results = query(tree, Base.Fix1(Extents.intersects, search_extent)) @test sort(results) == [1, 2] # Should match first two extents end +end +function test_dual_query_functionality(TreeType) @testset "Dual query functionality" begin - # Create two trees for dual query testing - tree1 = FlatNoTree([ + # Create two trees with overlapping extents + tree1 = TreeType([ Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)) ]) - tree2 = FlatNoTree([ + tree2 = TreeType([ Extents.Extent(X=(0.5, 1.5), Y=(0.5, 1.5)), Extents.Extent(X=(1.5, 2.5), Y=(1.5, 2.5)) ]) @@ -65,7 +83,7 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree all_pred = (x, y) -> true results = Tuple{Int, Int}[] do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) - @test sort(results) == [(1,1), (1,2), (2,1), (2,2)] + @test length(results) == 4 # 2 points in tree1 * 2 points in tree2 # Test dual query with a specific predicate intersects_pred = (x, y) -> Extents.intersects(x, y) @@ -75,3 +93,62 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree end end +function test_geometry_support(TreeType) + @testset "Geometry support" begin + # Create a tree with 100 points + points = tuple.(1:100, 1:100) + tree = TreeType(points) + + # Test basic interface + @test isleaf(tree) + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + + # Test child indices and extents + indices_extents = collect(child_indices_extents(tree)) + @test length(indices_extents) == 100 + + # Check first and last points + first_idx, first_extent = indices_extents[1] + last_idx, last_extent = indices_extents[100] + + @test first_idx == 1 + @test last_idx == 100 + @test first_extent == Extents.Extent(X=(1.0, 1.0), Y=(1.0, 1.0)) + @test last_extent == Extents.Extent(X=(100.0, 100.0), Y=(100.0, 100.0)) + + # Test query functionality + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == collect(1:100) + + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == collect(45:55) + end +end + +# Test FlatNoTree implementation +@testset "FlatNoTree" begin + test_basic_interface(FlatNoTree) + test_child_indices_extents(FlatNoTree) + test_query_functionality(FlatNoTree) + test_dual_query_functionality(FlatNoTree) + test_geometry_support(FlatNoTree) +end + +# Test STRtree implementation +@testset "STRtree" begin + test_basic_interface(STRtree) + test_child_indices_extents(STRtree) + test_query_functionality(STRtree) + test_dual_query_functionality(STRtree) + test_geometry_support(STRtree) +end + + + From 9efde48e7c88b84a4493f9cee2dcc32e6b744b63 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:44:54 -0400 Subject: [PATCH 10/30] actually include tests --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 1 + test/runtests.jl | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index d1a08c756f..ac93a7482a 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -8,6 +8,7 @@ import AbstractTrees # public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent export query, do_query +export FlatNoTree # ## Interface # Interface definition for spatial tree types. diff --git a/test/runtests.jl b/test/runtests.jl index 3d5465c96f..829b2e545c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,11 @@ end @safetestset "Types" begin include("types.jl") end @safetestset "Primitives" begin include("primitives.jl") end + +# Utils +@safetestset "LoopStateMachine" begin include("utils/LoopStateMachine.jl") end +@safetestset "SpatialTreeInterface" begin include("utils/SpatialTreeInterface.jl") end + # Methods @safetestset "Angles" begin include("methods/angles.jl") end @safetestset "Area" begin include("methods/area.jl") end From e5ad556c80a059ac97853af19801a30bbbcd2b7e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:49:18 -0400 Subject: [PATCH 11/30] implement `isspatialtree` for STR types --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index ac93a7482a..74913eec79 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -200,18 +200,20 @@ end using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode +isspatialtree(::Type{<: STRtree}) = true nchild(tree::STRtree) = nchild(tree.rootnode) getchild(tree::STRtree) = getchild(tree.rootnode) getchild(tree::STRtree, i) = getchild(tree.rootnode, i) isleaf(tree::STRtree) = isleaf(tree.rootnode) child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) - +isspatialtree(::Type{<: STRNode}) = true nchild(node::STRNode) = length(node.children) getchild(node::STRNode) = node.children getchild(node::STRNode, i) = node.children[i] isleaf(node::STRNode) = false # STRNodes are not leaves by definition +isspatialtree(::Type{<: STRLeafNode}) = true isleaf(node::STRLeafNode) = true child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) From 4ccf67bd790eeb3fd0781c46fdcec9924bfad7d2 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 07:51:30 -0400 Subject: [PATCH 12/30] fix tests --- test/utils/SpatialTreeInterface.jl | 17 +---------------- 1 file changed, 1 insertion(+), 16 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 33bde5a49d..50b5f25d23 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -14,8 +14,7 @@ function test_basic_interface(TreeType) # Create a simple tree with one extent extents = [Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0))] tree = TreeType(extents) - - @test isleaf(tree) + @test isspatialtree(tree) @test isspatialtree(typeof(tree)) end @@ -100,23 +99,9 @@ function test_geometry_support(TreeType) tree = TreeType(points) # Test basic interface - @test isleaf(tree) @test isspatialtree(tree) @test isspatialtree(typeof(tree)) - # Test child indices and extents - indices_extents = collect(child_indices_extents(tree)) - @test length(indices_extents) == 100 - - # Check first and last points - first_idx, first_extent = indices_extents[1] - last_idx, last_extent = indices_extents[100] - - @test first_idx == 1 - @test last_idx == 100 - @test first_extent == Extents.Extent(X=(1.0, 1.0), Y=(1.0, 1.0)) - @test last_extent == Extents.Extent(X=(100.0, 100.0), Y=(100.0, 100.0)) - # Test query functionality all_pred = x -> true results = query(tree, all_pred) From 0900db2a1f9f87ffc4f38fd6271e225094d5a7dc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 09:29:33 -0400 Subject: [PATCH 13/30] Apply suggestions from code review Co-authored-by: Rafael Schouten --- src/GeometryOps.jl | 6 ++---- src/utils/utils.jl | 5 ----- 2 files changed, 2 insertions(+), 9 deletions(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index dd7d4fcba4..db7efd2ee4 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -39,12 +39,10 @@ include("primitives.jl") include("not_implemented_yet.jl") include("utils/utils.jl") - include("utils/LoopStateMachine/LoopStateMachine.jl") -using .LoopStateMachine - include("utils/SpatialTreeInterface/SpatialTreeInterface.jl") -using .SpatialTreeInterface + +using .LoopStateMachine, .SpatialTreeInterface include("methods/angles.jl") diff --git a/src/utils/utils.jl b/src/utils/utils.jl index 9ea887bf74..fd882768c1 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -145,25 +145,20 @@ Returns some iterator, which yields tuples of points. Each tuple is an edge. It goes `(p1, p2), (p2, p3), (p3, p4), ...` etc. """ eachedge(geom) = eachedge(GI.trait(geom), geom, Float64) - function eachedge(geom, ::Type{T}) where T eachedge(GI.trait(geom), geom, T) end - # implementation for LineString and LinearRing function eachedge(trait::GI.AbstractCurveTrait, geom, ::Type{T}) where T return (_tuple_point.((GI.getpoint(geom, i), GI.getpoint(geom, i+1)), T) for i in 1:GI.npoint(geom)-1) end - # implementation for Polygon, MultiPolygon, MultiLineString, GeometryCollection function eachedge(trait::GI.AbstractGeometryTrait, geom, ::Type{T}) where T return Iterators.flatten((eachedge(r, T) for r in flatten(GI.AbstractCurveTrait, geom))) end - function eachedge(trait::GI.PointTrait, geom, ::Type{T}) where T return ArgumentError("Can't get edges from points, $geom was a PointTrait.") end - function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T return ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.") end From 06c76a93b831ebf83cb87c1cc1e8215ce108a1cc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 10:29:24 -0400 Subject: [PATCH 14/30] Implement and test a custom exception type for unknown actions in LSM --- src/utils/LoopStateMachine/LoopStateMachine.jl | 18 +++++++++++++++++- test/utils/LoopStateMachine.jl | 15 +++++++++++++++ 2 files changed, 32 insertions(+), 1 deletion(-) diff --git a/src/utils/LoopStateMachine/LoopStateMachine.jl b/src/utils/LoopStateMachine/LoopStateMachine.jl index e01eea06b0..77b62d4501 100644 --- a/src/utils/LoopStateMachine/LoopStateMachine.jl +++ b/src/utils/LoopStateMachine/LoopStateMachine.jl @@ -65,6 +65,21 @@ function Base.show(io::IO, action::Action{T}) where T end end +struct UnrecognizedActionException <: Base.Exception + name::Symbol +end + +function Base.showerror(io::IO, e::UnrecognizedActionException) + print(io, "Unrecognized action: ") + printstyled(io, e.name; color = :red, bold = true) + println(io, ".") + println(io, "Valid actions are:") + println(io, ALL_ACTION_DESCRIPTIONS) +end + +# We exclude the macro definition from code coverage computations, +# because I know it's tested but Codecov doesn't seem to think so. +# COV_EXCL_START """ @controlflow f(...) @@ -94,13 +109,14 @@ macro controlflow(expr) elseif $varname.name == :full_return return $varname else - throw(ArgumentError("Unknown action: $($varname.name)")) + throw(UnrecognizedActionException($varname.name)) end else $varname end end end +# COV_EXCL_STOP # You can define more actions as you desire. diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 99cdd99a7e..30d4979c81 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -1,4 +1,5 @@ using Test +import GeometryOps.LoopStateMachine using GeometryOps.LoopStateMachine: @controlflow, Action @testset "Continue action" begin @@ -73,3 +74,17 @@ end @test sprint(print, Action(:x)) == "Action(:x)" end +@testset "Unknown action" begin + f(i) = Action(:unknown_action, i) + @test_throws LoopStateMachine.UnrecognizedActionException begin + for i in 1:3 + @controlflow f(i) + end + end + + @test_throws ":continue" begin + for i in 1:3 + @controlflow f(i) + end + end +end From 728befcf3649f43dc51506081240e4efda41cede Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 10:46:02 -0400 Subject: [PATCH 15/30] Factor out SpatialTreeInterface into separate files Maybe this is a bad idea, single file was nice for readability. But we can keep going. --- .../SpatialTreeInterface.jl | 211 +----------------- .../depth_first_search.jl | 29 +++ .../dual_depth_first_search.jl | 56 +++++ .../SpatialTreeInterface/implementations.jl | 55 +++++ src/utils/SpatialTreeInterface/interface.jl | 85 +++++++ 5 files changed, 236 insertions(+), 200 deletions(-) create mode 100644 src/utils/SpatialTreeInterface/depth_first_search.jl create mode 100644 src/utils/SpatialTreeInterface/dual_depth_first_search.jl create mode 100644 src/utils/SpatialTreeInterface/implementations.jl create mode 100644 src/utils/SpatialTreeInterface/interface.jl diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 74913eec79..2cc15b6c6f 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -10,118 +10,19 @@ import AbstractTrees export query, do_query export FlatNoTree -# ## Interface -# Interface definition for spatial tree types. -# There is no abstract supertype here since it's impossible to enforce, -# but we do have a few methods that are common to all spatial tree types. +# The spatial tree interface and its implementations are defined here. +include("interface.jl") +include("implementations.jl") -""" - isspatialtree(tree)::Bool - -Return true if the object is a spatial tree, false otherwise. - -## Implementation notes - -For type stability, if your spatial tree type is `MyTree`, you should define -`isspatialtree(::Type{MyTree}) = true`, and `isspatialtree(::MyTree)` will forward -to that method automatically. -""" -isspatialtree(::T) where T = isspatialtree(T) -isspatialtree(::Type{<: Any}) = false - - -""" - getchild(node) - -Return an iterator over all the children of a node. -This may be materialized if necessary or available, -but can also be lazy (like a generator). -""" -getchild(node) = AbstractTrees.children(node) - -""" - getchild(node, i) - -Return the `i`-th child of a node. -""" -getchild(node, i) = getchild(node)[i] - -""" - nchild(node) +# Here we have some algorithms that use the spatial tree interface. +# The first file holds a single depth-first search, i.e., a single-tree query. +include("depth_first_search.jl") -Return the number of children of a node. -""" -nchild(node) = length(getchild(node)) - -""" - isleaf(node) - -Return true if the node is a leaf node, i.e., there are no "children" below it. -[`getchild`](@ref) should still work on leaf nodes, though, returning an iterator over the extents stored in the node - and similarly for `getnodes.` -""" -isleaf(node) = error("isleaf is not implemented for node type $(typeof(node))") - -""" - child_indices_extents(node) - -Return an iterator over the indices and extents of the children of a node. - -Each value of the iterator should take the form `(i, extent)`. - -This can only be invoked on leaf nodes! -""" -function child_indices_extents(node) - return zip(1:nchild(node), getchild(node)) -end - -""" - node_extent(node) - -Return the extent like object of the node. -Falls back to `GI.extent` by default, which falls back -to `Extents.extent`. - -Generally, defining `Extents.extent(node)` is sufficient here, and you -won't need to define this - -The reason we don't use that directly is to give users of this interface -a way to define bounding boxes that are not extents, like spherical caps -and other such things. -""" -node_extent(node) = GI.extent(node) - -# ## Query functions -# These are generic functions that work with any spatial tree type that implements the interface. - - -""" - do_query(f, predicate, tree) - -Call `f(i)` for each index `i` in the tree that satisfies `predicate(extent(i))`. - -This is generic to anything that implements the SpatialTreeInterface, particularly the methods -[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). -""" -function do_query(f::F, predicate::P, node::N) where {F, P, N} - if isleaf(node) - for (i, leaf_geometry_extent) in child_indices_extents(node) - if predicate(leaf_geometry_extent) - @controlflow f(i) - end - end - else - for child in getchild(node) - if predicate(node_extent(child)) - @controlflow do_query(f, predicate, child) - end - end - end -end -function do_query(predicate, node) - a = Int[] - do_query(Base.Fix1(push!, a), predicate, node) - return a -end +# The second file holds a dual depth-first search, i.e., a dual-tree query. +# This iterates over two trees simultaneously, and is substantially more efficient +# than two separate single-tree queries since it can prune branches in tandem as it +# descends into the trees. +include("dual_depth_first_search.jl") """ @@ -155,94 +56,4 @@ sanitize_predicate(::GI.AbstractTrait, pred) = sanitize_predicate(GI.extent(pred sanitize_predicate(pred::Extents.Extent) = Base.Fix1(Extents.intersects, pred) -""" - do_dual_query(f, predicate, node1, node2) - -Call `f(i1, i2)` for each index `i1` in `node1` and `i2` in `node2` that satisfies `predicate(extent(i1), extent(i2))`. - -This is generic to anything that implements the SpatialTreeInterface, particularly the methods -[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). -""" -function do_dual_query(f::F, predicate::P, node1::N1, node2::N2) where {F, P, N1, N2} - if isleaf(node1) && isleaf(node2) - # both nodes are leaves, so we can just iterate over the indices and extents - for (i1, extent1) in child_indices_extents(node1) - for (i2, extent2) in child_indices_extents(node2) - if predicate(extent1, extent2) - @controlflow f(i1, i2) - end - end - end - elseif isleaf(node1) # node2 is not a leaf, node1 is - recurse further into node2 - for child in getchild(node2) - if predicate(node_extent(node1), node_extent(child)) - @controlflow do_dual_query(f, predicate, node1, child) - end - end - elseif isleaf(node2) # node1 is not a leaf, node2 is - recurse further into node1 - for child in getchild(node1) - if predicate(node_extent(child), node_extent(node2)) - @controlflow do_dual_query(f, predicate, child, node2) - end - end - else # neither node is a leaf, recurse into both children - for child1 in getchild(node1) - for child2 in getchild(node2) - if predicate(node_extent(child1), node_extent(child2)) - @controlflow do_dual_query(f, predicate, child1, child2) - end - end - end - end -end - -# Finally, here's a sample implementation of the interface for STRtrees - -using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode - -isspatialtree(::Type{<: STRtree}) = true -nchild(tree::STRtree) = nchild(tree.rootnode) -getchild(tree::STRtree) = getchild(tree.rootnode) -getchild(tree::STRtree, i) = getchild(tree.rootnode, i) -isleaf(tree::STRtree) = isleaf(tree.rootnode) -child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) - -isspatialtree(::Type{<: STRNode}) = true -nchild(node::STRNode) = length(node.children) -getchild(node::STRNode) = node.children -getchild(node::STRNode, i) = node.children[i] -isleaf(node::STRNode) = false # STRNodes are not leaves by definition - -isspatialtree(::Type{<: STRLeafNode}) = true -isleaf(node::STRLeafNode) = true -child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) - - -""" - FlatNoTree(iterable_of_geoms_or_extents) - -Represents a flat collection with no tree structure, i.e., a brute force search. -This is cost free, so particularly useful when you don't want to build a tree! -""" -struct FlatNoTree{T} - geometries::T -end - -isspatialtree(::Type{<: FlatNoTree}) = true -isleaf(tree::FlatNoTree) = true - -# NOTE: use pairs instead of enumerate here, so that we can support -# iterators or collections that define custom `pairs` methods. -# This includes things like filtered extent lists, for example, -# so we can perform extent thinning with no allocations. -function child_indices_extents(tree::FlatNoTree{T}) where T - # This test only applies at compile time and should be optimized away in any case. - # And we can use multiple dispatch to override anyway, but it should be cost free I think. - if applicable(Base.keys, T) - return ((i, GI.extent(obj)) for (i, obj) in pairs(tree.geometries)) - else - return ((i, GI.extent(obj)) for (i, obj) in enumerate(tree.geometries)) - end -end - end # module SpatialTreeInterface \ No newline at end of file diff --git a/src/utils/SpatialTreeInterface/depth_first_search.jl b/src/utils/SpatialTreeInterface/depth_first_search.jl new file mode 100644 index 0000000000..dc1be59b3a --- /dev/null +++ b/src/utils/SpatialTreeInterface/depth_first_search.jl @@ -0,0 +1,29 @@ + +""" + depth_first_search(f, predicate, tree) + +Call `f(i)` for each index `i` in the tree that satisfies `predicate(extent(i))`. + +This is generic to anything that implements the SpatialTreeInterface, particularly the methods +[`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). +""" +function depth_first_search(f::F, predicate::P, node::N) where {F, P, N} + if isleaf(node) + for (i, leaf_geometry_extent) in child_indices_extents(node) + if predicate(leaf_geometry_extent) + @controlflow f(i) + end + end + else + for child in getchild(node) + if predicate(node_extent(child)) + @controlflow depth_first_search(f, predicate, child) + end + end + end +end +function depth_first_search(predicate, node) + a = Int[] + depth_first_search(Base.Fix1(push!, a), predicate, node) + return a +end diff --git a/src/utils/SpatialTreeInterface/dual_depth_first_search.jl b/src/utils/SpatialTreeInterface/dual_depth_first_search.jl new file mode 100644 index 0000000000..2c63d0aa11 --- /dev/null +++ b/src/utils/SpatialTreeInterface/dual_depth_first_search.jl @@ -0,0 +1,56 @@ +""" + dual_depth_first_search(f, predicate, tree1, tree2) + +Executes a dual depth-first search over two trees, descending into the children of +nodes `i` and `j` when `predicate(node_extent(i), node_extent(j))` is true, +and pruning that branch when `predicate(node_extent(i), node_extent(j))` is false. + +Finally, calls `f(i1, i2)` for each leaf-level index `i1::Int` in `tree1` and `i2::Int` in `tree2` +that satisfies `predicate(extent(i1), extent(i2))`. + +Here, `f(i1::Int, i2::Int)` may be any function that takes two integers as arguments. +It may optionally return an [`Action`](@ref LoopStateMachine.Action) to alter the control +flow of the `Action(:full_return, true)` to return `Action(:full_return, true)` from this +function and break out of the recursion. + +This is generic to anything that implements the SpatialTreeInterface, particularly the methods +[`isleaf`](@ref), [`getchild`](@ref), and [`child_indices_extents`](@ref). + +## Examples + +```julia +using NaturalEarth, +``` +""" +function dual_depth_first_search(f::F, predicate::P, node1::N1, node2::N2) where {F, P, N1, N2} + if isleaf(node1) && isleaf(node2) + # both nodes are leaves, so we can just iterate over the indices and extents + for (i1, extent1) in child_indices_extents(node1) + for (i2, extent2) in child_indices_extents(node2) + if predicate(extent1, extent2) + @controlflow f(i1, i2) + end + end + end + elseif isleaf(node1) # node2 is not a leaf, node1 is - recurse further into node2 + for child in getchild(node2) + if predicate(node_extent(node1), node_extent(child)) + @controlflow dual_depth_first_search(f, predicate, node1, child) + end + end + elseif isleaf(node2) # node1 is not a leaf, node2 is - recurse further into node1 + for child in getchild(node1) + if predicate(node_extent(child), node_extent(node2)) + @controlflow dual_depth_first_search(f, predicate, child, node2) + end + end + else # neither node is a leaf, recurse into both children + for child1 in getchild(node1) + for child2 in getchild(node2) + if predicate(node_extent(child1), node_extent(child2)) + @controlflow dual_depth_first_search(f, predicate, child1, child2) + end + end + end + end +end diff --git a/src/utils/SpatialTreeInterface/implementations.jl b/src/utils/SpatialTreeInterface/implementations.jl new file mode 100644 index 0000000000..a2df1e5dc1 --- /dev/null +++ b/src/utils/SpatialTreeInterface/implementations.jl @@ -0,0 +1,55 @@ +# # Interface implementations +# Below are some basic implementations of the interface, +# for STRTree and a "no-tree" implementation that is a flat list of extents. + +using SortTileRecursiveTree: STRtree, STRNode, STRLeafNode + +# ## SortTileRecursiveTree + +isspatialtree(::Type{<: STRtree}) = true +node_extent(tree::STRtree) = node_extent(tree.rootnode) +nchild(tree::STRtree) = nchild(tree.rootnode) +getchild(tree::STRtree) = getchild(tree.rootnode) +getchild(tree::STRtree, i) = getchild(tree.rootnode, i) +isleaf(tree::STRtree) = isleaf(tree.rootnode) +child_indices_extents(tree::STRtree) = child_indices_extents(tree.rootnode) + +isspatialtree(::Type{<: STRNode}) = true +node_extent(node::STRNode) = node.extent +nchild(node::STRNode) = length(node.children) +getchild(node::STRNode) = node.children +getchild(node::STRNode, i) = node.children[i] +isleaf(node::STRNode) = false # STRNodes are not leaves by definition + +isspatialtree(::Type{<: STRLeafNode}) = true +node_extent(node::STRLeafNode) = node.extent +isleaf(node::STRLeafNode) = true +child_indices_extents(node::STRLeafNode) = zip(node.indices, node.extents) + +# ## FlatNoTree +""" + FlatNoTree(iterable_of_geoms_or_extents) + +Represents a flat collection with no tree structure, i.e., a brute force search. +This is cost free, so particularly useful when you don't want to build a tree! +""" +struct FlatNoTree{T} + geometries::T +end + +isspatialtree(::Type{<: FlatNoTree}) = true +isleaf(tree::FlatNoTree) = true + +# NOTE: use pairs instead of enumerate here, so that we can support +# iterators or collections that define custom `pairs` methods. +# This includes things like filtered extent lists, for example, +# so we can perform extent thinning with no allocations. +function child_indices_extents(tree::FlatNoTree{T}) where T + # This test only applies at compile time and should be optimized away in any case. + # And we can use multiple dispatch to override anyway, but it should be cost free I think. + if applicable(Base.keys, T) + return ((i, GI.extent(obj)) for (i, obj) in pairs(tree.geometries)) + else + return ((i, GI.extent(obj)) for (i, obj) in enumerate(tree.geometries)) + end +end \ No newline at end of file diff --git a/src/utils/SpatialTreeInterface/interface.jl b/src/utils/SpatialTreeInterface/interface.jl new file mode 100644 index 0000000000..5e0bf08ab8 --- /dev/null +++ b/src/utils/SpatialTreeInterface/interface.jl @@ -0,0 +1,85 @@ +# # Interface +# Interface definition for spatial tree types. +# There is no abstract supertype here since it's impossible to enforce, +# but we do have a few methods that are common to all spatial tree types. + +""" + isspatialtree(tree)::Bool + +Return true if the object is a spatial tree, false otherwise. + +## Implementation notes + +For type stability, if your spatial tree type is `MyTree`, you should define +`isspatialtree(::Type{MyTree}) = true`, and `isspatialtree(::MyTree)` will forward +to that method automatically. +""" +isspatialtree(::T) where T = isspatialtree(T) +isspatialtree(::Type{<: Any}) = false + + +""" + getchild(node) + getchild(node, i) + +Accessor function to get the children of a node. + +If invoked as `getchild(node)`, return an iterator over all the children of a node. +This may be lazy, like a `Base.Generator`, or it may be materialized. + +If invoked as `getchild(node, i)`, return the `i`-th child of a node. +""" +function getchild end + +getchild(node) = AbstractTrees.children(node) + +""" + getchild(node, i) + +Return the `i`-th child of a node. +""" +getchild(node, i) = getchild(node)[i] + +""" + nchild(node) + +Return the number of children of a node. +""" +nchild(node) = length(getchild(node)) + +""" + isleaf(node) + +Return true if the node is a leaf node, i.e., there are no "children" below it. +[`getchild`](@ref) should still work on leaf nodes, though, returning an iterator over the extents stored in the node - and similarly for `getnodes.` +""" +isleaf(node) = error("isleaf is not implemented for node type $(typeof(node))") + +""" + child_indices_extents(node) + +Return an iterator over the indices and extents of the children of a node. + +Each value of the iterator should take the form `(i, extent)`. + +This can only be invoked on leaf nodes! +""" +function child_indices_extents(node) + return zip(1:nchild(node), getchild(node)) +end + +""" + node_extent(node) + +Return the extent like object of the node. +Falls back to `GI.extent` by default, which falls back +to `Extents.extent`. + +Generally, defining `Extents.extent(node)` is sufficient here, and you +won't need to define this + +The reason we don't use that directly is to give users of this interface +a way to define bounding boxes that are not extents, like spherical caps +and other such things. +""" +node_extent(node) = GI.extent(node) From bb7620723a60422e536370e166b8fd75766fa38a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 15:20:38 -0400 Subject: [PATCH 16/30] Fix tests --- test/utils/SpatialTreeInterface.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 50b5f25d23..94f09aaa02 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -2,7 +2,7 @@ using Test using GeometryOps.SpatialTreeInterface using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent -using GeometryOps.SpatialTreeInterface: query, do_query, do_dual_query +using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI @@ -81,13 +81,13 @@ function test_dual_query_functionality(TreeType) # Test dual query with a predicate that matches all all_pred = (x, y) -> true results = Tuple{Int, Int}[] - do_dual_query((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) + dual_depth_first_search((i, j) -> push!(results, (i, j)), all_pred, tree1, tree2) @test length(results) == 4 # 2 points in tree1 * 2 points in tree2 # Test dual query with a specific predicate intersects_pred = (x, y) -> Extents.intersects(x, y) results = Tuple{Int, Int}[] - do_dual_query((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) + dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end end From 8f4de31ca61c5dc3038596a64f16951e883e7f80 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 17 Apr 2025 15:21:03 -0400 Subject: [PATCH 17/30] no really --- src/utils/SpatialTreeInterface/SpatialTreeInterface.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl index 2cc15b6c6f..d306d75b45 100644 --- a/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl +++ b/src/utils/SpatialTreeInterface/SpatialTreeInterface.jl @@ -32,7 +32,7 @@ Return a sorted list of indices of the tree that satisfy the predicate. """ function query(tree, predicate) a = Int[] - do_query(Base.Fix1(push!, a), sanitize_predicate(predicate), tree) + depth_first_search(Base.Fix1(push!, a), sanitize_predicate(predicate), tree) return sort!(a) end From d5a49de75873f1c5f40fb149eb4ce384b308f208 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 20:01:37 -0400 Subject: [PATCH 18/30] add a test for multilevel dual tree query --- Project.toml | 3 +- test/utils/SpatialTreeInterface.jl | 89 ++++++++++++++++++++++++++---- 2 files changed, 80 insertions(+), 12 deletions(-) diff --git a/Project.toml b/Project.toml index d0eaaa7ff5..8aad7e7d4e 100644 --- a/Project.toml +++ b/Project.toml @@ -67,6 +67,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" +Polylabel = "49a44318-e865-4b63-9842-695152d634c1" Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689" @@ -76,4 +77,4 @@ TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "SafeTestsets", "Shapefile", "TGGeometry", "Test"] +test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "Polylabel", "SafeTestsets", "Shapefile", "TGGeometry", "Test"] diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index 94f09aaa02..dc7e4245d8 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -7,6 +7,8 @@ using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents using GeoInterface: GeoInterface as GI using SortTileRecursiveTree: STRtree +using NaturalEarth +using Polylabel # Generic test functions for spatial trees function test_basic_interface(TreeType) @@ -67,7 +69,7 @@ function test_query_functionality(TreeType) end function test_dual_query_functionality(TreeType) - @testset "Dual query functionality" begin + @testset "Dual query functionality - simple" begin # Create two trees with overlapping extents tree1 = TreeType([ Extents.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)), @@ -90,6 +92,65 @@ function test_dual_query_functionality(TreeType) dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end + + @testset "Dual query functionality - every country's polylabel against every country" begin + + # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced + # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) + # from Natural Earth, or by using GADM across many countries. + + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + all_adm0_a3 = all_countries.ADM0_A3 + all_geoms = all_countries.geometry + # US minor outlying islands - bug in Polylabel.jl + # A lot of small geoms have this issue, that there will be an error from the queue + # because the cell exists in the queue already. + # Not sure what to do about it, I don't want to check containment every time... + deleteat!(all_adm0_a3, 205) + deleteat!(all_geoms, 205) + + geom_tree = TreeType(all_geoms) + + polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] + polylabel_tree = TreeType(polylabels) + + found_countries = falses(length(polylabels)) + + dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + if i == j + found_countries[i] = true + end + end + + @test all(found_countries) + end +end + +function test_find_point_in_all_countries(TreeType) + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + tree = TreeType(all_countries.geometry) + + ber = (13.4050, 52.5200) # Berlin + nyc = (-74.0060, 40.7128) # New York City + sin = (103.8198, 1.3521) # Singapore + + @testset "locate points using query" begin + @testset let point = ber, name = "Berlin" + # Test Berlin (should be in Germany) + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "DEU", results) + end + @testset let point = nyc, name = "New York City" + # Test NYC (should be in USA) + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "USA", results) + end + @testset let point = sin, name = "Singapore" + # Test Singapore + results = query(tree, point) + @test any(i -> all_countries.ADM0_A3[i] == "SGP", results) + end + end end function test_geometry_support(TreeType) @@ -119,20 +180,26 @@ end # Test FlatNoTree implementation @testset "FlatNoTree" begin - test_basic_interface(FlatNoTree) - test_child_indices_extents(FlatNoTree) - test_query_functionality(FlatNoTree) - test_dual_query_functionality(FlatNoTree) - test_geometry_support(FlatNoTree) + @testset let TreeType = FlatNoTree + test_basic_interface(TreeType) + test_child_indices_extents(TreeType) + test_query_functionality(TreeType) + test_dual_query_functionality(TreeType) + test_geometry_support(TreeType) + test_find_point_in_all_countries(TreeType) + end end # Test STRtree implementation @testset "STRtree" begin - test_basic_interface(STRtree) - test_child_indices_extents(STRtree) - test_query_functionality(STRtree) - test_dual_query_functionality(STRtree) - test_geometry_support(STRtree) + @testset let TreeType = STRtree + test_basic_interface(TreeType) + test_child_indices_extents(TreeType) + test_query_functionality(TreeType) + test_dual_query_functionality(TreeType) + test_geometry_support(TreeType) + test_find_point_in_all_countries(TreeType) + end end From 06727f1bc23f6f916d3febcd7e689de9cc625357 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 20:32:31 -0400 Subject: [PATCH 19/30] Fix up the tests --- test/utils/SpatialTreeInterface.jl | 150 ++++++++++++++++------------- 1 file changed, 84 insertions(+), 66 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index dc7e4245d8..c2d38e37d1 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -92,37 +92,47 @@ function test_dual_query_functionality(TreeType) dual_depth_first_search((i, j) -> push!(results, (i, j)), intersects_pred, tree1, tree2) @test sort(results) == [(1,1), (2,1), (2,2)] end + @testset "Dual tree query with many boundingboxes" begin + xs = 1:100 + ys = 1:100 + extent_grid = [Extents.Extent(X=(x, x+1), Y=(y, y+1)) for x in xs, y in ys] |> vec + point_grid = [(x + 0.5, y + 0.5) for x in xs, y in ys] |> vec - @testset "Dual query functionality - every country's polylabel against every country" begin - - # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced - # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) - # from Natural Earth, or by using GADM across many countries. - - all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) - all_adm0_a3 = all_countries.ADM0_A3 - all_geoms = all_countries.geometry - # US minor outlying islands - bug in Polylabel.jl - # A lot of small geoms have this issue, that there will be an error from the queue - # because the cell exists in the queue already. - # Not sure what to do about it, I don't want to check containment every time... - deleteat!(all_adm0_a3, 205) - deleteat!(all_geoms, 205) - - geom_tree = TreeType(all_geoms) - - polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] - polylabel_tree = TreeType(polylabels) - - found_countries = falses(length(polylabels)) + extent_tree = TreeType(extent_grid) + point_tree = TreeType(point_grid) - dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + found_everything = falses(length(extent_grid)) + dual_depth_first_search(Extents.intersects, extent_tree, point_tree) do i, j if i == j - found_countries[i] = true + found_everything[i] = true end end + @test all(found_everything) + end +end - @test all(found_countries) +function test_geometry_support(TreeType) + @testset "Geometry support" begin + # Create a tree with 100 points + points = tuple.(1:100, 1:100) + tree = TreeType(points) + + # Test basic interface + @test isspatialtree(tree) + @test isspatialtree(typeof(tree)) + + # Test query functionality + all_pred = x -> true + results = query(tree, all_pred) + @test sort(results) == collect(1:100) + + none_pred = x -> false + results = query(tree, none_pred) + @test isempty(results) + + search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) + results = query(tree, Base.Fix1(Extents.intersects, search_extent)) + @test sort(results) == collect(45:55) end end @@ -153,54 +163,62 @@ function test_find_point_in_all_countries(TreeType) end end -function test_geometry_support(TreeType) - @testset "Geometry support" begin - # Create a tree with 100 points - points = tuple.(1:100, 1:100) - tree = TreeType(points) - - # Test basic interface - @test isspatialtree(tree) - @test isspatialtree(typeof(tree)) - - # Test query functionality - all_pred = x -> true - results = query(tree, all_pred) - @test sort(results) == collect(1:100) - - none_pred = x -> false - results = query(tree, none_pred) - @test isempty(results) - - search_extent = Extents.Extent(X=(45.0, 55.0), Y=(45.0, 55.0)) - results = query(tree, Base.Fix1(Extents.intersects, search_extent)) - @test sort(results) == collect(45:55) - end -end - # Test FlatNoTree implementation @testset "FlatNoTree" begin - @testset let TreeType = FlatNoTree - test_basic_interface(TreeType) - test_child_indices_extents(TreeType) - test_query_functionality(TreeType) - test_dual_query_functionality(TreeType) - test_geometry_support(TreeType) - test_find_point_in_all_countries(TreeType) - end + test_basic_interface(FlatNoTree) + test_child_indices_extents(FlatNoTree) + test_query_functionality(FlatNoTree) + test_dual_query_functionality(FlatNoTree) + test_geometry_support(FlatNoTree) + test_find_point_in_all_countries(FlatNoTree) end # Test STRtree implementation @testset "STRtree" begin - @testset let TreeType = STRtree - test_basic_interface(TreeType) - test_child_indices_extents(TreeType) - test_query_functionality(TreeType) - test_dual_query_functionality(TreeType) - test_geometry_support(TreeType) - test_find_point_in_all_countries(TreeType) - end + test_basic_interface(STRtree) + test_child_indices_extents(STRtree) + test_query_functionality(STRtree) + test_dual_query_functionality(STRtree) + test_geometry_support(STRtree) + test_find_point_in_all_countries(STRtree) end +# This testset is not used because Polylabel.jl has some issues. + +#= + + + @testset "Dual query functionality - every country's polylabel against every country" begin + + # Note that this is a perfectly balanced tree query - we don't yet have a test for unbalanced + # trees (but could easily add one, e.g. by getting polylabels of admin-1 or admin-2 regions) + # from Natural Earth, or by using GADM across many countries. + + all_countries = NaturalEarth.naturalearth("admin_0_countries", 10) + all_adm0_a3 = all_countries.ADM0_A3 + all_geoms = all_countries.geometry + # US minor outlying islands - bug in Polylabel.jl + # A lot of small geoms have this issue, that there will be an error from the queue + # because the cell exists in the queue already. + # Not sure what to do about it, I don't want to check containment every time... + deleteat!(all_adm0_a3, 205) + deleteat!(all_geoms, 205) + + geom_tree = TreeType(all_geoms) + + polylabels = [Polylabel.polylabel(geom; rtol = 0.019) for geom in all_geoms] + polylabel_tree = TreeType(polylabels) + + found_countries = falses(length(polylabels)) + + dual_depth_first_search(Extents.intersects, geom_tree, polylabel_tree) do i, j + if i == j + found_countries[i] = true + end + end + + @test all(found_countries) + end +=# \ No newline at end of file From 40246680097eddc9afc61caa794b84c055ad5ede Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:35:28 -0400 Subject: [PATCH 20/30] add tests for imbalanced dual tree --- test/utils/SpatialTreeInterface.jl | 39 ++++++++++++++++++++++++++++-- 1 file changed, 37 insertions(+), 2 deletions(-) diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl index c2d38e37d1..eb61f553b7 100644 --- a/test/utils/SpatialTreeInterface.jl +++ b/test/utils/SpatialTreeInterface.jl @@ -1,11 +1,10 @@ using Test - +import GeometryOps as GO, GeoInterface as GI using GeometryOps.SpatialTreeInterface using GeometryOps.SpatialTreeInterface: isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent using GeometryOps.SpatialTreeInterface: query, depth_first_search, dual_depth_first_search using GeometryOps.SpatialTreeInterface: FlatNoTree using Extents -using GeoInterface: GeoInterface as GI using SortTileRecursiveTree: STRtree using NaturalEarth using Polylabel @@ -109,6 +108,42 @@ function test_dual_query_functionality(TreeType) end @test all(found_everything) end + + @testset "Imbalanced dual query - tree 1 deeper than tree 2" begin + xs = 0:0.01:10 + ys = 0:0.01:10 + extent_grid = [Extents.Extent(X=(x, x+0.1), Y=(y, y+0.1)) for x in xs, y in ys] |> vec + point_grid = [(x + 0.5, y + 0.5) for x in 0:9, y in 0:9] |> vec + + extent_tree = TreeType(extent_grid) + point_tree = TreeType(point_grid) + + found_everything = falses(length(point_grid)) + dual_depth_first_search(Extents.intersects, extent_tree, point_tree) do i, j + if Extents.intersects(extent_grid[i], GI.extent(point_grid[j])) + found_everything[j] = true + end + end + @test all(found_everything) + end + + @testset "Imbalanced dual query - tree 2 deeper than tree 1" begin + xs = 0:0.01:10 + ys = 0:0.01:10 + extent_grid = [Extents.Extent(X=(x, x+0.1), Y=(y, y+0.1)) for x in xs, y in ys] |> vec + point_grid = [(x + 0.5, y + 0.5) for x in 0:9, y in 0:9] |> vec + + extent_tree = TreeType(extent_grid) + point_tree = TreeType(point_grid) + + found_everything = falses(length(point_grid)) + dual_depth_first_search(Extents.intersects, point_tree, extent_tree) do i, j + if Extents.intersects(GI.extent(point_grid[i]), extent_grid[j]) + found_everything[i] = true + end + end + @test all(found_everything) + end end function test_geometry_support(TreeType) From 3529ac24880640a0b6fc03eee6c3eb2ba50daebc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:35:58 -0400 Subject: [PATCH 21/30] add docs + animation to dfs --- .../depth_first_search.jl | 102 ++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/src/utils/SpatialTreeInterface/depth_first_search.jl b/src/utils/SpatialTreeInterface/depth_first_search.jl index dc1be59b3a..fafae8a42a 100644 --- a/src/utils/SpatialTreeInterface/depth_first_search.jl +++ b/src/utils/SpatialTreeInterface/depth_first_search.jl @@ -1,3 +1,80 @@ +#= +# Depth-first search + +```@meta +CollapsedDocStrings = true +``` + +```@docs; canonical=false +depth_first_search +``` + +This file houses a function that performs a depth-first search over a spatial tree, filtering +by `predicate` on the extents of the nodes, and calling `f` on every leaf node that matches +`predicate`. + +Here, `predicate` is a one-argument function that returns a boolean, and `f` is a function +that takes a single argument, the index of a geometry as stored in the leaf node. + +The `depth_first_search` function is generic to anything that implements the SpatialTreeInterface. + +## Example + +Here's an animated example of what `depth_first_search` is doing: + +![depth_first_search](./depth_first_search.mp4) + +```@example dfs +using Makie, GeoInterfaceMakie +using Extents, GeometryOps.SpatialTreeInterface, SortTileRecursiveTree +import GeoInterface as GI + +# Construct a tree of extents (in this case, just a grid) +xs, ys = 1.0:0.1:10.0, 1.0:0.1:10.0 +rects_vec = vec([GI.Polygon([[(x, y), (x+0.1, y), (x+0.1, y+0.1), (x, y+0.1), (x, y)]]) for x in xs, y in ys]) +ext2rect(ext) = Rect2f((ext.X[1], ext.Y[1]), (ext.X[2] - ext.X[1], ext.Y[2] - ext.Y[1])) + +# First, create a figure that holds some plots: +f, a, all_rects_plot = lines(rects_vec; linewidth = 0.5, axis = (; aspect = DataAspect())); +Legend(f[1, 2], [Makie.PolyElement(; color = :transparent, strokecolor = :blue, strokewidth = 0.5), Makie.PolyElement(; color = :green), Makie.PolyElement(; color = :red)], ["Grid", "Branch node", "Leaf node"]) + +target_extent = Extents.Extent(X = (5, 7), Y = (5, 7)) +target_rect_plot = poly!(a, ext2rect(target_extent); color = (:blue, 0.5)) + +current_rect = Observable(ext2rect(GI.extent(rects_vec[1]))) +current_color = Observable(:red) + +current_rect_plot = poly!(a, current_rect; color = current_color, alpha = 0.3) + + +record(f, "./depth_first_search.mp4"; framerate = 1) do io + + function leaf_f(i) + current_rect[] = ext2rect(GI.extent(rects_vec[i])) + current_color[] = :red + recordframe!(io) + return nothing + end + + function pred_f(ext, target) + current_rect[] = ext2rect(ext) + current_color[] = :green + recordframe!(io) + result = Extents.intersects(ext, target) + if result + current_color[] = :forestgreen + recordframe!(io) + end + return result + end + + SpatialTreeInterface.depth_first_search(leaf_f, Base.Fix2(pred_f, target_extent), extents_tree) +end + +```` + +=# + """ depth_first_search(f, predicate, tree) @@ -6,6 +83,31 @@ Call `f(i)` for each index `i` in the tree that satisfies `predicate(extent(i))` This is generic to anything that implements the SpatialTreeInterface, particularly the methods [`isleaf`](@ref), [`getchild`](@ref), and [`child_extents`](@ref). + +## Example + +```jldoctest +using Extents, GeometryOps.SpatialTreeInterface, SortTileRecursiveTree +# Construct a tree of extents (in this case, just a grid) +xs, ys = 1:10, 1:10 +extents_vec = vec([Extents.Extent(X = (x, x+1), Y = (y, y+1)) for x in xs, y in ys]) + +# Construct a tree of extents - this is an STRtree, +# but works on any SpatialTreeInterface-compatible tree. +extents_tree = STRtree(extents_vec) + +# Count the number of extents that intersect the extent (1,1) x (2,2) +target_extent = Extents.Extent(X = (1,2), Y = (1,2)) +count = 0 +SpatialTreeInterface.depth_first_search( + (i::Int -> global count += 1), # `f` + Base.Fix1(Extents.intersects, target_extent), + extents_tree +) +count +# output +4 +``` """ function depth_first_search(f::F, predicate::P, node::N) where {F, P, N} if isleaf(node) From b68fbba4f01b1f91194c38e5046d6f66a1a604e6 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:45:14 -0400 Subject: [PATCH 22/30] test full_return properly --- test/utils/LoopStateMachine.jl | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index 30d4979c81..eb5f53bd01 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -48,6 +48,31 @@ end @test f(1) == Action(:full_return, 1) @test f(2) == Action(:full_return, 2) @test f(3) == Action(:full_return, 3) + + function descend(i) + if i == 1 + k = 0 + for j in 1:4 + k = @controlflow h(i, j) + end + return k + else + descend(i-1) + end + end + + function h(i, j) + for _i in i:-1:1 + for _j in 1:j + @controlflow l(j, i) + end + end + end + + l(j, i) = i == 1 ? Action(:full_return, j) : i + + @test descend(4) == Action(:full_return, 1) + end @testset "Return value" begin From 86931a6b91947c8ad862d60186db26b8c78aaf8e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:50:24 -0400 Subject: [PATCH 23/30] implement and use node_extent for FlatNoTree --- src/utils/SpatialTreeInterface/implementations.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/utils/SpatialTreeInterface/implementations.jl b/src/utils/SpatialTreeInterface/implementations.jl index a2df1e5dc1..0762d503bf 100644 --- a/src/utils/SpatialTreeInterface/implementations.jl +++ b/src/utils/SpatialTreeInterface/implementations.jl @@ -39,6 +39,7 @@ end isspatialtree(::Type{<: FlatNoTree}) = true isleaf(tree::FlatNoTree) = true +node_extent(tree::FlatNoTree) = mapreduce(GI.extent, Extents.union, tree.geometries) # NOTE: use pairs instead of enumerate here, so that we can support # iterators or collections that define custom `pairs` methods. @@ -48,8 +49,8 @@ function child_indices_extents(tree::FlatNoTree{T}) where T # This test only applies at compile time and should be optimized away in any case. # And we can use multiple dispatch to override anyway, but it should be cost free I think. if applicable(Base.keys, T) - return ((i, GI.extent(obj)) for (i, obj) in pairs(tree.geometries)) + return ((i, node_extent(obj)) for (i, obj) in pairs(tree.geometries)) else - return ((i, GI.extent(obj)) for (i, obj) in enumerate(tree.geometries)) + return ((i, node_extent(obj)) for (i, obj) in enumerate(tree.geometries)) end end \ No newline at end of file From 8456cb7598bcf94215261e92ea0b55beee11c935 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:50:44 -0400 Subject: [PATCH 24/30] more robust implementation of `getchild` that may use `pairs` --- src/utils/SpatialTreeInterface/interface.jl | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/utils/SpatialTreeInterface/interface.jl b/src/utils/SpatialTreeInterface/interface.jl index 5e0bf08ab8..b6fcc20b87 100644 --- a/src/utils/SpatialTreeInterface/interface.jl +++ b/src/utils/SpatialTreeInterface/interface.jl @@ -65,7 +65,12 @@ Each value of the iterator should take the form `(i, extent)`. This can only be invoked on leaf nodes! """ function child_indices_extents(node) - return zip(1:nchild(node), getchild(node)) + children = getchild(node) + if applicable(Base.keys, typeof(children)) + return ((i, node_extent(obj)) for (i, obj) in pairs(children)) + else + return ((i, node_extent(obj)) for (i, obj) in enumerate(children)) + end end """ From a6b0c5247b67b0f5f34816abd43bcd505913b814 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:50:50 -0400 Subject: [PATCH 25/30] rename LSM test --- test/utils/LoopStateMachine.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl index eb5f53bd01..a734dbe05c 100644 --- a/test/utils/LoopStateMachine.jl +++ b/test/utils/LoopStateMachine.jl @@ -75,7 +75,7 @@ end end -@testset "Return value" begin +@testset "Return value without Action" begin results = Int[] for i in 1:3 val = @controlflow begin From ab1f5c8923716e96dea8399a563ab64a81d86671 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 22:52:21 -0400 Subject: [PATCH 26/30] add packages to docs Project --- docs/Project.toml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/Project.toml b/docs/Project.toml index f8e1c22448..3dee4fc743 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -16,6 +16,7 @@ DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365" DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78" Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" +Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37" GADM = "a8dd9ffe-31dc-4cf5-a379-ea69100a8233" GeoDataFrames = "62cb38b5-d8d2-4862-a48e-6a340996859f" @@ -42,6 +43,7 @@ Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d" TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b" [sources] From 73b67904778a78198f49f4a4b803ebc04f633bdc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 23:10:47 -0400 Subject: [PATCH 27/30] actually throw the error correctly --- src/utils/utils.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/utils/utils.jl b/src/utils/utils.jl index fd882768c1..a539f58592 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -157,10 +157,10 @@ function eachedge(trait::GI.AbstractGeometryTrait, geom, ::Type{T}) where T return Iterators.flatten((eachedge(r, T) for r in flatten(GI.AbstractCurveTrait, geom))) end function eachedge(trait::GI.PointTrait, geom, ::Type{T}) where T - return ArgumentError("Can't get edges from points, $geom was a PointTrait.") + throw(ArgumentError("Can't get edges from points, $geom was a PointTrait.")) end function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T - return ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.") + throw(ArgumentError("Can't get edges from MultiPoint, $geom was a MultiPointTrait.")) end """ From 8fad4934d2b362ca7ccd4da9ecc175e6ec2cf8a7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 23:10:58 -0400 Subject: [PATCH 28/30] add options + docstrings to utils --- src/utils/utils.jl | 37 +++++++++++++++++++++++++++++++------ 1 file changed, 31 insertions(+), 6 deletions(-) diff --git a/src/utils/utils.jl b/src/utils/utils.jl index a539f58592..481ac9ccce 100644 --- a/src/utils/utils.jl +++ b/src/utils/utils.jl @@ -168,9 +168,18 @@ end Convert a geometry into a vector of `GI.Line` objects with attached extents. """ -to_edgelist(geom, ::Type{T}) where T = +to_edgelist(geom, ::Type{T} = Float64) where T = [_lineedge(ps, T) for ps in eachedge(geom, T)] -function to_edgelist(ext::E, geom, ::Type{T}) where {E<:Extents.Extent,T} + +""" + to_edgelist(ext::E, geom, [::Type{T}])::(::Vector{GI.Line}, ::Vector{Int}) + +Filter the edges of `geom` for those that intersect +`ext`, and return: +- a vector of `GI.Line` objects with attached extents, +- a vector of indices into the original geometry. +""" +function to_edgelist(ext::E, geom, ::Type{T} = Float64) where {E<:Extents.Extent,T} edges_in = eachedge(geom, T) l1 = _lineedge(first(edges_in), T) edges_out = typeof(l1)[] @@ -186,22 +195,38 @@ function to_edgelist(ext::E, geom, ::Type{T}) where {E<:Extents.Extent,T} end function _lineedge(ps::Tuple, ::Type{T}) where T - l = GI.Line(SVector{2,NTuple{2, T}}(ps)) # TODO: make this flexible in dimension + l = GI.Line(StaticArrays.SVector{2,NTuple{2, T}}(ps)) # TODO: make this flexible in dimension e = GI.extent(l) return GI.Line(l.geom; extent=e) end -function lazy_edgelist(geom, ::Type{T}) where T +""" + lazy_edgelist(geom, [::Type{T}]) + +Return an iterator over `GI.Line` objects with attached extents. +""" +function lazy_edgelist(geom, ::Type{T} = Float64) where T (_lineedge(ps, T) for ps in eachedge(geom, T)) end -function edge_extents(geom) +""" + edge_extents(geom, [::Type{T}]) + +Return a vector of the extents of the edges (line segments) of `geom`. +""" +function edge_extents(geom, ::Type{T} = Float64) where T return [begin Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) end - for edge in eachedge(geom, Float64)] + for edge in eachedge(geom, T)] end +""" + lazy_edge_extents(geom) + +Return an iterator over the extents of the edges (line segments) of `geom`. +This is lazy but nonallocating. +""" function lazy_edge_extents(geom) return (begin Extents.Extent(X=extrema(GI.x, edge), Y=extrema(GI.y, edge)) From 1b86b5ec7699a80771bc977d526a4258a911dd87 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 23:11:06 -0400 Subject: [PATCH 29/30] add tests for new utils methods --- test/runtests.jl | 1 + test/utils/utils.jl | 192 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 193 insertions(+) create mode 100644 test/utils/utils.jl diff --git a/test/runtests.jl b/test/runtests.jl index 829b2e545c..7249981f24 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -13,6 +13,7 @@ end @safetestset "Primitives" begin include("primitives.jl") end # Utils +@safetestset "Utils" begin include("utils/utils.jl") end @safetestset "LoopStateMachine" begin include("utils/LoopStateMachine.jl") end @safetestset "SpatialTreeInterface" begin include("utils/SpatialTreeInterface.jl") end diff --git a/test/utils/utils.jl b/test/utils/utils.jl new file mode 100644 index 0000000000..79ddf4e301 --- /dev/null +++ b/test/utils/utils.jl @@ -0,0 +1,192 @@ +using Test + +import GeometryOps as GO, GeoInterface as GI +import Extents + +point = GI.Point(1.0, 1.0) +linestring = GI.LineString([(1.0, 1.0), (2.0, 2.0)]) +multilinestring = GI.MultiLineString([[(1.0, 1.0), (2.0, 2.0)], [(3.0, 3.0), (4.0, 4.0)]]) + +polygon = GI.Polygon([[(1.0, 1.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)]]) +multipolygon = GI.MultiPolygon([[[(1.0, 1.0), (2.0, 2.0), (2.0, 1.0), (1.0, 1.0)]], [[(3.0, 3.0), (4.0, 4.0), (4.0, 3.0), (3.0, 3.0)]]]) + +# Test eachedge +@testset "eachedge" begin + # Test LineString + edges = collect(GO.eachedge(linestring)) + @test length(edges) == 1 + @test edges[1] == ((1.0, 1.0), (2.0, 2.0)) + + # Test MultiLineString + edges = collect(GO.eachedge(multilinestring)) + @test length(edges) == 2 + @test edges[1] == ((1.0, 1.0), (2.0, 2.0)) + @test edges[2] == ((3.0, 3.0), (4.0, 4.0)) + + # Test Polygon + edges = collect(GO.eachedge(polygon)) + @test length(edges) == 3 + @test edges[1] == ((1.0, 1.0), (2.0, 2.0)) + @test edges[2] == ((2.0, 2.0), (2.0, 1.0)) + @test edges[3] == ((2.0, 1.0), (1.0, 1.0)) + + # Test MultiPolygon + edges = collect(GO.eachedge(multipolygon)) + @test length(edges) == 6 # 3 edges per polygon + @test edges[1] == ((1.0, 1.0), (2.0, 2.0)) + @test edges[2] == ((2.0, 2.0), (2.0, 1.0)) + @test edges[3] == ((2.0, 1.0), (1.0, 1.0)) + @test edges[4] == ((3.0, 3.0), (4.0, 4.0)) + @test edges[5] == ((4.0, 4.0), (4.0, 3.0)) + @test edges[6] == ((4.0, 3.0), (3.0, 3.0)) + + # Test error cases + @test_throws ArgumentError GO.eachedge(point) + @test_throws ArgumentError GO.eachedge(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)])) +end + +# Test to_edgelist +@testset "to_edgelist" begin + # Test LineString + edges = GO.to_edgelist(linestring) + @test length(edges) == 1 + @test GI.getpoint(edges[1], 1) == (1.0, 1.0) + @test GI.getpoint(edges[1], 2) == (2.0, 2.0) + + # Test MultiLineString + edges = GO.to_edgelist(multilinestring) + @test length(edges) == 2 + @test GI.getpoint(edges[1], 1) == (1.0, 1.0) + @test GI.getpoint(edges[1], 2) == (2.0, 2.0) + @test GI.getpoint(edges[2], 1) == (3.0, 3.0) + @test GI.getpoint(edges[2], 2) == (4.0, 4.0) + + # Test Polygon + edges = GO.to_edgelist(polygon) + @test length(edges) == 3 + @test GI.getpoint(edges[1], 1) == (1.0, 1.0) + @test GI.getpoint(edges[1], 2) == (2.0, 2.0) + @test GI.getpoint(edges[2], 1) == (2.0, 2.0) + @test GI.getpoint(edges[2], 2) == (2.0, 1.0) + @test GI.getpoint(edges[3], 1) == (2.0, 1.0) + @test GI.getpoint(edges[3], 2) == (1.0, 1.0) + + # Test MultiPolygon + edges = GO.to_edgelist(multipolygon) + @test length(edges) == 6 + + # Test with extent filtering + extent = Extents.Extent(X=(1.5, 2.5), Y=(1.5, 2.5)) + edges, indices = GO.to_edgelist(extent, linestring, Float64) + @test length(edges) == 1 + @test indices == [1] + + # Test with extent that doesn't intersect + extent = Extents.Extent(X=(3.0, 4.0), Y=(3.0, 4.0)) + edges, indices = GO.to_edgelist(extent, linestring, Float64) + @test isempty(edges) + @test isempty(indices) + + # Test error cases + @test_throws ArgumentError GO.to_edgelist(point) + @test_throws ArgumentError GO.to_edgelist(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)])) +end + +# Test lazy_edgelist +@testset "lazy_edgelist" begin + # Test LineString + edges = collect(GO.lazy_edgelist(linestring)) + @test length(edges) == 1 + @test GI.getpoint(first(edges), 1) == (1.0, 1.0) + @test GI.getpoint(first(edges), 2) == (2.0, 2.0) + + # Test MultiLineString + edges = collect(GO.lazy_edgelist(multilinestring)) + @test length(edges) == 2 + @test edges == GO.to_edgelist(multilinestring) + + # Test Polygon + edges = collect(GO.lazy_edgelist(polygon)) + @test length(edges) == 3 + @test edges == GO.to_edgelist(polygon) + + # Test MultiPolygon + edges = collect(GO.lazy_edgelist(multipolygon)) + @test length(edges) == 6 + @test edges == GO.to_edgelist(multipolygon) + + # Test error cases + @test_throws ArgumentError collect(GO.lazy_edgelist(point)) + @test_throws ArgumentError collect(GO.lazy_edgelist(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)]))) +end + +# Test edge_extents +@testset "edge_extents" begin + # Test LineString + extents = GO.edge_extents(linestring) + @test length(extents) == 1 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + + # Test MultiLineString + extents = GO.edge_extents(multilinestring) + @test length(extents) == 2 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + @test extents[2].X == (3.0, 4.0) + @test extents[2].Y == (3.0, 4.0) + + # Test Polygon + extents = GO.edge_extents(polygon) + @test length(extents) == 3 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + @test extents[2].X == (2.0, 2.0) + @test extents[2].Y == (1.0, 2.0) + @test extents[3].X == (1.0, 2.0) + @test extents[3].Y == (1.0, 1.0) + + # Test MultiPolygon + extents = GO.edge_extents(multipolygon) + @test length(extents) == 6 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + @test extents[4].X == (3.0, 4.0) + @test extents[4].Y == (3.0, 4.0) + + # Test error cases + @test_throws ArgumentError GO.edge_extents(point) + @test_throws ArgumentError GO.edge_extents(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)])) +end + +# Test lazy_edge_extents +@testset "lazy_edge_extents" begin + # Test LineString + extents = collect(GO.lazy_edge_extents(linestring)) + @test length(extents) == 1 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + + # Test MultiLineString + extents = collect(GO.lazy_edge_extents(multilinestring)) + @test length(extents) == 2 + @test extents[1].X == (1.0, 2.0) + @test extents[1].Y == (1.0, 2.0) + @test extents[2].X == (3.0, 4.0) + @test extents[2].Y == (3.0, 4.0) + + # Test Polygon + extents = collect(GO.lazy_edge_extents(polygon)) + @test length(extents) == 3 + @test extents == GO.edge_extents(polygon) + + # Test MultiPolygon + extents = collect(GO.lazy_edge_extents(multipolygon)) + @test length(extents) == 6 + @test extents == GO.edge_extents(multipolygon) + + # Test error cases + @test_throws ArgumentError collect(GO.lazy_edge_extents(point)) + @test_throws ArgumentError collect(GO.lazy_edge_extents(GI.MultiPoint([(1.0, 1.0), (2.0, 2.0)]))) +end + From b9fbec8472fcc5d0cca23179d214abaafe189ce9 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sat, 19 Apr 2025 23:32:56 -0400 Subject: [PATCH 30/30] don't run the gif I guess --- src/utils/SpatialTreeInterface/depth_first_search.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/utils/SpatialTreeInterface/depth_first_search.jl b/src/utils/SpatialTreeInterface/depth_first_search.jl index fafae8a42a..35705f1f1e 100644 --- a/src/utils/SpatialTreeInterface/depth_first_search.jl +++ b/src/utils/SpatialTreeInterface/depth_first_search.jl @@ -22,9 +22,8 @@ The `depth_first_search` function is generic to anything that implements the Spa Here's an animated example of what `depth_first_search` is doing: -![depth_first_search](./depth_first_search.mp4) -```@example dfs +```julia using Makie, GeoInterfaceMakie using Extents, GeometryOps.SpatialTreeInterface, SortTileRecursiveTree import GeoInterface as GI