Merge pull request #9 from Circuitscape/vl/patches

handle "habitat patches" in which a vertex takes up more than one pixel

Project.toml

@@ -1,7 +1,7 @@
 name = "SpatialGraphs"
 uuid = "c8109be6-b162-4503-8829-8b6b29b93f50"
 authors = ["Vincent A. Landau <[email protected]>"]
-version = "0.1.0"
+version = "0.1.1"
 
 [deps]
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"

docs/src/userguide.md

@@ -21,4 +21,18 @@ make_simple_raster_graph
 ```@docs
 weightedrastergraph
 make_weighted_raster_graph
+```
+
+### Making vertex rasters
+In SpatialGraphs.jl, vertex rasters serve as the spatial reference for the 
+vertices in a graph. SpatialGraphs.jl provides functions to generate these
+rasters. Often, it is done internally and the end user doesn't need to use these
+functions, but there are some cases where it will be necessary. For example, it 
+is often the case that a user may want a single graph vertex to occupy multiple
+pixels in space (e.g. when modeling habitat connectivity between two
+protected areas). SpatialGraphs.jl enables this by offering a method for the 
+`make_vertex_raster` function (below) that accepts a raster
+representing these patches as an argument.
+```@docs
+make_vertex_raster
 ```

src/SpatialGraphs.jl

@@ -14,6 +14,6 @@ export AbstractSpatialGraph, AbstractRasterGraph, RasterGraph,
 
 ## Raster graph
 export make_weighted_raster_graph, weightedrastergraph,
-    make_simple_raster_graph, rastergraph
+    make_simple_raster_graph, rastergraph, make_vertex_raster
 
 end # module

src/rastergraphs.jl

@@ -2,11 +2,10 @@
     make_vertex_raster(A::Raster)
 
 Constuct a vertex raster (a raster where the value of each pixel corresponds
-to its ID in a graph, and 0s correspond to NoData). Returns a `Raster``. This 
-function is recommended for internal use only.
+to its ID in a graph, and 0s correspond to NoData). Returns a `Raster`.
 
 ## Parameters
-`A`: The `Raster`` from which a graph will be built, which is used as the
+`A`: The `Raster` from which a graph will be built, which is used as the
 reference for building the vertex raster. Pixels with NoData (`A.missingval`)
 are skipped (no vertex is assigned). Pixels with NoData will get a value of 0 in
 the vertex raster.
@@ -25,6 +24,53 @@ function make_vertex_raster(A::Raster)
 end
 
 
+"""
+    make_vertex_raster(A::Raster, patches::Raster)
+
+Constuct a vertex raster (a raster where the value of each pixel corresponds
+to its ID in a graph, and 0s correspond to NoData). Returns a `Raster`.
+
+## Parameters
+`A`: The `Raster` from which a graph will be built, which is used as the
+reference for building the vertex raster. Pixels with NoData (`A.missingval`)
+are skipped (no vertex is assigned). Pixels with NoData will get a value of 0 in
+the vertex raster.
+
+`patches`: A `Raster` defining patches that may consist of multiple pixels
+but should be considered a single vertex in the graph. patches must have the 
+same dimensions a `A`. Each patch should have its own unique value, starting 
+with 1, and ending with n, where n is the total number of patches (e.g. if you 
+have 3 patches, pixels in the first patch should have a value of 1, the second
+patch a value of 2, and the third patch a value of 3). Non-patch pixels can 
+either be labeled with 0 of the raster's `missingval`.
+"""
+function make_vertex_raster(A::Raster, patches::Raster)
+    # Make an array of unique node identifiers 
+    nodemap = zeros(Int64, size(A.data))
+    is_node = (A.data .!= A.missingval) .&
+        ((!).(isnan.(A.data)))
+
+    # Start by filling in the patches
+
+    # Check for correctness of patch raster
+    patch_vals = patches.data[(patches.data .!= patches.missingval) .& (patches.data .!= 0)]
+    if sort(unique(patch_vals)) != 1:maximum(patch_vals)
+        throw(ArgumentError("Patches must be labeled with numbers from 1 to n, where n is the number of patches, with no numbers skipped."))
+    end
+
+    n_patch = maximum(patches)
+    for i in 1:n_patch
+        nodemap[is_node .& (patches.data .== i)] .= i
+    end
+
+    is_matrix_node = is_node .& (nodemap .== 0)
+    nodemap[is_matrix_node] .= (n_patch + 1):(n_patch + sum(is_matrix_node))
+    nodemap = Raster(nodemap, dims(A), missingval = 0)
+
+    nodemap
+end
+
+
 """
     weightedrastergraph(
         weight_raster::Raster;
@@ -117,6 +163,7 @@ function weightedrastergraph(
     return sg
 end
 
+
 """
     rastergraph(
         raster::Raster;
@@ -185,6 +232,7 @@ function rastergraph(
     return sg
 end
 
+
 """
     make_weighted_raster_graph(
         weight_raster::Raster,

test/patches.jl

@@ -0,0 +1,47 @@
+using Rasters, SpatialGraphs, Test, Graphs, SimpleWeightedGraphs
+
+resistance = [
+    1. 2.5 6 0.2;
+    5 6 12 1;
+    0.5 3 7 4;
+    11 0.5 9 -9999
+]
+
+res = Raster(
+    reshape(resistance, (size(resistance)..., 1)),
+    missingval=-9999,
+    dims=(X(1:4), Y(1:4), Band(1:1))
+)
+
+patches = [
+    1 1 0 3;
+    0 0 0 3;
+    0 0 0 0;
+    0 2 2 0
+]
+
+pat = Raster(
+    reshape(patches, (size(patches)..., 1)),
+    missingval=-9999,
+    dims=(X(1:4), Y(1:4), Band(1:1))
+)
+
+vertex_raster = make_vertex_raster(res, pat)
+
+# Check that patches have proper vertex ID
+for i in 1:maximum(pat)
+   @test all(vertex_raster[pat .== i] .== i)
+end
+
+@test sort(unique(vertex_raster[vertex_raster .> maximum(pat)])) ==
+    (maximum(pat) + 1):(maximum(pat) + length(res[(pat .== 0) .& (res .!= res.missingval)]))
+
+## Make sure that a graph can be successfully created from a the vertex raster
+g = WeightedRasterGraph(make_weighted_raster_graph(res, vertex_raster), vertex_raster)
+
+# Test that a few random weights are correct, assumes that defaults were used
+# in make_weighted_raster_graph above, particularly that combine = min and
+# connect_using_avg_weights = true. THESE ARE HARD CODED. Changes to res or pat
+# above will result in failed tests.
+@test get_weight(g.graph, 1, 4) == 3 # minimum of the mean resistances between both vertices labeled 1 and vertex labeled 4
+@test get_weight(g.graph, 7, 11) == SpatialGraphs.res_diagonal_avg(res[vertex_raster .== 7], res[vertex_raster .== 11])[1]

test/runtests.jl

@@ -20,6 +20,10 @@ end
     include("graph_interface.jl")
 end
 
+@testset "Graphs with patches" begin
+    include("patches.jl")
+end
+
 printstyled("Checking that Base.show works...\n", bold = true)
 
 A_array = Array{Float64}(undef, (3, 4, 1))