R/subset_nhdplus.R

In nhdplusTools: NHDPlus Tools

#' @title Subset NHDPlus
#' @description Saves a subset of the National Seamless database or other
#' nhdplusTools compatible data based on a specified collection of COMIDs.
#' This function uses \code{\link{get_nhdplus}} for the "download" data
#' source but returns data consistent with local data subsets in a subset
#' file.
#' @param comids integer vector of COMIDs to include.
#' @param output_file character path to save the output to defaults
#' to the directory of the nhdplus_data.
#' @param nhdplus_data character path to the .gpkg or .gdb containing
#' the national seamless database, a subset of NHDPlusHR,
#' or "download" to use a web service to download NHDPlusV2.1 data.
#' Not required if \code{\link{nhdplus_path}} has been set or the default
#' has been adopted. See details for more.
#' @param bbox object of class "bbox" as returned by sf::st_bbox in Latitude/Longitude.
#' If no CRS is present, will be assumed to be in WGS84 Latitude Longitude.
#' @param simplified boolean if TRUE (the default) the CatchmentSP layer
#' will be included. Not relevant to the "download" option or NHDPlusHR data.
#' @param overwrite boolean should the output file be overwritten
#' @param return_data boolean if FALSE path to output file is returned silently otherwise
#' data is returned in a list.
#' @param status boolean should the function print status messages
#' @param flowline_only boolean WARNING: experimental
#' if TRUE only the flowline network and attributes will be returned
#' @param streamorder integer only streams of order greater than or equal will be downloaded.
#' Not implemented for local data.
#' @param out_prj character override the default output CRS of NAD83 lat/lon (EPSG:4269)
#' @details
#'
#' This function relies on the National Seamless Geodatabase or Geopackage.
#' It can be downloaded
#' \href{https://www.epa.gov/waterdata/nhdplus-national-data}{here.}
#'
#' The "download" option of this function should be considered preliminary
#' and subject to revision. It does not include as many layers and may not
#' be available permanently.
#'
#' @return character path to the saved subset geopackage
#' @export

subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NULL, bbox = NULL,
                           simplified = TRUE, overwrite = FALSE, return_data = TRUE, status = TRUE,
                           flowline_only = NULL, streamorder = NULL, out_prj = 4269) {

  if(is.null(flowline_only)) {
    if(!is.null(nhdplus_data) && nhdplus_data == "download") {
      flowline_only <- TRUE
    } else {
      flowline_only <- FALSE
    }
  }

  if(!interactive()) status <- FALSE

  if(!is.null(comids) && length(comids) == 0) stop("comids must be NULL or non-empty")

  if (status) message("All intersections performed in latitude/longitude.")

  if(any(bbox > 180 | bbox < -180)) stop("invalid bbox entry")

  if(!is.null(bbox) && nhdplus_data == "download") {
    x_range <- bbox[3] - bbox[1]
    y_range <- bbox[4] - bbox[2]

    if((x_range * y_range) > 10) {
      warning("Large bounding box submitted for download. Performance may be slow or unstable.")
    }
  }

  if(!is.null(output_file)) {
    if (!grepl("*.gpkg$", output_file)) {
      stop("output_file must end in '.gpkg'")
    }

    if (file.exists(output_file) & !overwrite) {
      stop("output_file exists and overwrite is false.")
    } else if (file.exists(output_file) & overwrite) {
      unlink(output_file)
    }
  }

  if (is.null(nhdplus_data)) {
    nhdplus_data <- nhdplus_path()
  }

  check_nhd_data(nhdplus_data)

  if(is.null(bbox)) {
    if(is.null(comids)) stop("must provide comids or bounding box")

    comids <- round(comids)

    out_list <- c(get_flowline_subset(nhdplus_data, comids,
                                      output_file,
                                      status, out_prj))

    if(!flowline_only) {
      tryCatch({
        out_list <- c(out_list, get_catchment_subset(nhdplus_data, comids,
                                                     output_file, simplified,
                                                     status, out_prj))

        catch_layer <- get_catchment_layer_name(simplified, nhdplus_data)

        envelope <- sf::st_transform(sf::st_as_sfc(sf::st_bbox(out_list[[catch_layer]])),
                                     4326)

        intersection_names <- c("NHDArea", "NHDWaterbody", "NHDFlowline_NonNetwork")
      }, error = function(e) {
        warning(paste("error getting catchment from nhdplus_data\n", e))
      })

      if(!exists("intersection_names")) intersection_names <- c()

    } else {
      intersection_names <- c()
    }

  } else {
    out_list <- list()

    if(!is.null(comids)) warning("using bounding box rather than submitted comids")

    if(!is.null(attr(bbox, "crs"))) {
      envelope <- sf::st_transform(sf::st_as_sfc(bbox),
                                   4326)
    } else {
      if((length(bbox) != 4 | !is.numeric(bbox)) |
         (!(all(bbox >= -180) & all(bbox <= 180)))) stop("invalid bbox entry")
      names(bbox) <- c("xmin", "ymin", "xmax", "ymax")
      bbox <- sf::st_bbox(bbox, crs = sf::st_crs(4326))
      envelope <- sf::st_as_sfc(bbox)
    }

    intersection_names <- c(get_catchment_layer_name(simplified, nhdplus_data),
                            get_flowline_layer_name(nhdplus_data),
                            "NHDArea", "NHDWaterbody", "NHDFlowline_NonNetwork")

    if(flowline_only) intersection_names <- get_flowline_layer_name(nhdplus_data)
  }

  if (nhdplus_data == "download") {

    tryCatch({
    for (layer_name in intersection_names) {
      if(is.null(out_list[layer_name][[1]])) {
        layer <- sf::st_transform(envelope, 4326) %>%
          get_nhdplus_bybox(layer = tolower(layer_name), streamorder = streamorder)
      } else {
        layer <- out_list[layer_name][[1]]
      }

      if(!is.null(nrow(layer)) && nrow(layer) > 0) {

        layer <- hydroloom::check_valid(layer, out_prj)

        if(return_data) {
          out_list[layer_name] <- list(layer)
        }

        if(!is.null(output_file)) {
          sf::write_sf(clean_bbox(layer), output_file, layer_name)
        }
      }
    }
    }, error = function(e) {
      warning(e)
      return(NULL)
    })

  } else {
    if(!flowline_only) {
      if("Gage" %in% st_layers(nhdplus_data)$name) {
        intersection_names <- c(intersection_names, "Gage", "Sink")
      } else {
        intersection_names <- c(intersection_names, "NHDPlusSink")
        intersection_names <- intersection_names[which(intersection_names %in% st_layers(nhdplus_data)$name)]
      }
    }

    out_list <- c(out_list,
                  stats::setNames(lapply(intersection_names, intersection_write,
                                         data_path = nhdplus_data,
                                         envelope = envelope,
                                         output_file = output_file,
                                         status = status,
                                         out_prj = out_prj), intersection_names))
  }

  if(return_data) return(out_list)

  return(output_file)
}

intersection_write <- function(layer_name, data_path, envelope,
                               output_file, status, out_prj) {
  out_list <- list()

  if (status) message(paste("Reading", layer_name))
  layer <- sf::st_zm(sf::read_sf(data_path, layer_name))

  intersection_test <- c()

  try(intersection_test <- suppressMessages(sf::st_intersects(
    sf::st_transform(layer, 4326), envelope)), silent = TRUE)

  found <- lengths(intersection_test)

  if(length(found) > 0) {
    out <- dplyr::filter(layer, found > 0)
  } else {
    out <- data.frame()
  }

  if (nrow(out) > 0) {

    out <- hydroloom::check_valid(out, out_prj)

    if (status) message(paste("Writing", layer_name))
    if (is.null(output_file)) {
      return(out)
    } else {
      sf::write_sf(clean_bbox(out), output_file, layer_name)
      return(layer_name)
    }
  } else {
    if (status) message(paste("No features to write in", layer_name))
  }
}

#' @title Try to find staged NHDPlus data
#' @noRd
check_nhd_data <- function(nhdplus_data) {

  if(nhdplus_data == "download") return(invisible(TRUE))

  if(file.exists(nhdplus_data)) {
    return(invisible(TRUE))
  } else {
    stop("couldn't find nhdplus data")
  }
}

#' @title Get subset of flowline data later.
#' @noRd
get_flowline_subset <- function(nhdplus_data, comids, output_file,
                                status, out_prj) {

  layer_name <- get_flowline_layer_name(nhdplus_data)

  if (status) message(paste("Reading", layer_name))

  if (nhdplus_data == "download") {

    if (length(comids) > 3000) {
      warning("Download functionality not tested for this many comids")
    }

    fline <- get_nhdplus_byid(comids, tolower(layer_name))

    if(is.null(fline)) {
      return(NULL)
    }

  } else {

    if(!layer_name %in% st_layers(nhdplus_data)$name) {
      layer_name <- "NHDFlowline"
    }

    fline <- get_nhd_data(nhdplus_data,layer_name, comids, "COMID", status)

  }

  fline <- hydroloom::check_valid(fline, out_prj)

  if (status) message(paste("Writing", layer_name))

  if(!is.null(output_file)) {
    sf::write_sf(clean_bbox(fline), output_file, layer_name)
  }
  out <- list()
  out[layer_name] <- list(fline)

  return(out)
}

get_nhd_data <- function(nhdplus_data, layer_name, comids, id, status) {

  sets <- lapply(1:ceiling(length(comids) / 1000), function(x) {
    start <- 1000 * (x - 1) + 1

    end <- 1000 * x

    end <- ifelse(end > length(comids), length(comids), end)

    comids[start:end]
  })

  assign("cur_count", 0, envir = nhdplusTools_env)

  out <- lapply(sets, function(x, total) {
    if(status) {

      cur_count <-
        get("cur_count", envir = nhdplusTools_env) + length(x)

      assign("cur_count", cur_count, envir = nhdplusTools_env)

      message(paste(cur_count, "comids of", total))

    }

    align_nhdplus_names(
      sf::read_sf(nhdplus_data,
                  query = get_query(nhdplus_data, layer_name,
                                    id, x)))
  }, total = sum(lengths(sets)))

  do.call(rbind, out)

}

get_query <- function(nhdplus_data, layer_name, id, comids) {
  layer_atts <- sf::read_sf(nhdplus_data,
                            query = paste0("SELECT * FROM ",
                                           layer_name, " LIMIT 1"))

  update_atts <- align_nhdplus_names(layer_atts)

  id_att <- names(layer_atts)[which(names(update_atts) == id)]

  query <- paste0("SELECT * FROM ", layer_name,
                  " WHERE ", id_att, " IN (",
                  paste(comids, collapse = ", "), ")")
}

#' @title Get subset of catchment data layer.
#' @noRd
get_catchment_subset <- function(nhdplus_data, comids, output_file,
                                 simplified, status, out_prj) {

  layer_name <- get_catchment_layer_name(simplified, nhdplus_data)

  if (status) message(paste("Reading", layer_name))

  if (nhdplus_data == "download") {

    catchment <- get_nhdplus_byid(comids, tolower(layer_name))

    if(is.null(catchment)) {
      return(NULL)
    }

  } else {

    catchment <- get_nhd_data(nhdplus_data,layer_name, comids, "FEATUREID", status)

  }

  catchment <- hydroloom::check_valid(catchment, out_prj)

  if (status) message(paste("Writing", layer_name))

  if(!is.null(output_file)) {
    sf::write_sf(clean_bbox(catchment), output_file, layer_name)
  }
  out <- list()
  out[layer_name] <- list(catchment)
  return(out)
}

clean_bbox <- function(x) {
  if("bbox" %in% names(x) && inherits(x$bbox[1], "list")) {
    x$bbox <- sapply(x$bbox, paste, collapse = ",")
  }

  return(x)
}


get_catchment_layer_name <- function(simplified, nhdplus_data) {
  if(is.null(nhdplus_data) || nhdplus_data == "download") { # Only simplified via download
    layer_name <- "CatchmentSP"
  } else {
    if(simplified) { # Can get simplified from local data
      layer_name <- "CatchmentSP"
    } else { # Has to be full catchment
      layer_name <- "Catchment"
    }
    if(!layer_name %in% sf::st_layers(nhdplus_data)$name) # unless it's high res.
      layer_name <- "NHDPlusCatchment"
  }
  return(layer_name)
}

get_flowline_layer_name <- function(nhdplus_data) {
  layer_name <- "NHDFlowline_Network"
  if(nhdplus_data != "download" &&
     !is.null(nhdplus_data) &&
     !layer_name %in% sf::st_layers(nhdplus_data)$name) # unless it's high res.
    layer_name <- "NHDFlowline"
  layer_name
}

#' Subset by Vector Processing Unit
#' @description Calls \link{subset_rpu} for all raster processing units for the
#' requested vector processing unit.
#' @param fline sf data.frame NHD Flowlines with comid, pathlength, lengthkm,
#' hydroseq, levelpathi, rpuid, vpuid, and arbolatesu
#' (dnhydroseq is required if tocomid is not provided).
#' @param vpu character e.g. "01"
#' @param include_null_rpuid logical default TRUE. Note that there are some
#' flowlines that may have a NULL rpuid but be included in the vector
#' processing unit.
#' @param run_make_standalone logical default TRUE
#' should the run_make_standalone function be run on result?
#' @export
#' @importFrom dplyr filter select
#' @return data.frame containing subset network
subset_vpu <- function(fline, vpu,
                       include_null_rpuid = TRUE,
                       run_make_standalone = TRUE) {

  orig_names <- names(fline)

  fline <- check_names(fline, "subset_vpu", tolower = TRUE)

  all_rpuid <- unique(filter(st_drop_geometry(fline),
                             .data$vpuid == vpu)[["rpuid"]])

  all_rpuid <- all_rpuid[(!is.na(all_rpuid) & !is.null(all_rpuid))]

  all_vpu <- lapply(all_rpuid,
                    function(x, fline, run_ms) {
                      subset_rpu(fline, x, run_ms)
                    }, fline = fline, run_ms = run_make_standalone)

  all_vpu <- do.call(rbind, all_vpu)

  if(include_null_rpuid) {

    all_vpu <- rbind(all_vpu, filter(fline, .data$vpuid == vpu &
                                       (is.null(.data$rpuid) |
                                          is.na(.data$rpuid))))

  }

  return(recase_sf(all_vpu, orig_names))

}

recase_sf <- function(x, orig_names) {
  names(x) <- orig_names

  if(inherits(x, "sf")) {
    attr(x, "sf_column") <- orig_names[grepl(attr(x, "sf_column"),
                                             orig_names,
                                             ignore.case = TRUE)]
  }

  x
}


#' Subset by Raster Processing Unit
#' @description Given flowlines and an rpu_code, performs a network-safe subset such
#' that the result can be used in downstream processing. Has been tested to work
#' against the entire NHDPlusV2 domain and satisfies a number of edge cases.
#' @param fline sf data.frame NHD Flowlines with comid, pathlength, lengthkm,
#' hydroseq, levelpathi, rpuid, and arbolatesu (dnhydroseq
#' is required if tocomid is not provided).
#' @param rpu character e.g. "01a"
#' @param run_make_standalone logical default TRUE
#' should the run_make_standalone function be run on result?
#' @param strict logical if TRUE, paths that extend outside the RPU but
#' have no tributaries in the upstream RPU will be included in the output.
#' @export
#' @return data.frame containing subset network
#' @importFrom dplyr filter arrange summarize
#' @importFrom sf st_sf st_drop_geometry
subset_rpu <- function(fline, rpu, run_make_standalone = TRUE, strict = FALSE) {
  orig_names <- names(fline)

  fline <- check_names(fline, "subset_rpu", tolower = TRUE)

  if(!any(c("tocomid") %in% names(fline))) {
    # derive tocomid from dnhydroseq
    # this avoids using get_tocomid() which requires additional attributes.
    fline$tocomid <- fline$comid[match(fline$dnhydroseq, fline$hydroseq)]
    fline$tocomid[is.na(fline$tocomid)] <- 0
  }

  # Filter so only navigable flowlines are included.
  fline <- fline[fline$comid %in% get_all_navigable(fline, rpu), ]

  # For flowlines labeled as in the RPU, find the top and bottom of each
  # LevelPath. This was required for some unique network situations.
  fline_sub_in <- filter(st_drop_geometry(fline), .data$rpuid %in% rpu)

  fline_sub_in <- group_by(fline_sub_in, .data$levelpathi)

  fline_sub_in <- summarize(fline_sub_in,
                            lp_top = max(.data$hydroseq),
                            lp_bot = min(.data$hydroseq))

  if(!strict) {
    # find flowlines completely outside the RPU
    fline_sub_out <- filter(st_drop_geometry(fline), !.data$rpuid %in% rpu)

    # Need all paths in the domain outside the RPU and their tributary relations
    # used just below. This filter avoids flowlines along the same levelpath.
    lp_with_trib <-  filter(select(fline_sub_out, "levelpathi", "dnlevelpat"),
                            .data$levelpathi != .data$dnlevelpat)

    # join to the paths within the RPU so we have the lp_top to filter on.
    fline_sub_out <- left_join(fline_sub_out, ungroup(fline_sub_in),
                               by = "levelpathi")

    # filter so we have the stuff that is complete upstream of the rpu
    fline_sub_out <- filter(group_by(fline_sub_out, "levelpathi"),
                            .data$hydroseq > .data$lp_top)

    # want to keep anything left that does not have anything flowing to it.
    # to do this, we need fline_sub_out to include things we CAN remove later.

    # first select only levelpath and dnlevelpat and get rid of cruft
    fline_sub_out <- distinct(select(ungroup(fline_sub_out),
                                     "levelpathi", "dnlevelpat"))

    # levelpaths with tributaries that are outside the domain.
    lp_with_trib <- filter(lp_with_trib, .data$dnlevelpat %in% fline_sub_out$levelpathi)
    fline_sub_out <- filter(fline_sub_out, .data$levelpathi %in% lp_with_trib$dnlevelpat)
  }


  # Using the levelpath top and bottoms found above, filter the complete
  # domain to the hydrosequence of the levelpath top and bottoms instead
  # of trusting the rpuid to be useable.
  fline <- left_join(fline, fline_sub_in, by = "levelpathi")

  # first filter to levelpaths included in our domain.
  fline <- group_by(filter(fline, .data$levelpathi %in% fline_sub_in$levelpathi),
                    .data$levelpathi)

  if(strict) {

    # filter to the top and bottom that are required to connect things
    # fully within the rpu
    fline <- filter(fline, .data$hydroseq >= .data$lp_bot &
                      .data$hydroseq <= .data$lp_top)

  } else if(!nrow(fline_sub_out) == 0) {

    # fline_sub_out contains levelpaths that we can subject to the lp_top / lp_bot
    # filter. All other levelpaths should be left whole.

    fline <- filter(fline,
                    # for levelpaths in flinesub out, filter hydrosequence
                    (.data$levelpathi %in% fline_sub_out$levelpathi &
                      (.data$hydroseq >= .data$lp_bot &
                         .data$hydroseq <= .data$lp_top)) |
                      # keep everything not in fline_sub_out
                      !.data$levelpathi %in% fline_sub_out$levelpathi)

  }

  fline <- select(ungroup(fline), -"lp_top", -"lp_bot")

  if(run_make_standalone) {
    fline <- make_standalone(fline)
  }

  if(!any(grepl("tocomid", orig_names, ignore.case = TRUE))) {
    fline <- select(fline, -"toCOMID") # case if required for make_standalone
  }

  recase_sf(fline, orig_names)
}

#' @noRd
get_all_navigable <- function(fline, rpu) {

  # Find all outlets of current rpu and sort by size
  outlets <- filter(st_drop_geometry(fline), .data$rpuid %in% rpu)

  outlets <- filter(outlets, .data$hydroseq == .data$terminalpa |
                      !.data$tocomid %in% .data$comid)

  outlets <- arrange(outlets, desc(.data$arbolatesu))

  network <- get_sorted(fline[c("comid", "tocomid")],
                        split = FALSE, outlets = outlets$comid)

  network$comid
}

#' @noRd
get_nhdplus_byid <- function(comids, layer, streamorder = NULL) {

  if(layer == "nhdflowline_network"){
    query_usgs_oafeat(ids = comids, type = "nhd", filter = streamorder_filter_cql(streamorder))
  } else if(layer == "catchmentsp"){
    query_usgs_oafeat(ids = comids, type = "catchment")
  } else {
    stop("Layer must be one of catchmentsp, nhdflowline_network")
  }
}

#' @noRd
get_nhdplus_bybox <- function(box, layer, streamorder = NULL) {

  if(!layer %in% c("nhdarea", "nhdwaterbody", "nhdflowline_network",
                   "nhdflowline_nonnetwork", "catchmentsp")) {
    stop("Layer must be one of nhdarea, nhdwaterbody, nhdflowline_network, nhdflowline_nonnetwork, catchmentsp.")
  }

  type <- dplyr::filter(query_usgs_oafeat(),
                        .data$layer_name == layer)$user_call


  if(layer == "nhdflowline_network") {
    query_usgs_oafeat(AOI = box,
                         type = type,
                         filter = streamorder_filter_cql(streamorder))
  } else {
    query_usgs_oafeat(AOI = box,
                         type = type)
  }
}

#' @noRd
streamorder_filter_cql <- function(streamorder) {
  if(is.null(streamorder)){ return(NULL)}

  paste0("streamorde%20>%20", streamorder - 1)
}