Spatial tree interface

Project.toml

@@ -4,6 +4,7 @@ authors = ["Anshul Singhvi <anshulsinghvi@gmail.com>", "Rafael Schouten <rafaels
 version = "0.1.18"
 
 [deps]
+AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 AdaptivePredicates = "35492f91-a3bd-45ad-95db-fcad7dcfedb7"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
@@ -32,6 +33,7 @@ GeometryOpsProjExt = "Proj"
 GeometryOpsTGGeometryExt = "TGGeometry"
 
 [compat]
+AbstractTrees = "0.4"
 AdaptivePredicates = "1.2"
 CoordinateTransformations = "0.5, 0.6"
 DataAPI = "1"
@@ -65,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"
@@ -74,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"]

docs/.gitignore

@@ -2,4 +2,6 @@ node_modules/
 build/
 package-lock.json
 src/source/
-Manifest.toml
+Manifest.toml
+
+src/call_notes.md

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]

src/GeometryOps.jl

@@ -36,9 +36,15 @@ 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")
+include("utils/SpatialTreeInterface/SpatialTreeInterface.jl")
+
+using .LoopStateMachine, .SpatialTreeInterface
+
+
 include("methods/angles.jl")
 include("methods/area.jl")
 include("methods/barycentric.jl")

src/utils/LoopStateMachine/LoopStateMachine.jl

@@ -0,0 +1,123 @@
+"""
+    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{T}(:unnamed, x)
+Action(x::Symbol) = Action(x, nothing)
+
+function Base.show(io::IO, action::Action{T}) where T
+    print(io, "Action")
+    print(io, "(:$(action.name)")
+    if isnothing(action.x)
+        print(io, ")")
+    else
+        print(io, ", ",action.x, ")")
+    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(...)
+
+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(UnrecognizedActionException($varname.name))
+            end
+        else
+            $varname
+        end
+    end
+end
+# COV_EXCL_STOP
+
+# You can define more actions as you desire.
+
+end

src/utils/SpatialTreeInterface/README.md

@@ -0,0 +1,47 @@
+# SpatialTreeInterface.jl
+
+A simple interface for spatial tree types.
+
+## What is a spatial tree?
+
+- 2 dimensional extents
+- Parent nodes encompass all leaf nodes
+- Leaf nodes contain references to the geometries they represent as indices (or so we assume here)
+
+## Why is this useful?
+
+- It allows us to write algorithms that can work with any spatial tree type, without having to know the details of the tree type.
+    - for example, dual tree traversal / queries
+- It allows us to flexibly and easily swap out and use different tree types, depending on the problem at hand.
+
+This is also a zero cost interface if implemented correctly!  Verified implementations exist for "flat" trees like the "Natural Index" from `tg`, and "hierarchical" trees like the `STRtree` from `SortTileRecursiveTree.jl`.
+
+## Interface
+
+- `isspatialtree(tree)::Bool`
+- `isleaf(node)::Bool` - is the node a leaf node?  In this context, a leaf node is a node that does not have other nodes as its children, but stores a list of indices and extents (even if implicit).
+- `getchild(node)` - get the children of a node.  This may be materialized if necessary or available, but can also be lazy (like a generator).
+- `getchild(node, i)` - get the `i`-th child of a node.
+- `nchild(node)::Int` - the number of children of a node.
+- `child_indices_extents(node)` - an iterator over the indices and extents of the children of a **leaf** node.
+
+These are the only methods that are required to be implemented.  
+
+Optionally, one may define:
+- `node_extent(node)` - get the extent of a node.  This falls back to `GI.extent` but can potentially be overridden if you want to return a different but extent-like object.
+
+They enable the generic query functions described below:
+
+## Query functions
+
+- `do_query(f, predicate, node)` - call `f(i)` for each index `i` in `node` that satisfies `predicate(extent(i))`.
+- `do_dual_query(f, predicate, tree1, tree2)` - call `f(i1, i2)` for each index `i1` in `tree1` and `i2` in `tree2` that satisfies `predicate(extent(i1), extent(i2))`.
+
+These are both completely non-allocating, and will only call `f` for indices that satisfy the predicate.
+You can of course build a standard query interface on top of `do_query` if you want - that's simply:
+```julia
+a = Int[]
+do_query(Base.Fix1(push!, a), predicate, node)
+```
+where `predicate` might be `Base.Fix1(Extents.intersects, extent_to_query)`.
+

src/utils/SpatialTreeInterface/SpatialTreeInterface.jl

@@ -0,0 +1,59 @@
+module SpatialTreeInterface
+
+import ..LoopStateMachine: @controlflow
+
+import Extents
+import GeoInterface as GI
+import AbstractTrees
+
+# public isspatialtree, isleaf, getchild, nchild, child_indices_extents, node_extent
+export query, do_query
+export FlatNoTree
+
+# The spatial tree interface and its implementations are defined here.
+include("interface.jl")
+include("implementations.jl")
+
+# 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")
+
+# 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")
+
+
+"""
+    query(tree, predicate)
+
+Return a sorted list of indices of the tree that satisfy the predicate.
+"""
+function query(tree, predicate)
+    a = Int[]
+    depth_first_search(Base.Fix1(push!, a), sanitize_predicate(predicate), tree)
+    return sort!(a)
+end
+
+
+"""
+    sanitize_predicate(pred)
+
+Convert a predicate to a function that returns a Boolean.
+
+If `pred` is an Extent, convert it to a function that returns a Boolean by intersecting with the extent.
+If `pred` is a geometry, convert it to an extent first, then wrap in Extents.intersects.
+
+Otherwise, return the predicate unchanged.
+
+
+Users and developers may overload this function to provide custom behaviour when something is passed in.
+"""
+sanitize_predicate(pred) = sanitize_predicate(GI.trait(pred), pred)
+sanitize_predicate(::Nothing, pred) = pred
+sanitize_predicate(::GI.AbstractTrait, pred) = sanitize_predicate(GI.extent(pred))
+sanitize_predicate(pred::Extents.Extent) = Base.Fix1(Extents.intersects, pred)
+
+
+end # module SpatialTreeInterface