-
Notifications
You must be signed in to change notification settings - Fork 5
/
compute_weights.R
264 lines (236 loc) · 11 KB
/
compute_weights.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
#' Computes weights \code{HSweights} from \code{HS} and a river network
#'
#' The function computes weights either based on catchment areas (polygons),
#' or line segments which intersect a polygon with runoff information. Weights
#' are assigned so that the sum of weights for area or line features within a
#' polygon equal to 1.
#'
#' The river network needs to be a "clean", connected network. This means that
#' connected river segments share a node at the point of confluence, and all
#' linestrings are cut at segment intersections.
#'
#' Weights should be one of the following: "equal", "length", "strahler", "area",
#' or a numeric vector which specifies the weight for each river segment.
#' \itemize{
#' \item \code{equal} option assigns equal weights to all river segments within
#' a polygon.
#' \item \code{length} option weights river segments within a polygon based on
#' the length of the segment.
#' \item \code{strahler} option weights river segments based on the Strahler
#' number computed for the
#' supplied river network.
#' \item A numeric vector with length equal to the number of features in
#' \code{river}. Weights will be computed within a polygon: \eqn{x_i/sum(x)}
#' \item \code{area} option weights the river segments based on segment-specific
#' catchment area falling inside a polygon. See more details below.
#' }
#'
#' If line-based weights (equal, length, strahler, user specified vector) are
#' used, the river network is split at grid cell boundaries before determining
#' polygon-segment relationship.
#'
#' "area" weights can be used when catchment-based weighting is desired. If no
#' further data is supplied, the function computes a Voronoi diagram from the
#' river network (see \code{\link{river_voronoi}} for details). If the basins
#' are known, they can be supplied which allows skipping the delineation step
#' entirely.
#'
#'
#' @param HS A 'HS' object, obtained with \code{\link{raster_to_HS}}.
#' @param river An 'sf' linestring feature representing a river network.
#' @param basins An 'sf' polygon object. If weights are set to "area", providing
#' basins skips the delineation process. ID column must have the name as
#' \code{riverID} See details. Optional. Defaults to \code{NULL}.
#' @param dasymetric Column name in \code{river} or in \code{basins} to be used
#' as the ancillary information in dasymetric mapping. If \code{NULL}
#' (default), no dasymetric mapping is performed. See details.
#' @param pycnophylactic Column in \code{HS} to be used as basis in
#' pycnophylactic interpolation. If \code{NULL}, or if \code{is.null(basins)},
#' not performed. See details.
#' @param n Number of iterations when using pycnophylactic interpolation.
#' Default \code{25}.
#' @param intensive Whether the pycnophylactic variable is intensive (density,
#' like runoff in mm), or not (in which case it is extensive, or counts like
#' runoff in volume).
#' @param weights Name of a column in \code{river} to be used directly as
#' weights. Defaults to \code{NULL}. See Details.
#' @param aoi An area of interest ('sf' polygon object) used to intersect
#' \code{river}. Ignored, if basins provided (in which case, aoi is the union
#' of \code{basins}). Defaults to \code{NULL}.
#' @param riverID A character string which specifies the name of the column in
#' \code{river} containing unique river network identifiers. Defaults to
#' \code{"riverID"}.
#' @param verbose Whether or not print progress indicators.
#'
#'
#' @return Returns a list object with class 'HSweights', containing the following
#' elements:
#' \itemize{
#' \item \code{river}. The supplied river network with routing information.
#' See \code{\link{river_network}} for details.
#' \item \code{weights}. River network lines or catchment polygon network
#' which was used as the basis of weighting. See
#' \code{\link{compute_area_weights}} or \code{\link{compute_river_weights}}
#' for details.
#' \item \code{HS}. The input HS object containing runoff information.
#' }
#'
compute_HSweights <- function(HS,
river,
basins = NULL,
dasymetric = NULL,
pycnophylactic = NULL,
n = 20,
intensive = TRUE,
weights = NULL,
aoi=NULL,
riverID = "riverID",
verbose = FALSE) {
##############
# CHECK INPUTS
##############
test <- inherits(HS, "HS")
if(!test) {
stop("HS input should be of class 'HS', obtained with function
raster_to_HS() or create_HS()")
}
test <- inherits(river, "sf")
if(!test) {
stop("river input should be of class 'sf'")
}
### determine what to do
test <- is.null(basins)
if(test) {
type <- sf::st_geometry_type(river, by_geometry = FALSE)
if(type %in% c("LINESTRING", "MULTILINESTRING")) {
track <- "line"
if(verbose) message(paste0("No basins provided - downscaling using ",
"river lines"))
} else if(type %in% c("POLYGON", "MULTIPOLYGON")) {
track <- "area"
if(verbose) message(paste0("Basins provided - interpolating ",
"area -> area"))
basins <- river
} else {
stop("Geometry of input 'river' is not LINESTRING, MULTILINESTRING",
" POLYGON or MULTIPOLYGON. Consider using function ",
"sf::st_collection_extract() to extract the desired type.")
}
# track <- "line"
# if(verbose) message(paste0("No basins provided - downscaling using ",
# "river lines"))
} else {
track <- "area"
if(verbose) message(paste0("Basins provided - interpolating ",
"area -> area"))
}
#############
# AREA TRACK
#############
if (track == "area") {
# check column names
test <- hasName(river, riverID) && hasName(basins, riverID)
if(!test) stop("Both river and basins input must have identically ",
"named column",riverID,"specified by parameter riverID")
test <- riverID == "riverID"
if(!test) {
ind <- which(names(river) == riverID)
river <- tibble::add_column(river,
riverID = dplyr::pull(river, ind),
.before=1)
ind <- which(names(basins) == riverID)
basins <- tibble::add_column(basins,
riverID = dplyr::pull(basins, ind),
.before=1)
}
if(!is.null(weights)) {
test <- hasName(basins, weights)
if(!test) stop("No column ", weights, " found in basins input.")
}
# 1. Intersect river network with union of basins, and intersect HS
# with the union of basins - this is to make sure that the runoff
# units is consistent with the basins. If they are not, weights
# will not equal to 1 for every runoff unit.
select <- sf::st_intersects(river,
sf::st_geometry(sf::st_union(basins)),
sparse=FALSE)
river <- river[select,]
# 2. compute basin weights
if(!is.null(weights)) {
basins <- compute_area_weights(basins,
HS,
pycno = NULL,
dasy = NULL,
weights = weights)
} else if(!is.null(pycnophylactic)) {
basins <- compute_area_weights(basins,
HS,
pycno = pycnophylactic,
dasy = dasymetric,
intensive = intensive,
weights = NULL)
} else if(!is.null(dasymetric)) {
basins <- compute_area_weights(basins,
HS,
pycno = NULL,
dasy = dasymetric,
weights = NULL)
} else {
basins <- compute_area_weights(basins,
HS,
pycno = NULL,
dasy = NULL,
weights = NULL)
}
# create output
HSweights <- create_HSweights(target = river,
weights = basins,
source = HS)
}
############
# LINE TRACK
############
if (track == "line") {
# check column names
test <- riverID == "riverID"
if(!test) {
ind <- which(names(river) == riverID)
river <- tibble::add_column(river,
riverID = dplyr::pull(river, ind),
.before=1)
}
if(!is.null(weights)) {
test <- hasName(river, weights)
if(!test) stop("No column ", weights, " found in river input.")
}
# 1. crop river network to aoi if given, else intersect it with grid
if(!is.null(aoi)) {
select <- sf::st_intersects(river,
sf::st_geometry(sf::st_union(aoi)),
sparse=FALSE)
river <- river[select,]
} else {
select <- sf::st_intersects(river,
sf::st_geometry(sf::st_union(HS)),
sparse=FALSE)
river <- river[select,]
}
# 2. compute weights based on river lines
if(!is.null(dasymetric)) {
splitriver <- compute_river_weights(river,
HS,
seg_weights = dasymetric,
split=TRUE)
} else {
splitriver <- compute_river_weights(river,
HS,
seg_weights = NULL,
split=TRUE)
}
# create output
HSweights <- create_HSweights(target = river,
weights = splitriver,
source = HS)
}
return(HSweights)
}