-
Bock, Andy authoredBock, Andy authored
utils.R 15.82 KiB
hyfabric_source <- list.files(pattern = "hyfabric_.*.tar.gz")
hyfabric_version <- substr(hyfabric_source, 10, 14)
if(!(require(hyfabric, quietly = TRUE) && packageVersion("hyfabric") == hyfabric_version)) {
dir.create(path = Sys.getenv("R_LIBS_USER"), showWarnings = FALSE, recursive = TRUE)
.libPaths(Sys.getenv("R_LIBS_USER"))
install.packages(hyfabric_source, lib = Sys.getenv("R_LIBS_USER"),
repos = NULL, type = "source")
}
fix_headwaters <- function(nhd_rds, out_rds, new_atts = NULL,
nhdpath = nhdplusTools::nhdplus_path(),
par = 6) {
if(!file.exists(out_rds)) {
nhd <- readRDS(nhd_rds)
# Need to replace headwater flowlines with the trimmed burn line event.
# Some headwater flowlines are not fully within the catchment they go with.
ble <- sf::read_sf(nhdpath, "BurnLineEvent")
ble <- dplyr::left_join(dplyr::select(sf::st_drop_geometry(nhd), COMID),
ble, by = "COMID")
ble <- sf::st_zm(sf::st_as_sf(ble))
fline <- sf::st_zm(nhd)
sf::st_geometry(nhd)[nhd$StartFlag == 1] <- sf::st_geometry(ble)[nhd$StartFlag == 1]
if(!is.null(new_atts)) {
new_atts <- data.table::fread(new_atts,
integer64 = "character",
showProgress = FALSE)
if("weight" %in% names(new_atts)) new_atts <- rename(new_atts,
ArbolateSu = weight)
nhd <- select(nhd, COMID, VPUID, RPUID, FTYPE, FromNode, ToNode,
StartFlag, StreamOrde, StreamCalc, Divergence, WBAREACOMI,
VPUID, RPUID, DnMinorHyd) %>%
right_join(new_atts, by = c("COMID" = "comid")) %>%
nhdplusTools::align_nhdplus_names()
}
nhd <- sf::st_sf(nhd)
m_to_km <- (1/1000)
nhd$LENGTHKM <- as.numeric(sf::st_length(sf::st_transform(nhd, 5070))) * m_to_km
# is a headwater and for sure flows to something.
check <- !nhd$COMID %in% nhd$toCOMID & !(nhd$toCOMID == 0 | is.na(nhd$toCOMID) | !nhd$toCOMID %in% nhd$COMID)
check_direction <- dplyr::filter(nhd, check)$COMID
cl <- parallel::makeCluster(16)
# check and fix these.
new_geom <- pbapply::pblapply(check_direction, fix_flowdir,
network = nhd[nhd$COMID %in% check_direction, ], cl = cl)
parallel::stopCluster(cl)
# One is empty. Deal with it and remove from the replacement.
error_index <- sapply(new_geom, inherits, what = "try-error")
errors <- filter(nhd, COMID %in% check_direction[error_index])
check[which(nhd$COMID %in% errors$COMID)] <- FALSE
if(!all(sapply(sf::st_geometry(errors), sf::st_is_empty))) { stop("Errors other than empty geometry from fix flowdir")
}
new_geom <- new_geom[!error_index]
new_geom <- do.call(c, new_geom)
st_geometry(nhd)[check] <- new_geom
saveRDS(nhd, out_rds)
### ONLY RUN IF CHANGED ###
# sbtools::authenticate_sb()
# sbtools::item_replace_files("5dcd5f96e4b069579760aedb", "data/NHDPlusNationalData/nhdplus_flowline_update.rds")
}
return(out_rds)
}
#' Retrieves files from SB in facet file_structure
#' @param sbitem (character) ScienceBase item ID
#' @param out_dir (directory) directory where files are downloaded
#'
# Retrieve files from SB if files not listed under sbtools
retrieve_sb <- function(sbitem, out_dir){
item <- item_get(sbitem)
linked_files <-
lapply(item$facets, function(x)
{
list(name = x$name,
files = lapply(x$files, function(y, out_path)
{
out_file <- file.path(out_path, y$name)
download.file(y$downloadUri, out_file)
list(fname = y$name,
url = y$downloadUri)
},
out_path = file.path(data_dir, x$name)))
})
}
#' Documents why streamgage is excluded from POIs
#' @param source (character) data source
#' @param reason (character) reason gage is excluded
#' @return writes output to file
#'
gage_document <- function(data, source, drop_reason){
if(!file.exists("cache/dropped_gages.csv")){
write.csv(data %>% mutate(data_source = source, reason = drop_reason),
"cache/dropped_gages.csv", quote = F, row.names = F)
} else {
write.table(data %>%
mutate(data_source = source, reason = drop_reason),
file = "cache/dropped_gages.csv", append = T, sep = ",",
col.names = F, row.names = F, quote = F)
}
}
#' clean network
#' @param net data.frame nhdplus network
#' @param net data.frame new nhdplus network as from e2nhd
#' @param nwm data.frame with national water model new network
#' @return data.frame fully reconciled attributes
clean_net <- function(net, net_new, nwm) {
# NOTE: need to disconnected diverted paths as indicated in the tonode / fromnode attributes.
diverted_heads <- filter(net_new, divergence == 2)
net_new <- filter(net_new, !tocomid %in% diverted_heads$comid)
# NOTE: avoid fcode == 56600 (coastal)
avoid <- net$comid[net$fcode == 56600]
# NOTE: avoid large groups where things have been redirected in an un-natural way.
groups <- group_by(nwm, tocomid)
groups <- data.frame(tocomid = summarise(groups, .groups = "drop"), size = group_size(groups))
groups <- filter(groups, size > 10 & tocomid != 0)
# NOTE: also avoid all flowlines that go to coastal flowlines.
# This is due to unwanted great lakes connections in the NWM.
avoid <- c(avoid,
net_new$comid[net_new$tocomid %in% avoid],
net_new$comid[net_new$tocomid %in% groups$tocomid],
nwm$comid[nwm$tocomid %in% groups$tocomid])
avoid <- unique(avoid)
nwm <- filter(nwm, tocomid %in% c(net_new$comid, 0))
net_new <- left_join(net_new, select(nwm, comid,
nwm_tocomid = tocomid),
by = "comid") %>%
mutate(tocomid = ifelse(is.na(tocomid), 0, tocomid))
## NWM already uses 0 for outlets.
# NOTE: just being explicit about this for clarity.
candidate <- select(net_new, comid, tocomid, nwm_tocomid) %>%
filter(!comid %in% avoid)
mismatch <- candidate[candidate$tocomid != candidate$nwm_tocomid, ]
net_new <- left_join(net_new, select(mismatch, comid,
new_tocomid = nwm_tocomid),
by = "comid")
to_change <- filter(net_new, !is.na(new_tocomid))
net_new <- net_new %>%
mutate(divergence = ifelse(comid %in% to_change$new_tocomid & divergence == 2, 1,
ifelse(comid %in% to_change$tocomid & divergence == 1, 2,
divergence)),
tocomid = ifelse(comid %in% to_change$comid, new_tocomid, tocomid)) %>%
select(-nwm_tocomid, -new_tocomid)
net_new <- filter(net_new, fcode != 56600)
net_new <- mutate(net_new, tocomid = ifelse(!tocomid %in% comid, 0, tocomid))
return(net_new)
}
#' Merges geopackages together to create CONUs geopackage of features
#' @param feat (character) Type of feature to merge (POIs, segments) character
#' @param in_gpkg (character) geopackage to merge (GF, ND_, etc.)
#' @param out_gpkg (character) Name of output geopackage
#' @return output geopackage of CONUS for given features
#'
Merge_VPU <- function(feat, in_gpkg, out_gpkg){
if(needs_layer(out_gpkg, feat)) {
all_sf <- paste0(feat, "_CONUS")
#VPUs <- c("03N", "03S", "03W")
#VPUs <- c("10L", "10U") VPUs <-c(paste0("0", c(1:2)), "03N", "03S", "03W", paste0("0", c(4:9)),
as.character(11:18), "10U", "10L")
featCON <- do.call("rbind", lapply(VPUs, function(x) {
tryCatch({
layer <- ifelse(feat %in% c("nhd_flowline", "unassigned_gages", "unassigned_TE"),
feat, paste0(feat, "_", x))
read_sf(file.path("cache", paste0(in_gpkg, x, ".gpkg")),
layer)},
error = function(e) stop(paste("didn't find", x,
"\n Original error was", e)))
}))
write_sf(featCON, out_gpkg, feat)
}
}
#' Merges geopackages together to create CONUs geopackage of features
#' @param VPU VPU from NHDPlusV2
#' @param full_cats (sf data.frame) all catchments (sinks/flowline) within given VPU
#' @param divides (sf data.frame) divides layer for a given VPU
#' @param wbd (data.frame) HUC12 X-walk table, see get_data for full citation
#' @param wbd_SF (sf data.frame) HUC12 layer from NHDPlusv2.1 national geodatabase
#' @return (sf data.frame) intersection of nhd catchments and divides layer with HUC12
#'
get_HUC12_Xwalk <- function(VPU, full_cats = FALSE, divides = FALSE, wbd = FALSE, wbd_SF = FALSE){
# HUC02 includes some
if(VPU == "02"){
grep_exp <-"^02|^04"
} else if (VPU == "08") {
grep_exp <- "^03|^08"
} else {
grep_exp <- paste0("^", substr(VPU, start = 1, stop = 2))
}
# crosswalk to HU12, filter by VPU
nhd_to_HU12 <- read.csv(wbd, colClasses = c("character", "integer", "character")) %>%
filter(grepl(grep_exp, .data$HUC_12)) %>%
mutate(HUC_8 = substr(HUC_12, 1, 8))
# Read in the NHDPlus National GDB HUC12 layer
HUC12_lyr <- readRDS(file.path(wbd_SF)) %>%
#filter(grepl(paste0("^", VPU,".*"), .data$HUC_12)) %>%
filter(grepl(grep_exp, .data$HUC_12)) %>% #,
st_set_precision(10) %>%
st_transform(st_crs(full_cats)) %>%
st_make_valid() %>%
select(HUC_12, ACRES, HU_10_TYPE, HU_12_DS, HU_12_MOD, HU_12_TYPE)
reg_cats <- full_cats %>%
left_join(select(nhd_to_HU12, c(FEATUREID, HUC_12)), by = "FEATUREID")
# Intersect the HUC12s and catchments
cats_HUC12 <- sf::st_intersection(reg_cats, HUC12_lyr) %>%
# Convert area to squaker kilometers
mutate(intArea = as.numeric(st_area(.)) * 0.000001,
AreaHUC12 = ACRES * .00404686) %>%
rename(HUC_12_int = HUC_12.1)
divides_HUC12 <- sf::st_intersection(divides, HUC12_lyr) %>%
# Convert area to squaker kilometers
mutate(intArea = as.numeric(st_area(.)) * 0.000001,
AreaHUC12 = ACRES * .00404686) %>%
rename(HUC_12_int = HUC_12.1)
return(list(cats_HUC12 = cats_HUC12, divides_HUC12 = divides_HUC12))
}
#' Fix flow direction
#' @description If flowlines aren't digitized in the expected direction, #' this will reorder the nodes so they are.
#' @param x The COMID of the flowline to check
#' @param r The entire network to check from. Requires a "toCOMID" field.
#' @return a geometry for the feature that has been reversed if needed.
fix_flowdir <- function(comid, network) {
try({
library(sf)
f <- network[network$COMID == comid, ]
if(is.na(f$toCOMID)) {
check_line <- network[network$toCOMID == f$COMID, ][1, ]
check_position <- "start"
} else {
check_line <- network[network$COMID == f$toCOMID, ][1, ]
check_position <- "end"
}
suppressMessages(check_end <- sf::st_join(nhdplusTools::get_node(f, position = check_position),
dplyr::select(check_line, check_COMID = COMID)))
reverse <- is.na(check_end$check_COMID)
if(reverse) {
sf::st_geometry(f)[reverse] <- sf::st_reverse(sf::st_geometry(f)[reverse])
}
return(sf::st_geometry(f))
})
}
#' Capture all catchments (flowline, sink, etc.) in a given RPU
#' @param fcats (sf data.frame) NHDPlusv2.1 National geodatabase catchment layer
#' @param the_RPU (character) RPU to retrive full catchment set.
#' @param fl_rds character path to flowline rds file
#' @param nhd_gdb character path to NHD geodatabase
#' @return (sf data.frame) with FEATUREID, RPUID, FTYPE, TERMINAFL
cat_rpu <- function(fcats, nhd_gdb){
fl <- read_sf(nhd_gdb, "NHDFlowline_Network") %>%
align_nhdplus_names()
# read the CatchmentSP
cats <- readRDS(fcats) %>%
st_as_sf() %>%
st_drop_geometry() %>%
dplyr::select(FEATUREID, SOURCEFC)
# read the Flowlines
flowln_df <- fl %>%
st_as_sf() %>%
st_drop_geometry() %>%
select(COMID, FTYPE, TerminalFl, RPUID)%>%
rename(FEATUREID = COMID) %>%
rename(TERMINALFL = TerminalFl)%>%
filter(FEATUREID %in% cats$FEATUREID)
# Read in sinks
sink_df <- sf::st_read(nhd_gdb, layer = "Sink") %>%
st_drop_geometry() %>%
dplyr::select(SINKID, RPUID = InRPU) %>%
dplyr::rename(FEATUREID = SINKID) %>%
dplyr::mutate(FTYPE = "Sink", TERMINALFL = 0) %>% dplyr::filter(!FEATUREID %in% flowln_df$FEATUREID)
#combine all the RPUIDs
lkup_rpu <- rbind(flowln_df, sink_df)
# FEATUREID 10957920, 20131674, 24534636 - this are the ids of the missing, checked on map,
# they look like they should not be used
# add the records for the missing
missrec_df <- data.frame(FEATUREID=c(10957920, 20131674, 245346360),
RPUID = c("03a", "03d", "17b"),
TERMINALFL = c(0, 0, 0)) %>%
mutate(FTYPE = "")
lkup_rpu2 <- rbind(lkup_rpu, missrec_df)
return(lkup_rpu2)
}
### DEPRECATED moved to nhdplusTools as subset_vpu with different methods###
#' Subsets VPU for terminal outlets
#' @param nhd_rds path to staged nhd flowline rds file
#' @param VPU VPU from NHDPlusV2
#' @return sf data.frame containing subset VPU
#'
VPU_Subset <- function(nhd_rds, VPU){
# Read in COMIDs of outlets
regoutlets <- list(`05` = 1844789L,
`06` = 1861888L,
`07` = 5093446L,
`10U` = 17244390L,
`10L` = 6018266L,
`14` = c(20734041L, 18267749L),
`11` = c(14320629L, 941140164L, 15334434L))
nhd <- readRDS(nhd_rds)
if (!VPU %in% names(regoutlets[!names(regoutlets) %in% c("10L", "10U")])){
# For regions that terminate to the NHDPlus domain
# 1,2,3,4,8,9,12,15,17,18
# Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
nhd <- nhd %>%
filter(grepl(paste0("^", VPU, ".*"), .data$VPUID)) %>%
filter(FTYPE != "Coastline")
keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE, skip_toCOMID = TRUE)
nhd <- filter(nhd, COMID %in% keep$COMID)
} else if (VPU %in% c("05", "06", "07", "10U", "14")){
# Subset VPUs that are single outlet
outlet <- regoutlets[VPU]
nhd <- filter(nhd, COMID %in% unlist(get_UT(nhd, outlet)))
keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE, skip_toCOMID = TRUE)
nhd <- filter(nhd, COMID %in% keep$COMID, VPUID == VPU)
} else {
outlets <- unlist(regoutlets[VPU])
for (out in outlets){
nhd_US <- get_UT(nhd, out)
keep <- prepare_nhdplus(filter(nhd, COMID %in% unlist(nhd_US)), 0, 0, 0, FALSE,
skip_toCOMID = TRUE)
ifelse( exists("final") ,
final <- rbind(final, keep),
final <- keep)
}
nhd <- filter(nhd, COMID %in% final$COMID)
}
nhd <- select(nhd, -VPUID, -FromNode, -ToNode,
-StartFlag, -StreamOrde, -StreamCalc)
return(st_as_sf(nhd))
}
### DEPRECATED moved to nhdplusTools as is ####
# make sf1 compatible with sf2
st_compatibalize <- function(sf1, sf2) {
sf1 <- st_transform(sf1, st_crs(sf2))
g <- attr(sf1, "sf_column")
gp <- attr(sf2, "sf_column")
names(sf1)[names(sf1) == g] <- gp
attr(sf1, "sf_column") <- gp
sf1
}
### DEPRECATED moved to nhdplusTools as is ####
rename_geometry <- function(g, name){
current = attr(g, "sf_column")
names(g)[names(g)==current] = name
st_geometry(g) <- name
g
}
#' @param df data.frame
#' @param rp integer rows per chunk
#' @param nc integer number of rows
split_df <- function(df, rp = NULL, nc = NULL) {
nr <- nrow(df)
if(!is.null(rp)) {
n <- rp
} else if(!is.null(nc)) {
n <- ceiling(nr / nc)
} else {
stop("must provide rp or nc")
}
split(df,
rep(1:ceiling(nr / n),
length.out = nr,
each = n))
}