Skip to content

Commit 99d373d

Browse files
author
Vincent Landau
committed
start on constructor methods
1 parent 4800c72 commit 99d373d

File tree

3 files changed

+298
-1
lines changed

3 files changed

+298
-1
lines changed

‎docs/make.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@ const formats = Any[
55
assets = [
66
"assets/custom.css"
77
],
8-
edit_link = :commit,
8+
edit_link = :commit,
99
),
1010
]
1111

‎src/constructors.jl

Lines changed: 284 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,284 @@
1+
"""
2+
weightedrastergraph(
3+
A::GeoArray;
4+
directed::Bool = false,
5+
condition::Function = is_data,
6+
cardinal_neighbors_only::Bool = false,
7+
connect_using_avg_raster_val::Bool = true
8+
)
9+
10+
Construct a `WeightedRasterGraph` or `WeightedRasterDiGraph` (if
11+
`directed = true`) from a raster dataset.
12+
13+
## Parameters
14+
`weight_raster`: A GeoData.GeoArray contained values that, where applicable
15+
based on other arguments, determines which pixels to connect and the edge
16+
weights between pixels. Any pixel in `weight_raster` with a value not equal to
17+
`weight_raster.missingval` will be assigned a vertex in the graph (corresponding
18+
to its centroid).
19+
20+
##Arguments
21+
`directed`: A `Bool` determining whether the graph should be directed.
22+
23+
`condition_raster`: A raster with values that can be used to determine whether
24+
two neighboring pixels should be connected. For example, an elevation raster
25+
can be used to create a hydologic flow graph.
26+
27+
`condition`: A function that compares the values in `condition_raster` for two
28+
neighbors to determine if those neighbors should be connected. The function must
29+
compare two values and return either `true` or `false`. Useful functions to use
30+
here include `<`, `<=`, `==`, etc. The first argument to `condition` corresponds
31+
to the source vertex, and the second argument corresponds to the destination
32+
vertex. So, if you only want to connect sources to destinations with a lower
33+
value in `condition_raster` (e.g. in the case of developing a hydrologic flow
34+
graph based on elevation), then you would use `<` for `condition`. Defaults to
35+
`is_data`, which results in neighbors being connected as long as they are not
36+
NoData (`condition_raster.missingval`) in `condition_raster`. Note that if using
37+
an inequality function (or any function where the result depends on argument
38+
position), the graph must be directed. For undirected graphs, you can use either
39+
`is_data` or `==`, or any other custom function where argument position doesn't
40+
matter, e.g. a function that determines whether the values in `condition_raster`
41+
are within a certain distance of each other.
42+
43+
`cardinal_neighbors_only`: A `Bool` stating whether only cardinal neighbors
44+
should be connected. By default, both cardinal _and_ diagonal neighbors are
45+
connected. Note that when determining weights between diagonal neighbors, the
46+
increased distance between them (as compared to the distance between cardinal
47+
neighbors) is accounted for.
48+
49+
`connect_using_avg_raster_val`: `Bool`. This is intended to offer methods that
50+
complement those used in Circuitscape.jl and Omniscape.jl. In this context,
51+
weights (the values in `weight_raster`) are in units of electrical resistance.
52+
If `false`, the weight between two nodes with resistances R1 and R2 is
53+
calculated by converting resistance to conductances, taking the average, then
54+
taking the inverse of the result to convert back to resistance:
55+
`1 / ((1/R1 + 1/R2) / 2)`. `connect_using_avg_weights = false` correspondes to
56+
the default settings in Circuitscape. Defaults to `true', in which case the
57+
simple average of the weights (adjusted for distance in the case of diagonal
58+
neighbors) in `weight_raster` are used.
59+
"""
60+
function weightedrastergraph(
61+
weight_raster::GeoArray;
62+
directed::Bool = false,
63+
condition_raster::GeoArray = A,
64+
condition::Function = is_data,
65+
cardinal_neighbors_only::Bool = false,
66+
connect_using_avg_weights::Bool = true
67+
)
68+
v = make_vertex_raster(weight_raster)
69+
g = make_weighted_graph(
70+
weight_raster,
71+
vertex_raster,
72+
directed = directed,
73+
condition_raster = condition_raster,
74+
condition = condition,
75+
cardinal_neighbors_only = cardinal_neighbors_only,
76+
connect_using_avg_weights = connect_using_avg_weights
77+
)
78+
79+
if directed
80+
sg = WeightedRasterDiGraph(g, v)
81+
else
82+
sg = WeightedRasterGraph(g, v)
83+
end
84+
85+
return sg
86+
87+
end
88+
89+
function make_weighted_graph(
90+
weight_raster::GeoArray,
91+
vertex_raster::GeoArray;
92+
directed::Bool = false,
93+
condition_raster::GeoArray = weight_raster,
94+
condition::Function = is_data,
95+
cardinal_neighbors_only::Bool = false,
96+
connect_using_avg_weights::Bool = true
97+
)
98+
# Which averaging function to use
99+
card_avg = connect_using_avg_weights ? res_cardinal_avg : cond_cardinal_avg
100+
diag_avg = connect_using_avg_weights ? res_diagonal_avg : cond_diagonal_avg
101+
dims = size(cost)
102+
not_no_data = cost .!= no_data_val
103+
104+
if sum((cost .<= 0) .& not_no_data) != 0
105+
@error("cost surface contains 0 or negative values (aside " *
106+
"from the provided no_data_val), which is not supported")
107+
end
108+
109+
sources = Vector{Int64}()
110+
destinations = Vector{Int64}()
111+
node_weights = Vector{Float64}()
112+
113+
# Add the edges
114+
# Only need to do neighbors down or to the right for undirected graphs
115+
# because edge additions will be redundant.
116+
# Cardinal directions are in quotes since sometimes GeoArray are permuted
117+
# such that the top row doesn't necessarily correspond to the northern-most
118+
# pixels
119+
for row in 1:dims[1]
120+
for column in 1:dims[2]
121+
source_idx = CartesianIndex((row, column))
122+
if vertex_raster[source_idx] == 0 ||
123+
cost[source_idx] == weight_raster.missingval
124+
continue
125+
else
126+
# Cardinal destination indices needed for undirected graph
127+
east_idx = CartesianIndex((row, column + 1))
128+
south_idx = CartesianIndex((row + 1, column))
129+
130+
## Add cardinal neighbors
131+
# "East"
132+
if column != dims[2] && vertex_raster[east_idx] != 0 &&
133+
weight_raster[east_idx] != no_data_val
134+
if condition(condition_raster[source_idx],
135+
condition_raster[east_idx])
136+
res = card_avg(weight_raster[source_idx],
137+
weight_raster[east_idx])
138+
push!(sources, vertex_raster[source_idx])
139+
push!(destinations, vertex_raster[east_idx])
140+
push!(node_weights, res)
141+
end
142+
end
143+
# "South"
144+
if row != dims[1] && vertex_raster[south_idx] != 0 &&
145+
weight_raster[south_idx] != no_data_val
146+
if condition(condition_raster[source_idx],
147+
condition_raster[south_idx])
148+
res = card_avg(weight_raster[source_idx],
149+
weight_raster[south_idx])
150+
push!(sources, vertex_raster[source_idx])
151+
push!(destinations, vertex_raster[south_idx])
152+
push!(node_weights, res)
153+
end
154+
end
155+
156+
## Add diagonal neighbors if needed
157+
if !cardinal_neighbors_only
158+
# Diagonal destination indices needed for undirected graph
159+
ne_idx = CartesianIndex((row - 1, column + 1))
160+
se_idx = CartesianIndex((row + 1, column + 1))
161+
162+
# "Northeast"
163+
if column != dims[2] && row != 1 &&
164+
vertex_raster[ne_idx] != 0 &&
165+
weight_raster[ne_idx] != no_data_val
166+
if condition(condition_raster[source_idx],
167+
condition_raster[ne_idx])
168+
res = diag_avg(weight_raster[source_idx],
169+
weight_raster[ne_idx])
170+
push!(sources, vertex_raster[source_idx])
171+
push!(destinations, vertex_raster[ne_idx])
172+
push!(node_weights, res)
173+
end
174+
end
175+
# "Southeast"
176+
if row != dims[1] && column != dims[2] &&
177+
vertex_raster[se_idx] != 0 &&
178+
weight_raster[se_idx] != no_data_val
179+
if condition(condition_raster[source_idx],
180+
condition_raster[se_idx])
181+
res = diag_avg(weight_raster[source_idx],
182+
weight_raster[se_idx])
183+
push!(sources, vertex_raster[source_idx])
184+
push!(destinations, vertex_raster[se_idx])
185+
push!(node_weights, res)
186+
end
187+
end
188+
end
189+
190+
## If it's a directed graph we also need to account for west,
191+
## north, southwest, and northwest destinations
192+
if directed
193+
# Cardinal neighbor destination indices for DiGraph
194+
west_idx = CartesianIndex((row, column - 1))
195+
north_idx = CartesianIndex((row - 1, column))
196+
197+
# "West"
198+
if column != dims[2] && vertex_raster[west_idx] != 0 &&
199+
weight_raster[west_idx] != no_data_val
200+
if condition(condition_raster[source_idx],
201+
condition_raster[west_idx])
202+
res = card_avg(weight_raster[source_idx],
203+
weight_raster[west_idx])
204+
push!(sources, vertex_raster[source_idx])
205+
push!(destinations, vertex_raster[west_idx])
206+
push!(node_weights, res)
207+
end
208+
end
209+
210+
# "North"
211+
if column != dims[2] && vertex_raster[north_idx] != 0 &&
212+
weight_raster[north_idx] != no_data_val
213+
if condition(condition_raster[source_idx],
214+
condition_raster[north_idx])
215+
res = card_avg(weight_raster[source_idx],
216+
weight_raster[north_idx])
217+
push!(sources, vertex_raster[source_idx])
218+
push!(destinations, vertex_raster[north_idx])
219+
push!(node_weights, res)
220+
end
221+
end
222+
223+
## Diagonal neighbors for directed graphs
224+
if !cardinal_neighbors_only
225+
# Diagonal neighbor destination indices for DiGraph
226+
nw_idx = CartesianIndex((row - 1, column - 1))
227+
sw_idx = CartesianIndex((row + 1, column - 1))
228+
229+
# "Northwest"
230+
if column != dims[2] && row != 1 &&
231+
vertex_raster[nw_idx] != 0 &&
232+
weight_raster[nw_idx] != no_data_val
233+
if condition(condition_raster[source_idx],
234+
condition_raster[nw_idx])
235+
res = diag_avg(weight_raster[source_idx],
236+
weight_raster[nw_idx])
237+
push!(sources, vertex_raster[source_idx])
238+
push!(destinations, vertex_raster[nw_idx])
239+
push!(node_weights, res)
240+
end
241+
end
242+
243+
# "Southwest"
244+
if row != dims[1] && column != dims[2] &&
245+
vertex_raster[sw_idx] != 0 &&
246+
weight_raster[sw_idx] != no_data_val
247+
if condition(condition_raster[source_idx],
248+
condition_raster[sw_idx])
249+
res = diag_avg(weight_raster[source_idx],
250+
weight_raster[sw_idx])
251+
push!(sources, vertex_raster[source_idx])
252+
push!(destinations, vertex_raster[sw_idx])
253+
push!(node_weights, res)
254+
end
255+
end
256+
end
257+
end
258+
end
259+
end
260+
end
261+
262+
# push!'ing to vectors then combining is way faster than add_edge!
263+
# Sometimes, depending on the nodemap, there may be multiple edges with
264+
# different weights connecting the same two vertices (only when vertices
265+
# span multiple "pixels"). Default to using the min of all duplicates via
266+
# combine = min
267+
if directed
268+
g = SimpleWeightedDiGraph(
269+
sources,
270+
destinations,
271+
node_weights,
272+
combine = min
273+
)
274+
else
275+
g = SimpleWeightedGraph(
276+
sources,
277+
destinations,
278+
node_weights,
279+
combine = min
280+
)
281+
end
282+
283+
return g
284+
end

‎src/utils.jl

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,13 @@
1+
2+
"""
3+
is_data(a::Number, b::Number)
4+
5+
This function is the default used for deciding whether to connect two neighbors
6+
in a GeoArray when constructing a graph. It always returns `true`, so all
7+
neighbors will be connected
8+
9+
Returns `true`
10+
"""
11+
function is_data(a::Number, b::Number)
12+
return true
13+
end

0 commit comments

Comments
 (0)