In dblodgett-usgs/hygeo: Hydrologic Geospatial Fabric Creation Functions

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6,
  fig.height=4,
  eval=nzchar(Sys.getenv("BUILD_VIGNETTES")),
  cache=FALSE
)
options(scipen = 9999)
options("rgdal_show_exportToProj4_warnings"="none")

Source dependencies and set up a bunch of files.

library(nhdplusTools)
library(hyRefactor)
library(sf)
library(tidyr)
library(dplyr)
library(hygeo)

src_gpkg <- system.file("gpkg/sugar_creek_fort_mill.gpkg", package = "hygeo")
fdr_tif <- system.file("tiff/sugar_creek_fort_mill_fdr.tif", package = "hygeo")
fac_tif <- system.file("tiff/sugar_creek_fort_mill_fac.tif", package = "hygeo")
collapsed <- tempfile(fileext = ".gpkg")
refactored_small <- tempfile(fileext = ".gpkg")
reconciled <- tempfile(fileext = ".gpkg")
reconciled_small <- tempfile(fileext = ".gpkg")

unlink(c(collapsed, refactored_small, reconciled, reconciled_small), force = TRUE)

catchment_prefix <- "cat-"
waterbody_prefix <- "fp-"
nexus_prefix <- "nex-"

Now we will get some data and just plot up our area of interest.

options("rgdal_show_exportToProj4_warnings"="none")
nhd <- nhdplusTools::plot_nhdplus(list(9731454),
                                  gpkg = src_gpkg,
                                  overwrite = FALSE,
                                  nhdplus_data = src_gpkg)

For this demonstration we'll just use NWIS Sites as they exist. Here we go grab upstream with tributaries NWIS Sites and index them to the network. This is all placeholder for doing this using a more formal process for selection of network locations.

nwis <- nhdplusTools::navigate_nldi(list(featureSource = "comid", 
                                         featureID = 9731454),
                                    mode = "UT", 
                                    data_source = "nwissite", 
                                    distance_km = 9999)

if(is.list(nwis)) {
  nwis <- nwis$UT_nwissite
}

what_nwis_data <- dataRetrieval::whatNWISdata(siteNumber = gsub("USGS-", "", nwis$identifier))

nwis_sites <- filter(what_nwis_data, parm_cd == "00060" & data_type_cd == "uv") %>%
  st_as_sf(coords = c("dec_long_va", "dec_lat_va"), 
                       crs = 4269) %>%
  st_transform(5070)

nwis_sites <-  bind_cols(nwis_sites,
                         left_join(data.frame(id = seq_len(nrow(nwis_sites))),
                                   nhdplusTools::get_flowline_index(
                                     st_transform(nhd$flowline, 5070), 
                                     nwis_sites, search_radius = 50), by = "id")) %>%
  filter(!is.na(COMID)) %>%
  st_sf() %>%
  select(site_no, COMID, REACHCODE, REACH_meas, offset) %>%
  left_join(select(st_drop_geometry(nhd$flowline), COMID, FromMeas, ToMeas), by = "COMID")

# Check if we have catchment polygons for all our gage outlets.
all(nwis_sites$COMID %in% nhd$catchment$FEATUREID)

# only chose sites that are 1/4 or more up the catchment.
split_sites <- nwis_sites %>%
  filter((100 * (REACH_meas - FromMeas) / (ToMeas - FromMeas)) > params$gage_tolerance)

nwis_sites$split <- nwis_sites$site_no %in% split_sites$site_no

nhdplusTools::plot_nhdplus(list(9731454),
                                  gpkg = src_gpkg,
                                  overwrite = FALSE,
                                  nhdplus_data = src_gpkg,
                                  actually_plot = TRUE)
plot(st_transform(nwis_sites$geometry, 3857), pch = 24, bg = "darkgrey", add = TRUE)

hyRefactor::refactor_nhdplus(nhd$flowline,
                             split_flines_meters = params$split_m,
                             split_flines_cores = 2,
                             collapse_flines_meters = params$collapse_m,
                             collapse_flines_main_meters = params$collapse_m,
                             out_refactored = collapsed,
                             out_reconciled = reconciled,
                             three_pass = TRUE,
                             purge_non_dendritic = FALSE, 
                             events = split_sites)

collapse <- sf::read_sf(collapsed)
reconcile <- sf::read_sf(reconciled)

slope <- select(st_drop_geometry(reconcile), ID, member_COMID) %>%
  mutate(member_COMID = strsplit(member_COMID, ",")) %>%
  unnest(cols = member_COMID) %>%
  mutate(member_COMID = floor(as.numeric(member_COMID))) %>%
  left_join(select(st_drop_geometry(nhd$flowline), COMID, slope), 
            by = c("member_COMID" = "COMID")) %>%
  group_by(ID) %>%
  summarise(slope = mean(slope))

reconcile <- left_join(reconcile, slope, by = "ID")

nwis_sites <- left_join(nwis_sites, 
                        select(st_drop_geometry(collapse), event_REACHCODE, split_COMID = COMID),
                        by = c("REACHCODE" = "event_REACHCODE"))

nwis_sites$local_id <- 
  sapply(seq_len(nrow(nwis_sites)), 
         function(x, nwis_sites) {

           if(!is.na(nwis_sites$split_COMID[x])) {
             checker <- nwis_sites$split_COMID[x]
           } else {
             checker <- as.character(nwis_sites$COMID[x])
           }

           id <- reconcile[grepl(checker, 
                                 reconcile$member_COMID), ]
           if(nrow(id) > 1) {
             id <- filter(id, Hydroseq == min(Hydroseq))
           }

           paste0(waterbody_prefix, id$ID)
         }, nwis_sites = nwis_sites)


nwis_sites <- select(nwis_sites, -split, -split_COMID)

fdr <- raster::raster(fdr_tif)
fac <- raster::raster(fac_tif)

crs <- raster::crs(fdr)

nhd$catchment <- sf::st_transform(nhd$catchment, crs)
reconcile <- sf::st_transform(sf::st_sf(reconcile), crs)
collapse <- sf::st_transform(sf::st_sf(collapse), crs)

sf::st_precision(nhd$catchment) <- 30

reconcile_divides <- hyRefactor::reconcile_catchment_divides(nhd$catchment,
                                                             fdr = fdr,
                                                             fac = fac,
                                                             fline_ref = collapse,
                                                             fline_rec = reconcile,
                                                             para = 1)

network_order <- dplyr::select(sf::st_drop_geometry(nhd$flowline), COMID, Hydroseq)

nwis_lookup <- dplyr::select(sf::st_drop_geometry(nwis_sites),
                            site_no, local_id) %>%
  dplyr::mutate(local_id = gsub(waterbody_prefix, catchment_prefix, local_id))

nhd_crosswalk <- c(get_nhd_crosswalk(reconcile, catchment_prefix, network_order, nwis_lookup),
                   get_nhd_crosswalk(reconcile, waterbody_prefix, network_order, nwis_lookup))

write_sf(reconcile, "../tests/testthat/data/sugar_creek_hyRefactor.gpkg", "reconcile")
write_sf(reconcile_divides, "../tests/testthat/data/sugar_creek_hyRefactor.gpkg", "reconcile_divides")

We can now pass reconciled divides into the hygeo package functions to generate catchment areas, waterbodies, and nexuses.

nexus <- get_nexus(reconcile, 
                   nexus_prefix = nexus_prefix)

catchment_edge_list <- get_catchment_edges(reconcile,
                                           catchment_prefix = catchment_prefix,
                                           nexus_prefix = nexus_prefix)
waterbody_edge_list <- get_waterbody_edge_list(reconcile,
                                               waterbody_prefix = waterbody_prefix)
sqkm_per_sqm <- 1 / 1000^2

reconcile_divides$area_sqkm <- as.numeric(st_area(st_transform(reconcile_divides, 5070))) * sqkm_per_sqm

catchment_data <- get_catchment_data(reconcile_divides,
                                     catchment_edge_list,
                                     catchment_prefix = catchment_prefix)

reconcile <- dplyr::rename(reconcile, length_km = LENGTHKM)

flowpath_data <- get_flowpath_data(reconcile,
                                     waterbody_edge_list,
                                     catchment_prefix = waterbody_prefix) %>%
      mutate(realized_catchment = gsub(waterbody_prefix, 
                                     catchment_prefix, ID))

nexus_data <- get_nexus_data(nexus,
                             catchment_edge_list)

mapview::mapview(catchment_data, layer.name = "Catchment Area", col.regions = "tan") +
  mapview::mapview(flowpath_data, layer.name = "Flowpaths", color = "blue") +
  mapview::mapview(nexus_data, layer.name = "Nexuses", cex = 3, color = "yellow4", col.regions = "yellow4") + 
  mapview::mapview(nwis_sites, layer.name = "NWIS Sites", cex = 6, color = "darkgrey", col.regions = "grey")

The outputs can be rendered into csv or json:

hygeo_list <- list(catchment = catchment_data, 
                   flowpath = flowpath_data, 
                   nexus = nexus_data,
                   catchment_edges = catchment_edge_list,
                   waterbody_edges = waterbody_edge_list)

class(hygeo_list) <- "hygeo"

out_path <- params$out_path
dir.create(out_path, recursive = TRUE, showWarnings = FALSE)

out_path <- write_hygeo(hygeo_list, out_path = out_path, edge_map = TRUE, overwrite = TRUE)

nhd_crosswalk <- lapply(nhd_crosswalk, function(x) {
  out <- list(COMID = x$COMID, outlet_COMID = jsonlite::unbox(x$outlet_COMID))

  if(!is.null(x$site_no)) {
    out <- c(out, list(site_no = jsonlite::unbox(x$site_no)))
  }
  out
})

jsonlite::write_json(nhd_crosswalk, file.path(out_path, "crosswalk.json"), 
                     pretty = TRUE, auto_unbox = FALSE)

(hygeo_list_read <- read_hygeo(out_path))

nhd_crosswalk

If needed, we can further break up the flowline network ensuring the outlets in the catchment network are included.

outlet_comids <- select(st_drop_geometry(reconcile), ID, member_COMID) %>%
  mutate(member_COMID = strsplit(member_COMID, ",")) %>%
  unnest(cols = member_COMID) %>%
  mutate(member_COMID = as.integer(member_COMID)) %>%
  left_join(select(st_drop_geometry(nhd$flowline), COMID, Hydroseq),
            by = c("member_COMID" = "COMID")) %>%
  group_by(ID) %>%
  filter(Hydroseq == min(Hydroseq)) %>%
  ungroup()

hyRefactor::refactor_nhdplus(nhd$flowline,
                             split_flines_meters = 500,
                             split_flines_cores = 2,
                             collapse_flines_meters = 400,
                             collapse_flines_main_meters = 400,
                             exclude_cats = outlet_comids$member_COMID,
                             out_refactored = refactored_small,
                             out_reconciled = reconciled_small,
                             three_pass = TRUE,
                             purge_non_dendritic = FALSE)