diff --git a/.Rbuildignore b/.Rbuildignore index 3a5914c..9b36f49 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -14,3 +14,4 @@ data-raw ^\.ccache$ ^\.github$ ^tic\.R$ +Dockerfile \ No newline at end of file diff --git a/.github/workflows/R-CMD-check.yml b/.github/workflows/R-CMD-check.yml index 172f0eb..cb09333 100644 --- a/.github/workflows/R-CMD-check.yml +++ b/.github/workflows/R-CMD-check.yml @@ -1,16 +1,16 @@ on: push: branches: - - master + - main pull_request: branches: - - master + - main name: R-CMD-check jobs: R-CMD-check: - runs-on: ubuntu-18.04 + runs-on: ubuntu-20.04 env: R_REMOTES_NO_ERRORS_FROM_WARNINGS: true GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} diff --git a/DESCRIPTION b/DESCRIPTION index fe32f83..5f12aa8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,27 +9,24 @@ URL: https://github.com/dblodgett-usgs/hyrefactor BugReports: https://github.com/dblodgett-usgs/hyrefactor/issues Depends: R (>= 3.5.0) -Imports: - dplyr, - sf, - lwgeom, - units, - magrittr, - data.table, - raster, - nhdplusTools, - igraph, - tidyr, - pbapply, - xml2, - rvest, - httr, - stringr, - methods, - rlang, - rgeos, - rmapshaper -Suggests: testthat, knitr, rmarkdown, mapview, rgdal, snow +Imports: + dplyr, + sf, + lwgeom, + units, + data.table, + terra, + nhdplusTools, + tidyr, + pbapply, + rvest, + httr, + methods, + rlang, + rgeos, + rmapshaper, + utils +Suggests: testthat, knitr, rmarkdown Remotes: dblodgett-usgs/nhdplusTools License: CC0 Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index a2ea50b..63503dc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -9,6 +9,7 @@ export(collapse_flowlines) export(download_elev) export(download_fdr_fac) export(flowpaths_to_linestrings) +export(get_minimal_network) export(map_outlet_ids) export(reconcile_catchment_divides) export(reconcile_collapsed_flowlines) @@ -18,6 +19,7 @@ export(split_flowlines) export(union_linestrings_geos) export(union_polygons_geos) importFrom(data.table,rbindlist) +importFrom(dplyr,`%>%`) importFrom(dplyr,arrange) importFrom(dplyr,bind_rows) importFrom(dplyr,case_when) @@ -31,6 +33,7 @@ importFrom(dplyr,mutate) importFrom(dplyr,n) importFrom(dplyr,rename) importFrom(dplyr,right_join) +importFrom(dplyr,row_number) importFrom(dplyr,select) importFrom(dplyr,slice_max) importFrom(dplyr,slice_min) @@ -40,28 +43,15 @@ importFrom(dplyr,ungroup) importFrom(httr,RETRY) importFrom(httr,progress) importFrom(httr,write_disk) -importFrom(igraph,V) -importFrom(igraph,bfs) -importFrom(igraph,graph_from_data_frame) -importFrom(igraph,head_of) -importFrom(igraph,incident_edges) -importFrom(igraph,shortest_paths) -importFrom(igraph,topo_sort) -importFrom(magrittr,"%>%") -importFrom(methods,is) importFrom(methods,slot) +importFrom(nhdplusTools,add_plus_network_attributes) importFrom(nhdplusTools,get_node) importFrom(nhdplusTools,get_vaa) importFrom(nhdplusTools,prepare_nhdplus) importFrom(nhdplusTools,rename_geometry) -importFrom(raster,as.matrix) -importFrom(raster,cellFromXY) -importFrom(raster,crop) -importFrom(raster,crs) -importFrom(raster,mask) -importFrom(raster,raster) -importFrom(raster,rasterToPolygons) -importFrom(raster,rowColFromCell) +importFrom(parallel,makeCluster) +importFrom(parallel,stopCluster) +importFrom(pbapply,pblapply) importFrom(rgeos,gLineMerge) importFrom(rgeos,gUnaryUnion) importFrom(rlang,":=") @@ -70,7 +60,8 @@ importFrom(rmapshaper,ms_explode) importFrom(rmapshaper,ms_simplify) importFrom(rvest,html_attr) importFrom(rvest,html_nodes) -importFrom(sf,"st_crs<-") +importFrom(rvest,read_html) +importFrom(sf,`st_crs<-`) importFrom(sf,as_Spatial) importFrom(sf,st_area) importFrom(sf,st_as_sf) @@ -91,10 +82,12 @@ importFrom(sf,st_geometry_type) importFrom(sf,st_intersection) importFrom(sf,st_intersects) importFrom(sf,st_is_empty) +importFrom(sf,st_is_longlat) importFrom(sf,st_length) importFrom(sf,st_line_merge) importFrom(sf,st_linestring) importFrom(sf,st_make_valid) +importFrom(sf,st_multipolygon) importFrom(sf,st_nearest_feature) importFrom(sf,st_point) importFrom(sf,st_precision) @@ -107,8 +100,23 @@ importFrom(sf,st_touches) importFrom(sf,st_transform) importFrom(sf,st_union) importFrom(sf,st_within) -importFrom(tidyr,unnest) -importFrom(tidyr,unnest_longer) +importFrom(terra,as.matrix) +importFrom(terra,as.polygons) +importFrom(terra,cellFromRowCol) +importFrom(terra,cellFromXY) +importFrom(terra,crop) +importFrom(terra,crs) +importFrom(terra,mask) +importFrom(terra,rast) +importFrom(terra,res) +importFrom(terra,rowColFromCell) +importFrom(terra,setValues) +importFrom(terra,vect) +importFrom(terra,xFromCol) +importFrom(terra,xyFromCell) +importFrom(terra,yFromRow) +importFrom(tidyr,separate_rows) importFrom(units,drop_units) importFrom(units,set_units) -importFrom(xml2,read_html) +importFrom(utils,capture.output) +importFrom(utils,str) diff --git a/R/aggregate_catchments.R b/R/aggregate_catchments.R index 4198d14..4fe1a17 100644 --- a/R/aggregate_catchments.R +++ b/R/aggregate_catchments.R @@ -24,7 +24,9 @@ #' outlets <- data.frame(ID = c(31, 3, 5, 1, 45, 92), #' type = c("outlet", "outlet", "outlet", "terminal", "outlet", "outlet"), #' stringsAsFactors = FALSE) -#' aggregated <- aggregate_catchments(walker_fline_rec, walker_catchment_rec, outlets) +#' aggregated <- aggregate_catchments(flowpath = walker_fline_rec, +#' divide = walker_catchment_rec, +#' outlets = outlets) #' plot(aggregated$cat_sets$geom, lwd = 3, border = "red") #' plot(walker_catchment_rec$geom, lwd = 1.5, border = "green", col = NA, add = TRUE) #' plot(walker_catchment$geom, lwd = 1, add = TRUE) diff --git a/R/aggregate_network.R b/R/aggregate_network.R index a1950fe..638320a 100644 --- a/R/aggregate_network.R +++ b/R/aggregate_network.R @@ -39,10 +39,8 @@ #' nexus and that all catchments are well-connected. #' #' @export -#' @importFrom igraph graph_from_data_frame topo_sort incident_edges V bfs head_of shortest_paths #' @importFrom sf st_is_empty st_drop_geometry #' @importFrom dplyr filter mutate left_join select distinct case_when bind_rows -#' @importFrom tidyr unnest_longer #' @examples #' source(system.file("extdata", "walker_data.R", package = "nhdplusTools")) #' @@ -58,6 +56,8 @@ #' #' aggregated <- aggregate_network(fline, outlets) #' +#' aggregated <- aggregate_network(fline, outlets) +#' #' outlets <- dplyr::filter(fline, ID %in% outlets$ID) #' #' outlets <- nhdplusTools::get_node(outlets) @@ -66,185 +66,30 @@ #' plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE) #' plot(outlets$geometry, add = TRUE) #' +#' aggregate_network <- function(flowpath, outlets, - da_thresh = NA, only_larger = FALSE, - post_mortem_file = NA) { - - flowpath <- check_names(flowpath, "aggregate_network") - - geo <- inherits(flowpath, "sf") - - if(geo){flowpath_geo = flowpath; flowpath = sf::st_drop_geometry(flowpath)} - - if (any(!outlets$ID %in% flowpath$ID)) stop("Outlet IDs must all be in flowpaths.") - - flowpath$toID[flowpath$toID == 0] <- NA - - term <- flowpath$toID[flowpath$ID %in% outlets[outlets$type == "terminal", ]$ID] - - if (any(!is.na(term))) { - if(!is.na(post_mortem_file)) save(list = ls(), file = post_mortem_file) - stop("Terminal paths must have an NA or 0 toID") - } - - lps <- get_lps(flowpath) - - outlets <- make_outlets_valid(outlets, flowpath, lps,da_thresh = da_thresh, only_larger = only_larger) %>% - distinct() - - lps <- get_lps(flowpath) - - outlets <- mutate(outlets, ID = paste0("cat-", .data$ID)) - - # Build catchment and nexus data.frames to preserve sanity in graph traversal. - catchment <- flowpath %>% - mutate(toID = ifelse(is.na(.data$toID), -.data$ID, .data$toID)) - - # Join id to toID use ID as from nexus ID since we are assuming dendritic. - nexus <- left_join(select(drop_geometry(catchment), toID = .data$ID), - select(drop_geometry(catchment), fromID = .data$ID, .data$toID), - by = "toID") %>% - mutate(nexID = paste0("nex-", .data$toID), - fromID = paste0("cat-", .data$fromID), - toID = paste0("cat-", .data$toID)) %>% - select(.data$nexID, .data$fromID, .data$toID) - - # get fromID and toID straight with "nex-" prefix - catchment <- mutate(drop_geometry(catchment), - fromID = paste0("nex-", .data$ID), - toID = paste0("nex-", .data$toID), - ID = paste0("cat-", .data$ID)) %>% - select(.data$fromID, .data$toID, cat_ID = .data$ID) - - # Convert the catchment network to a directed graph - cat_graph <- graph_from_data_frame(d = catchment, - directed = TRUE) - - outlets <- outlets %>% - left_join(select(catchment, nexID_stem = .data$fromID, ID = .data$cat_ID), by = "ID") %>% - left_join(select(catchment, nexID_terminal = .data$toID, ID = .data$cat_ID), by = "ID") %>% - mutate( - nexID = case_when( - type == "outlet" ~ .data$nexID_stem, - type == "terminal" ~ .data$nexID_terminal - ) - ) %>% - select(.data$ID, .data$type, .data$nexID) %>% - distinct() - - # sorted version of the graph to re-order outlets in upstream-downstream order. - cat_graph_sort_verts <- topo_sort(cat_graph) - outlet_verts <- cat_graph_sort_verts[names(cat_graph_sort_verts) %in% outlets$nexID] - outlets <- outlets[match(names(outlet_verts), outlets$nexID), ] - - cat_sets <- data.frame(ID = outlets$ID, - nexID = outlets$nexID, - set = I(rep(list(list()), nrow(outlets)))) - - fline_sets <- data.frame(ID = outlets$nexID, - set = I(rep(list(list()), nrow(outlets)))) - - verts <- V(cat_graph) - - us_verts <- c() - - for (cat in seq_len(nrow(cat_sets))) { - - if(cat %% 10 == 0) message(paste(cat, "of", nrow(cat_sets))) - - outlet <- filter(outlets, .data$ID == cat_sets$ID[cat]) - - ut <- bfs(graph = cat_graph, - root = cat_sets$nexID[cat], - neimode = "in", - order = TRUE, - unreachable = FALSE, - restricted = verts) - - outlet_id <- as.integer(gsub("^cat-", "", outlets$ID[cat])) - head_id <- lps$head_ID[which(lps$ID == outlet_id)] - head_id <- head_id[head_id != outlet_id] - - if (length(head_id) == 0) head_id <- outlet_id # Then the outlet is a headwater. - - abort_code <- tryCatch( - { - um <- find_um(us_verts, cat_graph, - cat_id = cat_sets$nexID[cat], - head_id = head_id, - outlet_type = filter(outlets, .data$nexID == cat_sets$nexID[cat])$type) - abort_code <- "" - }, error = function(e) e - ) - - if(abort_code != "") { - if(!is.na(post_mortem_file)) { - save(list = ls(), file = post_mortem_file) - stop(paste("Upstream Main error, post mortem file:", post_mortem_file)) - } else { - stop(paste("Upstream Main error:", abort_code)) - } - } - - if (length(us_verts) > 0) { - vert <- us_verts[which(um %in% us_verts)] - - if (length(vert) == 1) { - # Since longest path is used in find_um if multiple paths are found - # there is a chance that multiple us verts get consumed in a single step? - us_verts <- us_verts[!us_verts %in% vert] - } - } - - um <- as.numeric(gsub("^nex-", "", um)) - - if(0 %in% um) um[um == 0] <- as.numeric(gsub("^cat-", "", outlets[outlets$nexID == "nex-0", ]$ID)) - - abort_code <- tryCatch( - { - fline_sets$set[[cat]] <- um - - # Excludes the search node to avoid grabbing the downstream catchment. - ut_verts <- ut$order[!is.na(ut$order) & names(ut$order) != cat_sets$nexID[cat]] - cat_sets$set[[cat]] <- filter(catchment, .data$fromID %in% names(ut_verts))$cat_ID - remove <- head_of(cat_graph, unlist(incident_edges(cat_graph, ut_verts, "in"))) - verts <- verts[!verts %in% remove] - - if (outlet$type == "outlet") { - cat_sets$set[[cat]] <- c(cat_sets$set[[cat]], outlet$ID) - verts <- verts[!names(verts) == outlet$nexID] - } - - abort_code <- "" - }, error = function(e) e - ) - - if(abort_code != "") { - if(!is.na(post_mortem_file)) { - save(list = ls(), file = post_mortem_file) - stop(paste("error getting geometry type or with union, post mortem file:", post_mortem_file)) - } else { - stop(paste("error getting geometry type or with union for line merge:", abort_code)) - } - } - - cat_sets$set[[cat]] <- as.numeric(gsub("^cat-", "", cat_sets$set[[cat]])) - - # us_verts are where we need to stop while stepping upstream. - us_verts <- c(us_verts, cat_sets$nexID[cat]) - } - - cat_sets <- select(cat_sets, -.data$nexID) - - cat_sets[["ID"]] <- as.numeric(gsub("^cat-", "", outlets$ID)) - - fline_sets[["ID"]] <- as.numeric(gsub("^cat-", "", outlets$ID)) + da_thresh = NA, only_larger = FALSE, + post_mortem_file = NA) { + + flowpath <- validate_flowpath(flowpath, outlets, post_mortem_file) + + # makes sure outlets connect + outlets <- make_outlets_valid(outlets, drop_geometry(flowpath), + da_thresh = da_thresh, + only_larger = only_larger) %>% + sort_outlets(flowpath) + + fline_sets <- get_catchment_sets(flowpath, outlets) + cat_sets <- fline_sets[[2]] + + fline_sets <- fline_sets[[1]] + # create long form ID to set member list - sets <- tidyr::unnest_longer(drop_geometry(fline_sets), col = "set") - + sets <- unnest_flines(fline_sets) + # Figure out what the ID of the downstream catchment is. - next_id <- left_join(sets, + next_id <- left_join(sets, select(drop_geometry(flowpath), .data$ID, .data$toID), by = c("set" = "ID")) %>% # first find the ID downstream of the outlet of each catchment. @@ -257,35 +102,54 @@ aggregate_network <- function(flowpath, outlets, left_join(select(sets, set_toID = .data$ID, .data$set), by = c("toID" = "set")) %>% select(.data$ID, toID = .data$set_toID) - - if(geo){ - fline_sets = data.frame( - setID = unlist(fline_sets$set), - ID = rep(fline_sets$ID, times = lengths(fline_sets$set))) %>% - left_join(dplyr::select(flowpath_geo, ID), by = c("setID" = "ID")) %>% - st_as_sf() %>% - filter(!sf::st_is_empty(.)) %>% - union_linestrings_geos(ID = "ID") %>% + + if(inherits(flowpath, "sf")) { + fline_sets <- + data.frame( + setID = unlist(fline_sets$set), + ID = rep(fline_sets$ID, times = lengths(fline_sets$set))) %>% + left_join(select(flowpath, .data$ID), by = c("setID" = "ID")) %>% + st_as_sf() %>% + filter(!sf::st_is_empty(.)) %>% + union_linestrings_geos(ID = "ID") %>% left_join(fline_sets, by = "ID") } - + fline_sets <- left_join(fline_sets, next_id, by = "ID") cat_sets <- left_join(cat_sets, next_id, by = "ID") - + return(list(cat_sets = cat_sets, fline_sets = fline_sets)) } +validate_flowpath <- function(flowpath, outlets, post_mortem_file) { + + flowpath <- check_names(flowpath, "aggregate_network") + + if (any(!outlets$ID %in% flowpath$ID)) stop("Outlet IDs must all be in flowpaths.") + + flowpath$toID[is.na(flowpath$toID)] <- 0 + + if (any(flowpath$toID[flowpath$ID %in% outlets[outlets$type == "terminal", ]$ID] != 0)) { + if(!is.na(post_mortem_file)) save(list = ls(), file = post_mortem_file) + stop("Terminal paths must have an NA or 0 toID") + } + + flowpath$toID <- methods::as(flowpath$toID, class(flowpath$ID)) + + return(flowpath) +} + get_lps <- function(flowpath) { flowpath <- drop_geometry(flowpath) %>% select(.data$ID, .data$LevelPathID, .data$Hydroseq) %>% group_by(.data$LevelPathID) - + headwaters <- filter(flowpath, .data$Hydroseq == max(.data$Hydroseq)) %>% ungroup() - + outlets <- filter(flowpath, .data$Hydroseq == min(.data$Hydroseq)) %>% ungroup() - + left_join(ungroup(flowpath), select(headwaters, head_ID = .data$ID, .data$LevelPathID), by = "LevelPathID") %>% left_join(select(outlets, tail_ID = .data$ID, .data$LevelPathID), @@ -302,25 +166,40 @@ get_outlets <- function(outlets, lps) { # Adds everything that contributes the same recieving catchment as a given tail id. fix_nexus <- function(flowpath, tail_id, da_thresh = NA, only_larger = FALSE) { tail <- filter(flowpath, .data$ID == tail_id) - + add <- filter(flowpath, .data$toID == tail$toID) - + if (only_larger) { add <- filter(add, .data$LevelPathID <= tail$LevelPathID) } - + if (!is.na(da_thresh)) { add <- filter(add, .data$TotDASqKM > da_thresh) } # Can add functionality here to filter which Add IDs to include. add_ids <- add$ID - + data.frame(ID = add_ids, type = rep("outlet", length(add_ids)), stringsAsFactors = FALSE) } +apply_fix_nexus <- function(bad_id, flowpath, da_thresh, only_larger) { + bind_rows( + lapply(bad_id, + function(bad_tail_id, + flowpath, + da_thresh, + only_larger) { + fix_nexus(flowpath, bad_tail_id, + da_thresh, only_larger) + }, + flowpath = flowpath, + da_thresh = da_thresh, + only_larger = only_larger)) +} + fix_tail <- function(flowpath, outlets, toid_tail_id, da_thresh = NA, only_larger = FALSE) { potential_add <- fix_nexus(flowpath, toid_tail_id, da_thresh, only_larger) new <- !potential_add$ID %in% outlets$ID @@ -332,44 +211,58 @@ fix_tail <- function(flowpath, outlets, toid_tail_id, da_thresh = NA, only_large #' catchment doesn't get orphaned in the middle of a larger catchment. #' @param outlets the outlet list of gages, etc. #' @param flowpath the reconciled flowline network -#' @param lps level paths #' @importFrom dplyr filter distinct select left_join group_by mutate #' @noRd -make_outlets_valid <- function(outlets, flowpath, lps, +make_outlets_valid <- function(outlets, flowpath, da_thresh = NA, only_larger = FALSE) { - + + outlets$ID <- methods::as(outlets$ID, class(flowpath$ID)) + + # Finds levelpaths and their unique head and outlet + lps <- get_lps(drop_geometry(flowpath)) + + # Need to check if outlets are on levelpath tails and add + # required additional outlets. + outlets <- bind_rows( + outlets, + apply_fix_nexus(filter(get_outlets(outlets, lps), + .data$ID == .data$tail_ID)$ID, + flowpath, da_thresh, only_larger) + ) + + # deduplicate outlets. outlets <- distinct(outlets) %>% group_by(.data$ID) %>% - filter(!(n() > 1 & .data$type == "outlet")) %>% + filter(!(n() > 1 & # this removes outlets that duplicate terminals. + # they can be added in the above outlets check. + .data$type == "outlet")) %>% ungroup() - + otl <- get_outlets(outlets, lps) - + count_while <- 0 - + while (!all(otl$tail_ID %in% otl$ID)) { - + bad_tail <- otl$tail_ID[which(!otl$tail_ID %in% otl$ID)] - + message(paste("Fixing", length(bad_tail), "missing outlets.")) - + outlets <- dplyr::bind_rows( outlets, - dplyr::bind_rows(lapply( - bad_tail, function(bad_tail_id, flowpath, da_thresh, only_larger) { - fix_nexus(flowpath, bad_tail_id, da_thresh, only_larger) - }, flowpath = flowpath, da_thresh = da_thresh, only_larger = only_larger)) + apply_fix_nexus(bad_tail, + flowpath, da_thresh, only_larger) ) - + otl <- get_outlets(outlets, lps) - + count_while <- count_while + 1 - + if (count_while > 20) { stop("Stuck in a while loop trying to fix disconnected outlets. Reduce drainage area threshold?") } } - + # Need to check that a "next down tributary" in the outlet set has a break along the # main stem that each outlet contributes to. otl <- left_join(otl, select(drop_geometry(flowpath), @@ -379,7 +272,7 @@ make_outlets_valid <- function(outlets, flowpath, lps, toID_tail_ID = .data$tail_ID, toID_LevelpathID = .data$LevelPathID), by = c("toID" = "ID")) - + # this grabs the most downstream if duplicates were generated. otl <- group_by(otl, .data$ID) %>% filter(.data$toID_hydroseq == min(.data$toID_hydroseq)) %>% @@ -388,52 +281,52 @@ make_outlets_valid <- function(outlets, flowpath, lps, group_by(.data$toID_tail_ID) %>% filter(n() > 1) %>% ungroup() - + # This grabs all the inflows to the nexus that each of these outlets is at. otl <- left_join(otl, select(drop_geometry(flowpath), toID_fromID = .data$ID, .data$toID), by = "toID") %>% mutate(type = "add_outlet") - + if (!is.na(da_thresh)) { otl <- left_join(otl, select(drop_geometry(flowpath), .data$ID, toID_fromID_TotDASqKM = .data$TotDASqKM), by = c("toID_fromID" = "ID")) %>% filter(.data$toID_fromID_TotDASqKM > da_thresh) } - + if (only_larger) { otl <- left_join(otl, select(drop_geometry(flowpath), .data$ID, toID_fromID_lp = .data$LevelPathID), by = c("toID_fromID" = "ID")) %>% filter(.data$toID_fromID_lp <= .data$toID_LevelpathID) } - + otl <- otl %>% select(ID = .data$toID_fromID, .data$type, LevelPathID = .data$toID_LevelpathID, .data$tail_ID) %>% distinct() - + # Need to verify that all ID == tail_ID instances are connected. # They might have been missed above. if (any(grepl("add_outlet", otl$type))) { otl$type <- "outlet" outlets <- distinct(rbind(outlets, otl[, c("ID", "type")])) } - + otl <- get_outlets(outlets, lps) tail_outlets <- which(otl$ID == otl$tail_ID & otl$type != "terminal") - + for (check in seq_along(tail_outlets)) { outlets <- rbind(outlets, fix_tail(flowpath, outlets, otl$ID[check], da_thresh = da_thresh, only_larger = only_larger)) - + connected <- FALSE while (!connected) { toid <- filter(flowpath, .data$ID == otl[["ID"]][check])$toID toid_tail_id <- filter(lps, .data$ID == toid)[["tail_ID"]] - + if (!all(toid_tail_id %in% otl$tail_ID)) { outlets <- rbind(outlets, fix_tail(flowpath, outlets, unique(toid_tail_id), @@ -444,48 +337,284 @@ make_outlets_valid <- function(outlets, flowpath, lps, otl <- get_outlets(outlets, lps) } } - + return(outlets) } -#' @description Finds upstream nearest upstream main stem shortest path to -#' either the list of verticies that might be upstream or the up main headwater. -#' @param us_verts upstream verticies that should be used as stop points. -#' @param cat_graph the catchment network to be searched. -#' @param cat_id the catchment identifier to start at. -#' @param head_id the upstream mainstem top catchment identifier. -#' @param outlet_type the type of outlet chosen from "outlet", or "terminal". -#' @noRd -find_um <- function(us_verts, cat_graph, cat_id, head_id, outlet_type) { - um <- c() +sort_outlets <- function(outlets, flowpath) { + # set toID to negative the ID for terminals that are not 0. + # This is important to make sure that all outlet nexues have a unique ID + # as opposed to 0 or na. This is meangful when sorting the graph + flowpath$toNexID <- ifelse(flowpath$toID == 0, -flowpath$ID, flowpath$toID) + + # this is a little awkward and nuanced -- read comments carefully. + outlets <- outlets %>% + # the nexID_stem is just the catchment identifier. Conceptually, + # this is the identifier for the nearest upstream nexus from an outlet. + # We have to use this identifier for an outlet along the flowpath + # of a catchment. It can be the same as the catchment because we are + # working with a dendritic assumption so each nexus has one and only + # one downstream catchment. Otherwise, we would have to treat things a + # bit differently. + mutate(nexID_stem = .data$ID) %>% + # the nexID_terminal is the to nexus of the terminal flowpath. + # We need to treat this differently because an outlet of type + # "terminal" must be downstream of all other outlets. Without a + # unique nexus ID like this, we can not be gaurunteed a sort + # order that will work correctly. + left_join(select(drop_geometry(flowpath), .data$ID, + nexID_terminal = .data$toNexID), by = "ID") %>% + # Now we can set the nexID (which is used to sort the outlets) appropriately. + mutate( + nexID = case_when( + type == "outlet" ~ .data$nexID_stem, + type == "terminal" ~ .data$nexID_terminal + ) + ) %>% + select(.data$ID, .data$type, .data$nexID) %>% + distinct() + + fp_sort <- nhdplusTools::get_sorted(drop_geometry(select(flowpath, .data$ID, .data$toNexID)))$ID + + fp_sort <- c(fp_sort, flowpath$toNexID[flowpath$toNexID < 0]) + + fp_sort <- fp_sort[fp_sort %in% outlets$nexID] + + left_join(data.frame(nexID = as.numeric(fp_sort)), outlets, by = "nexID") %>% + select(-.data$nexID) +} + +prep_flowpath <- function(flowpath) { + # flowpath <- drop_geometry(flowpath) + + # Make sure flowpath is sorted correctly + flowpath <- arrange(flowpath, desc(Hydroseq)) + + # do a little indirection with index ids + flowpath$toid <- match(flowpath$toID, flowpath$ID) + flowpath$id <- seq(1, nrow(flowpath)) + + flowpath +} + +prep_outlets <- function(outlets, flowpath) { + outlets <- left_join(outlets, select(flowpath, .data$ID, .data$id), by = "ID") + outlets$set <- 1:nrow(outlets) + outlets +} + +get_heads <- function(flowpath) { + flowpath$id[!flowpath$id %in% flowpath$toid] +} + +# need a little network walker function. +get_dwn <- function(ID, toid) { + next_dn <- toid[ID] + if(is.na(next_dn)) { + return(ID) + } else { + return(c(ID, get_dwn(next_dn, toid))) + } +} + +my_combine <- function(old, new) { + # can optimize later + if(identical(old, list())) { + new + } else { + c(old, new) + } +} + +unnest_flines <- function(x, col = "set") { - if (length(us_verts) > 0) { - suppressWarnings(paths <- shortest_paths(cat_graph, - from = cat_id, - to = us_verts, - mode = "in")) - - if (length(paths$vpath) > 1) { - path_lengths <- lengths(sapply(paths$vpath, names)) + times <- lengths(x[[col]]) + base_names <- names(x)[!names(x) == col] + + out <- as.data.frame(cbind(sapply(base_names, function(n) rep(x[[n]], times = times)))) + + names(out) <- base_names + + out[[col]] <- unlist(x[[col]]) + + out +} + +get_catchment_sets <- function(flowpath, outlets) { + + fline_sets <- data.frame(ID = outlets$ID, + set = I(rep(list(list()), nrow(outlets)))) + + cat_sets <- data.frame(ID = outlets$ID, + set = I(rep(list(list()), nrow(outlets)))) + + flowpath <- prep_flowpath(flowpath) + + outlets <- prep_outlets(outlets, flowpath) + + # We'll start at all these + heads <- get_heads(flowpath) + + # We can stop once these conditions have been met. + outlet_count <- nrow(outlets) + headwater_count <- length(heads) + + message(paste("Running", headwater_count, "headwaters for", + outlet_count, "outlets.")) + + # Set up counters for the while loop + o_c <- 1 + h_c <- 1 + + while(h_c <= headwater_count) { + head <- heads[h_c] + + path <- get_dwn(head, flowpath$toid) + + sets <- outlets$id %in% path + + breaks <- outlets[sets, ] + nr_breaks <- nrow(breaks) + + # If this path only goes to one outlet, it is within an aggregate + ## do nothing. Or maybe we can use this for aggregation later? + # If this path goes through more than one outlet, it contains flowpaths. + if(nr_breaks > 1) { + paths <- split(path, + cut(path, + breaks = c(0, breaks$id), + labels = c(breaks$set))) + + cat_sets$set[as.integer(names(paths))] <- + lapply(names(paths), function(x) { + my_combine(cat_sets$set[as.integer(x)][[1]], + flowpath$ID[paths[x][[1]]]) + }) + + path_outlets <- sapply(paths, function(x) x[length(x)]) + + flowpath$toid[path_outlets] <- NA + + # the top one isn't useful for mainstems. + paths <- paths[2:length(paths)] + + fline_sets$set[as.integer(names(paths))] <- + lapply(paths, function(x) flowpath$ID[x]) + + o_c <- o_c + (nr_breaks - 1) + } else { - path_lengths <- length(names(paths$vpath[[1]])) - } - - if (any(path_lengths > 0)) { - path_lengths[path_lengths == 0] <- NA - um <- names(paths$vpath[which(path_lengths == min(path_lengths, na.rm = TRUE))][[1]]) - um <- um[!um %in% us_verts] # Path includes the "to" but don't want to include it! + + outlet <- outlets[outlets$id == path[length(path)], ] + + # I think this can be added tothe flowpath finding for free + cat_sets$set[outlet$set][[1]] <- my_combine(cat_sets$set[outlet$set][[1]], + flowpath$ID[path]) + } + + h_c <- h_c + 1 + + if(h_c %% 1000 == 0) message(paste(h_c, "of", headwater_count)) + } - - # then look for headwater. - if (length(um) == 0) { - um <- names(shortest_paths(cat_graph, - from = cat_id, - to = paste0("nex-", head_id), - mode = "in")$vpath[[1]]) - } - - return(um) + + set_length <- lengths(fline_sets$set) + + head_outlets <- outlets[set_length == 0, ] %>% + left_join(distinct(select(flowpath, LevelPathID, ID, Hydroseq)), by = "ID") %>% + select(LevelPathID, set, head_out_Hydroseq = Hydroseq) + + head_paths <- filter(flowpath, LevelPathID %in% head_outlets$LevelPathID) %>% + left_join(head_outlets, by = "LevelPathID") %>% + filter(Hydroseq >= .data$head_out_Hydroseq) %>% + select(ID, set) %>% + group_by(set) %>% + summarise(ID = list(ID)) + + fline_sets$set[head_paths$set] <- head_paths$ID + + cat_sets$set <- lapply(cat_sets$set, function(x) unique(x)) + + list(fline_sets, cat_sets) } +#' Get Minimal Network +#' @description Given a set of outlets, will generate a minimal network by +#' calling \code{\link{aggregate_network}} and adding nhdplus attributes to the result. +#' +#' If geometry is included with the network, it will be merged and returned. +#' +#' @inheritParams aggregate_network +#' @return a data.frame (potentially including an sfc list column) with +#' attributes generated by \code{\link[nhdplusTools]{add_plus_network_attributes}} +#' and a list column "set" containing members of each output flowpath. +#' @importFrom nhdplusTools add_plus_network_attributes +#' @export +#' @examples +#' source(system.file("extdata", "walker_data.R", package = "nhdplusTools")) +#' fline <- walker_flowline +#' +#' outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329435, 5329817), +#' type = c("outlet", "outlet", "outlet", "outlet", "outlet")) +#' +#' #' Add toCOMID +#' fline <- nhdplusTools::get_tocomid(fline, add = TRUE) +#' +#' # get attributes set +#' fline <- dplyr::select(fline, ID = comid, toID = tocomid, +#' levelpathid = levelpathi, hydroseq = hydroseq, +#' areasqkm = areasqkm, lengthkm = lengthkm) +#' +#' min_net <- get_minimal_network(fline, outlets) +#' +#' plot(sf::st_geometry(fline), col = "blue") +#' plot(sf::st_geometry(min_net), lwd = 2, add = TRUE) +#' plot(sf::st_geometry(nhdplusTools::get_node(min_net)), add = TRUE) +#' +get_minimal_network <- function(flowpath, outlets) { + + flowpath_sort <- left_join( + data.frame(ID = flowpath$ID), + nhdplusTools::get_sorted(flowpath[, c("ID", "toID"), drop = TRUE], + split = TRUE), by = "ID") + + terminal_paths <- unique(flowpath_sort$terminalID[flowpath_sort$ID %in% outlets$ID]) + + # Grab terminal paths that matter and combine with outlets. + outlets <- rbind(outlets, + data.frame(ID = terminal_paths, + type = "terminal")) + + flowpath <- flowpath[flowpath_sort$terminalID %in% terminal_paths, ] + + minimal <- hyRefactor::aggregate_network( + flowpath, dplyr::filter(outlets, .data$ID %in% flowpath$ID), + da_thresh = NA, only_larger = TRUE) + + min_net <- unnest_flines(drop_geometry(minimal$fline_sets)) %>% + left_join(select(flowpath, .data$ID, .data$lengthkm, + .data$areasqkm, .data$levelpathid), + by = c("set" = "ID")) %>% + group_by(ID) %>% + summarise(toID = .data$toID[1], + lengthkm = sum(.data$lengthkm), + areasqkm = sum(.data$areasqkm), + outlet_levelpath = min(.data$levelpathid)) %>% + mutate(toID = ifelse(is.na(.data$toID), 0, .data$toID)) %>% + rename(comid = .data$ID, + tocomid = .data$toID, + nameID = .data$outlet_levelpath)%>% + add_plus_network_attributes() %>% + rename(ID = .data$comid, toID = .data$tocomid, + outlet_nhdpv2_levelpath = .data$nameID, + arbolate_sum = .data$weight) %>% + left_join(select(minimal$fline_sets, .data$ID, .data$set), + by = "ID") + + if(inherits(minimal$fline_sets, "sf")) { + sf::st_sf(min_net) + } else { + min_net + } +} diff --git a/R/catchment_geometry.R b/R/catchment_geometry.R index d098d01..3369ad6 100644 --- a/R/catchment_geometry.R +++ b/R/catchment_geometry.R @@ -11,7 +11,7 @@ add_areasqkm = function(x){ drop_units(set_units(st_area(x), "km2")) -} +} #' Fast POLYGON Union #' @description This is significantly faster then sf::st_union or summarize @@ -20,44 +20,44 @@ add_areasqkm = function(x){ #' @return sf object #' @export #' @importFrom sf as_Spatial st_as_sf st_cast st_make_valid -#' @importFrom dplyr mutate +#' @importFrom dplyr mutate #' @importFrom rgeos gUnaryUnion #' @importFrom methods slot #' @importFrom rlang := union_polygons_geos = function(poly, ID){ - + if(any((types <- sf::st_geometry_type(poly, by_geometry = TRUE) == "GEOMETRYCOLLECTION"))) { poly <- sf::st_cast(poly, "POLYGON") } - + SPDF = as_Spatial(poly) - + rownames(SPDF@data) <- sapply(slot(SPDF, "polygons"), function(x) slot(x, "ID")) - + tmp <- tryCatch({ - gUnaryUnion(spgeom = SPDF, id = poly[[ID]], checkValidity = 0) + gUnaryUnion(spgeom = SPDF, id = poly[[ID]], checkValidity = 0) }, error = function(x){ - gUnaryUnion(spgeom = SPDF, id = poly[[ID]], checkValidity = 2) + gUnaryUnion(spgeom = SPDF, id = poly[[ID]], checkValidity = 2) }) - - + + ids <- as.numeric(sapply(slot(tmp, "polygons"), function(x) slot(x, "ID"))) - + suppressWarnings({ tt = st_as_sf(tmp) %>% mutate("{ID}" := ids) %>% - mutate(areasqkm = add_areasqkm(.)) %>% - st_cast("POLYGON") %>% + mutate(areasqkm = add_areasqkm(.)) %>% + st_cast("POLYGON") %>% st_make_valid() }) - + if(any(grepl("COLLECTION", st_geometry_type(tt)))){ tt = st_collection_extract(tt, "POLYGON") } - + return(tt) - + } #' Fast LINESTRING union @@ -73,16 +73,16 @@ union_polygons_geos = function(poly, ID){ union_linestrings_geos = function(lines, ID){ SPDF = as_Spatial(lines) - + rownames(SPDF@data) <- sapply(slot(SPDF, "lines"), function(x) slot(x, "ID")) - + tmp <- rgeos::gLineMerge(SPDF, byid = TRUE, id = lines[[ID]]) - + ids <- as.numeric(sapply(slot(tmp, "lines"), function(x) slot(x, "ID"))) - - st_as_sf(tmp) %>% - mutate("{ID}" := ids) %>% - flowpaths_to_linestrings() + + st_as_sf(tmp) %>% + mutate("{ID}" := ids) %>% + flowpaths_to_linestrings() } #' Convert MULITLINESTINGS to LINESTRINGS @@ -99,25 +99,25 @@ flowpaths_to_linestrings = function(flowpaths){ sf::st_geometry(multis) = st_line_merge(sf::st_geometry(multis)) } singles = flowpaths[!bool, ] - + bind_rows(multis, singles) } #' Clean Catchment Geometry -#' @description Fixes geometry issues present in catchments that originate in the -#' CatchmentSP layers, or from the reconcile_catchments hyRefactor preocess. -#' These include, but are not limited to disjoint polygon fragments, artifacts +#' @description Fixes geometry issues present in catchments that originate in the +#' CatchmentSP layers, or from the reconcile_catchments hyRefactor preocess. +#' These include, but are not limited to disjoint polygon fragments, artifacts #' from the 30m DEM used to generate the catchments, and non-valid geometry topolgies. -#' A goal of this functions is also to provide means to reduce the data column -#' of the catchments by offering a topology preserving simplification -#' through \code{\link[rmapshaper]{ms_simplify}}. -#' Generally a "keep" parameter of .9 seems appropriate for the resolution of +#' A goal of this functions is also to provide means to reduce the data column +#' of the catchments by offering a topology preserving simplification +#' through \code{\link[rmapshaper]{ms_simplify}}. +#' Generally a "keep" parameter of .9 seems appropriate for the resolution of #' the data but can be modified in function #' @param catchments catchments geometries to fix #' @param ID name of uniquely identifying column -#' @param keep proportion of points to retain in geometry simplification -#' (0-1; default 0.05). See \code{\link[rmapshaper]{ms_simplify}}. +#' @param keep proportion of points to retain in geometry simplification +#' (0-1; default 0.05). See \code{\link[rmapshaper]{ms_simplify}}. #' If NULL, then no simplification will be executed. #' @param crs integer or object compatible with sf::st_crs coordinate reference. #' Should be a projection that supports area-calculations. @@ -129,48 +129,48 @@ flowpaths_to_linestrings = function(flowpaths){ #' @importFrom nhdplusTools rename_geometry #' @importFrom rlang := -clean_geometry = function(catchments, +clean_geometry <- function(catchments, ID = "ID", keep = .9, crs = 5070) { in_crs = st_crs(catchments) - + in_cat <- suppressWarnings({ catchments %>% dplyr::select(ID = !!ID) %>% st_transform(crs) %>% - mutate(areasqkm = add_areasqkm(.)) %>% + mutate(areasqkm = add_areasqkm(.)) %>% ms_explode() %>% filter(!duplicated(.)) %>% mutate(area = add_areasqkm(.)) }) - + ids <- filter(in_cat, duplicated(.data$ID))$ID - + if (length(ids) != 0) { - + # single geometry catchments cat_no_problem <- filter(in_cat, !.data$ID %in% ids) - + # multi geometry catchments challenges = filter(in_cat, .data$ID %in% ids) %>% - mutate(tmpID = 1:n()) %>% + mutate(tmpID = 1:n()) %>% filter(!st_is_empty(.)) - + # base cats: combo of biggest multi geom and single geom catchments base_cats = challenges %>% group_by(.data$ID) %>% slice_max(.data$area) %>% - ungroup() %>% - bind_rows(cat_no_problem) - + ungroup() %>% + bind_rows(cat_no_problem) + # Frags are polygon slivers not in the base catchments - frags = filter(challenges, !.data$tmpID %in% base_cats$tmpID) - + frags = filter(challenges, !.data$tmpID %in% base_cats$tmpID) + # frags smaller than a square meter can be dropped. frags <- filter(frags, area > 0.000001) - + # If fragments exist, then fix, else return empty data.frame if(nrow(frags) > 0){ frags = frags %>% @@ -178,20 +178,20 @@ clean_geometry = function(catchments, ms_explode() %>% mutate(area = as.numeric(st_area(.))) %>% st_make_valid() - + out = tryCatch({ suppressWarnings({ st_intersection(frags, st_make_valid(base_cats)) %>% - st_collection_extract("LINESTRING") + st_collection_extract("LINESTRING") }) }, error = function(e) { NULL }) - + } else { out = data.frame() } - + # message(prettyNum(nrow(frags), big.mark = ",", scientific = FALSE), " fragments to clean...") - + if(!is.null(out) & nrow(out) != 0){ # ints are the LINSTRINGS of intersection between base_cats and frags ints = out %>% @@ -199,19 +199,19 @@ clean_geometry = function(catchments, group_by(.data$rmapshaperid) %>% slice_max(.data$l, with_ties = FALSE) %>% ungroup() - + tj = right_join(frags, - dplyr::select(st_drop_geometry(ints), .data$ID, + dplyr::select(st_drop_geometry(ints), .data$ID, .data$rmapshaperid), by = "rmapshaperid") %>% bind_rows(base_cats) %>% dplyr::select(-.data$rmapshaperid, -.data$areasqkm, -.data$tmpID) %>% group_by(.data$ID) %>% mutate(n = n()) %>% - ungroup() %>% + ungroup() %>% nhdplusTools::rename_geometry('geometry') %>% ungroup() - + if(nrow(tj) == 0){ in_cat = base_cats } else { @@ -221,125 +221,127 @@ clean_geometry = function(catchments, mutate(tmpID = 1:n()) }) } - + } else { in_cat = base_cats } - + # More then one polygon represents a FP dups = in_cat$ID[duplicated(in_cat$ID)] - + if (length(dups) > 0) { - + for (i in 1:length(dups)) { - here = filter(in_cat, .data$ID == dups[i]) - + here = filter(in_cat, .data$ID == dups[i]) + dissolve = slice_min(here, .data$areasqkm, with_ties = FALSE) - + tmap = st_filter(in_cat, dissolve, .predicate = st_touches) - - suppressWarnings({ + + suppressWarnings({ opt <- st_intersection(dissolve, filter(tmap, !.data$tmpID %in% dissolve$tmpID)) - + # If the intersection is only multipoint geometry # Note that geometry_type will be longer than one if more than one feature # was returned from the st_intersection above. gt <- sf::st_geometry_type(opt) if(length(gt) == 1 && grepl("POINT", gt)) { opt <- sf::st_cast(opt, "LINESTRING") - } - + } + opt <- opt %>% st_collection_extract("LINESTRING") %>% mutate(l = st_length(.)) %>% slice_max(.data$l, with_ties = FALSE) - + }) - + ind = which(in_cat$tmpID == opt$tmpID.1) - + in_cat$geometry[ind] = st_union(dissolve$geometry, in_cat$geometry[ind]) in_cat = filter(in_cat, !.data$tmpID %in% dissolve$tmpID) } } } - + if (!is.null(keep)) { # message("Simplifying catchment boundaries: keep = ", keep) in_cat = ms_simplify(in_cat, keep = keep, keep_shapes = TRUE) } - - # since areasqkm will be added based on the new geometries, + + # since areasqkm will be added based on the new geometries, # we need to drop that column if it exists. We want to retain # any others that that came with the OG catchments - + if('areasqkm' %in% names(catchments)){ - catchments = select(catchments, -areasqkm) + catchments = select(catchments, -.data$areasqkm) } - + in_cat %>% mutate(areasqkm = add_areasqkm(.), tmpID = NULL) %>% st_transform(in_crs) %>% - dplyr::select("{ID}" := ID, .data$areasqkm) %>% - left_join(st_drop_geometry(catchments), by = ID) %>% + dplyr::select("{ID}" := ID, .data$areasqkm) %>% + left_join(st_drop_geometry(catchments), by = ID) %>% filter(!duplicated(.)) - + } #' Add Length Map to Refactored Network -#' @description This function replicates the member_COMID column of a refactored -#' network but adds a new notation Following each COMID is '.' which is proceeded -#' by the fraction of that COMID making up the new flowpath. For example 101.1 -#' would indicate 100% of COMID 101 is in the new ID. +#' @description This function replicates the member_COMID column of a refactored +#' network but adds a new notation Following each COMID is '.' which is proceeded +#' by the fraction of that COMID making up the new flowpath. For example 101.1 +#' would indicate 100% of COMID 101 is in the new ID. #' Equally 101.05 would indicate 50% of COMID 101 is present in the new ID'd flowpath #' @param flowpaths a refactored flowpath network containing an member_COMID column -#' @param length_table a table of NHDPlus COMIDs and LENGTH to use as weights. +#' @param length_table a table of NHDPlus COMIDs and LENGTH to use as weights. #' Can be found with \code{nhdplusTools::get_vaa("lengthkm")} #' @return sf object #' @export #' @examples #' \dontrun{ #' path <- system.file("extdata/walker_reconcile.gpkg", package = "hyRefactor") -#' fps <- add_lengthmap(flowpaths = sf::read_sf(path), length_table = nhdplusTools::get_vaa("lengthkm")) +#' fps <- add_lengthmap(flowpaths = sf::read_sf(path), +#' length_table = nhdplusTools::get_vaa("lengthkm")) #' } #'@importFrom dplyr select mutate filter left_join right_join arrange group_by summarize -#'@importFrom tidyr unnest #'@importFrom sf st_drop_geometry st_length st_as_sf #'@importFrom nhdplusTools get_vaa add_lengthmap = function(flowpaths, length_table){ - unnested = dplyr::select(st_drop_geometry(flowpaths), - .data$ID, + tmp = dplyr::select(st_drop_geometry(flowpaths), + .data$ID, COMID = .data$member_COMID) %>% - mutate(COMID = strsplit(.data$COMID, ",")) %>% - tidyr::unnest(cols = .data$COMID) - - unnested2 = filter(unnested, grepl("\\.", COMID)) %>% + mutate(COMID = strsplit(.data$COMID, ",")) + + unnested = data.frame(ID = rep(tmp$ID, times = lengths(tmp$COMID)), + COMID = as.character(unlist(tmp$COMID))) + + unnested2 = filter(unnested, grepl("\\.", COMID)) %>% mutate(baseCOMID = floor(as.numeric(COMID))) + + lengthm_fp = length_table %>% + select(baseCOMID = .data$comid, .data$lengthkm) %>% + mutate(lengthm_fp = .data$lengthkm * 1000, lengthkm = NULL) - lengthm_fp = length_table %>% - select(baseCOMID = .data$comid, .data$lengthkm) %>% - mutate(lengthm_fp = .data$lengthkm * 1000, lengthkm = NULL) - - lengthm_id = flowpaths %>% - mutate(lengthm_id = as.numeric(st_length(.))) %>% - select(.data$ID, .data$lengthm_id) %>% + lengthm_id = flowpaths %>% + mutate(lengthm_id = as.numeric(st_length(.))) %>% + select(.data$ID, .data$lengthm_id) %>% st_drop_geometry() - - map = left_join(left_join(unnested2, lengthm_fp, "baseCOMID"), lengthm_id, by = "ID") %>% - mutate(perLength = round(.data$lengthm_id/.data$lengthm_fp, 3)/10) %>% - select(.data$ID, .data$COMID, .data$perLength) %>% - right_join(unnested, by = c("ID", "COMID")) %>% - mutate(perLength = ifelse(is.na(.data$perLength), 1, - as.character(.data$perLength))) %>% - arrange(.data$ID) %>% - mutate(new = paste0(floor(as.numeric(.data$COMID)), ".", - gsub("0\\.", "", .data$perLength))) %>% - group_by(.data$ID) %>% - summarize(lengthMap = paste(.data$new, collapse = ",")) %>% - right_join(flowpaths) %>% + + map = left_join(left_join(unnested2, lengthm_fp, "baseCOMID"), lengthm_id, by = "ID") %>% + mutate(perLength = round(.data$lengthm_id/.data$lengthm_fp, 3)/10) %>% + select(.data$ID, .data$COMID, .data$perLength) %>% + right_join(unnested, by = c("ID", "COMID")) %>% + mutate(perLength = ifelse(is.na(.data$perLength), 1, + as.character(.data$perLength))) %>% + arrange(.data$ID) %>% + mutate(new = paste0(floor(as.numeric(.data$COMID)), ".", + gsub("0\\.", "", .data$perLength))) %>% + group_by(.data$ID) %>% + summarize(lengthMap = paste(.data$new, collapse = ",")) %>% + right_join(flowpaths) %>% st_as_sf() - + return(map) } diff --git a/R/collapse.R b/R/collapse.R index 0f9fd99..27e61a0 100644 --- a/R/collapse.R +++ b/R/collapse.R @@ -26,6 +26,8 @@ collapse_flowlines <- function(flines, thresh, add_category = FALSE, check_names(flines, "collapse_flowlines") + flines[["toCOMID"]][flines[["toCOMID"]] == 0] <- NA + # very large thresh if (is.null(mainstem_thresh)) { mainstem_thresh_use <- max(flines$LENGTHKM) diff --git a/R/download_fdr_fac.R b/R/download_fdr_fac.R index c97a92d..5bcd86c 100644 --- a/R/download_fdr_fac.R +++ b/R/download_fdr_fac.R @@ -1,13 +1,22 @@ +# Either export from nhdplusTools or leave this +check7z = function (){ + tryCatch({ + system("7z", intern = TRUE) + }, error = function(e) { + stop(simpleError("Please Install 7zip (Windows) or p7zip (MacOS/Unix). Choose accordingly:\n Windows: https://www.7-zip.org/download.html\n Mac: 'brew install p7zip' or 'sudo port install p7zip'\n Linux: https://sourceforge.net/projects/p7zip/")) + }) +} + #' Download Elevation and Derivatives #' @param product character DEM, hydroDEM, or FDRFAC. #' @param out_dir path to directory to store output. #' @param regions character vector of two digit hydrologic -#' @importFrom rvest html_nodes html_attr -#' @importFrom xml2 read_html +#' @importFrom rvest read_html html_nodes html_attr #' @importFrom httr RETRY write_disk progress #' @export download_elev <- function(product, out_dir, regions = NULL) { - dev_null <- nhdplusTools:::check7z() + + dev_null <- check7z() allowable <- c("DEM" = ".*NEDSnapshot.*", "hydroDEM" = ".*HydroDem.*", "FDRFAC" = ".*FdrFac.*") @@ -25,7 +34,7 @@ download_elev <- function(product, out_dir, regions = NULL) { while(i <= length(check)) { url <- check[i] - pg <- read_html(url) + pg <- rvest::read_html(url) pg_links <- paste0(url, html_attr(html_nodes(pg, "a"), "href")) diff --git a/R/reconcile.R b/R/reconcile.R index 68b67f0..afc0e78 100644 --- a/R/reconcile.R +++ b/R/reconcile.R @@ -100,9 +100,12 @@ reconcile_collapsed_flowlines <- function(flines, geom = NULL, id = "COMID") { if (is.null(geom_column)) stop("geom must contain an sf geometry column") - new_flines <- right_join(select(geom, member_COMID = id, geom_column), new_flines, - by = "member_COMID") %>% - sf::st_as_sf() %>% + new_flines <- right_join(select(geom, member_COMID = id, geom_column), + new_flines, + by = "member_COMID") + + new_flines <- new_flines %>% + drop_geometry() %>% group_by(ID) %>% summarise(toID = toID[1], LENGTHKM = LENGTHKM[1], @@ -110,9 +113,9 @@ reconcile_collapsed_flowlines <- function(flines, geom = NULL, id = "COMID") { LevelPathID = LevelPathID[1], Hydroseq = Hydroseq[1], member_COMID = list(unique(member_COMID))) %>% - sf::st_cast("MULTILINESTRING") %>% ungroup() %>% - sf::st_line_merge() + cbind(union_linestrings_geos(select(new_flines, .data$ID), "ID")) %>% + sf::st_as_sf() } return(new_flines) } @@ -126,8 +129,8 @@ reconcile_collapsed_flowlines <- function(flines, geom = NULL, id = "COMID") { #' \code{\link{reconcile_collapsed_flowlines}} #' @param catchment sf data.frame NHDPlus Catchment or CatchmentSP layers for #' included COMIDs -#' @param fdr raster D8 flow direction -#' @param fac raster flow accumulation +#' @param fdr SpatRast D8 flow direction +#' @param fac SpatRast flow accumulation #' @param para integer numer of cores to use for parallel execution #' @param cache path .rda to cache incremental outputs #' @param fix_catchments logical. should catchment geometries be rectified? @@ -145,61 +148,76 @@ reconcile_collapsed_flowlines <- function(flines, geom = NULL, id = "COMID") { #' for this function. #' @details Note that all inputs must be passed in the same projection. #' @export -#' @importFrom magrittr "%>%" -#' @importFrom sf st_drop_geometry st_dimension st_as_sf st_cast st_transform st_precision st_crs -#' @importFrom dplyr select filter mutate left_join -#' @importFrom data.table rbindlist -#' @importFrom methods is +#' @importFrom sf st_crs st_transform st_precision st_drop_geometry st_sf st_multipolygon st_as_sf st_cast st_dimension #' @importFrom nhdplusTools rename_geometry -#' -reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, - fdr = NULL, fac = NULL, para = 2, cache = NULL, - min_area_m = 800, snap_distance_m = 100, - simplify_tolerance_m = 40, vector_crs = NULL, - fix_catchments = TRUE, - keep = .9, crs = 5070) { +#' @importFrom terra crs res +#' @importFrom dplyr select distinct filter group_by summarise ungroup row_number mutate bind_rows +#' @importFrom tidyr separate_rows +#' @importFrom parallel makeCluster stopCluster +#' @importFrom pbapply pblapply +#' @importFrom data.table rbindlist +reconcile_catchment_divides <- function(catchment, + fline_ref, + fline_rec, + fdr = NULL, fac = NULL, + para = 2, cache = NULL, + min_area_m = 800, + snap_distance_m = 100, + simplify_tolerance_m = 40, + vector_crs = NULL, + fix_catchments = TRUE, + keep = .9, crs = 5070) { + in_crs <- st_crs(catchment) catchment <- rename_geometry(catchment, "geom") fline_ref <- rename_geometry(fline_ref, "geom") fline_rec <- rename_geometry(fline_rec, "geom") - # if(!is.null(fdr)){ - # catchment <- st_transform(catchment, st_crs(fdr)) - # st_precision(catchment) <- raster::res(fdr)[1] - # fline_ref <- st_transform(fline_ref, st_crs(fdr)) - # fline_rec <- st_transform(fline_rec, st_crs(fdr)) - # } + if(!is.null(fdr) & !is.null(fac)){ + + if(!inherits(fdr, "SpatRaster")){ + fdr = terra::rast(fdr) + } + + if(!inherits(fdr, "SpatRaster")){ + fac = terra::rast(fac) + } + + catchment <- st_transform(catchment, terra::crs(fdr)) + #st_precision(catchment) <- terra::res(fdr)[1] + fline_ref <- st_transform(fline_ref, terra::crs(fdr)) + fline_rec <- st_transform(fline_rec, terra::crs(fdr)) + } - check_proj(catchment, fline_ref, fdr) - reconciled <- st_drop_geometry(fline_rec) %>% dplyr::select(ID, member_COMID) - + rm(fline_rec) - + # Not all catchments have flowlines. Remove the flowlines without. comid <- fline_ref$COMID # Just being explicit here. featureid <- catchment$FEATUREID # type conversion below is annoying. # as.integer removes the .1, .2, semantic part but the response retains # the semantic component. If you don't know what this means, stop hacking. comid_with_catchment <- comid[as.integer(comid) %in% featureid] - + reconciled <- distinct(reconciled) %>% # had dups from prior steps. tidyr::separate_rows(member_COMID, sep = ",") %>% # Make long form dplyr::filter(member_COMID %in% comid_with_catchment) %>% # dplyr::group_by(ID) %>% - dplyr::summarise(member_COMID = stringr::str_c(member_COMID, collapse = ",")) + dplyr::summarise(member_COMID = paste(member_COMID, collapse = ",")) %>% + dplyr::ungroup() fline_ref <- fline_ref[as.integer(fline_ref$COMID) %in% catchment$FEATUREID, ] - + to_split_bool <- as.numeric(fline_ref$COMID) != as.integer(fline_ref$COMID) - + to_split_ids <- fline_ref$COMID[which(to_split_bool)] - + to_split_featureids <- unique(as.integer(to_split_ids)) - + cl <- NULL if(para > 1) { @@ -213,11 +231,13 @@ reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, } if(!exists("split_cats")) { - split_cats <- pbapply::pblapply(to_split_featureids, par_split_cat, + split_cats <- pbapply::pblapply(to_split_featureids, + par_split_cat, to_split_ids = to_split_ids, fline_ref = fline_ref, catchment = catchment, - fdr = fdr, fac = fac, + fdr = fdr, + fac = fac, min_area_m = min_area_m, snap_distance_m = snap_distance_m, simplify_tolerance_m = simplify_tolerance_m, @@ -225,15 +245,15 @@ reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, cl = cl) if(!is.null(cache)) save(split_cats, file = cache) } - + if(length(split_cats) == 0) { split_cats <- st_sf(FEATUREID = NA, geom = list(sf::st_multipolygon())) } else { - if(!is(split_cats, "sf")) { - split_cats <- st_as_sf(rbindlist(split_cats[!sapply(split_cats, is.null)])) + if(!inherits(split_cats, "sf")) { + split_cats <- sf::st_as_sf(data.table::rbindlist(split_cats[!sapply(split_cats, is.null)])) } - split_cats <- st_cast(split_cats, "MULTIPOLYGON") + split_cats <- sf::st_cast(split_cats, "MULTIPOLYGON") } split_cats <- dplyr::group_by(split_cats, FEATUREID) %>% @@ -243,11 +263,11 @@ reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, unsplit <- get_cat_unsplit(catchment, fline_ref, to_split_featureids) if(nrow(unsplit) > 0) { - split_cats <- st_as_sf(rbindlist(list( + split_cats <- sf::st_as_sf(data.table::rbindlist(list( unsplit, split_cats))) } - + combinations <- reconciled$member_COMID[grepl(",", reconciled$member_COMID)] unioned_cats <- pbapply::pblapply(combinations, @@ -255,22 +275,25 @@ reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, split_cats = split_cats, cache = cache, cl = cl) - - if(!is.null(cl)) + + + if(!is.null(cl)) { parallel::stopCluster(cl) + } if(length(unioned_cats) > 0) { - unioned_cats <- st_as_sf(rbindlist(unioned_cats)) - split_cats <- st_as_sf(rbindlist(list( - filter(split_cats, !FEATUREID %in% unioned_cats$FEATUREID), - unioned_cats))) + unioned_cats <- sf::st_as_sf(data.table::rbindlist(unioned_cats)) + split_cats <- sf::st_as_sf(data.table::rbindlist(list( + dplyr::filter(split_cats, !FEATUREID %in% unioned_cats$FEATUREID), + unioned_cats + ))) } - + out <- st_sf(right_join(dplyr::select(split_cats, member_COMID = FEATUREID), reconciled, - by = "member_COMID")) - + by = "member_COMID")) + missing <- is.na(st_dimension(out$geom)) - + if (any(missing)) { out_mp <- filter(out, !missing) %>% @@ -288,9 +311,9 @@ reconcile_catchment_divides <- function(catchment, fline_ref, fline_rec, if(fix_catchments){ # cat("Fixing Catchment Geometries...\n") clean_geometry(catchments = out, "ID", 0.9) %>% - st_transform(in_crs) + sf::st_transform(in_crs) } else { - st_transform(out, in_crs) + sf::st_transform(out, in_crs) } } @@ -308,7 +331,7 @@ par_split_cat <- function(fid, to_split_ids, fline_ref, catchment, fdr, fac, # nolint end split_set <- to_split_ids[which(grepl(paste0("^", as.character(fid)), to_split_ids))] to_split_flines <- dplyr::filter(fline_ref, COMID %in% split_set) - to_split_cat <- dplyr::filter(catchment, FEATUREID == fid) + to_split_cat <- dplyr::filter(catchment, FEATUREID == fid) split_cats <- split_catchment_divide(catchment = to_split_cat, fline = to_split_flines, @@ -329,6 +352,7 @@ par_split_cat <- function(fid, to_split_ids, fline_ref, catchment, fdr, fac, } get_split_cats <- function(cats, split_cats, cache = NULL) { + error_found <- TRUE error_file <- paste0(gsub(".rda", "", cache), "error.rda") @@ -336,8 +360,8 @@ get_split_cats <- function(cats, split_cats, cache = NULL) { cats_vec <- unlist(strsplit(cats, ",")) union_cats <- dplyr::filter(split_cats, FEATUREID %in% cats_vec) - - if (nrow(union_cats) != length(cats_vec)) { + + if (length(unique(union_cats$FEATUREID)) != length(cats_vec)) { if(!is.null(cache)) { save(list = ls(), file = error_file) @@ -345,7 +369,7 @@ get_split_cats <- function(cats, split_cats, cache = NULL) { cats_vec, "environment saved to:", cache)) } stop(paste("missing a split catchment for an expected flowline.", - cats_vec)) + paste(cats_vec, collapse = "\n"))) } diff --git a/R/refactor_nhdplus.R b/R/refactor_nhdplus.R index 88707ff..688d02b 100644 --- a/R/refactor_nhdplus.R +++ b/R/refactor_nhdplus.R @@ -42,9 +42,10 @@ #' package = "nhdplusTools")) #' #' nhdplus_flowlines <- sf::st_zm(sample_flines) +#' #' refactor_nhdplus(nhdplus_flines = nhdplus_flowlines, #' split_flines_meters = 2000, -#' split_flines_cores = 3, +#' split_flines_cores = 2, #' collapse_flines_meters = 500, #' collapse_flines_main_meters = 500, #' out_refactored = "temp.gpkg", @@ -52,9 +53,9 @@ #' three_pass = TRUE, #' purge_non_dendritic = FALSE, #' warn = FALSE) +#' #' unlink("temp.gpkg") #' unlink("temp_rec.gpkg") -#' refactor_nhdplus <- function(nhdplus_flines, split_flines_meters, diff --git a/R/split_catchment.R b/R/split_catchment.R index f063753..6765244 100644 --- a/R/split_catchment.R +++ b/R/split_catchment.R @@ -117,14 +117,29 @@ flowpath_filter <- function(us, fac_matrix) { return(us) } + + +#' Trace Upstream +#' @param start_point row col index +#' @param cat catchment +#' @param fdr flow direction grid +#' @param fac_matrix flow accumulation matrix +#' @param fdr_matrix flow direction matrix +#' @importFrom terra xFromCol yFromRow +#' @importFrom sf st_sfc st_linestring st_crs +#' @importFrom utils capture.output str +#' @return sfc + trace_upstream <- function(start_point, cat, fdr, fac_matrix, fdr_matrix) { s_p <- st_coordinates(start_point) - suppressWarnings(row_col <- get_row_col(fdr, c(s_p[1], s_p[2]), fac_matrix)) + suppressWarnings(row_col <- get_row_col(fdr, s_p, fac_matrix)) tryCatch({ + us_flowpath <- collect_upstream(row_col, fdr_matrix, fac_matrix, flowpath_mode = TRUE) + }, error = function(e) { stop(paste0("Error with: ", paste(capture.output(str(start_point, @@ -135,10 +150,10 @@ trace_upstream <- function(start_point, cat, fdr, fac_matrix, fdr_matrix) { us_flowpath <- us_flowpath[!is.na(us_flowpath[, 1]), ] - xy_nodes <- matrix(c(raster::xFromCol(fdr, col = us_flowpath[, 2]), - raster::yFromRow(fdr, row = us_flowpath[, 1])), + xy_nodes <- matrix(c(terra::xFromCol(fdr, col = us_flowpath[, 2]), + terra::yFromRow(fdr, row = us_flowpath[, 1])), ncol = 2) - + xy_nodes <- rbind(start_point[[1]], xy_nodes) st_sfc(st_linestring(xy_nodes[nrow(xy_nodes):1, ]), crs = st_crs(start_point)) @@ -151,29 +166,38 @@ trace_upstream <- function(start_point, cat, fdr, fac_matrix, fdr_matrix) { #' @param catchment sf data.frame with one catchment divide #' @param fline sf data.frame with one or more flowline segments in #' upstream downstream order. -#' @param fdr raster a flow direction raster that fully covers the catchment -#' @param fac raster a flow accumulation raster that fuller covers the catchment +#' @param fdr SpatRaster a flow direction SpatRaster that fully covers the catchment +#' @param fac SpatRaster a flow accumulation SpatRaster that fuller covers the catchment #' @param lr boolean should catchments be split along the left/right bank? #' @param min_area_m minimum area in m^2 to filter out slivers (caution, use with care!!) -#' @param snap_distance_m distance in meters to snap raster generated geometry to polygon geometry +#' @param snap_distance_m distance in meters to snap SpatRaster generated geometry to polygon geometry #' @param simplify_tolerance_m dTolerance in meters for simplification of grid-cell based polygons #' @param vector_crs any object compatible with sf::st_crs. Used for vector-based calculations in case that -#' raster projection is not suitable (e.g. lon/lat) -- must result in units of meters. +#' fdr projection is not suitable (e.g. lon/lat) -- must result in units of meters. #' @return Split catchment divides as an sfc geometry. -#' @importFrom raster raster crs crop mask rowColFromCell cellFromXY rasterToPolygons as.matrix +#' @importFrom terra rast setValues as.polygons #' @importFrom dplyr group_by ungroup filter select mutate lead n -#' @importFrom sf st_crs st_crs<- st_coordinates as_Spatial st_buffer st_combine +#' @importFrom sf st_crs st_coordinates as_Spatial st_buffer st_combine #' st_as_sf st_as_sfc st_sfc st_geometry st_simplify st_snap -#' st_difference st_cast st_sf st_area st_distance st_within st_point +#' st_difference st_cast st_sf st_area st_distance st_within st_point `st_crs<-` #' st_make_valid st_segmentize st_nearest_feature st_linestring #' @importFrom nhdplusTools get_node +#' @importFrom dplyr `%>%` #' @export #' split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, min_area_m = 800, snap_distance_m = 100, simplify_tolerance_m = 40, vector_crs = NULL) { - check_proj(catchment, fline, fdr) + #check_proj(catchment, fline, fdr) + + if(!inherits(fdr, "SpatRaster")){ + fdr = terra::rast(fdr) + } + + if(!inherits(fac, "SpatRaster")){ + fac = terra::rast(fac) + } cat_crs <- sf::st_crs(catchment) @@ -184,6 +208,7 @@ split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, ungroup() suppressWarnings(fdr_matrix <- prep_cat_fdr_fac(catchment, fdr, fac)) + fdr <- fdr_matrix$fdr fac <- fdr_matrix$fac fac_matrix <- fdr_matrix$fac_matrix @@ -196,12 +221,13 @@ split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, } for (cat in seq_len(nrow(outlets) - 1)) { - in_out <- st_within(st_sfc(st_point(c(outlets$X[cat], - outlets$Y[cat])), + + in_out <- st_within(st_sfc(st_point(c(outlets$X[cat], outlets$Y[cat])), crs = st_crs(fline)), catchment, prepared = FALSE)[[1]] + if (length(in_out) > 0 && in_out == 1) { - suppressWarnings(row_col <- get_row_col(fdr, c(outlets$X[cat], outlets$Y[cat]), fac_matrix)) + suppressWarnings(row_col <- get_row_col(fdr, start = cbind(outlets$X[cat], outlets$Y[cat]), fac_matrix)) tryCatch({ us_cells <- collect_upstream(row_col, fdr_matrix) @@ -213,21 +239,20 @@ split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, collapse = "\n"), " original error was: \n", e)) }) - out <- matrix(0, nrow = nrow(fdr_matrix), ncol = ncol(fdr_matrix)) - - out[us_cells] <- 1 + vals <- matrix(0, nrow = nrow(fdr_matrix), ncol = ncol(fdr_matrix)) - out <- suppressWarnings(raster::raster(out, template = fdr)) + vals[us_cells] <- 1 + out = terra::setValues(fdr, vals) + names(out) = "cats" + raster_function <- function(x) x == 1 - suppressWarnings(out <- st_as_sf( - raster::rasterToPolygons(out, - fun = raster_function, - dissolve = TRUE))) + out = st_as_sf(terra::as.polygons(out)) |> + filter(.data$cats == 1) smaller_than_one_pixel <- units::set_units(min_area_m, "m^2") - snap_distance <- units::set_units(snap_distance_m, "m") + snap_distance <- units::set_units(snap_distance_m, "m") if(!is.null(vector_crs)) { out <- sf::st_transform(out, vector_crs) @@ -278,7 +303,7 @@ split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, out <- st_as_sfc(out, crs = st_crs(catchment)) if(lr) { - return(split_lr(out, fline, fdr, fac_matrix, fdr_matrix)) + return(split_lr(cat = out, fline, fdr, fac_matrix, fdr_matrix)) } return(out) @@ -287,6 +312,7 @@ split_catchment_divide <- function(catchment, fline, fdr, fac, lr = FALSE, split_lr <- function(cat, fline, fdr, fac_matrix, fdr_matrix) { out <- lapply(c(1:nrow(fline)), function(x, cat, fline, fdr, fac_matrix, fdr_matrix) { + line <- st_cast(st_geometry(fline[x, ]), "LINESTRING") cat <- st_sfc(st_segmentize(cat[[x]], dfMaxLength = 100), crs = st_crs(cat)) @@ -323,9 +349,17 @@ split_lr <- function(cat, fline, fdr, fac_matrix, fdr_matrix) { do.call(c, out) } + + +#' Get Row and Column +#' @param fdr flow direction grid +#' @param start matrix (row, col) +#' @param fac_matrix flow accumulation matrix +#' @importFrom terra cellFromXY rowColFromCell + get_row_col <- function(fdr, start, fac_matrix) { - cell <- raster::cellFromXY(fdr, start) - row_col <- raster::rowColFromCell(fdr, cell) + cell <- terra::cellFromXY(fdr, start) + row_col <- terra::rowColFromCell(fdr, cell) neighbors <- get_neighbor_df(row_col[1], row_col[2], nrow(fac_matrix), ncol(fac_matrix)) @@ -344,22 +378,29 @@ get_row_col <- function(fdr, start, fac_matrix) { return(row_col) } +#' Prep catchment with FDR/FAC +#' @param cat catchment (sf object) +#' @param fdr flow direction grid +#' @param fac flow accumulation grid +#' @importFrom terra vect crop mask as.matrix +#' @importFrom sf st_is_longlat st_buffer + prep_cat_fdr_fac <- function(cat, fdr, fac) { - sp_cat_buffer <- as_Spatial(if(sf::st_is_longlat(cat)) { - st_buffer(cat, 0.002) + + sp_cat_buffer <- terra::vect( + if(sf::st_is_longlat(cat)) { + c = st_make_valid(st_buffer(cat, 0.002)) } else { - st_buffer(cat, 200) + st_make_valid(st_buffer(cat, 200)) }) + + fdr <- terra::crop(fdr, sp_cat_buffer, snap = "out") + fdr <- terra::mask(fdr, sp_cat_buffer) + fac <- terra::crop(fac, sp_cat_buffer, snap = "out") + fac <- terra::mask(fac, sp_cat_buffer) - fdr <- raster::crop(fdr, sp_cat_buffer, - snap = "out") - fdr <- raster::mask(fdr, sp_cat_buffer) - fac <- raster::crop(fac, sp_cat_buffer, - snap = "out") - fac <- raster::mask(fac, sp_cat_buffer) - - fdr_matrix <- raster::as.matrix(fdr) - fac_matrix <- raster::as.matrix(fac) + fdr_matrix <- terra::as.matrix(fdr, wide=TRUE) + fac_matrix <- terra::as.matrix(fac, wide=TRUE) return(list(fdr_matrix = fdr_matrix, fac_matrix = fac_matrix, @@ -368,15 +409,15 @@ prep_cat_fdr_fac <- function(cat, fdr, fac) { } check_proj <- function(catchment, fline, fdr = NULL) { - + er <- "All inputs must have the same projection." - + if(st_crs(fline) != st_crs(catchment)) { stop(er) } - + if(!is.null(fdr)) { - + proj <- st_crs(fdr) if (st_crs(catchment) != st_crs(proj) | st_crs(fline) != st_crs(proj)) { @@ -393,28 +434,34 @@ check_proj <- function(catchment, fline, fdr = NULL) { #' Set edge of fdr and river path to NA to control stop. #' @noRd #' @param start_point sfc point where trace should start -#' @param fdr raster flow direction grid +#' @param fdr SpatRaster flow direction grid #' @param distance numeric max distance in number of grid cells -#' @return sfc linestring path traced along flow direction in crs of rasterr +#' @return sfc linestring path traced along flow direction in crs of fdr grid +#' @importFrom terra rast crs rowColFromCell cellFromXY as.matrix xyFromCell cellFromRowCol +#' @importFrom sf st_transform st_coordinates st_linestring st_sfc st_crs #' @examples +#' #' source(system.file("extdata", "walker_data.R", package = "hyRefactor")) #' #' start_point <- sf::st_sfc(sf::st_point(c(-122.7, 38.126)), crs = 4326) -#' fdr <- walker_fdr #' distance <- 100 #' -#' line <- sf::st_transform(dplyr::filter(walker_flowline, -#' COMID == 5329435), -#' sf::st_crs(fdr)) +#' line <- sf::st_transform(dplyr::filter(walker_flowline, +#' COMID == 5329435), +#' sf::st_crs(walker_fdr)) #' -#' fdr <- raster::mask(fdr, line, inverse = TRUE) +#' fdr <- terra::mask(walker_fdr, line, inverse = TRUE) #' #' xy <- trace_downstream(start_point, fdr, distance) -#' -#' mapview::mapview(list(line, xy)) trace_downstream <- function(start_point, fdr, distance = 10000) { + + if(!inherits(fdr, "SpatRaster")){ + fdr = terra::rast(fdr) + } + + lookup_rowcol <- rep(list(list()), 128) lookup_rowcol[[1]] <- c(0, 1) lookup_rowcol[[2]] <- c(1, 1) @@ -426,12 +473,12 @@ trace_downstream <- function(start_point, fdr, distance = 10000) { lookup_rowcol[[128]] <- c(-1, 1) # start point in raster projection - stp <- sf::st_transform(start_point, sf::st_crs(fdr)) + stp <- sf::st_transform(start_point, terra::crs(fdr)) # 1X2 matrix giving start row and col in raster - sti <- raster::rowColFromCell(fdr, raster::cellFromXY(fdr, sf::as_Spatial(stp))) + sti <- terra::rowColFromCell(fdr, terra::cellFromXY(fdr, st_coordinates(stp))) - fdr_matrix <- raster::as.matrix(fdr) + fdr_matrix <- terra::as.matrix(fdr, wide = TRUE) # empty max distance X 2 matrix track <- matrix(nrow = distance, ncol = 2) @@ -461,7 +508,7 @@ trace_downstream <- function(start_point, fdr, distance = 10000) { track <- track[!is.na(track[,1]), ] # get xy for all row/col from track and fddr - xy <- raster::xyFromCell(fdr, raster::cellFromRowCol(fdr, track[, 1], track[, 2])) + xy <- terra::xyFromCell(fdr, terra::cellFromRowCol(fdr, track[, 1], track[, 2])) # return sfc linestring in projection of fdr sf::st_sfc(sf::st_linestring(xy), crs = sf::st_crs(fdr)) diff --git a/R/split_flowline.R b/R/split_flowline.R index 5dd318c..e9e6c0b 100644 --- a/R/split_flowline.R +++ b/R/split_flowline.R @@ -27,11 +27,7 @@ #' split <- split_flowlines(suppressWarnings(sf::st_cast(sf::st_transform( #' new_hope_flowline, 5070), "LINESTRING")), #' max_length = 2000, events = new_hope_events) -#' -#' mapview::mapview(list(new_hope_events, new_hope_flowline)) -#' -#' mapview::mapview(list(new_hope_events, split)) -#' + split_flowlines <- function(flines, max_length = NULL, events = NULL, para = 0, avoid = NA) { diff --git a/README.md b/README.md index 019d7fe..1dc6520 100644 --- a/README.md +++ b/README.md @@ -7,9 +7,9 @@ ### Installation: ``` -install.packages("devtools") +install.packages("remotes") install.packages("rgeos", repos="http://R-Forge.R-project.org", type="source") -devtools::install_github("dblodgett-usgs/hyRefactor") +remotes::install_github("dblodgett-usgs/hyRefactor") ``` This package is based around the same concepts as [nhdplusTools](https://usgs-r.github.io/nhdplusTools/) and uses its utilities extensively. diff --git a/inst/extdata/gage_01013500_data.R b/inst/extdata/gage_01013500_data.R index f2f04ac..b52fcf4 100644 --- a/inst/extdata/gage_01013500_data.R +++ b/inst/extdata/gage_01013500_data.R @@ -1,54 +1,54 @@ # nolint start -library(raster) -library(dplyr) -library(sf) -library(nhdplusTools) - -extdata <- system.file("extdata", package = "hyRefactor") - -# gage_01013500.gpkg turned into pre-processed sample data. -base = "/Volumes/Transcend/ngen/refactor-tests/base-data/" -# RDS NHD files (dont keep local for space) -cats = readRDS(file.path(base, "catchments_all.rds")) -fps = readRDS(file.path(base, "nhdplus_flowline_update_sb.rds")) - -UT_COMIDs = get_UT(st_drop_geometry(fps), 724696) -flowpaths = make_standalone(filter(fps, COMID %in% UT_COMIDs)) -catchments = filter(cats, FEATUREID %in% UT_COMIDs) - -fac <- raster::raster(file.path(base, "fdrfac/01a_fac.tif")) - -fac = crop(fac, st_transform(catchments, crs(fac)), snap = "out") - -fdr = raster(file.path(base, "fdrfac/01a_fdr.tif")) %>% - crop(st_transform(catchments, crs(fac)), snap = "out") - - -refactor_nhdplus(nhdplus_flines = flowpaths, - split_flines_meters = 10000, - collapse_flines_meters = 1000, - collapse_flines_main_meters = 1000, - split_flines_cores = 2, - out_refactored = "refactor.gpkg", - out_reconciled = "reconcile.gpkg", - three_pass = TRUE, - purge_non_dendritic = FALSE, - warn = FALSE) - -fline_ref <- st_transform(read_sf("refactor.gpkg"), st_crs(fac)) -fline_rec <- st_transform(read_sf("reconcile.gpkg"), st_crs(fac)) - -catchments = st_transform(catchments, st_crs(fac)) -st_precision(catchments) <- raster::res(fdr)[1] - -raw <- reconcile_catchment_divides(catchment = catchments, - fline_ref = fline_ref, - fline_rec = fline_rec, - fdr = fdr, - fac = fac, - para = 1, - cache = NULL, - fix_catchments = FALSE) - -write_sf(raw, file.path(extdata, "gage_01013500.gpkg"), "raw-divides") -unlink(list("refactor.gpkg", "reconcile.gpkg")) +# library(terra) +# library(dplyr) +# library(sf) +# library(nhdplusTools) +# +# extdata <- system.file("extdata", package = "hyRefactor") +# +# # gage_01013500.gpkg turned into pre-processed sample data. +# base = "/Volumes/Transcend/ngen/refactor-tests/base-data/" +# # RDS NHD files (dont keep local for space) +# cats = readRDS(file.path(base, "catchments_all.rds")) +# fps = readRDS(file.path(base, "nhdplus_flowline_update_sb.rds")) +# +# UT_COMIDs = get_UT(st_drop_geometry(fps), 724696) +# flowpaths = make_standalone(filter(fps, COMID %in% UT_COMIDs)) +# catchments = filter(cats, FEATUREID %in% UT_COMIDs) +# +# fac <- terra::rast(file.path(base, "fdrfac/01a_fac.tif")) +# +# fac = crop(fac, st_transform(catchments, crs(fac)), snap = "out") +# +# fdr = terra::rast(file.path(base, "fdrfac/01a_fdr.tif")) %>% +# crop(st_transform(catchments, crs(fac)), snap = "out") +# +# +# refactor_nhdplus(nhdplus_flines = flowpaths, +# split_flines_meters = 10000, +# collapse_flines_meters = 1000, +# collapse_flines_main_meters = 1000, +# split_flines_cores = 2, +# out_refactored = "refactor.gpkg", +# out_reconciled = "reconcile.gpkg", +# three_pass = TRUE, +# purge_non_dendritic = FALSE, +# warn = FALSE) +# +# fline_ref <- st_transform(read_sf("refactor.gpkg"), st_crs(fac)) +# fline_rec <- st_transform(read_sf("reconcile.gpkg"), st_crs(fac)) +# +# catchments = st_transform(catchments, st_crs(fac)) +# st_precision(catchments) <- terra::res(fdr)[1] +# +# raw <- reconcile_catchment_divides(catchment = catchments, +# fline_ref = fline_ref, +# fline_rec = fline_rec, +# fdr = fdr, +# fac = fac, +# para = 1, +# cache = NULL, +# fix_catchments = FALSE) +# +# write_sf(raw, file.path(extdata, "gage_01013500.gpkg"), "raw-divides") +# unlink(list("refactor.gpkg", "reconcile.gpkg")) diff --git a/inst/extdata/new_hope_data.R b/inst/extdata/new_hope_data.R index 1dfe32f..17ced65 100644 --- a/inst/extdata/new_hope_data.R +++ b/inst/extdata/new_hope_data.R @@ -1,13 +1,11 @@ # nolint start -options("rgdal_show_exportToProj4_warnings"="none") -library(rgdal) -library(raster) +library(terra) extdata <- system.file("extdata", package = "hyRefactor") -new_hope_fac <- suppressWarnings(raster::raster(file.path(extdata, "new_hope_fac.tif"))) -new_hope_fdr <- suppressWarnings(raster::raster(file.path(extdata, "new_hope_fdr.tif"))) +new_hope_fac <- suppressWarnings(terra::rast(file.path(extdata, "new_hope_fac.tif"))) +new_hope_fdr <- suppressWarnings(terra::rast(file.path(extdata, "new_hope_fdr.tif"))) -proj <- as.character(raster::crs(new_hope_fdr)) +proj <- as.character(terra::crs(new_hope_fdr)) nhpgpkg <- tempfile(fileext = ".gpkg") @@ -48,9 +46,9 @@ new_hope_events <- sf::read_sf(nhpgpkgev) # new_hope_subset <- subset_nhdplus(UT, "new_hope.gpkg", overwrite = TRUE) # flowline <- read_sf(new_hope_subset, "NHDFlowline_Network") -# fac <- fasterize::raster("~/Documents/Projects/NWM/4_data/nhdplus_raster/fac/NHDPlusSA/NHDPlus03N/NHDPlusFdrFac03a/fac.tif") -# fdr <- fasterize::raster("~/Documents/Projects/NWM/4_data/nhdplus_raster/fdr/NHDPlusSA/NHDPlus03N/NHDPlusFdrFac03a/fdr.tif") -# proj <- as.character(raster::crs(fdr)) +# fac <- terra::rast("~/Documents/Projects/NWM/4_data/nhdplus_raster/fac/NHDPlusSA/NHDPlus03N/NHDPlusFdrFac03a/fac.tif") +# fdr <- terra::rast("~/Documents/Projects/NWM/4_data/nhdplus_raster/fdr/NHDPlusSA/NHDPlus03N/NHDPlusFdrFac03a/fdr.tif") +# proj <- terra::crs(fdr) # # catchment <- read_sf(new_hope_subset, "CatchmentSP") %>% # st_transform(proj) @@ -60,10 +58,10 @@ new_hope_events <- sf::read_sf(nhpgpkgev) # st_buffer(1000) %>% # as_Spatial() # -# sub_fac <- raster::crop(fac, cropper) -# sub_fdr <- raster::crop(fdr, cropper) -# raster::writeRaster(sub_fac, "new_hope_fac.tif", overwrite = TRUE) -# raster::writeRaster(sub_fdr, "new_hope_fdr.tif", overwrite = TRUE) +# sub_fac <- terra::crop(fac, cropper) +# sub_fdr <- terra::crop(fdr, cropper) +# terra::writeRaster(sub_fac, "new_hope_fac.tif", overwrite = TRUE) +# terra::writeRaster(sub_fdr, "new_hope_fdr.tif", overwrite = TRUE) # ##### # # flowline <- sf::read_sf("new_hope.gpkg", "NHDFlowline_Network") diff --git a/inst/extdata/walker_data.R b/inst/extdata/walker_data.R index cc79ce8..5e9cb88 100644 --- a/inst/extdata/walker_data.R +++ b/inst/extdata/walker_data.R @@ -1,20 +1,17 @@ # nolint start -# options("rgdal_show_exportToProj4_warnings"="none") -library(rgdal) -library(raster) +suppressWarnings(library(terra)) extdata <- system.file("extdata", package = "hyRefactor") -walker_fac <- suppressWarnings(raster::raster(file.path(extdata, "walker_fac.tif"))) -walker_fdr <- suppressWarnings(raster::raster(file.path(extdata, "walker_fdr.tif"))) -proj <- as.character(raster::crs(walker_fdr)) +walker_fac <- suppressWarnings(terra::rast(file.path(extdata, "walker_fac.tif"))) +walker_fdr <- suppressWarnings(terra::rast(file.path(extdata, "walker_fdr.tif"))) wgpkg <- tempfile(fileext = ".gpkg") file.copy(file.path(extdata, "walker.gpkg"), wgpkg) walker_catchment <- sf::read_sf(wgpkg, "CatchmentSP") -walker_catchment <- sf::st_transform(walker_catchment, proj) -walker_flowline <- sf::read_sf(wgpkg, "NHDFlowline_Network") -walker_flowline <- sf::st_transform(walker_flowline, proj) +walker_catchment <- sf::st_transform(walker_catchment, terra::crs(walker_fdr)) +walker_flowline <- sf::read_sf(wgpkg, "NHDFlowline_Network") +walker_flowline <- sf::st_transform(walker_flowline, terra::crs(walker_fdr)) # walker.gpkg turned into pre-processed sample data. # run the above then: @@ -37,17 +34,17 @@ walker_flowline <- sf::st_transform(walker_flowline, proj) # r <- fasterize::raster("NHDPlusCA/fdr.tif") # # cropper <- catchment %>% -# st_transform(as.character(raster::crs(r))) %>% +# st_transform(terra::crs(r)) %>% # st_union() %>% # st_buffer(1000) %>% # as_Spatial() # # fac <- fasterize::raster("NHDPlusCA/fac.tif") -# sub_fac <- raster::crop(fac, cropper) -# sub_r <- raster::crop(r, cropper) -# raster::writeRaster(sub_fac, "data-raw/walker_fac.tif", overwrite = TRUE) -# raster::writeRaster(sub_r, "data-raw/walker_fdr.tif", overwrite = TRUE) +# sub_fac <- terra::crop(fac, cropper) +# sub_r <- terra::crop(r, cropper) +# terra::writeRaster(sub_fac, "data-raw/walker_fac.tif", overwrite = TRUE) +# terra::writeRaster(sub_r, "data-raw/walker_fdr.tif", overwrite = TRUE) # nolint end -walker_fline_ref <- sf::read_sf(file.path(extdata, "walker_refactor.gpkg")) -walker_fline_rec <- sf::read_sf(file.path(extdata, "walker_reconcile.gpkg")) +walker_fline_ref <- sf::read_sf(file.path(extdata, "walker_refactor.gpkg")) +walker_fline_rec <- sf::read_sf(file.path(extdata, "walker_reconcile.gpkg")) walker_catchment_rec <- sf::read_sf(file.path(extdata, "walker_cat_rec.gpkg")) diff --git a/man/add_lengthmap.Rd b/man/add_lengthmap.Rd index d1e5447..71427b9 100644 --- a/man/add_lengthmap.Rd +++ b/man/add_lengthmap.Rd @@ -9,22 +9,23 @@ add_lengthmap(flowpaths, length_table) \arguments{ \item{flowpaths}{a refactored flowpath network containing an member_COMID column} -\item{length_table}{a table of NHDPlus COMIDs and LENGTH to use as weights. +\item{length_table}{a table of NHDPlus COMIDs and LENGTH to use as weights. Can be found with \code{nhdplusTools::get_vaa("lengthkm")}} } \value{ sf object } \description{ -This function replicates the member_COMID column of a refactored -network but adds a new notation Following each COMID is '.' which is proceeded -by the fraction of that COMID making up the new flowpath. For example 101.1 -would indicate 100% of COMID 101 is in the new ID. +This function replicates the member_COMID column of a refactored +network but adds a new notation Following each COMID is '.' which is proceeded +by the fraction of that COMID making up the new flowpath. For example 101.1 +would indicate 100% of COMID 101 is in the new ID. Equally 101.05 would indicate 50% of COMID 101 is present in the new ID'd flowpath } \examples{ \dontrun{ path <- system.file("extdata/walker_reconcile.gpkg", package = "hyRefactor") -fps <- add_lengthmap(flowpaths = sf::read_sf(path), length_table = nhdplusTools::get_vaa("lengthkm")) +fps <- add_lengthmap(flowpaths = sf::read_sf(path), +length_table = nhdplusTools::get_vaa("lengthkm")) } } diff --git a/man/aggregate_catchments.Rd b/man/aggregate_catchments.Rd index faafd0f..98172b3 100644 --- a/man/aggregate_catchments.Rd +++ b/man/aggregate_catchments.Rd @@ -48,7 +48,9 @@ source(system.file("extdata", "walker_data.R", package = "hyRefactor")) outlets <- data.frame(ID = c(31, 3, 5, 1, 45, 92), type = c("outlet", "outlet", "outlet", "terminal", "outlet", "outlet"), stringsAsFactors = FALSE) -aggregated <- aggregate_catchments(walker_fline_rec, walker_catchment_rec, outlets) +aggregated <- aggregate_catchments(flowpath = walker_fline_rec, + divide = walker_catchment_rec, + outlets = outlets) plot(aggregated$cat_sets$geom, lwd = 3, border = "red") plot(walker_catchment_rec$geom, lwd = 1.5, border = "green", col = NA, add = TRUE) plot(walker_catchment$geom, lwd = 1, add = TRUE) diff --git a/man/aggregate_network.Rd b/man/aggregate_network.Rd index 5d9a2fd..efbc605 100644 --- a/man/aggregate_network.Rd +++ b/man/aggregate_network.Rd @@ -69,6 +69,8 @@ outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329303, 5329435, 532981 aggregated <- aggregate_network(fline, outlets) +aggregated <- aggregate_network(fline, outlets) + outlets <- dplyr::filter(fline, ID \%in\% outlets$ID) outlets <- nhdplusTools::get_node(outlets) @@ -77,4 +79,5 @@ plot(aggregated$fline_sets$geom, lwd = 3, col = "red") plot(walker_flowline$geom, lwd = .7, col = "blue", add = TRUE) plot(outlets$geometry, add = TRUE) + } diff --git a/man/clean_geometry.Rd b/man/clean_geometry.Rd index 31cc628..487e3dc 100644 --- a/man/clean_geometry.Rd +++ b/man/clean_geometry.Rd @@ -11,8 +11,8 @@ clean_geometry(catchments, ID = "ID", keep = 0.9, crs = 5070) \item{ID}{name of uniquely identifying column} -\item{keep}{proportion of points to retain in geometry simplification -(0-1; default 0.05). See \code{\link[rmapshaper]{ms_simplify}}. +\item{keep}{proportion of points to retain in geometry simplification +(0-1; default 0.05). See \code{\link[rmapshaper]{ms_simplify}}. If NULL, then no simplification will be executed.} \item{crs}{integer or object compatible with sf::st_crs coordinate reference. @@ -22,13 +22,13 @@ Should be a projection that supports area-calculations.} sf object } \description{ -Fixes geometry issues present in catchments that originate in the -CatchmentSP layers, or from the reconcile_catchments hyRefactor preocess. -These include, but are not limited to disjoint polygon fragments, artifacts +Fixes geometry issues present in catchments that originate in the +CatchmentSP layers, or from the reconcile_catchments hyRefactor preocess. +These include, but are not limited to disjoint polygon fragments, artifacts from the 30m DEM used to generate the catchments, and non-valid geometry topolgies. -A goal of this functions is also to provide means to reduce the data column -of the catchments by offering a topology preserving simplification -through \code{\link[rmapshaper]{ms_simplify}}. -Generally a "keep" parameter of .9 seems appropriate for the resolution of +A goal of this functions is also to provide means to reduce the data column +of the catchments by offering a topology preserving simplification +through \code{\link[rmapshaper]{ms_simplify}}. +Generally a "keep" parameter of .9 seems appropriate for the resolution of the data but can be modified in function } diff --git a/man/get_minimal_network.Rd b/man/get_minimal_network.Rd new file mode 100644 index 0000000..6c836cb --- /dev/null +++ b/man/get_minimal_network.Rd @@ -0,0 +1,49 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/aggregate_network.R +\name{get_minimal_network} +\alias{get_minimal_network} +\title{Get Minimal Network} +\usage{ +get_minimal_network(flowpath, outlets) +} +\arguments{ +\item{flowpath}{sf data.frame Flowpaths with ID, toID, LevelPathID, and Hydroseq attributes.} + +\item{outlets}{data.frame with "ID" and "type" columns. "ID" must be identifiers from +fowpath and divide data.frames. "type" should be "outlet", or "terminal". +"outlet" will include the specified ID. +"terminal" will be treated as a terminal node with nothing downstream.} +} +\value{ +a data.frame (potentially including an sfc list column) with +attributes generated by \code{\link[nhdplusTools]{add_plus_network_attributes}} +and a list column "set" containing members of each output flowpath. +} +\description{ +Given a set of outlets, will generate a minimal network by +calling \code{\link{aggregate_network}} and adding nhdplus attributes to the result. + +If geometry is included with the network, it will be merged and returned. +} +\examples{ +source(system.file("extdata", "walker_data.R", package = "nhdplusTools")) +fline <- walker_flowline + +outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329435, 5329817), + type = c("outlet", "outlet", "outlet", "outlet", "outlet")) + +#' Add toCOMID +fline <- nhdplusTools::get_tocomid(fline, add = TRUE) + +# get attributes set +fline <- dplyr::select(fline, ID = comid, toID = tocomid, + levelpathid = levelpathi, hydroseq = hydroseq, + areasqkm = areasqkm, lengthkm = lengthkm) + +min_net <- get_minimal_network(fline, outlets) + +plot(sf::st_geometry(fline), col = "blue") +plot(sf::st_geometry(min_net), lwd = 2, add = TRUE) +plot(sf::st_geometry(nhdplusTools::get_node(min_net)), add = TRUE) + +} diff --git a/man/get_row_col.Rd b/man/get_row_col.Rd new file mode 100644 index 0000000..7a454e8 --- /dev/null +++ b/man/get_row_col.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/split_catchment.R +\name{get_row_col} +\alias{get_row_col} +\title{Get Row and Column} +\usage{ +get_row_col(fdr, start, fac_matrix) +} +\arguments{ +\item{fdr}{flow direction grid} + +\item{start}{matrix (row, col)} + +\item{fac_matrix}{flow accumulation matrix} +} +\description{ +Get Row and Column +} diff --git a/man/prep_cat_fdr_fac.Rd b/man/prep_cat_fdr_fac.Rd new file mode 100644 index 0000000..064a333 --- /dev/null +++ b/man/prep_cat_fdr_fac.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/split_catchment.R +\name{prep_cat_fdr_fac} +\alias{prep_cat_fdr_fac} +\title{Prep catchment with FDR/FAC} +\usage{ +prep_cat_fdr_fac(cat, fdr, fac) +} +\arguments{ +\item{cat}{catchment (sf object)} + +\item{fdr}{flow direction grid} + +\item{fac}{flow accumulation grid} +} +\description{ +Prep catchment with FDR/FAC +} diff --git a/man/reconcile_catchment_divides.Rd b/man/reconcile_catchment_divides.Rd index 722f500..d301a73 100644 --- a/man/reconcile_catchment_divides.Rd +++ b/man/reconcile_catchment_divides.Rd @@ -31,9 +31,9 @@ included COMIDs} \item{fline_rec}{sf data.frame flowpaths as returned by \code{\link{reconcile_collapsed_flowlines}}} -\item{fdr}{raster D8 flow direction} +\item{fdr}{SpatRast D8 flow direction} -\item{fac}{raster flow accumulation} +\item{fac}{SpatRast flow accumulation} \item{para}{integer numer of cores to use for parallel execution} @@ -41,12 +41,12 @@ included COMIDs} \item{min_area_m}{minimum area in m^2 to filter out slivers (caution, use with care!!)} -\item{snap_distance_m}{distance in meters to snap raster generated geometry to polygon geometry} +\item{snap_distance_m}{distance in meters to snap SpatRaster generated geometry to polygon geometry} \item{simplify_tolerance_m}{dTolerance in meters for simplification of grid-cell based polygons} \item{vector_crs}{any object compatible with sf::st_crs. Used for vector-based calculations in case that -raster projection is not suitable (e.g. lon/lat) -- must result in units of meters.} +fdr projection is not suitable (e.g. lon/lat) -- must result in units of meters.} \item{fix_catchments}{logical. should catchment geometries be rectified?} diff --git a/man/refactor_nhdplus.Rd b/man/refactor_nhdplus.Rd index 8b5b118..3a0a34e 100644 --- a/man/refactor_nhdplus.Rd +++ b/man/refactor_nhdplus.Rd @@ -69,9 +69,10 @@ source(system.file("extdata", package = "nhdplusTools")) nhdplus_flowlines <- sf::st_zm(sample_flines) + refactor_nhdplus(nhdplus_flines = nhdplus_flowlines, split_flines_meters = 2000, - split_flines_cores = 3, + split_flines_cores = 2, collapse_flines_meters = 500, collapse_flines_main_meters = 500, out_refactored = "temp.gpkg", @@ -79,9 +80,9 @@ refactor_nhdplus(nhdplus_flines = nhdplus_flowlines, three_pass = TRUE, purge_non_dendritic = FALSE, warn = FALSE) + unlink("temp.gpkg") unlink("temp_rec.gpkg") - } \seealso{ In addition to `prepare_nhdplus` from the nhdplusTools package, diff --git a/man/split_catchment_divide.Rd b/man/split_catchment_divide.Rd index 7825f42..7b20d4a 100644 --- a/man/split_catchment_divide.Rd +++ b/man/split_catchment_divide.Rd @@ -22,20 +22,20 @@ split_catchment_divide( \item{fline}{sf data.frame with one or more flowline segments in upstream downstream order.} -\item{fdr}{raster a flow direction raster that fully covers the catchment} +\item{fdr}{SpatRaster a flow direction SpatRaster that fully covers the catchment} -\item{fac}{raster a flow accumulation raster that fuller covers the catchment} +\item{fac}{SpatRaster a flow accumulation SpatRaster that fuller covers the catchment} \item{lr}{boolean should catchments be split along the left/right bank?} \item{min_area_m}{minimum area in m^2 to filter out slivers (caution, use with care!!)} -\item{snap_distance_m}{distance in meters to snap raster generated geometry to polygon geometry} +\item{snap_distance_m}{distance in meters to snap SpatRaster generated geometry to polygon geometry} \item{simplify_tolerance_m}{dTolerance in meters for simplification of grid-cell based polygons} \item{vector_crs}{any object compatible with sf::st_crs. Used for vector-based calculations in case that -raster projection is not suitable (e.g. lon/lat) -- must result in units of meters.} +fdr projection is not suitable (e.g. lon/lat) -- must result in units of meters.} } \value{ Split catchment divides as an sfc geometry. diff --git a/man/split_flowlines.Rd b/man/split_flowlines.Rd index 0b0bb02..ebbe9cf 100644 --- a/man/split_flowlines.Rd +++ b/man/split_flowlines.Rd @@ -39,11 +39,6 @@ new_hope_flowline <- split <- split_flowlines(suppressWarnings(sf::st_cast(sf::st_transform( new_hope_flowline, 5070), "LINESTRING")), max_length = 2000, events = new_hope_events) - -mapview::mapview(list(new_hope_events, new_hope_flowline)) - -mapview::mapview(list(new_hope_events, split)) - } \seealso{ The \code{\link{refactor_nhdplus}} function implements a complete diff --git a/man/trace_upstream.Rd b/man/trace_upstream.Rd new file mode 100644 index 0000000..ee28481 --- /dev/null +++ b/man/trace_upstream.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/split_catchment.R +\name{trace_upstream} +\alias{trace_upstream} +\title{Trace Upstream} +\usage{ +trace_upstream(start_point, cat, fdr, fac_matrix, fdr_matrix) +} +\arguments{ +\item{start_point}{row col index} + +\item{cat}{catchment} + +\item{fdr}{flow direction grid} + +\item{fac_matrix}{flow accumulation matrix} + +\item{fdr_matrix}{flow direction matrix} +} +\value{ +sfc +} +\description{ +Trace Upstream +} diff --git a/tests/testthat.R b/tests/testthat.R index 0d975af..dbbceb8 100644 --- a/tests/testthat.R +++ b/tests/testthat.R @@ -1,7 +1,5 @@ -options("rgdal_show_exportToProj4_warnings"="none") - -library("testthat") -library("sf") -library("dplyr") -library("nhdplusTools") +library(testthat) +library(sf) +library(dplyr) +library(nhdplusTools) test_check("hyRefactor") diff --git a/tests/testthat/data/split_cat_test_file.R b/tests/testthat/data/split_cat_test_file.R index b74ce75..9171ba1 100644 --- a/tests/testthat/data/split_cat_test_file.R +++ b/tests/testthat/data/split_cat_test_file.R @@ -1,8 +1,8 @@ library(sf) library(hyRefactor) -fdr <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusNE/NHDPlus01/NHDPlusFdrFac01a/fdr/") -fac <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusNE/NHDPlus01/NHDPlusFdrFac01a/fac/") +fdr <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusNE/NHDPlus01/NHDPlusFdrFac01a/fdr/") +fac <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusNE/NHDPlus01/NHDPlusFdrFac01a/fac/") test_cat <- sf::read_sf("~/Documents/active_code/gfv2/workspace/cache/test.gpkg", "nhd_catchment") test_fline <- sf::read_sf("~/Documents/active_code/gfv2/workspace/cache/test.gpkg", "nhd_flowline") @@ -23,7 +23,7 @@ refactor_nhdplus(nhdplus_flines = test_fline, test_fline_ref <- sf::read_sf("~/temp/collapsed.gpkg", "collapsed") test_fline_rec <- sf::read_sf("~/temp/reconciled.gpkg", "reconciled") -crs <- raster::crs(fdr) +crs <- terra::crs(fdr) test_fline_rec <- st_transform(test_fline_rec, crs) test_fline_ref <- st_transform(test_fline_ref, crs) @@ -39,8 +39,8 @@ reconciled_cats <- reconcile_catchment_divides(test_cat, test_fline_ref, catchment <- read_sf("data/split_cat.gpkg", "catchment") fline <- read_sf("data/split_cat.gpkg", "fline") -fdr <- raster::raster("data/split_cat_fdr.tif") -fac <- raster::raster("data/split_cat_fac.tif") +fdr <- terra::rast("data/split_cat_fdr.tif") +fac <- terra::rast("data/split_cat_fac.tif") check <- split_catchment_divide(catchment, fline, fdr, fac) diff --git a/tests/testthat/data/test_outlets.rds b/tests/testthat/data/test_outlets.rds new file mode 100644 index 0000000..412e912 Binary files /dev/null and b/tests/testthat/data/test_outlets.rds differ diff --git a/tests/testthat/error.rda b/tests/testthat/error.rda new file mode 100644 index 0000000..1a61e4d Binary files /dev/null and b/tests/testthat/error.rda differ diff --git a/tests/testthat/test_aggregate_catchments.R b/tests/testthat/test_aggregate_catchments.R index 8dc198a..b51ffcd 100644 --- a/tests/testthat/test_aggregate_catchments.R +++ b/tests/testthat/test_aggregate_catchments.R @@ -1,7 +1,5 @@ context("aggregate catchment") -options("rgdal_show_exportToProj4_warnings"="none") - test_that("walker aggregate runs", { source(system.file("extdata", "walker_data.R", package = "hyRefactor")) @@ -24,23 +22,26 @@ aggregated_cat <- aggregated$cat_sets expect_true(all(aggregated_cat$ID %in% get_id(c("5329385", "5329843", "5329339.1", "5329303")))) expect_equal(sort(aggregated_fline$ID), sort(get_id(c("5329385", "5329843", "5329339.1", "5329303")))) expect_true(aggregated_cat$ID[1] %in% aggregated_cat$set[[1]], "outlet ids should be in the result") -expect_true(length(aggregated_cat$set[[3]]) == 5, "got the wrong number in catchment set") -expect_true(!5 %in% aggregated_cat$set[[3]], "an upstream outlet should not be in another set") +expect_true(all(aggregated_cat$set[aggregated_cat$ID == 31][[1]] %in% + c(29, 30, 85, 86, 31))) +expect_equal(length(aggregated_cat$set[aggregated_cat$ID == 31][[1]]), + 5) +expect_true(!5 %in% aggregated_cat$set[aggregated_cat$ID == 31][[1]], + "an upstream outlet should not be in another set") -expect_true(length(filter(aggregated_fline, ID == 31)$set[[1]]) == 2, "got the wrong number of flowpaths") +expect_true(length(dplyr::filter(aggregated_fline, ID == 31)$set[[1]]) == 2, + "got the wrong number of flowpaths") aggregate_lookup_fline <- dplyr::select(sf::st_drop_geometry(aggregated$fline_sets), ID, set) %>% - tidyr::unnest_longer(col = set) %>% + hyRefactor:::unnest_flines() %>% dplyr::rename(aggregated_ID = ID, reconciled_ID = set) expect_true(!all(walker_fline_rec$ID %in% aggregate_lookup_fline$reconciled_ID), "all input ids should not be in flowline output") -aggregate_lookup_cat <- dplyr::select(sf::st_drop_geometry(aggregated$cat_sets), ID, set) %>% - tidyr::unnest_longer(col = set) %>% - dplyr::rename(aggregated_ID = ID, reconciled_ID = set) +aggregate_lookup_cat <- dplyr::select(sf::st_drop_geometry(aggregated$cat_sets), ID, set) -expect_true(all(walker_fline_rec$ID %in% aggregate_lookup_cat$reconciled_ID), +expect_true(all(walker_fline_rec$ID %in% unlist(aggregate_lookup_cat$set)), "all input ids should be in catchment output") expect_equal(aggregated_cat$toID, get_id(c(NA, "5329843", "5329339.1", "5329303")), info = "Expect these toIDs") @@ -54,9 +55,13 @@ aggregated_fline <- st_transform(aggregated_fline, crs) aggregated_cat <- aggregated_cat[match(aggregated_fline$ID, aggregated_cat$ID), ] -new_geom <- do.call(c, lapply(c(1:nrow(aggregated_cat)), function(g, ac, af, fdr, fac) { - split_catchment_divide(ac[g, ], af[g, ], fdr, fac, lr = TRUE) -}, ac = aggregated_cat, af = aggregated_fline, fdr = walker_fdr, fac = walker_fac)) +new_geom <- do.call(c, lapply(c(1:nrow(aggregated_cat)), function(g) { + split_catchment_divide(catchment = aggregated_cat[g, ], + fline = aggregated_fline[g, ], + fdr = walker_fdr, + fac = walker_fac, + lr = TRUE) +})) expect_true(all(lengths(new_geom) == 2)) diff --git a/tests/testthat/test_aggregate_network.R b/tests/testthat/test_aggregate_network.R index 15c0af6..84ac2ef 100644 --- a/tests/testthat/test_aggregate_network.R +++ b/tests/testthat/test_aggregate_network.R @@ -4,7 +4,7 @@ test_that("example runs", { source(system.file("extdata", "walker_data.R", package = "nhdplusTools")) fline <- dplyr::right_join(dplyr::select(walker_flowline, COMID), - suppressWarnings(prepare_nhdplus(walker_flowline, 0, 0, 0, FALSE))) + suppressWarnings(nhdplusTools::prepare_nhdplus(walker_flowline, 0, 0, 0, FALSE))) fline <- dplyr::select(fline, ID = COMID, toID = toCOMID, LevelPathID = LevelPathI, Hydroseq) @@ -28,3 +28,44 @@ test_that("example runs", { expect_error(aggregate_network(fline, outlets), "Terminal paths must have an NA or 0 toID") }) + +test_that("minimal network", { + source(system.file("extdata", "walker_data.R", package = "nhdplusTools")) + fline <- walker_flowline + + outlets <- data.frame(ID = c(5329357, 5329317, 5329365, 5329435, 5329817), + type = c("outlet", "outlet", "outlet", "outlet", "outlet")) + + #' Add toCOMID + fline <- nhdplusTools::get_tocomid(fline, add = TRUE) + + fline <- dplyr::select(fline, ID = comid, toID = tocomid, + levelpathid = levelpathi, hydroseq = hydroseq, + areasqkm = areasqkm, lengthkm = lengthkm) + + min_net <- get_minimal_network(fline, outlets) + + expect_equal(nrow(min_net), 8) + + expect_s3_class(min_net, "sf") + + min_net <- get_minimal_network(sf::st_drop_geometry(fline), outlets) + + expect_s3_class(min_net, c("tbl_df","tbl","data.frame"), exact = TRUE) + + expect_true(all(outlets$ID %in% min_net$ID)) +}) + +test_that("missing outlet", { + outlets <- as.data.frame(list(ID = c(496338, 21125133, 21047474, 249354, + 21124683, 21124865, 21046242, 255614), + type = c("outlet", "outlet", "outlet", "outlet", + "outlet", "outlet", "outlet", "outlet"))) + + net <- readRDS(list.files(pattern = "test_outlets.rds", full.names = TRUE, recursive = TRUE)) + + outlets <- hyRefactor:::make_outlets_valid(outlets, net) + + expect_true(21047070 %in% outlets$ID) +}) + diff --git a/tests/testthat/test_coastal.R b/tests/testthat/test_coastal.R index 1bd4db8..ac03654 100644 --- a/tests/testthat/test_coastal.R +++ b/tests/testthat/test_coastal.R @@ -36,15 +36,15 @@ source_fdr <- list.files(pattern = "coastal_fdr.tif$", recursive = TRUE, full.na source_fac <- list.files(pattern = "coastal_fac.tif$", recursive = TRUE, full.names = TRUE) min_da_km <- 10 -# fdr <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fdr/") -# fac <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fac/") -# crs <- raster::crs(fac) +# fdr <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fdr/") +# fac <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusCA/NHDPlus18/NHDPlusFdrFac18c/fac/") +# crs <- terra::crs(fac) # -# fdr <- raster::crop(fdr, as_Spatial(st_transform(cats, crs))) -# fac <- raster::crop(fac, as_Spatial(st_transform(cats, crs))) +# fdr <- terra::crop(fdr, vect(st_transform(cats, crs))) +# fac <- terra::crop(fac, vect(st_transform(cats, crs))) # -# raster::writeRaster(fdr, "tests/testthat/data/coastal_fdr.tif") -# raster::writeRaster(fac, "tests/testthat/data/coastal_fac.tif") +# terra::writeRaster(fdr, "tests/testthat/data/coastal_fdr.tif") +# terra::writeRaster(fac, "tests/testthat/data/coastal_fac.tif") test_that("basic coastal aggregation", { nhd <- sf::read_sf(source_gpkg, "NHDFlowline_Network") %>% @@ -101,9 +101,9 @@ test_that("basic coastal aggregation", { exclude_cats = unique(c(avoid$COMID, outlets$COMID, little_terminal$COMID)), warn = FALSE) - fdr <- suppressWarnings(raster::raster(source_fdr)) - fac <- suppressWarnings(raster::raster(source_fac)) - crs <- raster::crs(fac) + fdr <- suppressWarnings(terra::rast(source_fdr)) + fac <- suppressWarnings(terra::rast(source_fac)) + crs <- terra::crs(fac) refactored <- sf::st_transform(sf::read_sf(tf), crs) reconciled <- sf::st_transform(sf::read_sf(tr), crs) diff --git a/tests/testthat/test_collapse.R b/tests/testthat/test_collapse.R index d291119..2704354 100644 --- a/tests/testthat/test_collapse.R +++ b/tests/testthat/test_collapse.R @@ -1,7 +1,8 @@ context("collapse_flowlines") test_that("collapse flowlines works as expected", { - flines <- readRDS("data/petapsco_network.rds") + flines <- readRDS(list.files(pattern = "petapsco_network.rds", + full.names = TRUE, recursive = TRUE)) flines <- sf::st_set_geometry(flines, NULL) flines <- suppressWarnings(nhdplusTools::prepare_nhdplus(flines, 20, 1)) flines_out <- collapse_flowlines(flines, 1) @@ -250,7 +251,7 @@ test_that("repeat collapse doesn't leave orphans", { where = "package:lwgeom", mode = "function")) { - flines <- split_flowlines(flines, 2000, para = 3) + flines <- split_flowlines(flines, 2000, para = 2) flines <- collapse_flowlines(sf::st_set_geometry(flines, NULL), (0.125), TRUE, (0.125)) diff --git a/tests/testthat/test_erie_canal.R b/tests/testthat/test_erie_canal.R index 20a2283..8c79db9 100644 --- a/tests/testthat/test_erie_canal.R +++ b/tests/testthat/test_erie_canal.R @@ -3,60 +3,59 @@ # nhd <- stage_national_data() # atts <- readRDS(nhd$attributes) # UT <- c(get_UT(atts, 21976351), get_UT(atts, 21972746)) -# -# +# +# # nhd <- nhdplusTools::subset_nhdplus(UT, # output_file = "data/erie.gpkg", # return_data = TRUE) -# +# # download_fdr_fac("data", regions = "04a") -# -# fac <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusGL/NHDPlus04/NHDPlusFdrFac04a/fac/") -# fdr <- raster::raster("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusGL/NHDPlus04/NHDPlusFdrFac04a/fdr/") -# proj <- as.character(raster::crs(fdr)) -# +# +# fac <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusGL/NHDPlus04/NHDPlusFdrFac04a/fac/") +# fdr <- terra::rast("~/Documents/active_code/gfv2/workspace/data/fdrfac/NHDPlusGL/NHDPlus04/NHDPlusFdrFac04a/fdr/") +# proj <- as.character(terra::crs(fdr)) +# # catchment <- nhd$CatchmentSP %>% # st_transform(proj) -# +# # cropper <- catchment %>% # st_union() %>% # st_buffer(1000) %>% # as_Spatial() -# -# sub_fac <- raster::crop(fac, cropper) -# sub_fdr <- raster::crop(fdr, cropper) -# raster::writeRaster(sub_fac, "data/erie_fac.tif", overwrite = TRUE) -# raster::writeRaster(sub_fdr, "data/erie_fdr.tif", overwrite = TRUE) +# +# sub_fac <- terra::crop(fac, cropper) +# sub_fdr <- terra::crop(fdr, cropper) +# terra::writeRaster(sub_fac, "data/erie_fac.tif", overwrite = TRUE) +# terra::writeRaster(sub_fdr, "data/erie_fdr.tif", overwrite = TRUE) ##### context("erie") test_that("Erie Canal area works", { -fac <- raster::raster(list.files(pattern = "erie_fac.tif$", full.names = TRUE, recursive = TRUE)) -fdr <- raster::raster(list.files(pattern = "erie_fdr.tif$", full.names = TRUE, recursive = TRUE)) -proj <- as.character(raster::crs(fdr)) +fac <- terra::rast(list.files(pattern = "erie_fac.tif$", full.names = TRUE, recursive = TRUE)) +fdr <- terra::rast(list.files(pattern = "erie_fdr.tif$", full.names = TRUE, recursive = TRUE)) gpkg <- list.files(pattern = "erie.gpkg$", full.names = TRUE, recursive = TRUE) flowline <- sf::read_sf(gpkg, "NHDFlowline_Network") %>% - sf::st_transform(proj) + sf::st_transform(terra::crs(fac)) catchment <- sf::read_sf(gpkg, "CatchmentSP") %>% - sf::st_transform(proj) + sf::st_transform(terra::crs(fac)) catchment <- sf::st_make_valid(catchment) out_refactor <- tempfile(fileext = ".gpkg") out_rec <- tempfile(fileext = ".gpkg") -flowline <- right_join(dplyr::select(flowline, COMID), - prepare_nhdplus(flowline, 0, 0, 0, FALSE, warn = FALSE), +flowline <- dplyr::right_join(dplyr::select(flowline, COMID), + nhdplusTools::prepare_nhdplus(flowline, 0, 0, 0, FALSE, warn = FALSE), by = "COMID") # library(RSQLite) # con <- dbConnect(SQLite(), dbname = "../../../gfv2/workspace/cache/04a.gpkg") # outlets <- RSQLite::dbReadTable(con, "outlets") # dbDisconnect(con) -# +# # outlets <- filter(outlets, ID %in% flowline$COMID) -# +# # saveRDS(outlets, "data/erie_outlets.rds") outlets <- readRDS(list.files(pattern = "erie_outlets.rds$", full.names = TRUE, recursive = TRUE)) @@ -74,20 +73,25 @@ refactor_nhdplus(nhdplus_flines = flowline, fline_ref <- sf::read_sf(out_refactor) %>% - sf::st_transform(proj) + sf::st_transform(terra::crs(fac)) + fline_rec <- sf::read_sf(out_rec) %>% - sf::st_transform(proj) + sf::st_transform(terra::crs(fac)) + +suppressWarnings(cat_rec <- reconcile_catchment_divides(catchment, + fline_ref, + fline_rec, + fdr, + fac, + para = 1)) -suppressWarnings(cat_rec <- reconcile_catchment_divides(catchment, fline_ref, - fline_rec, fdr, fac, - para = 1)) fline_rec$member_COMID <- strsplit(fline_rec$member_COMID, ",") fline_rec$member_COMID <- lapply(fline_rec$member_COMID, as.integer) -ids <- sapply(outlets$ID, +ids <- sapply(outlets$ID, function(x, fline) { - fline$ID[sapply(fline$member_COMID, + fline$ID[sapply(fline$member_COMID, function(l, y) y %in% l, y = x)] }, fline = fline_rec) @@ -99,16 +103,15 @@ ids <- lapply(ids, function(x) x[[1]]) outlets$ID <- as.integer(ids) -outlets_sub <- filter(outlets, !is.na(ID) & ID %in% cat_rec$ID) +outlets_sub <- dplyr::filter(outlets, !is.na(ID) & ID %in% cat_rec$ID) -cat_agg <- aggregate_catchments(fline_rec, cat_rec, - outlets_sub, - da_thresh = 1, +cat_agg <- aggregate_catchments(fline_rec, + cat_rec, + outlets_sub, + da_thresh = 1, only_larger = TRUE) -# mapview::mapview(cat_agg) - -expect_equal(nrow(cat_agg$cat_sets), 513) +expect_equal(nrow(cat_agg$cat_sets), 517) -expect_equal(nrow(cat_agg$fline_sets), 513) +expect_equal(nrow(cat_agg$fline_sets), 517) }) diff --git a/tests/testthat/test_raindrop.R b/tests/testthat/test_raindrop.R index 70d927e..8903e3d 100644 --- a/tests/testthat/test_raindrop.R +++ b/tests/testthat/test_raindrop.R @@ -3,15 +3,13 @@ test_that("basic raindrop", { source(system.file("extdata", "walker_data.R", package = "hyRefactor")) start_point <- sf::st_sfc(sf::st_point(c(-122.7, 38.126)), crs = 4326) - fdr <- walker_fdr distance <- 100 - line <- sf::st_transform(dplyr::filter(walker_flowline, - COMID == 5329435), - sf::st_crs(fdr)) - - fdr <- raster::mask(fdr, line, inverse = TRUE) + line <- sf::st_transform(dplyr::filter(walker_flowline, COMID == 5329435), + sf::st_crs(walker_fdr)) + fdr <- mask(walker_fdr, vect(line), inverse = TRUE) + xy <- hyRefactor:::trace_downstream(start_point, fdr, distance) expect_equal(nrow(sf::st_coordinates(xy)), 19) diff --git a/tests/testthat/test_reconcile_flowline.R b/tests/testthat/test_reconcile_flowline.R index e9c1c1c..5c81865 100644 --- a/tests/testthat/test_reconcile_flowline.R +++ b/tests/testthat/test_reconcile_flowline.R @@ -50,7 +50,7 @@ test_that("collapse works on a double pass", { nhdplus_flines <- readRDS("data/oswego_network.rds") split_flines_meters <- 2000 - split_flines_cores <- 3 + split_flines_cores <- 2 collapse_flines_meters <- 500 collapse_flines_main_meters <- 500 diff --git a/tests/testthat/test_refactor_nhdplus.R b/tests/testthat/test_refactor_nhdplus.R index 502cd50..823b995 100644 --- a/tests/testthat/test_refactor_nhdplus.R +++ b/tests/testthat/test_refactor_nhdplus.R @@ -6,10 +6,14 @@ test_that("refactor_nhdplus works as expected with three pass mode", { where = "package:lwgeom", mode = "function")) { - nhdplus_flines <- sf::st_zm(readRDS("data/north_network.rds")) - + + + nhdplus_flines <- sf::st_zm(readRDS(list.files(pattern = "north_network.rds$", + full.names = TRUE, + recursive = TRUE))) + split_flines_meters <- 2000 - split_flines_cores <- 3 + split_flines_cores <- 2 collapse_flines_meters <- collapse_flines_main_meters <- 1000 out_refactored <- "nhdplus_collapsed.gpkg" out_reconciled <- "nhdplus_reconciled.gpkg" @@ -73,7 +77,7 @@ test_that("The refactor_nhdplus function runs as expected", { m <- suppressWarnings(# Known warnings -- don't want. capture_messages(refactor_nhdplus(nhdplus_flines = nhdplus_flowlines, split_flines_meters = 2000, - split_flines_cores = 3, + split_flines_cores = 2, collapse_flines_meters = 500, collapse_flines_main_meters = 500, out_refactored = "temp.gpkg", @@ -93,7 +97,7 @@ test_that("The refactor_nhdplus function runs as expected", { refactor_nhdplus(nhdplus_flines = nhdplus_flowlines, split_flines_meters = 2000, - split_flines_cores = 3, + split_flines_cores = 2, collapse_flines_meters = 500, collapse_flines_main_meters = 500, out_refactored = "temp.gpkg", diff --git a/tests/testthat/test_split_flowlines.R b/tests/testthat/test_split_flowlines.R index bc90b25..068016e 100644 --- a/tests/testthat/test_split_flowlines.R +++ b/tests/testthat/test_split_flowlines.R @@ -83,7 +83,7 @@ test_that("split lines works", { sf::st_as_sf() %>% sf::st_cast("LINESTRING") %>% sf::st_transform(5070) %>% - split_flowlines(2000, para = 3)) + split_flowlines(2000, para = 2)) expect_true(length(which(grepl("1623361", as.character(flines$COMID)))) == 10) diff --git a/tests/testthat/test_split_reconcile_catchment.R b/tests/testthat/test_split_reconcile_catchment.R index d0946d8..ab5dbe6 100644 --- a/tests/testthat/test_split_reconcile_catchment.R +++ b/tests/testthat/test_split_reconcile_catchment.R @@ -1,7 +1,5 @@ context("split_catchment_divide") -options("rgdal_show_exportToProj4_warnings"="none") - test_that("split_catchment_divide works", { tf <- file.path(tempfile(fileext = ".gpkg")) tr <- file.path(tempfile(fileext = ".gpkg")) @@ -61,8 +59,12 @@ test_that("split_catchment_divide works", { test_fline_ref = st_transform(test_fline_ref, st_crs(walker_fdr)) test_fline_rec = st_transform(test_fline_rec, st_crs(walker_fdr)) - suppressWarnings(reconciled_cats <- reconcile_catchment_divides(test_cat, test_fline_ref, test_fline_rec, - walker_fdr, walker_fac, para = 1)) + suppressWarnings(reconciled_cats <- reconcile_catchment_divides(catchment = test_cat, + fline_ref = test_fline_ref, + fline_rec = test_fline_rec, + fdr = walker_fdr, + fac = walker_fac, + para = 1)) expect_true(nrow(reconciled_cats) == nrow(test_fline_ref), "got the wrong number of split catchments") expect_true(all(reconciled_cats$member_COMID %in% test_fline_ref$COMID)) @@ -156,16 +158,17 @@ test_that("reconcile catchments works with reconciled flowline from split", { # "166755072,8866562.2" # "8833300.1", "8833300.2" - fdr <- suppressWarnings(raster::raster( + fdr <- suppressWarnings(terra::rast( list.files(pattern = "reconcile_test_fdr.tif$", recursive = TRUE, full.names = TRUE))) - fac <- suppressWarnings(raster::raster( + + fac <- suppressWarnings(terra::rast( list.files(pattern = "reconcile_test_fac.tif$", recursive = TRUE, full.names = TRUE))) raster_proj <- sf::st_crs(fdr) - rec_gpkg <- list.files(pattern = "reconcile_test.gpkg", + rec_gpkg <- list.files(pattern = "reconcile_test.gpkg", recursive = TRUE, full.names = TRUE) test_fline_ref <- st_transform( @@ -225,23 +228,22 @@ test_that("doing nothing does nothing", { }) test_that("too small split", { - + # Tests the snap_distance_m parameter. - - catchment <- read_sf("data/split_bug.gpkg", "catchment") - fline <- read_sf("data/split_bug.gpkg", "fline") - fdr <- raster("data/split_bug_fdr.tif") - fac <- raster("data/split_bug_fac.tif") - - expect_error(split_catchment_divide(catchment, fline, fdr, fac, + catchment <- read_sf(list.files(pattern = "split_bug.gpkg$", full.names = TRUE, recursive = TRUE), "catchment") + fline <- read_sf(list.files(pattern = "split_bug.gpkg$", full.names = TRUE, recursive = TRUE), "fline") + + fdr <- rast( list.files(pattern = "split_bug_fdr.tif$", full.names = TRUE, recursive = TRUE)) + fac <- rast( list.files(pattern = "split_bug_fac.tif$", full.names = TRUE, recursive = TRUE)) + + expect_error(split_catchment_divide(catchment, fline, fdr, fac, min_area_m = 800, snap_distance_m = 100), "Nothing left over. Split too small?") - - out <- split_catchment_divide(catchment, fline, fdr, fac, + + out <- split_catchment_divide(catchment, fline, fdr, fac, min_area_m = 800, snap_distance_m = 40) - + expect_equal(length(out), 2) - }) test_that("no fdr", { @@ -267,6 +269,7 @@ test_that("no fdr", { fline_ref <- sf::read_sf(tf) %>% dplyr::arrange(COMID) + fline_rec <- sf::read_sf(tr) cat <- sf::st_transform(walker_catchment, sf::st_crs(fline_rec)) @@ -280,26 +283,28 @@ test_that("no fdr", { }) test_that("merrit dem", { - fdr <- raster::raster("data/ak_fdr.tif") - fac <- raster::raster("data/ak_fac.tif") - cat <- sf::read_sf("data/ak_vector.gpkg", "catchment") - ref <- sf::read_sf("data/ak_vector.gpkg", "fline_ref") - rec <- sf::read_sf("data/ak_vector.gpkg", "fline_rec") + fdr <- terra::rast(list.files(pattern = "ak_fdr.tif$", full.names = TRUE, recursive = TRUE)) + fac <- terra::rast(list.files(pattern = "ak_fac.tif$", full.names = TRUE, recursive = TRUE)) + cat <- sf::read_sf(list.files(pattern = "ak_vector.gpkg$", full.names = TRUE, recursive = TRUE), "catchment") + ref <- sf::read_sf(list.files(pattern = "ak_vector.gpkg$", full.names = TRUE, recursive = TRUE), "fline_ref") + rec <- sf::read_sf(list.files(pattern = "ak_vector.gpkg$", full.names = TRUE, recursive = TRUE), "fline_rec") + rec_cat <- suppressWarnings(reconcile_catchment_divides( - catchment = cat, + catchment = cat, fline_ref = ref, - fline_rec = rec, - fdr, - fac, - para = 1, min_area_m = 10000, - snap_distance_m = 5, - simplify_tolerance_m = 5, - vector_crs = 3338)) - + fline_rec = rec, + fdr, + fac, + para = 1, + min_area_m = 10000, + snap_distance_m = 5, + simplify_tolerance_m = 5, + vector_crs = 3338)) + expect_equal(rec_cat$member_COMID, c("81000012.1", "81000012.2")) - + # plot(rec_cat$geom[2]) # plot(rec_cat$geom[1], add = TRUE, col = "red") - + }) diff --git a/tests/testthat/test_sw.R b/tests/testthat/test_sw.R index 79fa288..814fc28 100644 --- a/tests/testthat/test_sw.R +++ b/tests/testthat/test_sw.R @@ -1,12 +1,24 @@ context("sw catchment split") test_that("long thin sw catchment", { -fac <- raster::raster("data/sw_fac.tif") -fdr <- raster::raster("data/sw_fdr.tif") -proj <- as.character(raster::crs(fdr)) + +fac <- terra::rast(list.files(pattern = "sw_fac.tif$", + full.names = TRUE, + recursive = TRUE)) -flowline <- read_sf("data/sw.gpkg", "fline") %>% +fdr <- terra::rast(list.files(pattern = "sw_fdr.tif$", + full.names = TRUE, + recursive = TRUE)) + +proj <- terra::crs(fdr) + +flowline <- sf::read_sf(list.files(pattern = "sw.gpkg$", + full.names = TRUE, + recursive = TRUE), "fline") %>% st_transform(proj) -catchment <- read_sf("data/sw.gpkg", "catchment") %>% + +catchment <- sf::read_sf(list.files(pattern = "sw.gpkg$", + full.names = TRUE, + recursive = TRUE), "catchment") %>% st_transform(proj) split <- hyRefactor::split_catchment_divide(catchment, flowline, fdr, fac) diff --git a/vignettes/refactor_catchment.Rmd b/vignettes/refactor_catchment.Rmd index 36afba6..669a0fa 100644 --- a/vignettes/refactor_catchment.Rmd +++ b/vignettes/refactor_catchment.Rmd @@ -191,16 +191,16 @@ Given these caveats, we can build a complete lookup table from source catchment refactor_lookup <- dplyr::select(st_drop_geometry(flowline_rec), ID, member_COMID) %>% dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>% - tidyr::unnest_longer(col = member_COMID) %>% + hyRefactor:::unnest_flines(col = "member_COMID") %>% dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates dplyr::rename(reconciled_ID = ID) aggregate_lookup_fline <- dplyr::select(st_drop_geometry(aggregated$fline_sets), ID, set) %>% - tidyr::unnest_longer(col = set) %>% + hyRefactor:::unnest_flines() %>% dplyr::rename(aggregated_flowline_ID = ID, reconciled_ID = set) aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(aggregated$cat_sets), ID, set) %>% - tidyr::unnest_longer(col = set) %>% + hyRefactor:::unnest_flines() %>% dplyr::rename(aggregated_catchment_ID = ID, reconciled_ID = set) (lookup_table <- tibble::tibble(NHDPlusV2_COMID = input_ids) %>% diff --git a/vignettes/refactor_nhdplus.Rmd b/vignettes/refactor_nhdplus.Rmd index 2c0d0f7..d721831 100644 --- a/vignettes/refactor_nhdplus.Rmd +++ b/vignettes/refactor_nhdplus.Rmd @@ -88,7 +88,7 @@ The complete refactor workflow has been packaged into a single function that can ## Not Run ## refactor_nhdplus(nhdplus_flines = nhdplus_flines, split_flines_meters = 2000, - split_flines_cores = 3, + split_flines_cores = 2, collapse_flines_meters = 500, collapse_flines_main_meters = 500, out_refactored = "nhdplus_collapsed.gpkg",