-
Notifications
You must be signed in to change notification settings - Fork 2
/
nsink_generate_flowpath.R
249 lines (232 loc) · 11.1 KB
/
nsink_generate_flowpath.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
#' Generate and clean a flowpath for N-Sink
#'
#' This function takes an XY location as a starting point and generates a
#' flowpath for use in the N-Sink nitrogen removal analysis. The flowpath is a
#' combination of a flow direction derived flowpath on land plus NHDPlus derived
#' stream-reach flowpath.
#'
#' @param starting_location An \code{\link{sf}} point location as a starting
#' point for the flowpath. Projection must match
#' projection in input_data.
#' @param input_data A list of input data with (at least) "fdr", "streams",
#' "tot", and "raster_template". These may be generated with
#' \code{\link{nsink_prep_data}}.
#' @return An \code{\link{sf}} LINESTRING object of the flowpath that starts at
#' the \code{starting_location} and ends at the ouflow of the HUC.
#' @export
#' @import sf
#' @examples
#' \dontrun{
#' library(nsink)
#' niantic_huc <- nsink_get_huc_id("Niantic River")$huc_12
#' niantic_data <- nsink_get_data(niantic_huc, data_dir = "nsink_data")
#' aea <- 5072
#' niantic_nsink_data <- nsink_prep_data(niantic_huc, projection = aea,
#' data_dir = "nsink_niantic_data")
#' pt <- c(1948121, 2295822)
#' start_loc <- st_sf(st_sfc(st_point(c(pt)), crs = aea))
#' fp <- nsink_generate_flowpath(start_loc, niantic_nsink_data)
#' }
nsink_generate_flowpath <- function(starting_location, input_data){
if(st_crs(starting_location) != st_crs(input_data$streams)){
stop(paste0("The coordinate reference systems for your starting location and the input data do not match. Re-project to a common reference system."))
}
starting_location <- st_transform(starting_location, crs = st_crs(input_data$fdr))
# flowPath code borrowed directly from raster, modified to work with terra
fp <- flowPath(input_data$fdr, st_coordinates(starting_location))
if(is.null(fp)){
warning("Flowpath returned NULL.")
return(list(flowpath_ends = NULL, flowpath_network = NULL))
} else {
fp <- terra::xyFromCell(input_data$fdr, fp)
}
# Fixes cases with a single point flowpath: rare but annoying
if(nrow(fp) == 1){
dist <- units::set_units(1, "m")
dist <- units::set_units(dist, st_crs(input_data$streams,
parameters = TRUE)$ud_unit,
mode = "standard")
fp <- rbind(fp, fp + as.numeric(dist))
}
fp <- st_sfc(st_linestring(fp), crs = st_crs(input_data$fdr))
fp <- st_transform(fp, st_crs(input_data$streams))
fp_ends <- nsink_get_flowpath_ends(fp, input_data$streams, input_data$tot)
# This is for cases where flowpath doesn't intersect existing flowlines
dist <- units::set_units(0.001, "m")
dist <- units::set_units(dist, st_crs(input_data$streams,
parameters = TRUE)$ud_unit,
mode = "standard")
if(nrow(st_as_sf(fp_ends)) >= 2){
if(length(unlist(st_is_within_distance(fp_ends, input_data$streams,
as.numeric(dist)))) > 0 &
any(fp_ends[1,] != fp_ends[2,])){
fp_flowlines <- nsink_get_flowline(fp_ends, input_data$streams, input_data$tot)
} else {
fp_flowlines <- NULL
}
} else {
fp_flowlines <- NULL
}
fp_ends <- st_sfc(st_geometry(fp_ends))
fp_ends <- st_sf(fp_ends, crs = st_crs(input_data$streams))
if(!is.null(fp_flowlines)){
fp_flowlines <- nsink_split_flowline(fp_ends, fp_flowlines)
}
list(flowpath_ends = fp_ends, flowpath_network = fp_flowlines)
}
#' Get flowpath beginning and ends
#'
#' Flowpath from land is only portion that needs to be generated from the flow
#' direction grid. This function extracts those portions of the generated
#' flowpath. There may be issues with off network waterbodies... May need to
#' find all sections without connected flowpaths...
#'
#' @param flowpath An \code{sf} LINESTRING of the flowpath, generated with \code{\link{nsink_generate_flowpath}}
#' @param streams NHDPlus streams from \code{\link{nsink_prep_data}}
#' @param tot NHDPlus time of travel from \code{\link{nsink_prep_data}} which
#' provides the from and to nodes.
#' @return An \code{sf} object of the portions of the flowpath that are not
#' represented by the NHDPlus flowlines
#' @import sf
#' @importFrom dplyr filter
#' @keywords internal
nsink_get_flowpath_ends <- function(flowpath, streams, tot){
#drops streams without traced network from and to nodes
#Removal from these is dealt with via off_network removal
if(nrow(streams) > 0){
streams <- suppressMessages(left_join(streams, tot))
streams <- filter(streams, !is.na(.data$fromnode))
streams <- filter(streams, !is.na(.data$tonode))
streams <- st_difference(st_combine(streams), st_combine(flowpath))
splits <- lwgeom::st_split(flowpath, st_combine(streams))
splits <- suppressWarnings(st_collection_extract(splits, "LINESTRING"))
splits <- nsink_generate_from_to_nodes(
st_sf(splits, data = data.frame(id = seq_along(splits))))
first <- which(!splits$fromnode %in% splits$tonode)
last <- which(!splits$tonode %in% splits$fromnode)
ends <- splits[c(first, last),]
} else {ends <- flowpath}
ends
}
#' Get flowlines that intersect with a flowpath
#'
#' Extract flowlines that intersect with flowpath ends. This uses the actual
#' flowlines as a part of the flowpath instead of simply using the raster
#' derived flowpaths which do not follow the flowlines exactly.
#'
#'
#' @param flowpath_ends An \code{sf} LINESTRING of the flowpath ends, generated
#' with \code{\link{nsink_get_flowpath_ends}}
#' @param streams NHDPlus streams from \code{\link{nsink_prep_data}}
#' @param tot NHDPlus time of travel from \code{\link{nsink_prep_data}} which
#' provides the from and to nodes.
#' @return An \code{sf} object of the NHDPlus flowlines that occur after a
#' raster flowpath intersects the stream network.
#' @import sf dplyr
#' @importFrom utils tail
#' @importFrom igraph graph_from_data_frame shortest_paths edge_attr
#' @importFrom rlang .data
#' @keywords internal
nsink_get_flowline <- function(flowpath_ends, streams, tot){
#filtering out streams without network from and to nodes
#These are dealt with via off_network removal
streams_tot <- suppressMessages(left_join(streams, tot))
streams_tot <- filter(streams_tot, !is.na(.data$fromnode))
streams_tot <- filter(streams_tot, !is.na(.data$tonode))
streams_df <- select(streams_tot, "fromnode", "tonode", "stream_comid")
st_geometry(streams_df) <- NULL
streams_df <- mutate_all(streams_df, as.character)
streams_g <- graph_from_data_frame(streams_df, directed = TRUE)
dist <- units::set_units(0.001, "m")
dist <- units::set_units(dist, st_crs(streams, parameters = TRUE)$ud_unit,
mode = "standard")
from_nd_idx <- unlist(st_is_within_distance(flowpath_ends[1,], streams_tot, dist))
to_nd_idx <- unlist(st_is_within_distance(flowpath_ends[2,], streams_tot, dist))
from_nd <- streams_df[from_nd_idx,]$fromnode
to_nd <- streams_df[to_nd_idx,]$tonode
#to_nd <- filter(streams_df, !tonode %in% fromnode)
#to_nd <- unique(pull(to_nd, tonode))
idx <- shortest_paths(streams_g, from_nd, to_nd, output = "epath",
mode = "out")$epath[[1]]
if(length(idx) == 0){
idx <- shortest_paths(streams_g, from_nd, output = "epath",
mode = "out")$epath
idx_lng <- unlist(lapply(idx, length))
idx_idx <- which(idx_lng == max(idx_lng))
idx <- shortest_paths(streams_g, from_nd, output = "epath",
mode = "out")$epath[[idx_idx]]
}
fl_comids <- edge_attr(streams_g, "stream_comid", idx)
st_agr(flowpath_ends) <- "constant"
fp_end_pt <- tail(st_cast(flowpath_ends[1,], "POINT"), 1)
tol1 <- units::set_units(1, "m")
tol1 <- units::set_units(tol1, st_crs(streams, parameters = TRUE)$ud_unit,
mode = "standard")
tol01 <- units::set_units(0.1, "m")
tol01 <- units::set_units(tol01, st_crs(streams, parameters = TRUE)$ud_unit,
mode = "standard")
fp_flowlines <- slice(streams_tot, match(fl_comids, streams_tot$stream_comid))
fp_flowlines <- st_snap(fp_flowlines, fp_end_pt, tolerance = tol1)
fp_flowlines <- lwgeom::st_split(fp_flowlines, st_combine(fp_end_pt))
fp_flowlines <- suppressWarnings(st_collection_extract(fp_flowlines,
"LINESTRING"))
#fp_flowlines1 <- filter(fp_flowlines, !st_overlaps(st_snap(fp_flowlines,
# flowpath_ends[1,],
# tol01),
# flowpath_ends[1,],
# sparse = FALSE))
fp_flowlines
}
#' Split flowlines where they intersect with a flowpath
#'
#' Takes a flowpath input
#'
#'
#' @param flowpath_ends The ends of the flowpath that are not a part of the
#' network
#' @param flowpath_network The flowpath network
#' @return An \code{sf} object of the NHDPlus flowlines split where the fp ends
#' intersect with the flowlines.
#' @import sf dplyr
#' @importFrom rlang .data
#' @importFrom lwgeom st_split
#' @importFrom sf st_snap st_crs st_collection_extract st_intersects
#' @keywords internal
nsink_split_flowline <- function(flowpath_ends, flowpath_network){
tol1 <- units::set_units(1, "m")
tol1 <- units::set_units(tol1, st_crs(flowpath_network, parameters = TRUE)$ud_unit,
mode = "standard")
splits <- st_collection_extract(lwgeom::st_split(flowpath_network,
flowpath_ends[1,]),
"LINESTRING")
if(nrow(splits) == nrow(flowpath_network)){
flowpath_network <- st_snap(flowpath_network, flowpath_ends[1,],
tolerance = 0)
splits <- st_collection_extract(lwgeom::st_split(flowpath_network,
flowpath_ends[1,]),
"LINESTRING")
}
if(nrow(flowpath_network) == 1){
return(splits)
}
splits <- mutate(splits, split_id = seq_along(.data$stream_comid))
split_reach_comid <- splits$stream_comid[duplicated(splits$stream_comid)]
if(length(split_reach_comid) > 0){
split_reach <- filter(splits, .data$stream_comid == split_reach_comid)
split_reach_tonode <- unique(split_reach$tonode)
next_reach <- filter(flowpath_network, .data$fromnode == split_reach_tonode)
split_reach <- st_snap(split_reach, flowpath_ends, tolerance = tol1)
split_idx <- !st_intersects(split_reach, next_reach, sparse = FALSE)
} else {
split_reach <- splits[st_is_within_distance(flowpath_ends[1,], splits, tol1,
sparse = FALSE),]
next_reach <- filter(split_reach, .data$tonode %in% .data$fromnode)
split_idx <- next_reach$split_id
}
if(nrow(next_reach) == 0){
return(splits)
}
split_ditch <- split_reach$split_id[split_idx]
splits <- filter(splits, .data$split_id != split_ditch)
splits
}