add GeoInterface/Extents compat methods (#22)

* add GeoInterface/Extents compat methods

* docs tweaks

* bump julia version to 1.6

* minor additions

- add Manifest to .gitignore
- use Aqua quality assurance
- sort Project.toml
- remove debugging statement

---------

Co-authored-by: Martijn Visser <[email protected]>

.github/workflows/CI.yml

@@ -13,7 +13,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
+          - '1.6'
           - '1'
           - 'nightly'
         os:

.gitignore

@@ -1,5 +1,4 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
-/deps/deps.jl
-/deps/usr/
+Manifest.toml

Project.toml

@@ -4,14 +4,18 @@ license = "MIT"
 version = "0.2.0"
 
 [deps]
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 LibSpatialIndex_jll = "00e98e2a-4326-5239-88cb-15dcbe1c18d0"
 
 [compat]
-julia = "1.3"
+Aqua = "0.7"
+GeoInterface = "1"
 LibSpatialIndex_jll = "1.8"
+julia = "1.6"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["Aqua", "Test"]

README.md

@@ -8,53 +8,82 @@ LibSpatialIndex.jl is a julia wrapper around the C API of [libspatialindex](http
 # Quick Guide
 
 A new RTree with 2 dimensions can be created using this package as follows:
+
 ```julia
 import LibSpatialIndex
 rtree = LibSpatialIndex.RTree(2)
 ```
 
 ## Insertion
-Items can be inserted using the `insert!` method, where
+
+Items can be inserted using the `insert!` method, where:
+
 ```julia
-LibSpatialIndex.insert!(rtree, 1, [0.,0.], [1.,1.])
-LibSpatialIndex.insert!(rtree, 2, [0.,0.], [2.,2.])
+LibSpatialIndex.insert!(rtree, 1, [0.0, 0.0], [1.0, 1.0])
+LibSpatialIndex.insert!(rtree, 2, Extexts.Extent(X=(0.0, 2.0), (0.0, 2.0)))
 ```
+
 inserts two items,
 
 - the first with `id` 1, associated with the box specified by `[xmin=0.0,ymin=0.0]` and `[xmax=1.0,ymax=1.0]`.
 - the second with `id` 2, associated with the box specified by `[xmin=0.0,ymin=0.0]` and `[xmax=2.0,ymax=2.0]`.
 
 ## Queries
-Thereafter, you can perform queries on the `rtree` using either (i) `intersects(rtree, minvalues, maxvalues)` for all items intersecting the box specified by `minvalues` and `maxvalues`, or (ii) `knn(rtree, minvalues, maxvalues, k)` for the `k` nearest items in `rtree` to the box specified by `minvalues` and `maxvalues`.
+
+Thereafter, you can perform queries on the `rtree` using either (i) `intersects(rtree, minvalues, maxvalues)`
+for all items intersecting the box specified by `minvalues` and `maxvalues`, or (ii)
+`knn(rtree, minvalues, maxvalues, k)` for the `k` nearest items in `rtree` to the box
+specified by `minvalues` and `maxvalues`.
 
 ### Intersection
-So for instance,
+
+So for instance:
+
 ```julia
-LibSpatialIndex.intersects(rtree, [0.,0.],[1.,1.])
+LibSpatialIndex.intersects(rtree, [0.0, 0.0], [1.0, 1.0])
+LibSpatialIndex.intersects(rtree, Extents.extent(X=(0.0, 1.0), Y=(0.0, 1.0)))
 ```
-will return the vector `[1,2]` on the `rtree` constructed earlier, to indicate that items with ids `1` and `2` intersects the box specified by `[xmin=0.0,ymin=0.0]` and `[xmax=1.0,ymax=1.0]`.
 
-You can also perform queries on an individual point, so 
+will return the vector `[1, 2]` on the `rtree` constructed earlier, to indicate that items
+with ids `1` and `2` intersects the box specified by `[xmin=0.0, ymin=0.0]` and `[xmax=1.0, ymax=1.0]`.
+
+Any GeoInterface.jl or Extents.jl compatible object can be used directly instead
+of defining the `Extent` manually - the extent will either be detected or calculated.
+
+
+You can also perform queries on any individual GeoInterface.jl compatible point, so: 
+
 ```julia
-LibSpatialIndex.intersects(rtree, [1.,1.])
+LibSpatialIndex.intersects(rtree, (1.0, 1.0))
 ```
-will return the ids `[1,2]` in the `rtree` constructed earlier, and
+
+will return the ids `[1, 2]` in the `rtree` constructed earlier, and:
+
 ```julia
-LibSpatialIndex.intersects(rtree, [2.,2.])
+LibSpatialIndex.intersects(rtree, [2.0, 2.0])
 ```
-will only return the vector `[2]`, because item 1 does not contain the point `[2,2]`.
+
+will only return the vector `[2]`, because item 1 does not contain the point `[2.0, 2.0]`.
 
 ### k Nearest Neighbors
-For `knn` queries,
+
+For `knn` queries:
+
 ```julia
-LibSpatialIndex.knn(rtree, [2.,2.], 1)
+LibSpatialIndex.knn(rtree, [2.0, 2.0], 1)
 ```
+
 returns the vector `[2]` because the item with id `2` is closest to the point `[2.0, 2.0]`, and
+
 ```julia
-sort(LibSpatialIndex.knn(rtree, [2.,2.], 2))
+sort(LibSpatialIndex.knn(rtree, [2.0, 2.0], 2))
 ```
-returns the vector `[1,2]`. If the value of `k` exceeds the number of items in the `rtree`, then fewer than `k` items will be returned, so
+
+returns the vector `[1, 2]`. If the value of `k` exceeds the number of items in the `rtree`,
+then fewer than `k` items will be returned, so:
+
 ```julia
-sort(SI.knn(rtree, [2.,2.], 3))
+sort(SI.knn(rtree, [2.0, 2.0], 3))
 ```
-will return the vector `[1,2]`.
+
+will return the vector `[1, 2]`.

src/LibSpatialIndex.jl

@@ -1,5 +1,7 @@
 module LibSpatialIndex
 
+    import GeoInterface as GI
+
     include("capi.jl")
 
     version() = unsafe_string(C.SIDX_Version())
@@ -11,6 +13,8 @@ module LibSpatialIndex
     end
 
     """
+        RTree(ndim::Integer; kw...)
+
     The RTree index [guttman84] is a balanced tree structure that consists of
     index nodes, leaf nodes and data.
 
@@ -22,9 +26,12 @@ module LibSpatialIndex
     They cannot be empty though. A `fillfactor` specifies the minimum number of
     entries allowed in any node. The fill factor is usually close to `70%`.
 
-    # Options
+    ## Arguments
+
+    `ndim`: Dimensionality of the data that will be inserted.
+
+    ## Keywords
 
-    * `ndim`: Dimensionality of the data that will be inserted.
     * `indextype`: one of `RT_RTree` (default), `RT_MVRTree`, or `RT_TPRTree`.
     * `variant`: one of `RT_Linear`, `RT_Quadratic`, or `RT_Star` (default).
     * `storage`: one of `RT_Memory` (default), `RT_Disk`, or `RT_Custom`.
@@ -39,7 +46,7 @@ module LibSpatialIndex
     * `fillfactor`: The fill factor. Default is `0.7`.
     * `splitdistributionfactor`: Default is `0.4`.
     * `reinsertfactor`: Default is `0.3`.
-    
+
     # Performance
 
     Dataset size, data density, etc. have nothing to do with capacity and page
@@ -152,12 +159,22 @@ module LibSpatialIndex
     end
 
     """
+        insert!(rtree::RTree, id::Integer, minvalues::Vector{Float64}, maxvalues::Vector{Float64})
+        insert!(rtree::RTree, id::Integer, extent::Extent)
+        insert!(rtree::RTree, id::Integer, obj)
+
     Inserts an item into the `rtree` with given `id` and boundingbox specified
     by `minvalues` and `maxvalues`, where the item lies within the interval
-    `[minvalues[i],maxvalues[i]]` for each axis `i` in 1, ..., `ndim`.
+    `[minvalues[i], maxvalues[i]]` for each axis `i` in 1, ..., `ndim`,
+    or similar for the `Extent` of `obj`.
+
+    If `obj` is passed it will be detected as a `GeoInterface.PointTrait`
+    and used as a point, or otherwise `GeoInterface.extent` will be called to
+    detect or calculate the objects `Extent`, falling back to `Extents.extent`.
+
+    In these cases `minvalues` and `maxvalues` are taken from the point or extent.
     """
-    function insert!(
-            rtree::RTree,
+    function insert!(rtree::RTree,
             id::Integer,
             minvalues::Vector{Float64},
             maxvalues::Vector{Float64}
@@ -166,13 +183,30 @@ module LibSpatialIndex
             pointer(maxvalues), UInt32(length(minvalues)), Ptr{UInt8}(0), Cint(0)
         )
     end
+    function insert!(rtree::RTree, id::Integer, extent::GI.Extent)
+        insert!(rtree, id, _ext2vecs(extent)...)
+    end
+    insert!(rtree::RTree, id::Integer, extent::Nothing) = _not_point_or_ext_error()
+    function insert!(rtree::RTree, id::Integer, obj)
+        insert!(rtree, id, GI.extent(obj))
+    end
 
     """
+        intersects(rtree::RTree, minvalues::Vector{Float64}, maxvalues::Vector{Float64})
+        intersects(rtree::RTree, extent::Extent)
+        intersects(rtree::RTree, obj)
+
     Returns a vector of `id`s corresponding to items in `rtree` that intersects
     the box specified by `minvalues` and `maxvalues`.
 
-    Each item intersects the interval `[minvalues[i],maxvalues[i]]` for each
-    axis `i` in 1, ..., `ndim`.
+    Each item intersects the interval `[minvalues[i], maxvalues[i]]` for each
+    axis `i` in 1, ..., `ndim`, or similar for the `Extent` of `obj`.
+
+    If `obj` is passed it will be detected as a `GeoInterface.PointTrait`
+    and used as a point, or otherwise `GeoInterface.extent` will be called to
+    detect or calculate the objects `Extent`, falling back to `Extents.extent`.
+
+    In these cases `minvalues` and `maxvalues` are taken from the point or extent.
     """
     function intersects(
             rtree::RTree,
@@ -187,21 +221,38 @@ module LibSpatialIndex
         _checkresult(result, "Index_Intersects_id: Failed to evaluate")
         unsafe_wrap(Array, items[], nresults[])
     end
-
-    """
-    Returns a vector of `id`s corresponding to items in `rtree` that intersects
-    the coordinates specified by `point`.
-    """
     intersects(rtree::RTree, point::Vector{Float64}) = intersects(rtree, point, point)
+    function intersects(rtree::RTree, obj)
+        if GI.trait(obj) isa GI.PointTrait
+            intersects(rtree, _point2vec(obj))
+        else
+            intersects(rtree, GI.extent(obj))
+        end
+    end
+    function intersects(rtree::RTree, extent::GI.Extent)
+        intersects(rtree::RTree, _ext2vecs(extent)...)
+    end
+    intersects(rtree::RTree, ::Nothing) = _not_point_or_ext_error()
 
     """
+        knn(rtree::RTree, minvalues::Vector{Float64}, maxvalues::Vector{Float64}, k::Integer)
+        knn(rtree::RTree, extent::Extent, k::Integer)
+        knn(rtree::RTree, obj, k::Integer)
+
     Returns a vector of `id`s corresponding to the `k` items in `rtree`
     that are nearest to the box specified by `minvalues` and `maxvalues`.
+    Each item intersects the interval `[minvalues[i], maxvalues[i]]` for each
+    axis `i` in 1, ..., `ndim`, or similar for the `Extent` of `obj`.
 
-    Each item intersects the interval `[minvalues[i],maxvalues[i]]` for each
-    axis `i` in 1, ..., `ndim`. If there are fewer than `k` items in `rtree`,
+    If there are fewer than `k` items in `rtree`,
     it will return less than `k` items. On the other hand, if there are ties
     between some of the items, it might return more than `k` items.
+
+    If `obj` is passed it will be detected as a `GeoInterface.PointTrait`
+    and used as a point, or otherwise `GeoInterface.extent` will be called to
+    detect or calculate the objects `Extent`, falling back to `Extents.extent`.
+
+    In these cases `minvalues` and `maxvalues` are taken from the point or extent.
     """
     function knn(
             rtree::RTree,
@@ -216,15 +267,38 @@ module LibSpatialIndex
         _checkresult(result, "Index_NearestNeighbors_id: Failed to evaluate")
         unsafe_wrap(Array, items[], nresults[])
     end
+    knn(rtree::RTree, point::Vector{Float64}, k::Integer) = knn(rtree, point, point, k)
+    knn(rtree::RTree, extent::GI.Extent, k::Integer) = knn(rtree::RTree, _ext2vecs(extent)..., k)
+    knn(rtree::RTree, extent::Nothing, k::Integer) = _not_point_or_ext_error()
+    function knn(rtree::RTree, obj, k::Integer)
+        if GI.trait(obj) isa GI.PointTrait
+            knn(rtree, _point2vec(obj), k)
+        else
+            knn(rtree, GI.extent(obj), k)
+        end
+    end
 
-    """
-    Returns a vector of `id`s corresponding to the `k` items in `rtree`
-    that are nearest to the box specified by `minvalues` and `maxvalues`.
+    # Utils
+    function _ext2vecs(ex::GI.Extent)
+        haskey(ex, :X) && haskey(ex, :Y) || throw(ArgumentError("Extent does not have X and Y keys"))
+
+        min, max = if haskey(ex, :Z)
+            Float64[ex.X[1], ex.Y[1], ex.Z[1]], Float64[ex.X[2], ex.Y[2], ex.Z[1]]
+        else
+            Float64[ex.X[1], ex.Y[1]], Float64[ex.X[2], ex.Y[2]]
+        end
+
+        return min, max
+    end
+
+    function _point2vec(p)
+        if GI.is3d(p)
+            Float64[GI.x(p), GI.y(p), GI.z(p)]
+        else
+            Float64[GI.x(p), GI.y(p)]
+        end
+    end
+
+    _not_point_or_ext_error() = throw(ArgumentError("object is not a point, and does not have an extent"))
 
-    If there are fewer than `k` items in `rtree`, it will return less than `k`
-    items. On the other hand, if there are ties between some of the items,
-    it might return more than `k` items.
-    """
-    knn(rtree::RTree, point::Vector{Float64}, k::Integer) = knn(rtree, point, point, k)
-
 end # module

test/runtests.jl

@@ -1,6 +1,8 @@
 using LibSpatialIndex
 using Test
-const SI = LibSpatialIndex
+import GeoInterface as GI
+import LibSpatialIndex as SI
+import Aqua
 
 @testset "Simple Tutorial" begin
     # based on https://github.com/libspatialindex/libspatialindex/wiki/Simple-Tutorial
@@ -51,3 +53,23 @@ end
     @test sort(SI.knn(rtree, [2.,2.], 1)) == [2]
     @test sort(SI.knn(rtree, [2.,2.], 2)) == [1,2]
 end
+
+@testset "GeoInterface/Extents Operations" begin
+    rtree = SI.RTree(2)
+    result = SI.insert!(rtree, 1, GI.Extent(X=(0.0, 1.0), Y=(0.0, 1.0)))
+    @test result == SI.C.RT_None
+
+    polygon = GI.Polygon([GI.LinearRing([(0.0, 0.0), (0.5, 0.0), (2.0, 0.5), (0.0, 2.0), (0.0, 0.0)])])
+    result = SI.insert!(rtree, 2, polygon)
+
+    @test result == SI.C.RT_None
+    @test SI.intersects(rtree, GI.LineString([(0.0, 0.0), (1.0, 1.0)])) == [1, 2]
+    @test SI.intersects(rtree, GI.Point(0.0, 0.0)) == [1, 2]
+    @test SI.intersects(rtree, (X=2.0, Y=2.0)) == [2]
+    @test sort(SI.knn(rtree, GI.Extent(X=(2.0, 2.0), Y=(2.0, 2.0)), 1)) == [2]
+    @test sort(SI.knn(rtree, (2.0, 2.0), 1)) == [2]
+end
+
+@testset "Aqua" begin
+    Aqua.test_all(LibSpatialIndex)
+end