diff --git a/Project.toml b/Project.toml index 8ee10bf44..8aad7e7d4 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Anshul Singhvi ", "Rafael Schouten 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) + 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 000000000..2c63d0aa1 --- /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 000000000..0762d503b --- /dev/null +++ b/src/utils/SpatialTreeInterface/implementations.jl @@ -0,0 +1,56 @@ +# # 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 +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. +# 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, node_extent(obj)) for (i, obj) in pairs(tree.geometries)) + else + return ((i, node_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 000000000..b6fcc20b8 --- /dev/null +++ b/src/utils/SpatialTreeInterface/interface.jl @@ -0,0 +1,90 @@ +# # 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) + 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 + +""" + 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) diff --git a/src/utils.jl b/src/utils/utils.jl similarity index 55% rename from src/utils.jl rename to src/utils/utils.jl index 540a7211e..481ac9ccc 100644 --- a/src/utils.jl +++ b/src/utils/utils.jl @@ -126,3 +126,110 @@ 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 + throw(ArgumentError("Can't get edges from points, $geom was a PointTrait.")) +end +function eachedge(trait::GI.MultiPointTrait, geom, ::Type{T}) where T + throw(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} = Float64) where T = + [_lineedge(ps, T) for ps in eachedge(geom, 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)[] + 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(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 + +""" + 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 + +""" + 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, 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)) + end + for edge in eachedge(geom, Float64)) +end diff --git a/test/runtests.jl b/test/runtests.jl index 3d5465c96..7249981f2 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,6 +11,12 @@ end @safetestset "Types" begin include("types.jl") 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 + # Methods @safetestset "Angles" begin include("methods/angles.jl") end @safetestset "Area" begin include("methods/area.jl") end diff --git a/test/utils/LoopStateMachine.jl b/test/utils/LoopStateMachine.jl new file mode 100644 index 000000000..a734dbe05 --- /dev/null +++ b/test/utils/LoopStateMachine.jl @@ -0,0 +1,115 @@ +using Test +import GeometryOps.LoopStateMachine +using GeometryOps.LoopStateMachine: @controlflow, Action + +@testset "Continue action" begin + count = 0 + f(i) = begin + count += 1 + if i == 3 + return Action(: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 Action(: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 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) + + 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 without Action" begin + results = Int[] + for i in 1:3 + val = @controlflow begin + i * 2 + end + push!(results, val) + 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 + +@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 diff --git a/test/utils/SpatialTreeInterface.jl b/test/utils/SpatialTreeInterface.jl new file mode 100644 index 000000000..eb61f553b --- /dev/null +++ b/test/utils/SpatialTreeInterface.jl @@ -0,0 +1,259 @@ +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 SortTileRecursiveTree: STRtree +using NaturalEarth +using Polylabel + +# 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 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 + @test indices_extents[1] == (1, extents[1]) + @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) + @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 +end + +function test_dual_query_functionality(TreeType) + @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)), + Extents.Extent(X=(1.0, 2.0), Y=(1.0, 2.0)) + ]) + 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)) + ]) + + # Test dual query with a predicate that matches all + all_pred = (x, y) -> true + results = Tuple{Int, Int}[] + 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}[] + 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 + + extent_tree = TreeType(extent_grid) + point_tree = TreeType(point_grid) + + found_everything = falses(length(extent_grid)) + dual_depth_first_search(Extents.intersects, extent_tree, point_tree) do i, j + if i == j + found_everything[i] = true + end + 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) + @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 + +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 + +# 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) + test_find_point_in_all_countries(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) + 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 diff --git a/test/utils/utils.jl b/test/utils/utils.jl new file mode 100644 index 000000000..79ddf4e30 --- /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 +