diff --git a/.Rbuildignore b/.Rbuildignore index 400e4bd8..09964125 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,6 +1,6 @@ ^.*\.Rproj$ ^\.Rproj\.user$ -^.*cache +^.*cache$ LICENSE .travis.yml .lintr diff --git a/DESCRIPTION b/DESCRIPTION index 0597cb3a..b564dfcd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nhdplusTools Type: Package Title: NHDPlus Tools -Version: 1.0.1 +Version: 1.1.0 Authors@R: c(person(given = "David", family = "Blodgett", role = c("aut", "cre"), @@ -28,7 +28,7 @@ URL: https://doi-usgs.github.io/nhdplusTools/ https://github.com/doi-usgs/nhdplu BugReports: https://github.com/doi-usgs/nhdplusTools/issues/ Depends: R (>= 4.0) -Imports: hydroloom, dataRetrieval, dplyr, sf, units, magrittr, jsonlite, httr, xml2, R.utils, utils, tidyr, methods, maptiles, mapsf, fst, arrow, tools, zip, pbapply +Imports: hydroloom, dataRetrieval, dplyr, sf, units, magrittr, jsonlite, httr, xml2, R.utils, utils, tidyr, methods, maptiles, mapsf, fst, arrow, tools, zip, pbapply, memoise, digest Suggests: testthat, knitr, rmarkdown, ggmap, ggplot2, lwgeom, gifski, leaflet, httptest, future, future.apply License: CC0 Encoding: UTF-8 diff --git a/NAMESPACE b/NAMESPACE index a8201e73..554b8292 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,6 +14,7 @@ export(download_rf1) export(download_vaa) export(download_wbd) export(fix_flowdir) +export(get_3dhp) export(get_DD) export(get_DM) export(get_UM) @@ -26,8 +27,6 @@ export(get_flowline_index) export(get_gagesII) export(get_hr_data) export(get_huc) -export(get_huc12) -export(get_huc8) export(get_hydro_location) export(get_levelpaths) export(get_nhdarea) @@ -64,6 +63,7 @@ export(make_standalone) export(map_nhdplus) export(navigate_network) export(navigate_nldi) +export(nhdplusTools_cache_settings) export(nhdplusTools_data_dir) export(nhdplus_path) export(plot_nhdplus) @@ -77,6 +77,7 @@ export(subset_rpu) export(subset_vpu) importFrom(arrow,open_dataset) importFrom(arrow,s3_bucket) +importFrom(digest,digest) importFrom(dplyr,across) importFrom(dplyr,all_of) importFrom(dplyr,any_of) @@ -103,6 +104,7 @@ importFrom(fst,read.fst) importFrom(httr,GET) importFrom(httr,POST) importFrom(httr,RETRY) +importFrom(httr,content) importFrom(httr,progress) importFrom(httr,write_disk) importFrom(hydroloom,accumulate_downstream) @@ -130,6 +132,9 @@ importFrom(hydroloom,sort_network) importFrom(hydroloom,st_compatibalize) importFrom(jsonlite,fromJSON) importFrom(magrittr,`%>%`) +importFrom(memoise,cache_filesystem) +importFrom(memoise,cache_memory) +importFrom(memoise,memoise) importFrom(methods,as) importFrom(methods,is) importFrom(pbapply,pbapply) @@ -137,6 +142,7 @@ importFrom(pbapply,pblapply) importFrom(pbapply,pboptions) importFrom(sf,"st_geometry<-") importFrom(sf,read_sf) +importFrom(sf,sf_use_s2) importFrom(sf,st_as_sf) importFrom(sf,st_as_sfc) importFrom(sf,st_bbox) @@ -149,6 +155,7 @@ importFrom(sf,st_filter) importFrom(sf,st_geometry) importFrom(sf,st_geometry_type) importFrom(sf,st_layers) +importFrom(sf,st_make_valid) importFrom(sf,st_multilinestring) importFrom(sf,st_nearest_feature) importFrom(sf,st_sf) diff --git a/NEWS.md b/NEWS.md index 079d6fb6..0a31c693 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,22 @@ -nhdplusTools 1.0.1 +nhdplusTools 1.1.0 ========== +This release has a significant consolidation of code to work with web services. +It incorporates an ArcGIS REST service for the first time. It has very few +backward compatibility issues with previous versions and introduces new functions +and behavior. + +A key change in behavior is the introduction of caching with the `memoise` package. +See #366, pull request #364, and documentation in `nhdplusTools_cache_settings()` +for more. + - fixed nhd and nhdplushr urls #368 - fix issue with nhdplusTools_data_dir() #365 - fixed bug with large nhdplus downloads with empty tiles. #361 +- Added 3DHP_all service client. #363 +- Removed deprecated function `get_huc12` and `get_huc8` +- Updated documentation of `get_huc()` and other web service functions. +- added `nhdplusTools_cache_settings()` to control use of a `memoise` cache. #366 nhdplusTools 1.0.0 ========== diff --git a/R/A_nhdplusTools.R b/R/A_nhdplusTools.R index 7a3730c0..a4fa1cc2 100644 --- a/R/A_nhdplusTools.R +++ b/R/A_nhdplusTools.R @@ -93,6 +93,9 @@ assign("nhdplus_attributes", nhdplus_attributes, envir = nhdplusTools_env) assign("geoserver_root", "https://labs.waterdata.usgs.gov/geoserver/", envir = nhdplusTools_env) +assign("arcrest_root", "https://hydro.nationalmap.gov/arcgis/rest/services/", + envir = nhdplusTools_env) + assign("split_flowlines_attributes", c("COMID", "toCOMID", "LENGTHKM"), envir = nhdplusTools_env) @@ -350,6 +353,76 @@ nhdplus_path <- function(path = NULL, warn = FALSE) { } } +#' @title nhdplusTools cache settings +#' @description +#' Provides an interface to adjust nhdplusTools `memoise` cache. +#' +#' Mode and timeout can also be set using environment variables. +#' `NHDPLUSTOOLS_MEMOISE_CACHE` and `NHDPLUSTOOLS_MEMOISE_TIMEOUT` are +#' used unless overriden with this function. +#' +#' @param mode character 'memory' or 'filesystem' +#' @param timeout numeric number of seconds until caches invalidate +#' @return list containing settings at time of calling. If inputs are +#' NULL, current settings. If settings are altered, previous setting values. +#' @export +#' +nhdplusTools_cache_settings <- function(mode = NULL, timeout = NULL) { + current_mode <- get("nhdpt_mem_cache", envir = nhdplusTools_env) + current_timeout <- get("nhdpt_cache_timeout", envir = nhdplusTools_env) + + if(!is.null(mode) && mode %in% c("memory", "filesystem")) { + assign("nhdpt_mem_cache", mode, envir = nhdplusTools_env) + } + + if(!is.null(timeout) && is.numeric(timeout)) { + assign("nhdpt_cache_timeout", timeout, envir = nhdplusTools_env) + } + + return(invisible(list(mode = current_mode, timeout = current_timeout))) +} + +#' @importFrom memoise memoise cache_memory cache_filesystem +#' @importFrom digest digest +nhdplusTools_memoise_cache <- function() { + sys_memo_cache <- Sys.getenv("NHDPLUSTOOLS_MEMOISE_CACHE") + ses_memo_cache <- try(get("nhdpt_mem_cache", envir = nhdplusTools_env), silent = TRUE) + + # if it hasn't been set up yet, try to use the system env + if(!inherits(ses_memo_cache, "try-error")) { + return(ses_memo_cache) + } else { + if(sys_memo_cache == "memory") { + memoise::cache_memory() + } else { + dir.create(nhdplusTools_data_dir(), showWarnings = FALSE, recursive = TRUE) + + memoise::cache_filesystem(nhdplusTools_data_dir()) + } + + } + +} + +nhdplusTools_memoise_timeout <- function() { + sys_timeout <- Sys.getenv("NHDPLUSTOOLS_MEMOISE_TIMEOUT") + ses_timeout <- try(get("nhdpt_cache_timeout", envir = nhdplusTools_env), silent = TRUE) + + # if it hasn't been set up yet, try to use the system env + if(!inherits(ses_timeout, "try-error")) { + return(ses_timeout) + } else { + if(sys_timeout != "") { + as.numeric(sys_timeout) + } else { + # default to one day + oneday_seconds <- 60 * 60 * 24 + } + } +} + +assign("nhdpt_mem_cache", nhdplusTools_memoise_cache(), envir = nhdplusTools_env) +assign("nhdpt_cache_timeout", nhdplusTools_memoise_timeout(), envir = nhdplusTools_env) #' @title Align NHD Dataset Names #' @description this function takes any NHDPlus dataset and aligns the attribute names with those used in nhdplusTools. diff --git a/R/arcrest_tools.R b/R/arcrest_tools.R new file mode 100644 index 00000000..6386d871 --- /dev/null +++ b/R/arcrest_tools.R @@ -0,0 +1,219 @@ +get_3dhp_service_info <- memoise::memoise(function() { + + # TODO: support more services? + server <- "3DHP_all" + + url_base <- paste0(get("arcrest_root", envir = nhdplusTools_env), + server, + "/MapServer/") + + all_layers <- jsonlite::read_json(paste0(url_base, "?f=json")) + + list(url_base = url_base, all_layers = all_layers) + +}) + +# this function is intended to be a very general internal utility. It's initial +# implementation is quite specific but this should become more general over time. +#' @title Query USGS Hydro ESRI Rest Server +#' @description Query the USGS Hydro ESRI Rest Server for spatial data by location, +#' area, or ID. +#' @details The returned object(s) will have the same +#' Spatial Reference System (SRS) as the input AOI. If a individual or set of +#' IDs are used to query, then the default CRS of EPSG:4269 is +#' preserved. In all cases, a user-defined SRS can be passed to \code{t_srs} +#' which will override all previous SRS (either input or default). +#' All buffer and distance operations are handled internally using in +#' EPSG:5070 Albers Equal Area projection +#' @param AOI sf (MULTI)POINT or (MULTI)POLYGON. An 'area of interest' can +#' be provided as either a location (sf POINT) or area (sf POLYGON) +#' in any Spatial Reference System. +#' @param ids character. A set of identifier(s) from the data +#' type requested, for 3dhp, this is id3dhp. +#' @param type character. Type of feature to return +#' ("hydrolocation", "flowline", "waterbody", "drainage area", "catchment"). +#' If NULL (default) a data.frame of available resources is returned +#' @param where character. An where clause to pass to the server. +#' @param t_srs character (PROJ string or EPSG code) or numeric (EPSG code). +#' A user specified - target -Spatial Reference System (SRS/CRS) for returned objects. +#' Will default to the CRS of the input AOI if provided, and to 4326 for ID requests. +#' @param buffer numeric. The amount (in meters) to buffer a POINT AOI by for an +#' extended search. Default = 0.5 +#' @return a simple features (sf) object or valid types if no type supplied +#' @keywords internal +#' @importFrom sf st_crs st_geometry_type st_buffer st_transform st_zm read_sf st_bbox st_as_sfc +#' @importFrom httr RETRY content +#' @importFrom dplyr filter +#' @importFrom methods as +query_usgs_arcrest <- memoise::memoise(function(AOI = NULL, ids = NULL, + type = NULL, where = NULL, + t_srs = NULL, + buffer = 0.5){ + + si <- get_3dhp_service_info() + + source <- data.frame(user_call = sapply(si$all_layers$layers, \(x) tolower(x$name)), + layer = sapply(si$all_layers$layers, \(x) x$id)) + + group_layers <- si$all_layers$layers[sapply(si$all_layers$layers, + \(x) grepl("Group Layer", x$type))] + + if(is.null(type)) { + message("`type` input must be one of: \n\t\"", + paste(source$user_call, collapse = "\"\n\t\""), "\"") + return(source) + } + + if(!type %in% source$user_call) { + warning("\"", type, "\" not in `type` input. Must be one of: \n\t\"", + paste(source$user_call, collapse = "\"\n\t\""), "\"") + return(NULL) + } + + if(grepl(paste(sapply(group_layers, \(x) x$name), + collapse = "|"), + type, ignore.case = TRUE)) { + layer_id <- filter(source, .data$user_call == !!type)$layer + + group_layer <- group_layers[[sapply(group_layers, \(x) x$id == layer_id)]] + + need_layers <- as.integer(group_layer$subLayerIds) + } else { + need_layers <- as.integer(filter(source, .data$user_call == !!type)$layer) + } + + all_out <- rep(list(list()), length(need_layers)) + + for(l in seq_len(length(all_out))) { + + layer <- need_layers[l] + + type <- filter(source, .data$layer == !!layer)$user_call + + AOI <- check_query_params(AOI, ids, type, where, source, t_srs, buffer) + t_srs <- AOI$t_srs + AOI <- AOI$AOI + + here <- filter(source, .data$user_call == !!type) + + URL <- paste0(si$url_base, here$layer, "/query") + + spat_filter <- spatial_filter(AOI, tile = FALSE, format = "esri") + + if(!is.null(ids)) { + + if(!is.null(where)) stop("can't specify ids and where") + + where <- paste0("id3dhp IN ('", + paste(ids, collapse = "', '"), "')") + } + + post_body <- list(where = where, + f = "json", + returnIdsOnly = "true") + + if(length(spat_filter[[1]]) > 0) { + post_body <- c(list(geometry = spat_filter[[1]]), + post_body) + } + + tryCatch({ + if(nhdplus_debug()) { + message(paste(URL, "\n")) + message(post_body) + } + + all_ids <- content(RETRY("POST", + URL, + body = post_body, + encode = "form")) + all_ids <- unlist(all_ids$objectIds) + + }, error = function(e) { + warning("Something went wrong trying to access a service.") + out <- NULL + }) + + length_ids <- length(all_ids) + + if(is.null(all_ids) | length(all_ids) == 0) { + warning(paste("No", here$user_call, "features found in area of interest."), call. = FALSE) + out <- NULL + } else { + + chunk_size <- 2000 + all_ids <- split(all_ids, ceiling(seq_along(all_ids)/chunk_size)) + + out <- rep(list(list()), length(all_ids)) + + for(i in seq_along(all_ids)) { + + if(length_ids > chunk_size) { + top <- i * chunk_size + if(top > length_ids) top <- length_ids + message("Getting features ", (i - 1) * chunk_size, " to ", top, " of ", length_ids) + } + + post_body <- list(objectIds = paste(all_ids[[i]], collapse = ","), + outFields = "*", + f = "geojson") + + tryCatch({ + if(nhdplus_debug()) { + message(paste(URL, "\n")) + message(post_body) + } + + out[[i]] <- rawToChar(httr::RETRY("POST", + URL, + body = post_body, + encode = "form")$content) + }, error = function(e) { + warning("Something went wrong trying to access a service.") + out <- NULL + }) + + out[[i]] <- tryCatch({ + sf::st_zm(sf::read_sf(out[[i]]))}, + error = function(e) NULL) + } + + if(inherits(out[[1]], "data.frame")) { + out <- bind_rows(unify_types(out)) + + out <- check_valid(out[!duplicated(out[["id3dhp"]]), ], + out_prj = t_srs) + + } else { + + out <- NULL + + } + + if(any(is.null(out), nrow(out) == 0)) { + + out = NULL + + } else if(!is.null(AOI)){ + + out = sf::st_filter(sf::st_transform(out, t_srs), + sf::st_transform(AOI, t_srs)) + + if(nrow(out) == 0){ + warning(paste(here$user_call, "features found but were outside area of interest polygon.")) + out = NULL + + } + } + + if(!is.null(out)) { + out <- select(out, -any_of("OBJECTID")) + } else { + warning(paste("No", here$user_call, "features found"), call. = FALSE) + out <- NULL + } + } + all_out[[l]] <- out + } + tryCatch(sf::st_sf(data.table::rbindlist(all_out)), error = function(e) NULL) +}, ~memoise::timeout(nhdplusTools_memoise_timeout()), cache = nhdplusTools_memoise_cache()) diff --git a/R/geoserver_tools.R b/R/geoserver_tools.R index 39a98fcd..c70b1be5 100644 --- a/R/geoserver_tools.R +++ b/R/geoserver_tools.R @@ -25,19 +25,16 @@ #' @return a simple features (sf) object #' @keywords internal #' @importFrom sf st_crs st_geometry_type st_buffer st_transform st_zm read_sf st_bbox st_as_sfc +#' @importFrom sf sf_use_s2 #' @importFrom httr POST RETRY #' @importFrom dplyr filter #' @importFrom methods as +#' @importFrom xml2 read_xml -query_usgs_geoserver <- function(AOI = NULL, ids = NULL, +query_usgs_geoserver <- memoise::memoise(function(AOI = NULL, ids = NULL, type = NULL, filter = NULL, t_srs = NULL, - buffer = 0.5){ - - # If t_src is not provided set to AOI CRS - if(is.null(t_srs)){ t_srs <- sf::st_crs(AOI)} - # If AOI CRS is NA (e.g st_crs(NULL)) then set to 4326 - if(is.na(t_srs)) { t_srs <- 4326} + buffer = 0.5) { source <- data.frame(server = 'wmadata', user_call = c('huc02', 'huc04', 'huc06', @@ -80,44 +77,22 @@ query_usgs_geoserver <- function(AOI = NULL, ids = NULL, if(is.null(type)){ return(source) } - if(!is.null(AOI) & !is.null(ids)){ - # Check if AOI and IDs are both given - stop("Either IDs or a spatial AOI can be passed.", .call = FALSE) - } else if(is.null(AOI) & is.null(ids)){ - # Check if AOI and IDs are both NULL - stop("IDs or a spatial AOI must be passed.",.call = FALSE) - } else if(!(type %in% source$user_call)) { - # Check that "type" is valid - stop(paste("Type not available must be one of:", - paste(source$user_call, collapse = ", ")), - call. = FALSE) - } - - if(!is.null(AOI)){ - - if(length(sf::st_geometry(AOI)) > 1) { - stop("AOI must be one an only one feature.") - } + AOI <- check_query_params(AOI, ids, type, NULL, source, t_srs, buffer) + t_srs <- AOI$t_srs + AOI <- AOI$AOI - if(sf::st_geometry_type(AOI) == "POINT"){ - # If input is a POINT, buffer by 1/2 meter (in equal area projection) - AOI = sf::st_buffer(sf::st_transform(AOI, 5070), buffer) %>% - sf::st_bbox() %>% sf::st_as_sfc() %>% sf::st_make_valid() - } - } + here <- filter(source, .data$user_call == !!type) - here <- dplyr::filter(source, .data$user_call == !!type) + URL <- paste0(get("geoserver_root", envir = nhdplusTools_env), + here$server, + "/ows") - URL <- paste0(get("geoserver_root", envir = nhdplusTools_env), - here$server, - "/ows") + use_s2 <- sf_use_s2() + sf_use_s2(FALSE) - use_s2 <- sf::sf_use_s2() - sf::sf_use_s2(FALSE) + on.exit(sf_use_s2(use_s2), add = TRUE) - on.exit(sf::sf_use_s2(use_s2), add = TRUE) - - out <- spatial_filter(AOI, type = here$geoserver, tile = here$page, geom_name = here$geom_name) + out <- spatial_filter(AOI, tile = here$page, geom_name = here$geom_name) for(i in 1:length(out)) { @@ -156,19 +131,19 @@ query_usgs_geoserver <- function(AOI = NULL, ids = NULL, tryCatch({ if(nhdplus_debug()) { message(paste(URL, "\n")) - message(as.character(xml2::read_xml(filterXML))) + message(as.character(read_xml(filterXML))) } - out[[i]] <- rawToChar(httr::RETRY("POST", - URL, - body = filterXML)$content) + out[[i]] <- rawToChar(RETRY("POST", + URL, + body = filterXML)$content) }, error = function(e) { warning("Something went wrong trying to access a service.") return(NULL) }) out[[i]] <- tryCatch({ - sf::st_zm(sf::read_sf(out[[i]]))}, + st_zm(read_sf(out[[i]]))}, error = function(e) return(NULL)) } @@ -188,8 +163,8 @@ query_usgs_geoserver <- function(AOI = NULL, ids = NULL, } else if(!is.null(AOI)){ - out = sf::st_filter(sf::st_transform(out, t_srs), - sf::st_transform(AOI, t_srs)) + out = st_filter(st_transform(out, t_srs), + st_transform(AOI, t_srs)) if(nrow(out) == 0){ out = NULL } @@ -202,7 +177,7 @@ query_usgs_geoserver <- function(AOI = NULL, ids = NULL, NULL } -} +}, ~memoise::timeout(nhdplusTools_memoise_timeout()), cache = nhdplusTools_memoise_cache()) unify_types <- function(out) { all_class <- bind_rows(lapply(out, function(x) { @@ -250,36 +225,16 @@ assign("bb_break_size", value = 2, nhdplusTools_env) #' @noRd #' @importFrom sf st_geometry_type st_buffer st_transform st_bbox -spatial_filter <- function(AOI, type = 'catchmentsp', +spatial_filter <- function(AOI, break_size = get("bb_break_size", nhdplusTools_env), tile = TRUE, - geom_name = "the_geom"){ + geom_name = "the_geom", format = 'ogc') { if(is.null(AOI)) return(list(list())) - # "The current GeoServer implementation will return all the streets whose - # ENVELOPE overlaps the BBOX provided. This may includes some nearby - # streets that do not strictly fall within the bounding box. - # This is a fast response - if a more accurate but slower response is required, - # use Intersects or Within." - - # if (as.character(sf::st_geometry_type(AOI)) == "POINT") { - # - # AOI <- sf::st_coordinates(AOI) - # - # return(paste0("?service=WFS", - # "&version=1.0.0", - # "&request=GetFeature", - # "&typeName=wmadata:", type, - # "&outputFormat=application%2Fjson", - # "&srsName=EPSG:4326", - # "&CQL_FILTER=INTERSECTS%28the_geom,%20POINT%20%28", - # AOI[1], "%20", AOI[2], "%29%29")) - # } - - - bb <- sf::st_transform(AOI, 4326) %>% - sf::st_bbox() + + bb <- st_transform(AOI, 4326) %>% + st_bbox() if(tile) { x_breaks <- seq(bb$xmin, bb$xmax, length.out = (ceiling((bb$xmax - bb$xmin) / break_size) + 1)) @@ -293,15 +248,37 @@ spatial_filter <- function(AOI, type = 'catchmentsp', } else { bb_list <- list(bb) } - lapply(bb_list, function(bb) { - paste0('', - '', geom_name, '', - '', - '', bb$ymin, " ", bb$xmin, '', - '', bb$ymax, " ", bb$xmax, '', - '', - '') - }) + + if(format == "ogc") { + lapply(bb_list, function(bb) { + paste0('', + '', geom_name, '', + '', + '', bb$ymin, " ", bb$xmin, '', + '', bb$ymax, " ", bb$xmax, '', + '', + '') + }) + } else if(format == "esri") { + + # { + # "xmin": -109.55, + # "ymin": 25.76, + # "xmax": -86.39, + # "ymax": 49.94, + # "spatialReference": { + # "wkid": 4326 + # } + # } + + lapply(bb_list, function(bb) { + jsonlite::toJSON( + c(bb, list(spatialReference = list(wkid = 4326))), + auto_unbox = TRUE) + }) + } else { + stop("format must be ogc or esri") + } } #' @title Construct an attribute filter for geoservers @@ -382,7 +359,7 @@ tc <- function(x) { #' @importFrom httr RETRY GET #' @importFrom jsonlite fromJSON -extact_comid_nwis <- function(nwis){ +extact_comid_nwis <- memoise::memoise(function(nwis){ # We could export this from dataRetrieval dataRetrieval:::pkg.env$nldi_base #but currently its not... baseURL <- paste0(get_nldi_url(), "/linked-data/") @@ -390,4 +367,40 @@ extact_comid_nwis <- function(nwis){ c <- rawToChar(httr::RETRY("GET", url)$content) f.comid <- jsonlite::fromJSON(c, simplifyVector = TRUE) f.comid$features$properties$comid +}, ~memoise::timeout(nhdplusTools_memoise_timeout()), cache = nhdplusTools_memoise_cache()) + +#' @importFrom sf st_make_valid st_as_sfc st_bbox st_buffer st_transform +check_query_params <- function(AOI, ids, type, where, source, t_srs, buffer) { + # If t_src is not provided set to AOI CRS + if(is.null(t_srs)){ t_srs <- st_crs(AOI) } + # If AOI CRS is NA (e.g st_crs(NULL)) then set to 4326 + if(is.na(t_srs)) { t_srs <- 4326 } + + if(!is.null(AOI) & !is.null(ids)) { + # Check if AOI and IDs are both given + stop("Either IDs or a spatial AOI can be passed.", .call = FALSE) + } else if(is.null(AOI) & is.null(ids) & !(!is.null(where) && grepl("IN", where))) { + # Check if AOI and IDs are both NULL + stop("IDs or a spatial AOI must be passed.", .call = FALSE) + } else if(!(type %in% source$user_call)) { + # Check that "type" is valid + stop(paste("Type not available must be one of:", + paste(source$user_call, collapse = ", ")), + call. = FALSE) + } + + if(!is.null(AOI)){ + + if(length(st_geometry(AOI)) > 1) { + stop("AOI must be one an only one feature.") + } + + if(st_geometry_type(AOI) == "POINT"){ + # If input is a POINT, buffer by 1/2 meter (in equal area projection) + AOI = st_buffer(st_transform(AOI, 5070), buffer) %>% + st_bbox() %>% st_as_sfc() %>% st_make_valid() + } + } + + return(list(AOI = AOI, t_srs = t_srs)) } diff --git a/R/get_hydro.R b/R/get_hydro.R index 704ce821..18b9e418 100644 --- a/R/get_hydro.R +++ b/R/get_hydro.R @@ -1,17 +1,30 @@ -#' @title Find WBD HUC 12 unit subsets -#' @description Subsets the WBD level 12 features by location (POINT), -#' area (POLYGON), or set of IDs. -#' @inherit query_usgs_geoserver details return -#' @inheritParams query_usgs_geoserver +#' @title Find WBD HUC unit subsets +#' @description Subsets WBD features by location (POINT), +#' area (POLYGON), or set of HUC IDs. +#' +#' @inherit query_usgs_geoserver details return params #' @param id WBD HUC ID(s) #' @param type character. Type of feature to return -#' ('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12'). -#' See /link{download_nhdplusv2} for documentation of that dataset. +#' ('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12', 'huc12_nhdplusv2'). +#' +#' Pulls `huc02`-`huc12` from a web service that hosts a snapshot of the +#' Watershed Boundary Dataset from October, 2020. +#' +#' See for full source data. +#' +#' See https://labs.waterdata.usgs.gov/geoserver/web/ for the web service. +#' +#' `huc12_nhdplusv2` derives from a snapshot of the WBD available from the nhdplusv2. +#' See \link{download_nhdplusv2} for source data documentation. +#' +#' +#' #' @export #' get_huc <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, type = "huc12") { - allow_types <- c('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12') + allow_types <- c('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12', + 'huc12_nhdplusv2') if(!type %in% allow_types) { stop("type must be one of ", paste(allow_types, collapse = " ")) @@ -22,42 +35,9 @@ get_huc <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, type = "hu } -#' @title Find WBD HUC 08 unit subsets (DEPRECATED) -#' @description Subsets the WBD level 08 features by location (POINT), -#' area (POLYGON), or set of IDs. -#' @inherit query_usgs_geoserver details return -#' @inheritParams query_usgs_geoserver -#' @param id WBD HUC08 ID(s) -#' @export -get_huc8 <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ - - warning("this function is deprecated -- use get_huc(..., type = \"huc08\") instead") - - query_usgs_geoserver(AOI = AOI, ids = id, type = "huc08_legacy", - t_srs = t_srs, buffer = buffer) -} - -#' @title Find WBD HUC 12 unit subsets (DEPRECATED) -#' @description Subsets the WBD level 12 features by location (POINT), -#' area (POLYGON), or set of IDs. Derived from a static snapshot of -#' HUC 12s from: -#' @inherit query_usgs_geoserver details return -#' @inheritParams query_usgs_geoserver -#' @param id WBD HUC12 ID(s) -#' See /link{download_nhdplusv2} for documentation of that dataset. -#' @export -#' -get_huc12 <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ - - warning("this function is deprecated -- use get_huc(..., type = \"huc12\") instead") - - query_usgs_geoserver(AOI = AOI, ids = id, type = "huc12_nhdplusv2", - t_srs = t_srs, buffer = buffer) -} - -#' @title Find NHD Water Bodies -#' @description Subsets NHD waterbody features by location (POINT), -#' area (POLYGON), or set of IDs. +#' @title Find NHDPlusV2 Water Bodies +#' @description Subsets NHDPlusV2 waterbody features by location (POINT), +#' area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data documentation. #' @inherit query_usgs_geoserver details return #' @inheritParams query_usgs_geoserver #' @param id NHD Waterbody COMID(s) @@ -70,9 +50,9 @@ get_waterbodies <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ buffer = buffer) } -#' @title Find NHD Areas -#' @description Subsets NHD Area features by location (POINT), -#' area (POLYGON), or set of IDs. +#' @title Find NHDPlusV2 Areas +#' @description Subsets NHDPlusV2 Area features by location (POINT), +#' area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data documentation. #' @inherit query_usgs_geoserver details return #' @inheritParams query_usgs_geoserver #' @param id NHD Area COMID(s) @@ -86,7 +66,7 @@ get_nhdarea <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ #' @title Find gagesII Features #' @description Subsets the gagesII dataset by location (POINT), -#' area (POLYGON), or set of IDs. +#' area (POLYGON), or set of IDs. See for documentation of source data. #' @inherit query_usgs_geoserver details return #' @inheritParams query_usgs_geoserver #' @param id character NWIS Gage ID(s) @@ -197,3 +177,77 @@ get_nwis <- function(AOI = NULL, t_srs = NULL, buffer = 20000){ return(st_transform(sites_sf, t_srs)) } } + +#' Get 3DHP Data +#' @description +#' Calls the 3DHP_all web service and returns sf data.frames for the selected +#' layers. See https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer +#' for source data documentation. +#' +#' @inherit query_usgs_arcrest details return params +#' @param ids character vector of id3dhp ids or mainstem uris +#' @param universalreferenceid character vector of hydrolocation universal +#' reference ids such as reachcodes +#' @export +#' @examples +#' \donttest{ +#' AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816, +#' xmax = -89.24681, ymax = 43.17192), +#' crs = "+proj=longlat +datum=WGS84 +no_defs")) +#' +#' # get flowlines and hydrolocations +#' flowlines <- get_3dhp(AOI = AOI, type = "flowline") +#' hydrolocation <- get_3dhp(AOI = AOI, type = "hydrolocation") +#' waterbody <- get_3dhp(AOI = AOI, type = "waterbody") +#' +#' if(!is.null(waterbody) & !is.null(flowlines) & !is.null(hydrolocation)) { +#' plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey") +#' plot(sf::st_geometry(flowlines), col = "blue", add = TRUE) +#' plot(sf::st_geometry(hydrolocation), col = "grey", pch = "+", add = TRUE) } +#' +#' # given mainstem ids from any source, can query for them in ids. +#' +#' CO <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/29559", +#' type = "flowline") +#' +#' if(!is.null(CO)) +#' plot(sf::st_geometry(CO), col = "blue") +#' +#' # get all the waterbodies along the CO river +#' CO_wb <- get_3dhp(ids = unique(CO$waterbodyid3dhp), type = "waterbody") +#' +#' if(!is.null(CO_wb)) { +#' plot(sf::st_geometry(CO_wb[grepl("Powell", CO_wb$gnisidlabel),]), +#' col = "blue", border = "NA") } +#' +#' # given universalreferenceid (reachcodes), can query for them but only +#' # for hydrolocations. This is useful for looking up mainstem ids. +#' +#' get_3dhp(universalreferenceid = unique(hydrolocation$universalreferenceid), +#' type = "hydrolocation") +#'} +get_3dhp <- function(AOI = NULL, ids = NULL, type = NULL, + universalreferenceid = NULL, + t_srs = NULL, buffer = 0.5) { + + if(!is.null(universalreferenceid) & !grepl("outlet|reach|hydrolocation", type)) { + stop("universalereferenceid can only be specified for hydrolocation features") + } + + where <- NULL + if(!is.null(universalreferenceid)) { + where <- paste(paste0("universalreferenceid IN ('", + paste(universalreferenceid, collapse = "', '"), "')")) + if(!is.null(ids)) stop("can not specify both universalreferenceid and other ids") + } + + if(!is.null(ids) && grepl("^https://", ids[1])) { + where <- paste(paste0("mainstemid IN ('", + paste(ids, collapse = "', '"), "')")) + ids <- NULL + } + + query_usgs_arcrest(AOI, ids, type, where, t_srs, buffer) + +} + diff --git a/R/get_nldi.R b/R/get_nldi.R index 325a7bc0..56323a38 100644 --- a/R/get_nldi.R +++ b/R/get_nldi.R @@ -272,7 +272,7 @@ get_nldi_index <- function(location) { #' @importFrom httr GET #' @importFrom jsonlite fromJSON #' @noRd -query_nldi <- function(query, base_path = "/linked-data", parse_json = TRUE) { +query_nldi <- memoise::memoise(function(query, base_path = "/linked-data", parse_json = TRUE) { nldi_base_url <- paste0(get_nldi_url(), base_path) url <- paste(nldi_base_url, query, @@ -303,7 +303,7 @@ query_nldi <- function(query, base_path = "/linked-data", parse_json = TRUE) { warning("Something went wrong accessing the NLDI.\n", e) NULL }) -} +}, ~memoise::timeout(nhdplusTools_memoise_timeout()), cache = nhdplusTools_memoise_cache()) #' @noRd check_nldi_feature <- function(nldi_feature, convert = TRUE) { diff --git a/R/rescale_catchments.R b/R/rescale_catchments.R index 842c462c..1aad331a 100644 --- a/R/rescale_catchments.R +++ b/R/rescale_catchments.R @@ -198,6 +198,10 @@ rescale_catchment_characteristics <- function(vars, lookup_table, catchment_characteristics <- get_catchment_characteristics(varname = vars$characteristic_id, ids = unique(lookup_table$comid)) + if(is.null(catchment_characteristics)) { + return(NULL) + } + } var_names <- unique(catchment_characteristics$characteristic_id) diff --git a/_pkgdown.yml b/_pkgdown.yml index e4ad5a47..1fe168d6 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -37,8 +37,6 @@ reference: - '`get_nhdplus`' - '`get_gagesII`' - '`get_huc`' - - '`get_huc12`' - - '`get_huc8`' - '`get_nhdarea`' - '`get_waterbodies`' - '`get_nwis`' diff --git a/man/get_3dhp.Rd b/man/get_3dhp.Rd new file mode 100644 index 00000000..8548d437 --- /dev/null +++ b/man/get_3dhp.Rd @@ -0,0 +1,91 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_hydro.R +\name{get_3dhp} +\alias{get_3dhp} +\title{Get 3DHP Data} +\usage{ +get_3dhp( + AOI = NULL, + ids = NULL, + type = NULL, + universalreferenceid = NULL, + t_srs = NULL, + buffer = 0.5 +) +} +\arguments{ +\item{AOI}{sf (MULTI)POINT or (MULTI)POLYGON. An 'area of interest' can +be provided as either a location (sf POINT) or area (sf POLYGON) +in any Spatial Reference System.} + +\item{ids}{character vector of id3dhp ids or mainstem uris} + +\item{type}{character. Type of feature to return +("hydrolocation", "flowline", "waterbody", "drainage area", "catchment"). +If NULL (default) a data.frame of available resources is returned} + +\item{universalreferenceid}{character vector of hydrolocation universal +reference ids such as reachcodes} + +\item{t_srs}{character (PROJ string or EPSG code) or numeric (EPSG code). +A user specified - target -Spatial Reference System (SRS/CRS) for returned objects. +Will default to the CRS of the input AOI if provided, and to 4326 for ID requests.} + +\item{buffer}{numeric. The amount (in meters) to buffer a POINT AOI by for an +extended search. Default = 0.5} +} +\value{ +a simple features (sf) object or valid types if no type supplied +} +\description{ +Calls the 3DHP_all web service and returns sf data.frames for the selected +layers. See https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer +for source data documentation. +} +\details{ +The returned object(s) will have the same +Spatial Reference System (SRS) as the input AOI. If a individual or set of +IDs are used to query, then the default CRS of EPSG:4269 is +preserved. In all cases, a user-defined SRS can be passed to \code{t_srs} +which will override all previous SRS (either input or default). +All buffer and distance operations are handled internally using in +EPSG:5070 Albers Equal Area projection +} +\examples{ +\donttest{ +AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816, + xmax = -89.24681, ymax = 43.17192), + crs = "+proj=longlat +datum=WGS84 +no_defs")) + +# get flowlines and hydrolocations +flowlines <- get_3dhp(AOI = AOI, type = "flowline") +hydrolocation <- get_3dhp(AOI = AOI, type = "hydrolocation") +waterbody <- get_3dhp(AOI = AOI, type = "waterbody") + +if(!is.null(waterbody) & !is.null(flowlines) & !is.null(hydrolocation)) { +plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey") +plot(sf::st_geometry(flowlines), col = "blue", add = TRUE) +plot(sf::st_geometry(hydrolocation), col = "grey", pch = "+", add = TRUE) } + +# given mainstem ids from any source, can query for them in ids. + +CO <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/29559", + type = "flowline") + +if(!is.null(CO)) + plot(sf::st_geometry(CO), col = "blue") + +# get all the waterbodies along the CO river +CO_wb <- get_3dhp(ids = unique(CO$waterbodyid3dhp), type = "waterbody") + +if(!is.null(CO_wb)) { +plot(sf::st_geometry(CO_wb[grepl("Powell", CO_wb$gnisidlabel),]), + col = "blue", border = "NA") } + +# given universalreferenceid (reachcodes), can query for them but only +# for hydrolocations. This is useful for looking up mainstem ids. + +get_3dhp(universalreferenceid = unique(hydrolocation$universalreferenceid), + type = "hydrolocation") +} +} diff --git a/man/get_gagesII.Rd b/man/get_gagesII.Rd index 23fd1f6e..8e8b27a0 100644 --- a/man/get_gagesII.Rd +++ b/man/get_gagesII.Rd @@ -28,7 +28,7 @@ a simple features (sf) object } \description{ Subsets the gagesII dataset by location (POINT), -area (POLYGON), or set of IDs. +area (POLYGON), or set of IDs. See for documentation of source data. } \details{ The returned object(s) will have the same diff --git a/man/get_huc.Rd b/man/get_huc.Rd index 6733d11a..829ba9a0 100644 --- a/man/get_huc.Rd +++ b/man/get_huc.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/get_hydro.R \name{get_huc} \alias{get_huc} -\title{Find WBD HUC 12 unit subsets} +\title{Find WBD HUC unit subsets} \usage{ get_huc(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5, type = "huc12") } @@ -21,15 +21,24 @@ Will default to the CRS of the input AOI if provided, and to 4326 for ID request extended search. Default = 0.5} \item{type}{character. Type of feature to return -('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12'). -See /link{download_nhdplusv2} for documentation of that dataset.} +('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12', 'huc12_nhdplusv2'). + +Pulls `huc02`-`huc12` from a web service that hosts a snapshot of the +Watershed Boundary Dataset from October, 2020. + +See for full source data. + +See https://labs.waterdata.usgs.gov/geoserver/web/ for the web service. + +`huc12_nhdplusv2` derives from a snapshot of the WBD available from the nhdplusv2. +See \link{download_nhdplusv2} for source data documentation.} } \value{ a simple features (sf) object } \description{ -Subsets the WBD level 12 features by location (POINT), -area (POLYGON), or set of IDs. +Subsets WBD features by location (POINT), +area (POLYGON), or set of HUC IDs. } \details{ The returned object(s) will have the same diff --git a/man/get_huc8.Rd b/man/get_huc8.Rd deleted file mode 100644 index 5de67d1c..00000000 --- a/man/get_huc8.Rd +++ /dev/null @@ -1,38 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/get_hydro.R -\name{get_huc8} -\alias{get_huc8} -\title{Find WBD HUC 08 unit subsets (DEPRECATED)} -\usage{ -get_huc8(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5) -} -\arguments{ -\item{AOI}{sf (MULTI)POINT or (MULTI)POLYGON. An 'area of interest' can -be provided as either a location (sf POINT) or area (sf POLYGON) -in any Spatial Reference System.} - -\item{id}{WBD HUC08 ID(s)} - -\item{t_srs}{character (PROJ string or EPSG code) or numeric (EPSG code). -A user specified - target -Spatial Reference System (SRS/CRS) for returned objects. -Will default to the CRS of the input AOI if provided, and to 4326 for ID requests.} - -\item{buffer}{numeric. The amount (in meters) to buffer a POINT AOI by for an -extended search. Default = 0.5} -} -\value{ -a simple features (sf) object -} -\description{ -Subsets the WBD level 08 features by location (POINT), -area (POLYGON), or set of IDs. -} -\details{ -The returned object(s) will have the same -Spatial Reference System (SRS) as the input AOI. If a individual or set of -IDs are used to query, then the default geoserver CRS of EPSG:4326 is -preserved. In all cases, a user-defined SRS can be passed to \code{t_srs} -which will override all previous SRS (either input or default). -All buffer and distance operations are handled internally using in -EPSG:5070 Albers Equal Area projection -} diff --git a/man/get_nhdarea.Rd b/man/get_nhdarea.Rd index 6c4a5936..f0e67600 100644 --- a/man/get_nhdarea.Rd +++ b/man/get_nhdarea.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/get_hydro.R \name{get_nhdarea} \alias{get_nhdarea} -\title{Find NHD Areas} +\title{Find NHDPlusV2 Areas} \usage{ get_nhdarea(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5) } @@ -24,8 +24,8 @@ extended search. Default = 0.5} a simple features (sf) object } \description{ -Subsets NHD Area features by location (POINT), -area (POLYGON), or set of IDs. +Subsets NHDPlusV2 Area features by location (POINT), +area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data documentation. } \details{ The returned object(s) will have the same diff --git a/man/get_waterbodies.Rd b/man/get_waterbodies.Rd index bc6b2d34..b72ef34e 100644 --- a/man/get_waterbodies.Rd +++ b/man/get_waterbodies.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/get_hydro.R \name{get_waterbodies} \alias{get_waterbodies} -\title{Find NHD Water Bodies} +\title{Find NHDPlusV2 Water Bodies} \usage{ get_waterbodies(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5) } @@ -24,8 +24,8 @@ extended search. Default = 0.5} a simple features (sf) object } \description{ -Subsets NHD waterbody features by location (POINT), -area (POLYGON), or set of IDs. +Subsets NHDPlusV2 waterbody features by location (POINT), +area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data documentation. } \details{ The returned object(s) will have the same diff --git a/man/nhdplusTools_cache_settings.Rd b/man/nhdplusTools_cache_settings.Rd new file mode 100644 index 00000000..22597a96 --- /dev/null +++ b/man/nhdplusTools_cache_settings.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/A_nhdplusTools.R +\name{nhdplusTools_cache_settings} +\alias{nhdplusTools_cache_settings} +\title{nhdplusTools cache settings} +\usage{ +nhdplusTools_cache_settings(mode = NULL, timeout = NULL) +} +\arguments{ +\item{mode}{character 'memory' or 'filesystem'} + +\item{timeout}{numeric number of seconds until caches invalidate} +} +\value{ +list containing settings at time of calling. If inputs are +NULL, current settings. If settings are altered, previous setting values. +} +\description{ +Provides an interface to adjust nhdplusTools `memoise` cache. + +Mode and timeout can also be set using environment variables. +`NHDPLUSTOOLS_MEMOISE_CACHE` and `NHDPLUSTOOLS_MEMOISE_TIMEOUT` are +used unless overriden with this function. +} diff --git a/man/get_huc12.Rd b/man/query_usgs_arcrest.Rd similarity index 54% rename from man/get_huc12.Rd rename to man/query_usgs_arcrest.Rd index 306990d0..f8d26cc9 100644 --- a/man/get_huc12.Rd +++ b/man/query_usgs_arcrest.Rd @@ -1,18 +1,31 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/get_hydro.R -\name{get_huc12} -\alias{get_huc12} -\title{Find WBD HUC 12 unit subsets (DEPRECATED)} +% Please edit documentation in R/arcrest_tools.R +\name{query_usgs_arcrest} +\alias{query_usgs_arcrest} +\title{Query USGS Hydro ESRI Rest Server} \usage{ -get_huc12(AOI = NULL, id = NULL, t_srs = NULL, buffer = 0.5) +query_usgs_arcrest( + AOI = NULL, + ids = NULL, + type = NULL, + where = NULL, + t_srs = NULL, + buffer = 0.5 +) } \arguments{ \item{AOI}{sf (MULTI)POINT or (MULTI)POLYGON. An 'area of interest' can be provided as either a location (sf POINT) or area (sf POLYGON) in any Spatial Reference System.} -\item{id}{WBD HUC12 ID(s) -See /link{download_nhdplusv2} for documentation of that dataset.} +\item{ids}{character. A set of identifier(s) from the data +type requested, for 3dhp, this is id3dhp.} + +\item{type}{character. Type of feature to return +("hydrolocation", "flowline", "waterbody", "drainage area", "catchment"). +If NULL (default) a data.frame of available resources is returned} + +\item{where}{character. An where clause to pass to the server.} \item{t_srs}{character (PROJ string or EPSG code) or numeric (EPSG code). A user specified - target -Spatial Reference System (SRS/CRS) for returned objects. @@ -22,19 +35,19 @@ Will default to the CRS of the input AOI if provided, and to 4326 for ID request extended search. Default = 0.5} } \value{ -a simple features (sf) object +a simple features (sf) object or valid types if no type supplied } \description{ -Subsets the WBD level 12 features by location (POINT), -area (POLYGON), or set of IDs. Derived from a static snapshot of -HUC 12s from: +Query the USGS Hydro ESRI Rest Server for spatial data by location, +area, or ID. } \details{ The returned object(s) will have the same Spatial Reference System (SRS) as the input AOI. If a individual or set of -IDs are used to query, then the default geoserver CRS of EPSG:4326 is +IDs are used to query, then the default CRS of EPSG:4269 is preserved. In all cases, a user-defined SRS can be passed to \code{t_srs} which will override all previous SRS (either input or default). All buffer and distance operations are handled internally using in EPSG:5070 Albers Equal Area projection } +\keyword{internal} diff --git a/tests/testthat/helper.R b/tests/testthat/helper.R index eec009b7..8434ee27 100644 --- a/tests/testthat/helper.R +++ b/tests/testthat/helper.R @@ -3,6 +3,8 @@ options("rgdal_show_exportToProj4_warnings"="none") library("sf") library("dplyr") +nhdplusTools_cache_settings(mode = "memory", timeout = 1) + sf::sf_use_s2(TRUE) unlink(file.path(tempdir(check = TRUE), "*"), recursive = TRUE) diff --git a/tests/testthat/test_03_get_functions.R b/tests/testthat/test_03_get_functions.R index 00423fb0..8c49e642 100644 --- a/tests/testthat/test_03_get_functions.R +++ b/tests/testthat/test_03_get_functions.R @@ -47,28 +47,20 @@ test_that("query water labs...",{ test_that("huc8", { testthat::skip_on_cran() #Point - expect_warning({ - ptHUC8 = get_huc8(AOI = pt) - }) + ptHUC8 = get_huc(AOI = pt, type = "huc08") expect_equal(nrow(ptHUC8), 1) expect_equal(ptHUC8$huc8, "17010101") expect_equal(sf::st_crs(ptHUC8)$epsg, 4326) - expect_warning({ - expect_error(get_huc8(AOI = rbind(pt,pt2)), + expect_error(get_huc(AOI = rbind(pt, pt2), type = "huc08"), "AOI must be one an only one feature.") - }) #Area - expect_warning({ - areaHUC8 = get_huc8(AOI = area, t_srs = 5070) - }) + areaHUC8 = get_huc(AOI = area, t_srs = 5070, type = "huc08") expect_equal(sf::st_crs(areaHUC8)$epsg, 5070) expect_equal(nrow(areaHUC8), 1) #ID - expect_warning({ - ptHUC8id = get_huc8(id = "17010101") - }) + ptHUC8id = get_huc(id = "17010101", type = "huc08") expect_identical(ptHUC8$huc8, ptHUC8id$huc8) expect_true(sf::st_crs(ptHUC8) == sf::st_crs(ptHUC8id) ) }) @@ -77,27 +69,22 @@ test_that("huc8", { test_that("huc", { testthat::skip_on_cran() + + expect_error(get_huc(AOI = pt, type = "borked"), + "type must be one of huc02 huc04 huc06 huc08 huc10 huc12 huc12_nhdplusv2") #Point - expect_warning({ - ptHUC12 = get_huc12(AOI = pt) - }) + ptHUC12 = get_huc(AOI = pt, type = "huc12_nhdplusv2") expect_equal(nrow(ptHUC12), 1) expect_equal(ptHUC12$huc12, "170101010306") #Area - expect_warning({ - areaHUC12 = get_huc12(AOI = area) + areaHUC12 = get_huc(AOI = area, type = "huc12_nhdplusv2") expect_equal(nrow(areaHUC12), 2) - }) #ID - expect_warning({ - HUC12id = get_huc12(id = "170101010306") - }) + HUC12id = get_huc(id = "170101010306", type = "huc12_nhdplusv2") expect_identical(ptHUC12$huc12, HUC12id$huc12) # multi-id... only need to check once - expect_warning({ - HUC12id2 = get_huc12(id = areaHUC12$huc12) %>% + HUC12id2 = get_huc(id = areaHUC12$huc12, type = "huc12_nhdplusv2") %>% sf::st_transform(sf::st_crs(area)) - }) expect_identical(HUC12id2$geometry, areaHUC12$geometry) diff --git a/tests/testthat/test_arcrest.R b/tests/testthat/test_arcrest.R new file mode 100644 index 00000000..a078950f --- /dev/null +++ b/tests/testthat/test_arcrest.R @@ -0,0 +1,72 @@ +AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816, + xmax = -89.24681, ymax = 43.17192), + crs = "+proj=longlat +datum=WGS84 +no_defs")) + +test_that("spatial_filter", { + + expect_error(nhdplusTools:::spatial_filter(AOI, format = "test"), + "ogc or esri") + + expect_equal( + nhdplusTools:::spatial_filter(AOI, format = "ogc")[[1]], + "the_geom42.99816 -89.5668443.17192 -89.24681") + + expect_equal(as.character(nhdplusTools:::spatial_filter(AOI, format = "esri")[[1]]), + '{"xmin":-89.5668,"ymin":42.9982,"xmax":-89.2468,"ymax":43.1719,"spatialReference":{"wkid":4326}}') +}) + +test_that("check_query_params", { + + ids <- NULL + type <- "test" + source <- data.frame(user_call = c("test", "test2")) + t_srs <- NULL + buffer <- 0.5 + where <- NULL + + expect_equal( + nhdplusTools:::check_query_params(AOI, ids, type, where, source, t_srs, buffer), + list(AOI = AOI, t_srs = sf::st_crs(AOI))) + + expect_error(nhdplusTools:::check_query_params(AOI, c(1,2), type, where, source, t_srs, buffer), + "Either IDs or") + + expect_error(nhdplusTools:::check_query_params(NULL, NULL, type, where, source, t_srs, buffer), + "IDs or a spatial AOI must be passed.") + + expect_error(nhdplusTools:::check_query_params(AOI, ids, "test3", where, source, t_srs, buffer), + "test, test2") + + expect_error(nhdplusTools:::check_query_params(c(AOI, AOI), ids, type, where, source, t_srs, buffer), + "AOI must be one an only one feature.") + + AOI <- sf::st_sfc(sf::st_point(c(-89.56684, 42.99816)), crs = 4326) + + expect_equal( + nhdplusTools:::check_query_params(AOI, ids, type, where, source, t_srs, buffer)$AOI[[1]], + structure(list(structure(c(521279.507824339, 521279.507824339, + 521280.507824339, 521280.507824339, + 521279.507824339, 2240146.81449532, + 2240147.81449532, 2240147.81449532, + 2240146.81449532, 2240146.81449532), + dim = c(5L, 2L))), + class = c("XY", "POLYGON", "sfg")), tolerance = 0.1) + +}) + +test_that("basic 3dhp service requests", { + skip_on_cran() + + AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.5, ymin = 43.0, + xmax = -89.4, ymax = 43.1), + crs = "+proj=longlat +datum=WGS84 +no_defs")) + + expect_message(expect_s3_class(nhdplusTools:::query_usgs_arcrest(AOI), + "data.frame")) + + expect_warning(expect_warning(nhdplusTools:::query_usgs_arcrest(AOI, type = "hydrolocation"))) + + test_data <- nhdplusTools:::query_usgs_arcrest(AOI, type = "reach code, external connection") + + expect_s3_class(test_data, "sf") +}) diff --git a/tests/testthat/test_get_3dhp.R b/tests/testthat/test_get_3dhp.R new file mode 100644 index 00000000..b8cfb5c5 --- /dev/null +++ b/tests/testthat/test_get_3dhp.R @@ -0,0 +1,24 @@ +test_that("get_3dhp", { + + expect_error(get_3dhp(universalreferenceid = "01234", type = "test"), + "only be specified for hydrolocation") + + expect_error(get_3dhp(universalreferenceid = "01234", type = "hydrolocation", + ids = "1"), + "can not specify both") + + skip_on_cran() + + ms <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/377002", + type = "flowline") + + expect_equal(unique(ms$mainstemid), "https://geoconnex.us/ref/mainstems/377002") + + expect_s3_class(ms, "sf") + + suppressWarnings({ + hl <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/377002", + type = "hydrolocation") + }) + expect_equal(unique(hl$mainstemid), "https://geoconnex.us/ref/mainstems/377002") +}) diff --git a/tests/testthat/test_get_vaa.R b/tests/testthat/test_get_vaa.R index 69efb95b..e38bd1af 100644 --- a/tests/testthat/test_get_vaa.R +++ b/tests/testthat/test_get_vaa.R @@ -2,7 +2,8 @@ test_that("vaa examples", { skip_on_cran() - testthat::skip_on_os("linux") + skip_on_os("linux") + skip_on_os("mac") vaa_names <- get_vaa_names() @@ -36,7 +37,8 @@ test_that("vaa examples", { test_that("catchment chars", { skip_on_cran() - testthat::skip_on_os("linux") + skip_on_os("linux") + skip_on_os("mac") httptest::without_internet({ suppressMessages(expect_warning(w <- get_characteristics_metadata(cache = FALSE))) diff --git a/tests/testthat/test_rescale_catchments.R b/tests/testthat/test_rescale_catchments.R index 483f094b..bf3d7db3 100644 --- a/tests/testthat/test_rescale_catchments.R +++ b/tests/testthat/test_rescale_catchments.R @@ -1,6 +1,7 @@ test_that("rescale", { skip_on_cran() + skip_on_os("mac") vars <- data.frame(characteristic_id = c("CAT_EWT", "CAT_EWT", "CAT_EWT", "CAT_EWT", "CAT_BASIN_AREA"), summary_statistic = c("area_weighted_mean", "min", "sum", "max", "sum")) diff --git a/vignettes/get_3dhp_data.Rmd b/vignettes/get_3dhp_data.Rmd new file mode 100644 index 00000000..d6a05065 --- /dev/null +++ b/vignettes/get_3dhp_data.Rmd @@ -0,0 +1,114 @@ +--- +title: "3DHP Data Access" +author: "dblodgett@usgs.gov" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{3DHP Data Access} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + + +```{r setup, include = FALSE} +library(nhdplusTools) + +local <- (Sys.getenv("BUILD_VIGNETTES") == "TRUE") +if(local) { + cache_path <- file.path(nhdplusTools_data_dir(), "3dhp_v_cache") +} else { + cache_path <- tempdir() +} + +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width=6, + fig.height=4, + eval=local, + cache=local, + cache.path=(cache_path) +) + +oldoption <- options(scipen = 9999) + +``` + +This vignette shows how to work with a [new service based on data from the 3D Hydrography Program (3DHP)](https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer) released in 2024. 3DHP data uses "mainstem" identifiers as the primary river ID. We'll work with one for the Baraboo River in Wisconsin. + +We can find one to start with at a url like: + +`https://reference.geoconnex.us/collections/mainstems/items?filter=name_at_outlet ILIKE '%Baraboo%'` + +This url searches the reference mainstem collection for rivers with names like "Baraboo". Using that url query, we can find that the Baraboo mainstem ID is `https://geoconnex.us/ref/mainstems/359842` and that it flows to the Wisconsin and Mississippi `https://geoconnex.us/ref/mainstems/323742` and `https://geoconnex.us/ref/mainstems/312091` respectively. At the reference mainstem page, we can also find that the outlet NHDPlusV2 COMID is `https://geoconnex.us/nhdplusv2/comid/937070225`. + +The code block just below generates starting locations based on information derived from the reference mainstem for the Baraboo and a point location near the Baraboo's outlet. It queries the 3DHP service for a flowline within a 10m buffer of the point specified just to get us started. + +```{r} + + comid <- "937070225" + + point <- c(-89.441, 43.487) |> + sf::st_point() |> + sf::st_sfc(crs = 4326) |> sf::st_sf() + + dm <- c('https://geoconnex.us/ref/mainstems/323742', + 'https://geoconnex.us/ref/mainstems/312091') + + flowline <- get_3dhp(point, type = "flowline", buffer = 10) + + mapview::mapview(list(point, flowline)) +``` + +Now, we can use [dataRetrieval](https://cran.r-project.org/package=dataRetrieval) to retrieve a basin upstream of this outlet location. We specify the NHDPlusV2 comid here since we know it, but could have used the point location just as well. + +With the basin pulled from the [`findNLDI()`](https://doi-usgs.github.io/dataRetrieval/reference/findNLDI.html) function, we can pass it to `get_3dhp()` as a Area of Interest for flowlines, waterbodies, and hydrolocations. We can then use the mainstem ids of our two downstream mainstem rivers found above to get the downtream mainstem rivers. + +```{r} + basin <- dataRetrieval::findNLDI(comid = comid, find = "basin") + + mapview::mapview(basin$basin) + + network <- get_3dhp(basin$basin, type = "flowline") + water <- get_3dhp(basin$basin, type = "waterbody") + hydrolocation <- get_3dhp(basin$basin, type = "hydrolocation") + + down_mains <- get_3dhp(ids = dm, type = "flowline") + + map_data <- list(outlet = basin$origin, + basin = basin$basin, + flowline = network, + waterbody = water, + hydrolocation = hydrolocation, + down_mains = down_mains) + + mapview::mapview(map_data) +``` + +Neat. Now we have flowlines, waterbodies, hydrologic locations, a basin boundary, and major rivers downstream to work with. These could be saved to a local file for use in a GIS or used in some further analysis. + +Say you want to know how to relate NHD data to 3DHP data. 3DHP hydrolocations include "reachcode" tops and bottoms to hel with that. Here, we lookup the hydrolocation with reachcode `07070004002889`then grab the rest of the reachcode hydrolocations along the mainstem that the reachcode is part of. + +```{r} + + reachcode <- "07070004002889" + + hydrolocation <- get_3dhp(universalreferenceid = reachcode, + type = "reach code, external connection") + + hydrolocation + + mainstem_points <- get_3dhp(ids = hydrolocation$mainstemid, type = "reach code, external connection") + mainstem_lines <- get_3dhp(ids = hydrolocation$mainstemid, type = "flowline") + + mapview::mapview(list(mainstem_points, mainstem_lines, hydrolocation)) +``` +With this functionality, we have what we need to inter operate between older NHD data and newly released 3DHP data. This service is brand new and `nhdplusTools` is just getting going with functionality to support what's shown here. Please [submit issues](https://github.com/DOI-USGS/nhdplusTools/issues) if you find problems or have questions!! + + +```{r teardown, include=FALSE} +options(oldoption) + +if(Sys.getenv("BUILD_VIGNETTES") != "TRUE") { + unlink(work_dir, recursive = TRUE) +} +```