Skip to content
Snippets Groups Projects
Commit 3c7b612c authored by Bock, Andy's avatar Bock, Andy
Browse files

Non Dend Updates

parent 1d1c1003
No related branches found
No related tags found
1 merge request!161NonDend updates
Hide whitespace changes
Inline Side-by-side
workspace/07-2_NonDend.Rmd
+ 48
− 77
View file @ 3c7b612c
......@@ -175,7 +175,9 @@ if(needs_layer(ND_gpkg, divides_nd)){
# DEBUG:
# HRU layer
if(debug) {
protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
protoHRUs <- divides_lu %>%
group_by("POI_ID") %>%
summarize(do_union = T) %>%
sfheaders::sf_remove_holes(.) %>%
st_make_valid()
write_sf(protoHRUs, ND_gpkg, HRU_layer)
......@@ -213,7 +215,9 @@ if(!"HUC_12" %in% unique(divides_lu$aggStep)){
if(debug) {
# DEBUG:
# HRU layer
protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
protoHRUs <- divides_lu %>%
group_by(POI_ID) %>%
summarize(do_union = T) %>%
sfheaders::sf_remove_holes(.) %>%
st_make_valid()
write_sf(protoHRUs, ND_gpkg, HRU_layer)
......@@ -247,7 +251,9 @@ if("Coastline" %in% unique(full_nhd$FTYPE)){
if(debug) {
# DEBUG:
# HRU layer
protoHRUs <- union_polygons_geos(divides_lu, "POI_ID") %>%
protoHRUs <- divides_lu %>%
group_by(POI_ID) %>%
summarize(do_union = T) %>%
sfheaders::sf_remove_holes(.) %>%
st_make_valid()
# Write out updates
......@@ -277,89 +283,54 @@ This is where alot of the more complex steps happen. For most VPUs, there are no
if(needs_layer(ND_gpkg, "missing_cats")){
print(paste0("Currently there are ", sum(is.na(divides_lu$POI_ID)), " cats without a POI_ID"))
# Handle nhd sinks
# arrange divides_lu by member_COMID and populate with RowID for easy iteration within next two steps
divides_lu <- divides_lu %>%
arrange(member_COMID, POI_ID) %>%
mutate(row_ID = row_number()) %>%
arrange(row_ID) %>%
NHD_sinks(area_thresh = min_da_km/2,
if(sum(is.na(divides_lu$POI_ID)) > 0){
# Handle nhd sinks
# arrange divides_lu by member_COMID and populate with RowID for easy iteration within next two steps
divides_lu <- divides_lu %>%
arrange(member_COMID, POI_ID) %>%
mutate(row_ID = row_number()) %>%
arrange(row_ID)
HUC_sinks <- NHD_sinks(divides_lu, area_thresh = min_da_km/2,
HUC12_table =file.path(data_paths$nhdplus_dir, "HUC12.rds"),
NHD_sinks = read_sf(data_paths$nhdplus_gdb, "Sink"))
# identify terminal outlets and upstream areas not yet integrated
if(needs_layer(ND_gpkg, missing_terms)){
# Read in HUCs
vpu_WBD <- vpu_WBD %>%
st_transform(crs)
# join to HUC12s information
term_outlets_wPOI <- outlets_close(full_nhd, nhd, divides_lu) %>%
left_join(st_intersection(., vpu_WBD) %>%
st_drop_geometry(.) %>%
dplyr::select(outlet_COMID, HUC_12), by = "outlet_COMID") %>%
mutate(outlet_HUC12 = HUC_12,
outlet_COMID = as(outlet_COMID,
class(full_nhd$COMID))) %>%
dplyr::select(-HUC_12)
# get unassigned catchments that hold terminal outlets and their contributing
# catchments upstream
term_outlets_wPOI <- dplyr::select(term_outlets_wPOI, -c(POI_ID, HUC_12_int.y)) %>%
left_join(miss_term_assign(
term_outlets = term_outlets_wPOI,
divides_poi = divides_lu,
nhd = full_nhd, elev = elev),
by = "outlet_COMID")
write_sf(term_outlets_wPOI, ND_gpkg, missing_terms)
} else {
term_outlets_wPOI <- read_sf(ND_gpkg, missing_terms)
}
# Assign cats POI_ID based on geometric/topological relationships
# derive downstream or aggregated catchment to assign to
term_outlets_sol <- assignable_cats(
list(divides_lu, full_nhd, term_outlets_wPOI, xwalk_divides_wbd)
)
# write out missing cats
if (nrow(filter(term_outlets_sol$divides_lu, is.na(POI_ID))) > 0){
# get missing cats
# Get unassigned catchments with no associated terminal outlet/flowline
# Bring over to divides
divides_lu <-
left_join(term_outlets_sol$divides_lu,
miss_term_assign(
term_outlets = filter(term_outlets_sol$divides_lu,
is.na(POI_ID)) %>%
rename(outlet_COMID = member_COMID),
divides_poi = term_outlets_sol$divides_lu,
nhd = full_nhd, elev = elev) %>%
select(outlet_COMID, new_POI_ID = POI_ID),
by = c("member_COMID" = "outlet_COMID")) %>%
mutate(POI_ID = ifelse(!is.na(new_POI_ID), new_POI_ID, POI_ID),
aggStep = ifelse(!is.na(new_POI_ID), "boundary DEM", aggStep)) %>%
select(-new_POI_ID)
write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats)
if(length(HUC_sinks) == 2){
divides_lu <- HUC_sinks$divides_poi
sinks_table <- HUC_sinks$sink_cats_table
}
if(debug)
write_sf(term_outlets_sol$term_outlets, ND_gpkg, "unassigned_term")
if(sum(is.na(divides_lu$POI_ID)) > 0) {
divides_dem <- miss_term_assign(divides_lu, full_nhd, elev)
divides_lu <- divides_lu %>%
left_join(select(divides_dem, member_COMID, agg_ID),
by = "member_COMID") %>%
mutate(POI_ID = ifelse(!is.na(agg_ID), agg_ID, POI_ID),
aggStep = ifelse(!is.na(agg_ID), "boundary DEM", aggStep)) %>%
select(-agg_ID)
if(exists("sinks_table")){
sinks_table_fin <- filter(sinks_table, !member_COMID %in% divides_dem$member_COMID) %>%
rbind(divides_dem)
} else {
sinks_table_fin <- divides_dem
}
if(sum(is.na(divides_lu$POI_ID)) > 0){
write_sf(filter(divides_lu, is.na(POI_ID)), ND_gpkg, missing_cats)
}
write_sf(sinks_table_fin, ND_gpkg, ND_table)
}
} else {
print ("all unaggregated catchments assigned")
divides_lu <- term_outlets_sol$divides_lu
}
# Same issue as AK, for some MP HRUs, the wrong cat part is kept
all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>%
union_polygons_geos(., "POI_ID") %>%
sfheaders::sf_remove_holes(.) %>%
st_make_valid()
# Prob HRU - filter(all_hrus, POI_ID == 140402000209)
all_hrus <- filter(divides_lu, !is.na(POI_ID)) %>%
group_by(POI_ID) %>%
......
workspace/R/non_dend.R
+ 32
− 12
View file @ 3c7b612c
......@@ -453,6 +453,10 @@ outlets_close <- function(nhd_full, nhd_sub, divides_poi){
arrange(COMID) %>%
mutate(COMID = as.character(COMID))
if(nrow(little_outlets_fin) == 0){
return(NULL)
}
# some VPUs may not have these outlets
if (nrow(little_outlets_fin) == 0){
missing_outlets <- NA
......@@ -707,7 +711,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){
if(nrow(filter(sink_cats, is.na(POI_ID))) > 0){
# Aggretate sink cats to HUC12 and add area
sink_cats_agg <- filter(sink_cats, is.na(POI_ID)) %>%
group_by("HUC_12_int") %>%
group_by(HUC_12_int) %>%
summarize(do_union = T) %>%
mutate(agg_id = row_number())
......@@ -732,7 +736,8 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){
agg_cats_index <- sapply(cats_int, function(s) if (length(s) == 0) NA else s)
# Populate ID of HRUs where applicable
sink_cats$POI_ID_new <- agg_cats_large[agg_cats_index,]$HUC_12
sink_cats$POI_ID_new <- agg_cats_large[agg_cats_index,]$HUC_12_int
sink_cats <- mutate(sink_cats, aggStep = ifelse(!is.na(POI_ID_new), "HUC12_sinks, large", aggStep))
sink_cats_final <- dplyr::select(st_drop_geometry(sink_cats), member_COMID, POI_ID_new) %>%
filter(!is.na(POI_ID_new))
......@@ -781,11 +786,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){
group_by(agg_id) %>%
filter(n_distinct(bound_cat_POI_ID) == 1) %>%
ungroup()
sink_cats_assign <- dplyr::select(sink_cats_small, -POI_ID_new) %>%
inner_join(dplyr::select(small_sinks_POI, agg_id, POI_ID_new = bound_cat_POI_ID), by = "agg_id") %>%
dplyr::select(member_COMID, POI_ID_new)
sink_cats_assign <- dplyr::select(sink_cats_nhd_table, member_COMID, agg_id) %>%
inner_join(dplyr::select(small_sinks_POI, agg_id, POI_ID_new = bound_cat_POI_ID), by = "agg_id") %>%
distinct()
......@@ -799,8 +800,20 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){
}
}
large_sinks <- select(sink_cats, member_COMID, HUC_12 = HUC_12_int, agg_ID = POI_ID_new, aggStep) %>%
st_drop_geometry()
small_sinks <- select(sink_cats_nhd_table, member_COMID, HUC_12 = HUC_12_int, neighbor_COMID = bound_cat,
neighbor_HUC12 = bound_cat_HUC12, agg_ID = bound_cat_POI_ID, agg_id, aggStep) %>%
filter(member_COMID %in% sink_cats_small$member_COMID, !member_COMID %in% sink_cats_final$member_COMID) %>%
select(-agg_id) %>%
filter(HUC_12 == neighbor_HUC12) %>%
mutate(aggStep = "HUC12_sinks, small") %>%
st_drop_geometry()
sinks_table <- data.table::rbindlist(list(large_sinks, small_sinks), fill = TRUE)
}
return(divides_poi)
return(list(divides_poi = divides_poi, sink_cats_table = sinks_table))
}
#' Find and lump remaining areas where possible
......@@ -813,7 +826,7 @@ NHD_sinks <- function(divides_poi, area_thresh, HUC12_table, NHD_sinks){
#' @return (sf data frame) data frame of likely POI_ID assignments of reamining sinks
#' Keep the below link here until publishd
#' Topo data bins: https://www.sciencebase.gov/catalog/item/5f5154ba82ce4c3d12386a02
miss_term_assign <- function(term_outlets, divides_poi, nhd, elev){
miss_term_assign <- function(divides_poi, nhd, elev){
# comid <- single comid of a missing catchment (member_COMID in divides_lu)
#' Find and lump remaining areas where possible
......@@ -885,18 +898,25 @@ miss_term_assign <- function(term_outlets, divides_poi, nhd, elev){
return(buf_df2)
}
}
term_outlets <- filter(divides_poi, is.na(POI_ID) | aggStep == "HUC12_sinks, small") %>%
dplyr::rename(outlet_COMID = member_COMID)
out_df <- do.call(rbind,
pbapply::pblapply(unique(term_outlets$outlet_COMID),
assign_func, divides_poi, nhd, elev, cl = NULL))
assign_func, divides_poi, nhd, elev, cl = NULL)) %>%
rename(member_COMID = outlet_COMID)
# out_df <- lapply(unique(term_outlets$outlet_COMID),
# assign_func, divides_poi, nhd, elev)
out_df <- out_df %>%
left_join(dplyr::select(st_drop_geometry(divides_poi), HUC12_neighbor = HUC_12_int, member_COMID),
left_join(dplyr::select(st_drop_geometry(divides_poi), HUC_12 = HUC_12_int, member_COMID),
by = "member_COMID") %>%
left_join(dplyr::select(st_drop_geometry(divides_poi), neighbor_HUC12 = HUC_12_int, member_COMID),
by = c("neighbor_COMID" = "member_COMID")) %>%
dplyr::select(outlet_COMID, Elev, neighbor_COMID, HUC12_neighbor, POI_ID)
dplyr::select(member_COMID, neighbor_COMID, HUC_12, neighbor_HUC12, agg_ID = POI_ID) %>%
mutate(aggStep = "boundary DEM")
return(out_df)
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment