Newer
Older
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
#' create refactor input
#' @param flowline flowlines to be used in refactor
#' @param refactor_event_table table of locations refactor should split flowlines
#' @param pois all POIs with geometry
#' @param pois_data attributes for POIs
#' @return list with all events, outlets and comids to exclude
create_refactor_input <- function(flowline, refactor_event_table, pois, pois_data) {
events_refactor <- filter(refactor_event_table, hy_id %in% flowline$COMID) %>%
rename(COMID = hy_id) %>%
distinct(COMID, REACHCODE, REACH_meas,
poi_id, geom) %>%
arrange(poi_id) %>%
mutate(event_identifier = row_number())
POIs_ref <- pois |>
inner_join(select(st_drop_geometry(flowline), TotDASqKM, COMID, DnHydroseq),
by = c("hy_id" = "COMID"))
# Also need to avoid modification to flowlines immediately downstream of POIs
# This can cause some hydrologically-incorrect catchment aggregation
POI_downstream <- filter(flowline, Hydroseq %in% POIs_ref$DnHydroseq, AreaSqKM > 0)
# build final outlets set
term_poi <- filter(pois_data, hl_reference == "Type_Term")
outlets <- POIs_ref |>
mutate(type = ifelse(poi_id %in% term_poi$poi_id, "terminal", "outlet"))
# Need to avoid refactoring catchments that are long and thing and
# reasonably large. They cause problems in a downstream step.
avoid <- dplyr::filter(flowline, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
exclude <- c(outlets$COMID,
avoid$COMID,
POI_downstream$COMID)
outlets <- rename(outlets, COMID = hy_id) %>%
filter(!poi_id %in% events_refactor$poi_id)
list(events_refactor = events_refactor,
outlets = outlets,
exclude = exclude)
}
#' run refactor nhdplus
#' @param flowline flowlines to be used in refactor
#' @param split_meters distance to use to split flowlines
#' @param para_split_flines how many jobs to run at a time for splitting
#' @param combine_meters how long to use as combine target
#' @param events_refactor event locations where flowlines should be split
#' @param refactor_exclusion_list flowlines that should not be refactored
#' @param rpu_code raster processing unit
#' @param proj_crs projected coordinate reference system to put outputs in
run_refactor_nhdplus <- function(flowline, split_meters,
para_split_flines, combine_meters,
events_refactor, refactor_exclusion_list,
rpu_code, proj_crs) {
tf <- paste0("temp/refactored_", rpu_code,".gpkg")
tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
unlink(list(tf, tr))
refactor_nhdplus(nhdplus_flines = dplyr::select(flowline, -FTYPE), # Select controls whether prepare_nhdplus is used.
split_flines_meters = split_meters,
split_flines_cores = para_split_flines,
collapse_flines_meters = combine_meters,
collapse_flines_main_meters = combine_meters,
out_refactored = tf,
out_reconciled = tr,
three_pass = TRUE,
purge_non_dendritic = FALSE,
events = events_refactor,
exclude_cats = refactor_exclusion_list,
warn = TRUE)
out <- list(refactored = st_transform(read_sf(tf), proj_crs),
reconciled = st_transform(read_sf(tr), proj_crs))
unlink(list(tf, tr))
return(out)
}
#' create refactor lookup
#' @param refactored_network reconciled and refactored network as ouput from run_refactor_nhdplus
#' @param flowline flowlines as used in refactor
#' @param events_refactor event locations where flowlines should be split
#' @param outlets all outlets
#' @return list with lookup table, outlets, and duplicate checks
create_refactor_lookup <- function(refactored_network, flowline, events_refactor, outlets) {
# create lookup for ref flowlines to use in the non-dendritic steps
refactor_lookup <- st_drop_geometry(refactored_network$reconciled) %>%
dplyr::select(ID, member_COMID) %>%
dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
tidyr::unnest(cols = member_COMID) %>%
dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
dplyr::rename(reconciled_ID = ID)
if(is.character(refactor_lookup$reconciled_ID))
refactor_lookup$reconciled_ID <- as.integer(refactor_lookup$reconciled_ID)
lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID")
# readr::write_csv(lookup_table, lookup_table_file)
# write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
# Join refactored to original NHD
refactored <- refactored_network$refactored %>%
select(member_COMID = COMID, Hydroseq, event_identifier, event_REACHCODE) %>%
inner_join(select(st_drop_geometry(flowline), orig_COMID = COMID, Hydroseq), by = "Hydroseq")
if(nrow(events_refactor) > 0){
# Subset for events
refactored_events <- refactored %>%
filter(!is.na(event_REACHCODE), !is.na(event_identifier)) %>%
mutate(event_identifier = as.numeric(event_identifier))
event_outlets <- events_refactor %>%
inner_join(st_drop_geometry(refactored_events), by = "event_identifier") %>%
select(hy_id = COMID, event_id = event_identifier, poi_id, member_COMID)%>%
inner_join(select(lookup_table, -NHDPlusV2_COMID), by = "member_COMID")
# subset for refactored outlets (non-events)
refactored_outlets <- filter(refactored,
!member_COMID %in% event_outlets$member_COMID)
# get ref_COMId for other outlets
outlets_ref <- outlets %>%
filter(!poi_id %in% event_outlets$poi_id) %>%
inner_join(lookup_table,
by = c("COMID" = "NHDPlusV2_COMID")) %>%
group_by(COMID) %>%
filter(member_COMID == max(member_COMID)) %>%
select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, type)
outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets),
fill = TRUE) %>%
st_as_sf()
} else {
# get ref_COMId for other outlets
outlets_ref_COMID <- outlets %>%
inner_join(lookup_table,
by = c("COMID" = "NHDPlusV2_COMID")) %>%
group_by(COMID) %>%
filter(member_COMID == max(member_COMID)) %>%
select(hy_id = COMID, poi_id, member_COMID, reconciled_ID, type)
}
# write_sf(outlets_ref_COMID, out_refac_gpkg, outlets_layer)
check_dups_poi <- outlets_ref_COMID %>%
group_by(reconciled_ID) %>%
filter(n() > 1) %>%
ungroup()
# if(nrow(check_dups_poi) > 1){
# print("Double-check for double POIs")
# write_sf(check_dups_poi, out_refac_gpkg, paste0(dup_pois, "_", rpu_code))
# } else {
# print("no double POIs detected")
# }
out <- list(lookup = lookup_table, outlets = outlets_ref_COMID, dups = check_dups_poi)
}
reconcile_divides <- function(cats, refactored_network, fdr_fac, para_reconcile, cache_split) {
divides <- reconcile_catchment_divides(catchment = cats,
fline_ref = refactored_network$refactored,
fline_rec = refactored_network$reconciled,
fdr = fdr_fac$fdr,
fac = fdr_fac$fac,
para = para_reconcile,
cache = cache_split,
keep = NULL)
load(cache_split)
split_cats <- bind_rows(split_cats[lengths(split_cats) > 0])
split_cats$areasqkm <- as.numeric(units::set_units(sf::st_area(split_cats), "km^2"))
split_cats <- sf::st_transform(split_cats, proj_crs)
list(divides = divides, split_cats = split_cats)