reword filenames, add tests, update docs

Author: Vincent Landau

.gitignore

@@ -3,3 +3,4 @@ run-dev.sh
 Manifest.toml
 docs/build
 docs/make-local.jl
+local

Project.toml

@@ -12,5 +12,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
 ArchGDAL = "0.6"
-GeoData = "0.4.6"
 julia = "1.6"

docs/make.jl

@@ -15,7 +15,8 @@ makedocs(
     authors = "Vincent A. Landau",
     sitename = "SpatialGraphs.jl",
     pages = ["About" => "index.md",
-             "Graph Types" => "graphtypes.md"],
+             "Graph Types" => "graphtypes.md",
+             "User Guide" => "userguide.md"],
 )
 
 deploydocs(

docs/src/graphtypes.md

@@ -15,18 +15,19 @@ AbstractSpatialGraph
 The `AbstractRasterGraph` type is as subtype of `AbstractSpatialGraph`. All
 `AbstractRasterGraph` subtypes contain a field called `vertex_raster`, which is 
 a `GeoData.GeoArray` that describes the spatial locations for the vertices in 
-the graph. If your graph has 10 vertices, then the corresponding `GeoArray` 
+the graph. If your graph has 10 vertices, then the corresponding `vertex_raster` 
 must have 10 unique values, starting at 1 (which corresponds to the first vertex
 in the graph) and going up to 10 (which corresponds to the tenth vertex in the 
 graph). To find the location of a given graph vertex, _n_, you simply identify 
-the pixel(s) with a value of _n_.
+the pixel(s) with a value of _n_. For pixels/elements in `vertex_raster` for
+which there shouldn't be a graph vertex, use a value of 0.
 ```@docs
 AbstractRasterGraph
 ```
 
 ## AbstractRasterGraph Subtypes
 SpatialGraphs.jl has raster graph types for undirected, directed, unweighted,
-and weighted graphsm each detailed below.
+and weighted graphs, each detailed below.
 
 ```@docs
 SimpleRasterGraph

docs/src/userguide.md

@@ -0,0 +1,11 @@
+# User Guide
+
+## Building Graphs from Rasters
+
+SpatialGraphs.jl offers options methods for constructing graphs from raster
+data.
+
+```@docs
+weightedrastergraph
+make_raster_graph
+```

src/SpatialGraphs.jl

@@ -4,8 +4,13 @@ using LightGraphs, SimpleWeightedGraphs, ArchGDAL, GeoData
 
 include("structs.jl")
 include("graph_interface.jl")
+include("rastergraphs.jl")
+include("utils.jl")
 ## Types and Structs
 export AbstractSpatialGraph, AbstractRasterGraph, SimpleRasterGraph,
        SimpleRasterDiGraph, WeightedRasterGraph, WeightedRasterDiGraph
 
+## Raster graph
+export make_raster_graph, weightedrastergraph
+
 end # module

src/graph_interface.jl

@@ -1,5 +1,14 @@
 ## Methods for the LightGraphs ad SimpleWeightedGraphs interfaces
 ## With these methods defined, functions from LightGraphs should "just work"
+import LightGraphs: 
+    nv, ne, vertices, edges, eltype, edgetype, has_edge, has_vertex,
+    inneighbors, outneighbors, is_directed, add_edge!
+
+import SimpleWeightedGraphs:
+    add_edge!, get_weight
+
+import Base:
+    eltype, zero
 
 ### LightGraphs interface
 nv(g::AbstractSpatialGraph) = nv(g.graph)
@@ -14,3 +23,8 @@ inneighbors(g::AbstractSpatialGraph, v) = inneighbors(g.graph, v)
 outneighbors(g::AbstractSpatialGraph, v) = outneighbors(g.graph, v)
 is_directed(g::AbstractSpatialGraph) = is_directed(g.graph)
 Base.zero(g::AbstractSpatialGraph) = zero(g.graph)
+add_edge!(g::AbstractSpatialGraph, a::Integer, b::Integer, c::Number) = add_edge!(g.graph, a, b, c)
+
+### SimpleWeightedGraphs
+get_weight(g::WeightedRasterGraph, a::Integer, b::Integer) = get_weight(g.graph, a, b)
+get_weight(g::WeightedRasterDiGraph, a::Integer, b::Integer) = get_weight(g.graph, a, b)

src/constructors.jl renamed to src/rastergraphs.jl

@@ -1,14 +1,44 @@
+"""
+    make_vertex_raster(A::GeoArray)
+
+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 GeoArray. This 
+function is recommended for internal use only.
+
+## Parameters
+`A`: The GeoArray 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.
+"""
+function make_vertex_raster(A::GeoData.GeoArray)
+    # Make an array of unique node identifiers
+    nodemap = zeros(Int64, size(A.data))
+    is_node = (A.data .!= A.missingval) .&
+        ((!).(isnan.(A.data)))
+
+    nodemap[is_node] = 1:sum(is_node)
+
+    nodemap = GeoData.GeoArray(nodemap, dims(A))
+
+    nodemap
+end
+
+
 """
     weightedrastergraph(
-        A::GeoArray;
+        weight_raster::GeoArray;
         directed::Bool = false,
         condition::Function = is_data,
         cardinal_neighbors_only::Bool = false,
         connect_using_avg_raster_val::Bool = true
     )
 
 Construct a `WeightedRasterGraph` or `WeightedRasterDiGraph` (if 
-`directed = true`) from a raster dataset.
+`directed = true`) from a raster dataset. The weight raster,  denotes
+the edge weight correponding to each vertex. Since edges are between rather than
+on vertices, edge weights are calculated as the average of the weights for each
+vertex.
 
 ## Parameters
 `weight_raster`: A GeoData.GeoArray contained values that, where applicable 
@@ -53,27 +83,33 @@ If `false`, the weight between two nodes with resistances R1 and R2 is
 calculated by converting resistance to conductances, taking the average, then 
 taking the inverse of the result to convert back to resistance: 
 `1 / ((1/R1 + 1/R2) / 2)`. `connect_using_avg_weights = false` correspondes to 
-the default settings in Circuitscape. Defaults to `true', in which case the 
-simple average of the weights (adjusted for distance in the case of diagonal 
-neighbors) in `weight_raster` are used.
+the default settings in Circuitscape. Defaults to `true`, in which case the 
+simple average (adjusted for distance in the case of diagonal 
+neighbors) of the weights  in `weight_raster` is used.
+
+`combine`: In the case that there are multiple edges between a single pair
+of vertices, how should the weight be chosen? Defaults to `min`. See the docs
+for `SparseArrays.sparse()` for more information.
 """
 function weightedrastergraph(
         weight_raster::GeoArray;
         directed::Bool = false,
-        condition_raster::GeoArray = A,
+        condition_raster::GeoArray = weight_raster,
         condition::Function = is_data,
         cardinal_neighbors_only::Bool = false,
-        connect_using_avg_weights::Bool = true
+        connect_using_avg_weights::Bool = true,
+        combine::Function = min
     )
     v = make_vertex_raster(weight_raster)
-    g = make_weighted_graph(
+    g = make_raster_graph(
         weight_raster,
-        vertex_raster,
+        v,
         directed = directed,
         condition_raster = condition_raster,
         condition = condition,
         cardinal_neighbors_only = cardinal_neighbors_only,
-        connect_using_avg_weights = connect_using_avg_weights
+        connect_using_avg_weights = connect_using_avg_weights,
+        combine = combine
     )
 
     if directed
@@ -86,30 +122,106 @@ function weightedrastergraph(
 
 end
 
-function make_weighted_graph(
+"""
+    make_raster_graph(
+        weight_raster::GeoArray,
+        vertex_raster::GeoArray;
+        directed::Bool = false,
+        condition_raster::GeoArray = weight_raster,
+        condition::Function = is_data,
+        cardinal_neighbors_only::Bool = false,
+        connect_using_avg_weights::Bool = true,
+        combine::Function = min
+    )
+
+Construct a `SimpleWeightedGraph` or `SimpleWeightedDiGraph` (if 
+`directed = true`) based on vertex and weight rasters. This function is useful 
+if you already have a custom vertex raster and don't want SpatialGraphs.jl to
+make one for you. The vertex raster denotes the spatial locations of each vertex
+in the graph, and the weight raster denotes the edge weight correponding to each
+vertex. Since edges are between rather than on vertices, edge weights are
+calculated as the average of the weights for each vertex being connected.
+
+## Parameters
+`weight_raster`: A `GeoData.GeoArray` containing values that, where applicable 
+based on other arguments, determine which pixels to connect and the edge 
+weights between pixels. Any pixel in `weight_raster` with a value not equal to 
+`weight_raster.missingval` will be assigned a vertex in the graph (corresponding
+to its centroid). 
+
+`vertex_raster`: A `GeoData.GeoArray` with integer values ranging from 1:n, 
+where n is the number of unique vertices in the graph. 
+
+##Arguments
+`directed`: A `Bool` determining whether the graph should be directed.
+
+`condition_raster`: A raster with values that can be used to determine whether
+two neighboring pixels should be connected. For example, an elevation raster 
+can be used to create a hydologic flow graph.
+
+`condition`: A function that compares the values in `condition_raster` for two
+neighbors to determine if those neighbors should be connected. The function must
+compare two values and return either `true` or `false`. Useful functions to use 
+here include `<`, `<=`, `==`, etc. The first argument to `condition` corresponds
+to the source vertex, and the second argument corresponds to the destination 
+vertex. So, if you only want to connect sources to destinations with a lower
+value in `condition_raster` (e.g. in the case of developing a hydrologic flow 
+graph based on elevation), then you would use `<` for `condition`. Defaults to 
+`is_data`, which results in neighbors being connected as long as they are not 
+NoData (`condition_raster.missingval`) in `condition_raster`. Note that if using
+an inequality function (or any function where the result depends on argument 
+position), the graph must be directed. For undirected graphs, you can use either
+`is_data` or `==`, or any other custom function where argument position doesn't
+matter, e.g. a function that determines whether the values in `condition_raster`
+are within a certain distance of each other.
+
+`cardinal_neighbors_only`: A `Bool` stating whether only cardinal neighbors
+should be connected. By default, both cardinal _and_ diagonal neighbors are
+connected. Note that when determining weights between diagonal neighbors, the
+increased distance between them (as compared to the distance between cardinal
+neighbors) is accounted for.
+
+`connect_using_avg_raster_val`: `Bool`. This is intended to offer methods that
+complement those used in Circuitscape.jl and Omniscape.jl. In this context,
+weights (the values in `weight_raster`) are in units of electrical resistance. 
+If `false`, the weight between two nodes with resistances R1 and R2 is 
+calculated by converting resistance to conductances, taking the average, then 
+taking the inverse of the result to convert back to resistance: 
+`1 / ((1/R1 + 1/R2) / 2)`. `connect_using_avg_weights = false` correspondes to 
+the default settings in Circuitscape. Defaults to `true', in which case the 
+simple average of the weights (adjusted for distance in the case of diagonal 
+neighbors) in `weight_raster` are used.
+
+`combine`: In the case that there are multiple edges defined for a single pair
+of vertices, how should the weight be chosen? Defaults to `min`. See the docs
+for `SparseArrays.sparse()` for more information.
+"""
+function make_raster_graph(
         weight_raster::GeoArray,
         vertex_raster::GeoArray;
         directed::Bool = false,
         condition_raster::GeoArray = weight_raster,
         condition::Function = is_data,
         cardinal_neighbors_only::Bool = false,
-        connect_using_avg_weights::Bool = true
+        connect_using_avg_weights::Bool = true,
+        combine::Function = min
     )
     # Which averaging function to use
     card_avg = connect_using_avg_weights ? res_cardinal_avg : cond_cardinal_avg
     diag_avg = connect_using_avg_weights ? res_diagonal_avg : cond_diagonal_avg
-    dims = size(cost)
-    not_no_data = cost .!= no_data_val
+    dims = size(weight_raster)
+    no_data_val = weight_raster.missingval
+    not_no_data = weight_raster .!= no_data_val
 
-    if sum((cost .<= 0) .& not_no_data) != 0
-        @error("cost surface contains 0 or negative values (aside " * 
-                "from the provided no_data_val), which is not supported")
+    if sum((weight_raster .<= 0) .& not_no_data) != 0
+        @error("weight_raster contains 0 or negative values (aside " * 
+                "from weight_raster.missingval), which is not supported")
     end
 
     sources = Vector{Int64}()
     destinations = Vector{Int64}()
     node_weights = Vector{Float64}()
-
+    
     # Add the edges
     # Only need to do neighbors down or to the right for undirected graphs
     # because edge additions will be redundant.
@@ -120,7 +232,7 @@ function make_weighted_graph(
         for column in 1:dims[2]
             source_idx = CartesianIndex((row, column))
             if vertex_raster[source_idx] == 0 || 
-                    cost[source_idx] == weight_raster.missingval
+                    weight_raster[source_idx] == no_data_val
                 continue
             else
                 # Cardinal destination indices needed for undirected graph
@@ -269,16 +381,16 @@ function make_weighted_graph(
             sources,
             destinations,
             node_weights,
-            combine = min
+            combine = combine
         )
     else    
         g = SimpleWeightedGraph(
             sources,
             destinations,
             node_weights,
-            combine = min
+            combine = combine
         )
     end
 
     return g
-end
+end

src/utils.jl

@@ -11,3 +11,21 @@ Returns `true`
 function is_data(a::Number, b::Number)
     return true
 end
+
+
+# averaging functions, to operate on resistances
+# cond_* functions convert to conductance, calulate mean, then convert to resistance
+# res_* functions directly operate on the resistances
+cond_cardinal_avg(x, y) = 1 / ((1/x + 1/y) / 2)
+cond_diagonal_avg(x, y) = 1 / ((1/x + 1/y) / (2 * √2))
+res_cardinal_avg(x, y) = (x + y) / 2
+res_diagonal_avg(x, y) = ((x + y) * √2) / 2
+
+
+# function get_dims(A::GeoData.GeoArray)
+#     y_first = dims(mytif)[1] isa XDim
+#     first_dim_type = y_first ? YDim : XDim
+#     second_dim_type = y_first ? XDim : YDim
+
+#     (dims(A, first_dim_type), dims(A, second_dim_type), Band(1:1, Categorical(order = Ordered())))
+# end