From a88bd880a5aad480de0741fab68a4d2d7854ddb0 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Mon, 8 Sep 2025 14:50:30 -0500 Subject: [PATCH 01/10] reorganize to prepare for geoserver removal #436 --- R/A_nhdplusTools.R | 28 ++++++++++++ R/downloading_tools.R | 64 ++++++++++++++++++++++++++++ R/geoserver_tools.R | 99 ------------------------------------------- R/get_geoconnex.R | 85 ------------------------------------- R/get_nhdplus.R | 19 +++++++++ R/oafeat_tools.R | 74 ++++++++++++++++++++++++++++++++ 6 files changed, 185 insertions(+), 184 deletions(-) create mode 100644 R/oafeat_tools.R diff --git a/R/A_nhdplusTools.R b/R/A_nhdplusTools.R index eeb695cd..52adf6dd 100644 --- a/R/A_nhdplusTools.R +++ b/R/A_nhdplusTools.R @@ -506,6 +506,34 @@ align_nhdplus_names <- function(x){ } +filter_list_kvp <- \(l, key, val, type = NULL, n = NULL) { + ret <- l[vapply(l, \(x) x[[key]] == val, TRUE)] + + + if(!is.null(type)) { + ret <- ret[vapply(ret, \(x) x[["type"]] == type, TRUE)] + } + + if(!is.null(n)) { + ret <- ret[[n]] + } + + ret +} + +extract <- `[[` + +#' @title Trim and Cull NULLs +#' @description Remove NULL arguments from a list +#' @param x a list +#' @keywords internal +#' @return a list +#' @noRd + +tc <- function(x) { + Filter(Negate(is.null), x) +} + #' @importFrom hydroloom st_compatibalize #' @export hydroloom::st_compatibalize diff --git a/R/downloading_tools.R b/R/downloading_tools.R index fb04d9a6..6e258819 100644 --- a/R/downloading_tools.R +++ b/R/downloading_tools.R @@ -1,3 +1,7 @@ +##################################################################### +# File contains general downloading tools and API utility functions # +##################################################################### + #' Download NHDPlus HiRes #' @param nhd_dir character directory to save output into #' @param hu_list character vector of hydrologic region(s) to download. @@ -342,3 +346,63 @@ check7z <- function() { } +#' memoise get json +#' @description +#' attempts to get a url as JSON and return the content. +#' +#' Will return NULL if anything fails +#' +#' @param url character url to get +#' @return list containing parsed json on success, NULL otherwise +#' @noRd +mem_get_json <- memoise::memoise(\(url) { + tryCatch({ + retn <- httr::RETRY("GET", url, httr::accept_json()) + + if(retn$status_code == 200 & grepl("json", retn$headers$`content-type`)) { + return(httr::content(retn, simplifyVector = FALSE, type = "application/json")) + } else { + warning("Can't access json from ", url) + return(NULL) + } + }, error = function(e) { + warning("Error accessing ", url, "\n\n", e) + return(NULL) + }) +}) + +#' @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/geoserver_tools.R b/R/geoserver_tools.R index 74f48330..1cf85176 100644 --- a/R/geoserver_tools.R +++ b/R/geoserver_tools.R @@ -179,38 +179,6 @@ query_usgs_geoserver <- function(AOI = NULL, ids = NULL, } -unify_types <- function(out) { - all_class <- bind_rows(lapply(out, function(x) { - vapply(x, function(x) class(x)[1], character(1)) - })) - - set_type <- function(out, n, type) { - lapply(out, function(x) { - x[[n]] <- as(x[[n]], type) - x - }) - } - - for(n in names(all_class)) { - if(length(unique(all_class[[n]])) > 1) { - if("numeric" %in% all_class[[n]]) { # prefer numeric - out <- set_type(out, n, "numeric") - } else if("integer" %in% all_class[[n]]) { # then integer - out <- set_type(out, n, "integer") - } else if("cheracter" %in% all_class[[n]]) { - out <- set_type(out, n, "character") - } - } - } - - rows <- sapply(out, nrow) - - if(any(rows > 0)) - out <- out[rows > 0] - - out -} - assign("bb_break_size", value = 2, nhdplusTools_env) #' @title Construct a BBOX spatial filter for geoservers @@ -337,70 +305,3 @@ streamorder_filter <- function(streamorder){ '', streamorder - 1, '', '') } - -#' @title Trim and Cull NULLs -#' @description Remove NULL arguments from a list -#' @param x a list -#' @keywords internal -#' @return a list -#' @noRd - -tc <- function(x) { - Filter(Negate(is.null), x) -} - -#' @title Identify NHD features by collocated NWIS ID(s) -#' @description Use the NLDI to identify the COMIDs associated -#' with a given NWIS ID. -#' @param nwis character or numeric. A vector of USGS NWIS id(s) -#' @keywords internal -#' @return a vector of COMIDs -#' @noRd -#' @importFrom httr RETRY GET -#' @importFrom jsonlite fromJSON - -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/") - url <- paste0(baseURL, "nwissite/USGS-", nwis) - c <- rawToChar(httr::RETRY("GET", url)$content) - f.comid <- jsonlite::fromJSON(c, simplifyVector = TRUE) - f.comid$features$properties$comid -}) - -#' @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_geoconnex.R b/R/get_geoconnex.R index 979049f7..fc199138 100644 --- a/R/get_geoconnex.R +++ b/R/get_geoconnex.R @@ -1,45 +1,3 @@ -#' memoise get json -#' @description -#' attempts to get a url as JSON and return the content. -#' -#' Will return NULL if anything fails -#' -#' @param url character url to get -#' @return list containing parsed json on success, NULL otherwise -#' @noRd -mem_get_json <- memoise::memoise(\(url) { - tryCatch({ - retn <- httr::RETRY("GET", url, httr::accept_json()) - - if(retn$status_code == 200 & grepl("json", retn$headers$`content-type`)) { - return(httr::content(retn, simplifyVector = FALSE, type = "application/json")) - } else { - warning("Can't access json from ", url) - return(NULL) - } - }, error = function(e) { - warning("Error accessing ", url, "\n\n", e) - return(NULL) - }) -}) - -filter_list_kvp <- \(l, key, val, type = NULL, n = NULL) { - ret <- l[vapply(l, \(x) x[[key]] == val, TRUE)] - - - if(!is.null(type)) { - ret <- ret[vapply(ret, \(x) x[["type"]] == type, TRUE)] - } - - if(!is.null(n)) { - ret <- ret[[n]] - } - - ret -} - -extract <- `[[` - #' discover geoconnex reference feature layers #' @description #' Queries the geoconnex.us reference feature server for available layers and @@ -84,49 +42,6 @@ discover_geoconnex_reference <- function() { dplyr::left_join(collections_meta, q_ables, by = "id") } -get_features_paging <- function(base_call, limit = 1000, status = TRUE) { - - if(!grepl("\\?", base_call)) { - base_call <- paste0(base_call, "?") - } else { - base_call <- paste0(base_call, "&") - } - - offset <- 0 - - keep_going <- TRUE - - if(status) message("Starting download of first set of features.") - - out <- rep(list(list()), 1e6) - i <- 1 - - while(keep_going) { - req <- paste0(base_call, "limit=", limit, "&offset=", offset) - - out[[i]] <- try(read_sf(req)) - - if(inherits(out[[i]], "sf") & nrow(out[[i]]) == limit) { - offset <- offset + limit - } - - if(nrow(out[[i]]) < limit) keep_going <- FALSE - - if(!inherits(out[[i]], "sf")) { - warning("Something went wrong requesting data.") - keep_going <- FALSE - } - - if(status & keep_going) message("Starting next download from ", offset, ".") - - i <- i + 1 - } - - out <- out[1:(i - 1)] - - sf::st_sf(dplyr::bind_rows(unify_types(out))) -} - #' get geoconnex reference feature layers #' @description #' Queries the geoconnex reference feature server for features of interest. diff --git a/R/get_nhdplus.R b/R/get_nhdplus.R index cee5efcb..efb7e1e2 100644 --- a/R/get_nhdplus.R +++ b/R/get_nhdplus.R @@ -117,3 +117,22 @@ get_nhdplus <- function(AOI = NULL, return(geoms) } +#' @title Identify NHD features by collocated NWIS ID(s) +#' @description Use the NLDI to identify the COMIDs associated +#' with a given NWIS ID. +#' @param nwis character or numeric. A vector of USGS NWIS id(s) +#' @keywords internal +#' @return a vector of COMIDs +#' @noRd +#' @importFrom httr RETRY GET +#' @importFrom jsonlite fromJSON + +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/") + url <- paste0(baseURL, "nwissite/USGS-", nwis) + c <- rawToChar(httr::RETRY("GET", url)$content) + f.comid <- jsonlite::fromJSON(c, simplifyVector = TRUE) + f.comid$features$properties$comid +}) diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R new file mode 100644 index 00000000..8b4342b4 --- /dev/null +++ b/R/oafeat_tools.R @@ -0,0 +1,74 @@ +get_features_paging <- function(base_call, limit = 1000, status = TRUE) { + + if(!grepl("\\?", base_call)) { + base_call <- paste0(base_call, "?") + } else { + base_call <- paste0(base_call, "&") + } + + offset <- 0 + + keep_going <- TRUE + + if(status) message("Starting download of first set of features.") + + out <- rep(list(list()), 1e6) + i <- 1 + + while(keep_going) { + req <- paste0(base_call, "limit=", limit, "&offset=", offset) + + out[[i]] <- try(read_sf(req)) + + if(inherits(out[[i]], "sf") & nrow(out[[i]]) == limit) { + offset <- offset + limit + } + + if(nrow(out[[i]]) < limit) keep_going <- FALSE + + if(!inherits(out[[i]], "sf")) { + warning("Something went wrong requesting data.") + keep_going <- FALSE + } + + if(status & keep_going) message("Starting next download from ", offset, ".") + + i <- i + 1 + } + + out <- out[1:(i - 1)] + + sf::st_sf(dplyr::bind_rows(unify_types(out))) +} + +unify_types <- function(out) { + all_class <- bind_rows(lapply(out, function(x) { + vapply(x, function(x) class(x)[1], character(1)) + })) + + set_type <- function(out, n, type) { + lapply(out, function(x) { + x[[n]] <- as(x[[n]], type) + x + }) + } + + for(n in names(all_class)) { + if(length(unique(all_class[[n]])) > 1) { + if("numeric" %in% all_class[[n]]) { # prefer numeric + out <- set_type(out, n, "numeric") + } else if("integer" %in% all_class[[n]]) { # then integer + out <- set_type(out, n, "integer") + } else if("cheracter" %in% all_class[[n]]) { + out <- set_type(out, n, "character") + } + } + } + + rows <- sapply(out, nrow) + + if(any(rows > 0)) + out <- out[rows > 0] + + out +} From 6bc384815978c8bc066981e5c16f935eeaccca45 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Mon, 8 Sep 2025 15:06:50 -0500 Subject: [PATCH 02/10] discover_oafeat() internal function. #436 --- R/get_geoconnex.R | 28 +--------------------------- R/oafeat_tools.R | 41 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 27 deletions(-) diff --git a/R/get_geoconnex.R b/R/get_geoconnex.R index fc199138..9352300f 100644 --- a/R/get_geoconnex.R +++ b/R/get_geoconnex.R @@ -12,34 +12,8 @@ #' discover_geoconnex_reference <- function() { - landing <- mem_get_json(get("gocnx_ref_base_url", envir = nhdplusTools_env)) + discover_oafeat(get("gocnx_ref_base_url", envir = nhdplusTools_env)) - collections <- landing$links |> - filter_list_kvp("rel", "data", n = 1) |> - extract("href") |> - mem_get_json() - - collections_meta <- dplyr::bind_rows( - lapply(collections$collections, - \(x) c(x[c("id", "title", "description")], - list(url = filter_list_kvp(x$links, - "rel", "self", n = 1)$href)))) - - - q_ables <- dplyr::bind_rows(lapply(collections$collections, \(x) { - q <- filter_list_kvp(x$links, "rel", "http://www.opengis.net/def/rel/ogc/1.0/queryables", - type = "application/schema+json", n = 1)$href |> - mem_get_json() |> - (\(y) list(id = x$id, qs = y$properties))() - - q$qs <- q$qs[vapply(q$qs, \(x) all(c("title", "type") %in% names(x)), TRUE)] - - data.frame(id = rep(q$id, length(q$qs)), - attribute = vapply(q$qs, \(x) x$title, ""), - type = vapply(q$qs, \(x) x$type, ""), row.names = NULL) - })) - - dplyr::left_join(collections_meta, q_ables, by = "id") } #' get geoconnex reference feature layers diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index 8b4342b4..476f6910 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -1,3 +1,44 @@ +#' discover OGC API feature layers +#' @description +#' Queries an OGC API feature server for available layers and +#' attributes. +#' +#' @return data.frame containing layers available and fields that are available to query. +#' @param landing_url url of landing page of OGC API +#' @noRd +#' +discover_oafeat <- function(landing_url) { + + landing <- mem_get_json(landing_url) + + collections <- landing$links |> + filter_list_kvp("rel", "data", n = 1) |> + extract("href") |> + mem_get_json() + + collections_meta <- dplyr::bind_rows( + lapply(collections$collections, + \(x) c(x[c("id", "title", "description")], + list(url = filter_list_kvp(x$links, + "rel", "self", n = 1)$href)))) + + + q_ables <- dplyr::bind_rows(lapply(collections$collections, \(x) { + q <- filter_list_kvp(x$links, "rel", "http://www.opengis.net/def/rel/ogc/1.0/queryables", + type = "application/schema+json", n = 1)$href |> + mem_get_json() |> + (\(y) list(id = x$id, qs = y$properties))() + + q$qs <- q$qs[vapply(q$qs, \(x) all(c("title", "type") %in% names(x)), TRUE)] + + data.frame(id = rep(q$id, length(q$qs)), + attribute = vapply(q$qs, \(x) x$title, ""), + type = vapply(q$qs, \(x) x$type, ""), row.names = NULL) + })) + + dplyr::left_join(collections_meta, q_ables, by = "id") +} + get_features_paging <- function(base_call, limit = 1000, status = TRUE) { if(!grepl("\\?", base_call)) { From 55edeacf7966607b4dae73104892e73eceb3b5a9 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Mon, 8 Sep 2025 15:22:49 -0500 Subject: [PATCH 03/10] get_oafeat() for #436 --- R/get_geoconnex.R | 40 +------------------------------- R/oafeat_tools.R | 59 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 39 deletions(-) diff --git a/R/get_geoconnex.R b/R/get_geoconnex.R index 9352300f..ef3f71c6 100644 --- a/R/get_geoconnex.R +++ b/R/get_geoconnex.R @@ -52,47 +52,9 @@ get_geoconnex_reference <- function(AOI, buffer = 0.5, status = TRUE) { - avail <- discover_geoconnex_reference() - - if(is.null(type)) { - warning("type is required, returning choices.") - return(avail) - } - base <- get("gocnx_ref_base_url", envir = nhdplusTools_env) - if(!type %in% avail$id) stop("Type must be in available ids. ", - "Check discover_geoconnex_reference()") - - base_call <- paste0(base, "collections/", type, "/items") - - if(is.character(AOI)) { - - AOI <- try(sf::read_sf(AOI)) - - if(!inherits(AOI, "sf")) { - stop("AOI did not return an sf object when read") - } - - } - - if(!inherits(AOI, "bbox")) { - AOI <- st_bbox(AOI) - } else if(!inherits(AOI, "bbox") && - grepl("point", sf::st_geometry_type(AOI), ignore.case = TRUE)) { - AOI <- sf::st_buffer(AOI, units::as_units(buffer, "m")) - } - - # pull features with paging if necessary - - bbox <- paste(AOI, collapse = ",") - - base_call <- paste0(base_call, "?bbox=", bbox) - - out <- get_features_paging(base_call, status = status) - - if(!is.null(t_srs)) out <- st_transform(out, t_srs) + get_oafeat(base, AOI, type, t_srs, buffer, status) - out } diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index 476f6910..7ec97e66 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -39,6 +39,65 @@ discover_oafeat <- function(landing_url) { dplyr::left_join(collections_meta, q_ables, by = "id") } +#' get geoconnex reference feature layers +#' @description +#' Queries the geoconnex reference feature server for features of interest. +#' +#' @param AOI bbox, sf polygon or point, or a URL that will return an sf object when passed to +#' \link[sf]{read_sf} +#' @param type character the feature type desired +#' @inheritParams query_usgs_geoserver +#' @param status boolean print status or not +#' @return sf data.frame containing requested features +#' @noRd +get_oafeat <- function(base, + AOI, + type = NULL, + t_srs = NULL, + buffer = 0.5, + status = TRUE) { + + avail <- discover_oafeat(base) + + if(is.null(type)) { + warning("type is required, returning choices.") + return(avail) + } + + if(!type %in% avail$id) stop("Type must be in available ids: ", paste(avail$id, collapse = ", "), ".") + + base_call <- paste0(base, "collections/", type, "/items") + + if(is.character(AOI)) { + + AOI <- try(sf::read_sf(AOI)) + + if(!inherits(AOI, "sf")) { + stop("AOI did not return an sf object when read") + } + + } + + if(!inherits(AOI, "bbox")) { + AOI <- st_bbox(AOI) + } else if(!inherits(AOI, "bbox") && + grepl("point", sf::st_geometry_type(AOI), ignore.case = TRUE)) { + AOI <- sf::st_buffer(AOI, units::as_units(buffer, "m")) + } + + # pull features with paging if necessary + + bbox <- paste(AOI, collapse = ",") + + base_call <- paste0(base_call, "?bbox=", bbox) + + out <- get_features_paging(base_call, status = status) + + if(!is.null(t_srs)) out <- st_transform(out, t_srs) + + out +} + get_features_paging <- function(base_call, limit = 1000, status = TRUE) { if(!grepl("\\?", base_call)) { From 1d26c9a44ce7abe382e44ec9ccd19fe032f05a51 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Tue, 9 Sep 2025 14:22:11 -0500 Subject: [PATCH 04/10] initial work on query_usgs_oafeat #436 --- R/A_nhdplusTools.R | 3 + R/downloading_tools.R | 2 +- R/oafeat_tools.R | 88 ++++++++++++++++++++++++-- man/query_usgs_oafeat.Rd | 53 ++++++++++++++++ tests/testthat/test_03_get_functions.R | 18 ++++++ 5 files changed, 158 insertions(+), 6 deletions(-) create mode 100644 man/query_usgs_oafeat.Rd diff --git a/R/A_nhdplusTools.R b/R/A_nhdplusTools.R index 52adf6dd..6798788d 100644 --- a/R/A_nhdplusTools.R +++ b/R/A_nhdplusTools.R @@ -101,6 +101,9 @@ assign("arcrest_root", "https://hydro.nationalmap.gov/arcgis/rest/services/", assign("gocnx_ref_base_url", "https://reference.geoconnex.us/", envir = nhdplusTools_env) +assign("usgs_water_root", "https://nhgf.dev-wma.chs.usgs.gov/api/fabric/pygeoapi/", + envir = nhdplusTools_env) + assign("split_flowlines_attributes", c("COMID", "toCOMID", "LENGTHKM"), envir = nhdplusTools_env) diff --git a/R/downloading_tools.R b/R/downloading_tools.R index 6e258819..ee4f601f 100644 --- a/R/downloading_tools.R +++ b/R/downloading_tools.R @@ -357,7 +357,7 @@ check7z <- function() { #' @noRd mem_get_json <- memoise::memoise(\(url) { tryCatch({ - retn <- httr::RETRY("GET", url, httr::accept_json()) + retn <- httr::GET(url, httr::accept_json()) if(retn$status_code == 200 & grepl("json", retn$headers$`content-type`)) { return(httr::content(retn, simplifyVector = FALSE, type = "application/json")) diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index 7ec97e66..be04eff5 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -1,3 +1,76 @@ +#' @title Query USGS Water OGC API Features +#' @description Query the USGS Water OGC API 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 server 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 an +#' 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 or numeric. A set of identifier(s) from the data +#' type requested, for example if NHDPlusV2, then a set of COMID(s). +#' @param type character. Type of feature to return +#' ('huc08','huc12', 'nhd', 'catchment', 'waterbodies', 'gagesII'). +#' If NULL (default) a data.frame of available resources is returned +#' @param filter character. An filter to pass to the query +#' @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 +#' @keywords internal +query_usgs_oafeat <- function(AOI = NULL, ids = NULL, + type = NULL, filter = NULL, + t_srs = NULL, + buffer = 0.5) { + + base <- get("usgs_water_root", envir = nhdplusTools_env) + + source <- data.frame(server = 'usgs_oafeat', + user_call = c('huc08_legacy', "huc12_nhdplusv2", + 'nhd','catchment', 'nhdarea', + 'nonnetwork', + 'waterbodies', + 'gagesII', "gagesII-basin"), + layer_name = c("nhdplusv2-huc08", "nhdplusv2-huc12", + "nhdflowline_network", "catchmentsp", 'nhdarea', + "nhdflowline_nonnetwork", + "nhdwaterbody", + "gagesii", "gagesii-basins"), + geom_name = c("the_geom", "the_geom", + "the_geom", "the_geom", "the_geom", + "the_geom", + "the_geom", + "the_geom", "the_geom"), + ids = c("huc8", "huc12", + "comid", "featureid", "comid", + "comid", + "comid", + "staid", "gage_id"), + page = c(FALSE, FALSE, + TRUE, TRUE, TRUE, + TRUE, + TRUE, + FALSE, TRUE)) + + if(is.null(type)) { + warning("type is required, returning choices.") + return(source) + } + + AOI <- check_query_params(AOI, ids, type, NULL, source, t_srs, buffer) + t_srs <- AOI$t_srs + AOI <- AOI$AOI + + get_oafeat(base, AOI = AOI, type = type, ids = ids, t_srs = t_srs, buffer = buffer, status = FALSE) + +} + #' discover OGC API feature layers #' @description #' Queries an OGC API feature server for available layers and @@ -18,16 +91,20 @@ discover_oafeat <- function(landing_url) { collections_meta <- dplyr::bind_rows( lapply(collections$collections, - \(x) c(x[c("id", "title", "description")], + \(x) c(x[c("id", "title", "description", "itemType")], list(url = filter_list_kvp(x$links, - "rel", "self", n = 1)$href)))) + "rel", "self", n = 1)$href)))) |> + dplyr::filter(.data$itemType == "feature") |> + dplyr::select(-"itemType") q_ables <- dplyr::bind_rows(lapply(collections$collections, \(x) { + if(x$itemType != "feature") return(NULL) + q <- filter_list_kvp(x$links, "rel", "http://www.opengis.net/def/rel/ogc/1.0/queryables", type = "application/schema+json", n = 1)$href |> mem_get_json() |> - (\(y) list(id = x$id, qs = y$properties))() + (\(y) if(is.null(y)) y else list(id = x$id, qs = y$properties))() q$qs <- q$qs[vapply(q$qs, \(x) all(c("title", "type") %in% names(x)), TRUE)] @@ -37,6 +114,7 @@ discover_oafeat <- function(landing_url) { })) dplyr::left_join(collections_meta, q_ables, by = "id") + } #' get geoconnex reference feature layers @@ -51,7 +129,7 @@ discover_oafeat <- function(landing_url) { #' @return sf data.frame containing requested features #' @noRd get_oafeat <- function(base, - AOI, + AOI, ids = NULL, type = NULL, t_srs = NULL, buffer = 0.5, @@ -64,7 +142,7 @@ get_oafeat <- function(base, return(avail) } - if(!type %in% avail$id) stop("Type must be in available ids: ", paste(avail$id, collapse = ", "), ".") + if(!type %in% avail$id) stop("Type must be in available ids: ", paste(unique(avail$id), collapse = ", "), ".") base_call <- paste0(base, "collections/", type, "/items") diff --git a/man/query_usgs_oafeat.Rd b/man/query_usgs_oafeat.Rd new file mode 100644 index 00000000..3ff12069 --- /dev/null +++ b/man/query_usgs_oafeat.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/oafeat_tools.R +\name{query_usgs_oafeat} +\alias{query_usgs_oafeat} +\title{Query USGS Water OGC API Features} +\usage{ +query_usgs_oafeat( + AOI = NULL, + ids = NULL, + type = NULL, + filter = 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 or numeric. A set of identifier(s) from the data +type requested, for example if NHDPlusV2, then a set of COMID(s).} + +\item{type}{character. Type of feature to return +('huc08','huc12', 'nhd', 'catchment', 'waterbodies', 'gagesII'). +If NULL (default) a data.frame of available resources is returned} + +\item{filter}{character. An filter to pass to the query} + +\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{ +Query the USGS Water OGC API 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 server 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 an +EPSG:5070 Albers Equal Area projection +} +\keyword{internal} diff --git a/tests/testthat/test_03_get_functions.R b/tests/testthat/test_03_get_functions.R index d445651e..43558a34 100644 --- a/tests/testthat/test_03_get_functions.R +++ b/tests/testthat/test_03_get_functions.R @@ -32,6 +32,24 @@ test_that("query water geoserver...",{ expect_error(query_usgs_geoserver(AOI = pt, id = 17010101, type = 'huc8_legacy')) }) +test_that("query water oafeat...",{ + testthat::skip_on_cran() + #available? + expect_warning(df <- nhdplusTools:::query_usgs_oafeat(), "required") + expect_equal(ncol(df), 6) + + # errors + # Bad type request + expect_error(nhdplusTools:::query_usgs_oafeat(AOI = pt, type = 'wrong'), + "Type") + # Missing AOI and ID(s) + expect_error(nhdplusTools:::query_usgs_oafeat(AOI = NULL, id = NULL, type = 'catchmentsp'), + "IDs or a spatial AOI") + # Providing both an AOI and ID(s) + expect_error(nhdplusTools:::query_usgs_oafeat(AOI = pt, id = 17010101, type = 'huc8_legacy'), + "Either") +}) + # Walk our way through the 7 different offerings... # server user_call geoserver ids # 1 wmadata huc08 huc08 huc8 From 578155e84eceecd2c414949ca345e365539d31b8 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Fri, 19 Sep 2025 14:28:27 -0500 Subject: [PATCH 05/10] get_huc functioning #436 --- R/A_nhdplusTools.R | 2 +- R/downloading_tools.R | 10 +++-- R/get_hydro.R | 4 +- R/oafeat_tools.R | 51 ++++++++++++++++++-------- tests/testthat/test_03_get_functions.R | 2 +- 5 files changed, 46 insertions(+), 23 deletions(-) diff --git a/R/A_nhdplusTools.R b/R/A_nhdplusTools.R index 6798788d..5d4cee7f 100644 --- a/R/A_nhdplusTools.R +++ b/R/A_nhdplusTools.R @@ -101,7 +101,7 @@ assign("arcrest_root", "https://hydro.nationalmap.gov/arcgis/rest/services/", assign("gocnx_ref_base_url", "https://reference.geoconnex.us/", envir = nhdplusTools_env) -assign("usgs_water_root", "https://nhgf.dev-wma.chs.usgs.gov/api/fabric/pygeoapi/", +assign("usgs_water_root", "https://labs-beta.waterdata.usgs.gov/api/fabric/pygeoapi/", envir = nhdplusTools_env) assign("split_flowlines_attributes", diff --git a/R/downloading_tools.R b/R/downloading_tools.R index ee4f601f..930dc283 100644 --- a/R/downloading_tools.R +++ b/R/downloading_tools.R @@ -371,7 +371,7 @@ mem_get_json <- memoise::memoise(\(url) { }) }) -#' @importFrom sf st_make_valid st_as_sfc st_bbox st_buffer st_transform +#' @importFrom sf st_make_valid st_as_sfc st_bbox st_buffer st_transform st_crs 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) } @@ -399,8 +399,12 @@ check_query_params <- function(AOI, ids, type, where, source, t_srs, buffer) { 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() + AOI = st_transform(AOI, 5070) |> + st_buffer(buffer) |> + st_bbox() |> + st_as_sfc() |> + st_make_valid() |> + st_transform(st_crs(AOI)) } } diff --git a/R/get_hydro.R b/R/get_hydro.R index a84f69b2..6cadb588 100644 --- a/R/get_hydro.R +++ b/R/get_hydro.R @@ -30,8 +30,8 @@ get_huc <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, type = "hu stop("type must be one of ", paste(allow_types, collapse = " ")) } - query_usgs_geoserver(AOI = AOI, ids = id, type = type, - t_srs = t_srs, buffer = buffer) + query_usgs_oafeat(AOI = AOI, ids = id, type = type, + t_srs = t_srs, buffer = buffer) } diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index be04eff5..fbedbb29 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -32,7 +32,7 @@ query_usgs_oafeat <- function(AOI = NULL, ids = NULL, base <- get("usgs_water_root", envir = nhdplusTools_env) source <- data.frame(server = 'usgs_oafeat', - user_call = c('huc08_legacy', "huc12_nhdplusv2", + user_call = c('huc08', "huc12_nhdplusv2", 'nhd','catchment', 'nhdarea', 'nonnetwork', 'waterbodies', @@ -59,7 +59,6 @@ query_usgs_oafeat <- function(AOI = NULL, ids = NULL, FALSE, TRUE)) if(is.null(type)) { - warning("type is required, returning choices.") return(source) } @@ -67,6 +66,14 @@ query_usgs_oafeat <- function(AOI = NULL, ids = NULL, t_srs <- AOI$t_srs AOI <- AOI$AOI + here <- filter(source, .data$user_call == !!type) + + type <- here$layer_name + + if(!is.null(ids)) { + ids <- stats::setNames(list(ids), here$ids) + } + get_oafeat(base, AOI = AOI, type = type, ids = ids, t_srs = t_srs, buffer = buffer, status = FALSE) } @@ -144,30 +151,42 @@ get_oafeat <- function(base, if(!type %in% avail$id) stop("Type must be in available ids: ", paste(unique(avail$id), collapse = ", "), ".") + avail <- filter(avail, id == !!type) + base_call <- paste0(base, "collections/", type, "/items") - if(is.character(AOI)) { + if(!is.null(AOI)) { + if(is.character(AOI)) { + + AOI <- try(sf::read_sf(AOI)) + + if(!inherits(AOI, "sf")) { + stop("AOI did not return an sf object when read") + } - AOI <- try(sf::read_sf(AOI)) + } - if(!inherits(AOI, "sf")) { - stop("AOI did not return an sf object when read") + if(!inherits(AOI, "bbox")) { + AOI <- st_bbox(AOI) + } else if(!inherits(AOI, "bbox") && + grepl("point", sf::st_geometry_type(AOI), ignore.case = TRUE)) { + AOI <- sf::st_buffer(AOI, units::as_units(buffer, "m")) } - } + # pull features with paging if necessary - if(!inherits(AOI, "bbox")) { - AOI <- st_bbox(AOI) - } else if(!inherits(AOI, "bbox") && - grepl("point", sf::st_geometry_type(AOI), ignore.case = TRUE)) { - AOI <- sf::st_buffer(AOI, units::as_units(buffer, "m")) - } + bbox <- paste(AOI, collapse = ",") - # pull features with paging if necessary + base_call <- paste0(base_call, "?bbox=", bbox) + } else { + + id_attribute <- names(ids) - bbox <- paste(AOI, collapse = ",") + ids <- ids[[1]] - base_call <- paste0(base_call, "?bbox=", bbox) + base_call <- paste0(base_call, "?", id_attribute, "=", paste(ids, collapse = ",")) + + } out <- get_features_paging(base_call, status = status) diff --git a/tests/testthat/test_03_get_functions.R b/tests/testthat/test_03_get_functions.R index 43558a34..63b3b263 100644 --- a/tests/testthat/test_03_get_functions.R +++ b/tests/testthat/test_03_get_functions.R @@ -35,7 +35,7 @@ test_that("query water geoserver...",{ test_that("query water oafeat...",{ testthat::skip_on_cran() #available? - expect_warning(df <- nhdplusTools:::query_usgs_oafeat(), "required") + df <- nhdplusTools:::query_usgs_oafeat() expect_equal(ncol(df), 6) # errors From ad1aabcfcadaf1e47f2a183672a634a12b59d744 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Wed, 24 Sep 2025 13:55:39 -0500 Subject: [PATCH 06/10] migrate nhdplus to use oafeat #436 --- R/geoserver_tools.R | 8 +++- R/get_geoconnex.R | 4 +- R/get_hydro.R | 26 +++++----- R/get_nhdplus.R | 12 ++--- R/oafeat_tools.R | 66 ++++++++++++++++++++++---- R/subset_nhdplus.R | 22 +++++++-- man/get_gagesII.Rd | 4 +- man/get_huc.Rd | 4 +- man/get_nhdarea.Rd | 4 +- man/get_nhdplus.Rd | 4 +- man/get_nwis.Rd | 4 +- man/get_waterbodies.Rd | 4 +- tests/testthat/test_03_get_functions.R | 2 + 13 files changed, 115 insertions(+), 49 deletions(-) diff --git a/R/geoserver_tools.R b/R/geoserver_tools.R index 1cf85176..cdd62f22 100644 --- a/R/geoserver_tools.R +++ b/R/geoserver_tools.R @@ -252,7 +252,7 @@ spatial_filter <- function(AOI, #' @title Construct an attribute filter for geoservers #' @description From a user provided vector of IDs and an attribute name, #' generate a WMS filter to pass to a geoserver. -#' @inheritParams query_usgs_geoserver +#' @inheritParams query_usgs_oafeat #' @param name character. The name of the id attribute field in the desired dataset #' @return a character string #' @keywords internal @@ -305,3 +305,9 @@ streamorder_filter <- function(streamorder){ '', streamorder - 1, '', '') } + +streamorder_filter_cql <- function(streamorder) { + if(is.null(streamorder)){ return(NULL)} + + paste0("streamorde%20>%20", streamorder - 1) +} diff --git a/R/get_geoconnex.R b/R/get_geoconnex.R index ef3f71c6..c137b238 100644 --- a/R/get_geoconnex.R +++ b/R/get_geoconnex.R @@ -23,7 +23,7 @@ discover_geoconnex_reference <- function() { #' @param AOI bbox, sf polygon or point, or a URL that will return an sf object when passed to #' \link[sf]{read_sf} #' @param type character the feature type chosen from \link{discover_geoconnex_reference} -#' @inheritParams query_usgs_geoserver +#' @inheritParams query_usgs_oafeat #' @param status boolean print status or not #' @return sf data.frame containing requested reference features #' @export @@ -54,7 +54,7 @@ get_geoconnex_reference <- function(AOI, base <- get("gocnx_ref_base_url", envir = nhdplusTools_env) - get_oafeat(base, AOI, type, t_srs, buffer, status) + get_oafeat(base, AOI, type = type, t_srs = t_srs, buffer = buffer, status = status) } diff --git a/R/get_hydro.R b/R/get_hydro.R index 6cadb588..fa120a67 100644 --- a/R/get_hydro.R +++ b/R/get_hydro.R @@ -2,7 +2,7 @@ #' @description Subsets WBD features by location (POINT), #' area (POLYGON), or set of HUC IDs. #' -#' @inherit query_usgs_geoserver details return params +#' @inherit query_usgs_oafeat details return params #' @param id WBD HUC ID(s) #' @param type character. Type of feature to return #' ('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12', 'huc12_nhdplusv2'). @@ -38,13 +38,13 @@ get_huc <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, type = "hu #' @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 +#' @inherit query_usgs_oafeat details return +#' @inheritParams query_usgs_oafeat #' @param id NHD Waterbody COMID(s) #' @export get_waterbodies <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ - query_usgs_geoserver(AOI = AOI, ids = id, + query_usgs_oafeat(AOI = AOI, ids = id, type = "waterbodies", t_srs = t_srs, buffer = buffer) @@ -53,13 +53,13 @@ get_waterbodies <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ #' @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 +#' @inherit query_usgs_oafeat details return +#' @inheritParams query_usgs_oafeat #' @param id NHD Area COMID(s) #' @export get_nhdarea <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ - query_usgs_geoserver(AOI = AOI, ids = id, type = "nhdarea", + query_usgs_oafeat(AOI = AOI, ids = id, type = "nhdarea", t_srs = t_srs, buffer = buffer) } @@ -67,8 +67,8 @@ 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. See for documentation of source data. -#' @inherit query_usgs_geoserver details return -#' @inheritParams query_usgs_geoserver +#' @inherit query_usgs_oafeat details return +#' @inheritParams query_usgs_oafeat #' @param id character NWIS Gage ID(s) #' @param basin logical should the gagesII basin also be returned? If True, #' return value will be a list with "site" and "basin" elements. @@ -77,12 +77,12 @@ get_nhdarea <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5){ get_gagesII <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, basin = FALSE){ - out <- query_usgs_geoserver(AOI = AOI, ids = id, type = "gagesII", + out <- query_usgs_oafeat(AOI = AOI, ids = id, type = "gagesII", t_srs = t_srs, buffer = buffer) if(basin) { return(list(site = out, - basin = query_usgs_geoserver( + basin = query_usgs_oafeat( ids = out[["staid"]], type = "gagesII-basin", t_srs = t_srs, buffer = buffer))) } @@ -94,8 +94,8 @@ get_gagesII <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, #' @description Returns a POINT feature class of active, stream network, #' NWIS gages for an Area of Interest. If a POINT feature is used as an AOI, #' then the returned sites within the requested buffer, are sorted by distance (in meters) from that POINT. -#' @inherit query_usgs_geoserver details return -#' @inheritParams query_usgs_geoserver +#' @inherit query_usgs_oafeat details return +#' @inheritParams query_usgs_oafeat #' @param buffer numeric. The amount (in meters) to buffer a POINT AOI by #' for an extended search. Default = 20,000. Returned results are arrange #' by distance from POINT AOI diff --git a/R/get_nhdplus.R b/R/get_nhdplus.R index efb7e1e2..7083d83f 100644 --- a/R/get_nhdplus.R +++ b/R/get_nhdplus.R @@ -2,8 +2,8 @@ #' @description Subsets NHDPlusV2 features by location (POINT), area (POLYGON), #' or set of COMIDs. Multi realizations are supported allowing you to query #' for flowlines, catchments, or outlets. -#' @inherit query_usgs_geoserver details -#' @inheritParams query_usgs_geoserver +#' @inherit query_usgs_oafeat details +#' @inheritParams query_usgs_oafeat #' @param comid numeric or character. Search for NHD features by COMID(s) #' @param nwis numeric or character. Search for NHD features by #' collocated NWIS identifiers @@ -85,15 +85,15 @@ get_nhdplus <- function(AOI = NULL, } if("catchment" %in% realization){ - geoms$catchment <- query_usgs_geoserver(AOI = AOI, ids = comid, + geoms$catchment <- query_usgs_oafeat(AOI = AOI, ids = comid, type = "catchment", t_srs = t_srs) } if(any(c("flowline", "outlet") %in% realization)){ - geoms$flowline <- query_usgs_geoserver(AOI = AOI, ids = comid, type = 'nhd', - filter = streamorder_filter(streamorder), - t_srs = t_srs) + geoms$flowline <- query_usgs_oafeat(AOI = AOI, ids = comid, type = 'nhd', + filter = streamorder_filter_cql(streamorder), + t_srs = t_srs) if("outlet" %in% realization){ geoms$outlet <- geoms$flowline diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index fbedbb29..fe9ecfcc 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -74,7 +74,29 @@ query_usgs_oafeat <- function(AOI = NULL, ids = NULL, ids <- stats::setNames(list(ids), here$ids) } - get_oafeat(base, AOI = AOI, type = type, ids = ids, t_srs = t_srs, buffer = buffer, status = FALSE) + out <- get_oafeat(base, AOI = AOI, type = type, ids = ids, filter = filter, t_srs = t_srs, buffer = buffer, status = FALSE) + + out <- check_valid(out) + + if(any(is.null(out), nrow(out) == 0)) { + + out = NULL + + } else if(!is.null(AOI)){ + + out = st_filter(st_transform(out, t_srs), + st_transform(AOI, t_srs)) + if(nrow(out) == 0){ + out = NULL + } + } + + if(!is.null(out)) { + return(out) + } else { + warning(paste("No", here$user_call, "features found"), call. = FALSE) + NULL + } } @@ -131,13 +153,16 @@ discover_oafeat <- function(landing_url) { #' @param AOI bbox, sf polygon or point, or a URL that will return an sf object when passed to #' \link[sf]{read_sf} #' @param type character the feature type desired -#' @inheritParams query_usgs_geoserver +#' @param filter character CQL filter string +#' @inheritParams query_usgs_oafeat #' @param status boolean print status or not #' @return sf data.frame containing requested features #' @noRd get_oafeat <- function(base, - AOI, ids = NULL, + AOI, + ids = NULL, type = NULL, + filter = NULL, t_srs = NULL, buffer = 0.5, status = TRUE) { @@ -178,16 +203,32 @@ get_oafeat <- function(base, bbox <- paste(AOI, collapse = ",") base_call <- paste0(base_call, "?bbox=", bbox) - } else { + + } else if(!is.null(ids)) { + + if(!is.null(filter)) stop("filter not supported with ids") id_attribute <- names(ids) ids <- ids[[1]] - base_call <- paste0(base_call, "?", id_attribute, "=", paste(ids, collapse = ",")) + # filter=comid+IN+(101,13293454) + if(length(ids) == 1) { + filt = paste0("?", id_attribute, "=", ids) + } else if(is.character(ids)) { + filt <- paste0("?filter=", id_attribute, '+IN+("', paste(ids, collapse = '","'),'")') + } else { + filt <- paste0("?filter=", id_attribute, '+IN+(', paste(ids, collapse = ','),')') + } + + base_call <- paste0(base_call, filt) } + if(!is.null(filter)) { + base_call <- paste0(add_sep(base_call), "filter=", filter) + } + out <- get_features_paging(base_call, status = status) if(!is.null(t_srs)) out <- st_transform(out, t_srs) @@ -195,13 +236,18 @@ get_oafeat <- function(base, out } -get_features_paging <- function(base_call, limit = 1000, status = TRUE) { - - if(!grepl("\\?", base_call)) { - base_call <- paste0(base_call, "?") +add_sep <- function(bc) { + if(!grepl("\\?", bc)) { + bc <- paste0(bc, "?") } else { - base_call <- paste0(base_call, "&") + bc <- paste0(bc, "&") } + bc +} + +get_features_paging <- function(base_call, limit = 1000, status = TRUE) { + + base_call <- add_sep(base_call) offset <- 0 diff --git a/R/subset_nhdplus.R b/R/subset_nhdplus.R index 30ce0cb4..b5e37eb6 100644 --- a/R/subset_nhdplus.R +++ b/R/subset_nhdplus.R @@ -556,6 +556,18 @@ check_valid <- function(x, out_prj = sf::st_crs(x)) { })) } + if(as.character(sf::st_geometry_type(x, by_geometry = FALSE)) == "MULTIPOLYGON") { + if(all(sapply(sf::st_geometry(x), length) == 1)) { + try(suppressWarnings(x <- sf::st_cast(x, "POLYGON"))) + } + } + + if(as.character(sf::st_geometry_type(x, by_geometry = FALSE)) == "MULTILINESTRING") { + if(all(sapply(sf::st_geometry(x), length) == 1)) { + try(suppressWarnings(x <- sf::st_cast(x, "LINESTRING"))) + } + } + if (sf::st_crs(x) != sf::st_crs(out_prj)) { x <- sf::st_transform(x, out_prj) } @@ -829,9 +841,9 @@ get_all_navigable <- function(fline, rpu) { get_nhdplus_byid <- function(comids, layer, streamorder = NULL) { if(layer == "nhdflowline_network"){ - query_usgs_geoserver(ids = comids, type = "nhd", filter = streamorder_filter(streamorder)) + query_usgs_oafeat(ids = comids, type = "nhd", filter = streamorder_filter_cql(streamorder)) } else if(layer == "catchmentsp"){ - query_usgs_geoserver(ids = comids, type = "catchment") + query_usgs_oafeat(ids = comids, type = "catchment") } else { stop("Layer must be one of catchmentsp, nhdflowline_network") } @@ -850,11 +862,11 @@ get_nhdplus_bybox <- function(box, layer, streamorder = NULL) { if(layer == "nhdflowline_network") { - query_usgs_geoserver(AOI = box, + query_usgs_oafeat(AOI = box, type = type, - filter = streamorder_filter(streamorder)) + filter = streamorder_filter_cql(streamorder)) } else { - query_usgs_geoserver(AOI = box, + query_usgs_oafeat(AOI = box, type = type) } } diff --git a/man/get_gagesII.Rd b/man/get_gagesII.Rd index 8e8b27a0..fff0249a 100644 --- a/man/get_gagesII.Rd +++ b/man/get_gagesII.Rd @@ -33,9 +33,9 @@ area (POLYGON), or set of IDs. See for documentation of s \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 server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } diff --git a/man/get_huc.Rd b/man/get_huc.Rd index ddc25ad9..ffb8939e 100644 --- a/man/get_huc.Rd +++ b/man/get_huc.Rd @@ -43,9 +43,9 @@ area (POLYGON), or set of HUC 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 +IDs are used to query, then the default server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } diff --git a/man/get_nhdarea.Rd b/man/get_nhdarea.Rd index f0e67600..f284e67b 100644 --- a/man/get_nhdarea.Rd +++ b/man/get_nhdarea.Rd @@ -30,9 +30,9 @@ area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data doc \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 server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } diff --git a/man/get_nhdplus.Rd b/man/get_nhdplus.Rd index 6cb5c30b..43037cc2 100644 --- a/man/get_nhdplus.Rd +++ b/man/get_nhdplus.Rd @@ -47,10 +47,10 @@ for flowlines, catchments, or outlets. \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 server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } \examples{ diff --git a/man/get_nwis.Rd b/man/get_nwis.Rd index 66e82b1d..00b3d43c 100644 --- a/man/get_nwis.Rd +++ b/man/get_nwis.Rd @@ -30,9 +30,9 @@ then the returned sites within the requested buffer, are sorted by distance (in \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 server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } diff --git a/man/get_waterbodies.Rd b/man/get_waterbodies.Rd index b72ef34e..1669e59d 100644 --- a/man/get_waterbodies.Rd +++ b/man/get_waterbodies.Rd @@ -30,9 +30,9 @@ area (POLYGON), or set of IDs. See \link{download_nhdplusv2} for source data doc \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 server 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 +All buffer and distance operations are handled internally using an EPSG:5070 Albers Equal Area projection } diff --git a/tests/testthat/test_03_get_functions.R b/tests/testthat/test_03_get_functions.R index 63b3b263..a171656d 100644 --- a/tests/testthat/test_03_get_functions.R +++ b/tests/testthat/test_03_get_functions.R @@ -48,6 +48,8 @@ test_that("query water oafeat...",{ # Providing both an AOI and ID(s) expect_error(nhdplusTools:::query_usgs_oafeat(AOI = pt, id = 17010101, type = 'huc8_legacy'), "Either") + + # nhdplusTools:::query_usgs_oafeat(AOI = pt2, type = "huc08_legacy") }) # Walk our way through the 7 different offerings... From dfd3b3982bf66acf56b1974e385eb6050d950571 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Fri, 26 Sep 2025 09:50:25 -0500 Subject: [PATCH 07/10] add 2020 wbd layers and debug some test issues --- R/get_hydro.R | 8 +++++- R/oafeat_tools.R | 20 +++++++++---- R/subset_nhdplus.R | 4 +-- tests/testthat/test_03_get_functions.R | 25 +++++++++------- tests/testthat/test_arcrest.R | 40 +++++++++++++------------- 5 files changed, 58 insertions(+), 39 deletions(-) diff --git a/R/get_hydro.R b/R/get_hydro.R index fa120a67..a3787d8b 100644 --- a/R/get_hydro.R +++ b/R/get_hydro.R @@ -24,12 +24,18 @@ get_huc <- function(AOI = NULL, id = NULL, t_srs = NULL, buffer = .5, type = "huc12") { allow_types <- c('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12', - 'huc12_nhdplusv2') + 'huc02_2020', 'huc04_2020', 'huc06_2020', 'huc08_2020', 'huc10_2020', 'huc12_2020', + 'huc08_nhdplusv2', 'huc12_nhdplusv2') if(!type %in% allow_types) { stop("type must be one of ", paste(allow_types, collapse = " ")) } + if(type %in% c('huc02', 'huc04', 'huc06', 'huc08', 'huc10', 'huc12')) { + type <- paste0(type, "_2020") + message("defaulting to 2020 version of WBD") + } + query_usgs_oafeat(AOI = AOI, ids = id, type = type, t_srs = t_srs, buffer = buffer) diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index fe9ecfcc..d2d5a761 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -32,27 +32,37 @@ query_usgs_oafeat <- function(AOI = NULL, ids = NULL, base <- get("usgs_water_root", envir = nhdplusTools_env) source <- data.frame(server = 'usgs_oafeat', - user_call = c('huc08', "huc12_nhdplusv2", + user_call = c('huc02_2020', 'huc04_2020', 'huc06_2020', + 'huc08_2020', 'huc10_2020', 'huc12_2020', + 'huc08_nhdplusv2', "huc12_nhdplusv2", 'nhd','catchment', 'nhdarea', 'nonnetwork', 'waterbodies', 'gagesII', "gagesII-basin"), - layer_name = c("nhdplusv2-huc08", "nhdplusv2-huc12", + layer_name = c("wbd02_20201006", "wbd04_20201006", "wbd06_20201006", + "wbd08_20201006", "wbd10_20201006", "wbd12_20201006", + "nhdplusv2-huc08", "nhdplusv2-huc12", "nhdflowline_network", "catchmentsp", 'nhdarea', "nhdflowline_nonnetwork", "nhdwaterbody", "gagesii", "gagesii-basins"), - geom_name = c("the_geom", "the_geom", + geom_name = c("SHAPE", "SHAPE", "SHAPE", + "SHAPE", "SHAPE", "SHAPE", + "the_geom", "the_geom", "the_geom", "the_geom", "the_geom", "the_geom", "the_geom", "the_geom", "the_geom"), - ids = c("huc8", "huc12", + ids = c("huc2", "huc4", "huc6", + "huc8", "huc10", "huc12", + "huc8", "huc12", "comid", "featureid", "comid", "comid", "comid", "staid", "gage_id"), - page = c(FALSE, FALSE, + page = c(FALSE, FALSE, FALSE, + FALSE, FALSE, FALSE, + FALSE, FALSE, TRUE, TRUE, TRUE, TRUE, TRUE, diff --git a/R/subset_nhdplus.R b/R/subset_nhdplus.R index b5e37eb6..397a2886 100644 --- a/R/subset_nhdplus.R +++ b/R/subset_nhdplus.R @@ -857,8 +857,8 @@ get_nhdplus_bybox <- function(box, layer, streamorder = NULL) { stop("Layer must be one of nhdarea, nhdwaterbody, nhdflowline_network, nhdflowline_nonnetwork, catchmentsp.") } - type <- dplyr::filter(query_usgs_geoserver(), - .data$geoserver == layer)$user_call + type <- dplyr::filter(query_usgs_oafeat(), + .data$layer_name == layer)$user_call if(layer == "nhdflowline_network") { diff --git a/tests/testthat/test_03_get_functions.R b/tests/testthat/test_03_get_functions.R index a171656d..2f14a083 100644 --- a/tests/testthat/test_03_get_functions.R +++ b/tests/testthat/test_03_get_functions.R @@ -67,20 +67,20 @@ test_that("query water oafeat...",{ test_that("huc8", { testthat::skip_on_cran() #Point - ptHUC8 = get_huc(AOI = pt, type = "huc08") + ptHUC8 = get_huc(AOI = pt, type = "huc08_2020") expect_equal(nrow(ptHUC8), 1) expect_equal(ptHUC8$huc8, "17010101") expect_equal(sf::st_crs(ptHUC8)$epsg, 4326) - expect_error(get_huc(AOI = rbind(pt, pt2), type = "huc08"), + expect_error(get_huc(AOI = rbind(pt, pt2), type = "huc08_2020"), "AOI must be one an only one feature.") #Area - areaHUC8 = get_huc(AOI = area, t_srs = 5070, type = "huc08") + areaHUC8 = get_huc(AOI = area, t_srs = 5070, type = "huc08_2020") expect_equal(sf::st_crs(areaHUC8)$epsg, 5070) expect_equal(nrow(areaHUC8), 1) #ID - ptHUC8id = get_huc(id = "17010101", type = "huc08") + ptHUC8id = get_huc(id = "17010101", type = "huc08_2020") expect_identical(ptHUC8$huc8, ptHUC8id$huc8) expect_true(sf::st_crs(ptHUC8) == sf::st_crs(ptHUC8id) ) }) @@ -91,7 +91,7 @@ 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") + "type must be one of") #Point ptHUC12 = get_huc(AOI = pt, type = "huc12_nhdplusv2") expect_equal(nrow(ptHUC12), 1) @@ -108,27 +108,30 @@ test_that("huc", { expect_identical(HUC12id2$geometry, areaHUC12$geometry) - hu12 <- get_huc(AOI = pt, type = "huc12") + hu12 <- get_huc(AOI = pt, type = "huc12_2020") expect_equal(hu12$huc12, "170101010806") - hu10 <- get_huc(AOI = pt, type = "huc10") + expect_message(get_huc(AOI = pt, type = "huc12"), + "defaulting to 2020 version of WBD") # TODO default to final WBD + + hu10 <- get_huc(AOI = pt, type = "huc10_2020") expect_equal(hu10$huc10, "1701010108") - hu08 <- get_huc(AOI = pt, type = "huc08") + hu08 <- get_huc(AOI = pt, type = "huc08_2020") expect_equal(hu08$huc8, "17010101") - hu06 <- get_huc(AOI = pt, type = "huc06") + hu06 <- get_huc(AOI = pt, type = "huc06_2020") expect_equal(hu06$huc6, "170101") - hu04 <- get_huc(AOI = pt, type = "huc04") + hu04 <- get_huc(AOI = pt, type = "huc04_2020") expect_equal(hu04$huc4, "1701") - hu02 <- get_huc(AOI = pt, type = "huc02") + hu02 <- get_huc(AOI = pt, type = "huc02_2020") expect_equal(hu02$huc2, "17") diff --git a/tests/testthat/test_arcrest.R b/tests/testthat/test_arcrest.R index 222a30a5..dab84669 100644 --- a/tests/testthat/test_arcrest.R +++ b/tests/testthat/test_arcrest.R @@ -1,17 +1,17 @@ -AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816, +AOI_esri <- 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"), + expect_error(nhdplusTools:::spatial_filter(AOI_esri, format = "test"), "ogc or esri") expect_equal( - nhdplusTools:::spatial_filter(AOI, format = "ogc")[[1]], + nhdplusTools:::spatial_filter(AOI_esri, format = "ogc")[[1]], "the_geom42.99816 -89.5668443.17192 -89.24681") - expect_equal(as.character(nhdplusTools:::spatial_filter(AOI, format = "esri")[[1]]), + expect_equal(as.character(nhdplusTools:::spatial_filter(AOI_esri, format = "esri")[[1]]), '{"xmin":-89.5668,"ymin":42.9982,"xmax":-89.2468,"ymax":43.1719,"spatialReference":{"wkid":4326}}') }) @@ -25,30 +25,30 @@ test_that("check_query_params", { 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))) + nhdplusTools:::check_query_params(AOI_esri, ids, type, where, source, t_srs, buffer), + list(AOI = AOI_esri, t_srs = sf::st_crs(AOI_esri))) - expect_error(nhdplusTools:::check_query_params(AOI, c(1,2), type, where, source, t_srs, buffer), + expect_error(nhdplusTools:::check_query_params(AOI_esri, 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), + expect_error(nhdplusTools:::check_query_params(AOI_esri, 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), + expect_error(nhdplusTools:::check_query_params(c(AOI_esri, AOI_esri), 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) + AOI_new <- 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), + nhdplusTools:::check_query_params(AOI_new, ids, type, where, source, t_srs, buffer)$AOI[[1]], + structure(list(structure(c(-89.5668465687776, -89.5668457346386, + -89.5668334312216, -89.5668342653621, + -89.5668465687776, 42.9981558371707, + 42.9981647683392, 42.9981641628291, + 42.9981552316606, 42.9981558371707), dim = c(5L, 2L))), class = c("XY", "POLYGON", "sfg")), tolerance = 0.1) @@ -57,17 +57,17 @@ test_that("check_query_params", { test_that("basic 3dhp service requests", { skip_on_cran() - AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.5, ymin = 43.0, + AOI_new <- 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, service = "3DHP_all"), + expect_message(expect_s3_class(nhdplusTools:::query_usgs_arcrest(AOI_new, service = "3DHP_all"), "data.frame")) - expect_warning(nhdplusTools:::query_usgs_arcrest(AOI, service = "3DHP_all", + expect_warning(nhdplusTools:::query_usgs_arcrest(AOI_new, service = "3DHP_all", type = "hydrolocation")) - test_data <- nhdplusTools:::query_usgs_arcrest(AOI, , service = "3DHP_all", + test_data <- nhdplusTools:::query_usgs_arcrest(AOI_new, , service = "3DHP_all", type = "hydrolocation - reach code, external connection") expect_s3_class(test_data, "sf") From f91d5634f1409c1fda6c9c986ea508cf90fd7efc Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Fri, 26 Sep 2025 16:47:22 -0500 Subject: [PATCH 08/10] add ability to POST cql --- R/geoserver_tools.R | 7 +++++++ R/oafeat_tools.R | 49 +++++++++++++++++++++++++++++++++++---------- 2 files changed, 45 insertions(+), 11 deletions(-) diff --git a/R/geoserver_tools.R b/R/geoserver_tools.R index cdd62f22..30585dff 100644 --- a/R/geoserver_tools.R +++ b/R/geoserver_tools.R @@ -286,6 +286,13 @@ id_filter <- function(ids, name = "comid"){ } +id_filter_cql <- function(ids, name = "comid"){ + + jsonlite::toJSON(list(op = "in", args = list(list(property = name), c(ids))), + auto_unbox = TRUE) + +} + #' @title Stream Order Filter #' @description Generate a WMS filter to pass to a geoserver. That constrains #' the returned flowlines to only those with a stream order diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index d2d5a761..f01c0eb2 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -189,6 +189,7 @@ get_oafeat <- function(base, avail <- filter(avail, id == !!type) base_call <- paste0(base, "collections/", type, "/items") + post_body <- "" if(!is.null(AOI)) { if(is.character(AOI)) { @@ -222,24 +223,19 @@ get_oafeat <- function(base, ids <- ids[[1]] - # filter=comid+IN+(101,13293454) if(length(ids) == 1) { - filt = paste0("?", id_attribute, "=", ids) - } else if(is.character(ids)) { - filt <- paste0("?filter=", id_attribute, '+IN+("', paste(ids, collapse = '","'),'")') + base_call <- paste0(base_call, paste0("?", id_attribute, "=", ids)) } else { - filt <- paste0("?filter=", id_attribute, '+IN+(', paste(ids, collapse = ','),')') + post_body <- id_filter_cql(ids, id_attribute) } - base_call <- paste0(base_call, filt) - } if(!is.null(filter)) { base_call <- paste0(add_sep(base_call), "filter=", filter) } - out <- get_features_paging(base_call, status = status) + out <- get_features_paging(base_call, post_body, status = status) if(!is.null(t_srs)) out <- st_transform(out, t_srs) @@ -255,7 +251,38 @@ add_sep <- function(bc) { bc } -get_features_paging <- function(base_call, limit = 1000, status = TRUE) { +make_request <- function(req, body = "") { + tryCatch({ + if(nhdplus_debug()) { + message(paste(req, "\n")) + message(body) + } + + if(body != "") { + out <- rawToChar(RETRY("POST", + req, + body = body, + httr::content_type("application/query-cql-json"))$content) + } else { + + out <- rawToChar(RETRY("GET", req)$content) + + } + + + }, error = function(e) { + warning("Something went wrong trying to access a service.") + return(NULL) + }) + + out <- tryCatch({ + st_zm(read_sf(out))}, + error = function(e) return(NULL)) + + out +} + +get_features_paging <- function(base_call, post_body = "", limit = 1000, status = TRUE) { base_call <- add_sep(base_call) @@ -271,9 +298,9 @@ get_features_paging <- function(base_call, limit = 1000, status = TRUE) { while(keep_going) { req <- paste0(base_call, "limit=", limit, "&offset=", offset) - out[[i]] <- try(read_sf(req)) + out[[i]] <- make_request(req, post_body) - if(inherits(out[[i]], "sf") & nrow(out[[i]]) == limit) { + if(!is.null(out[[i]]) && inherits(out[[i]], "sf") & nrow(out[[i]]) == limit) { offset <- offset + limit } From 4a76802e0e89aeffb8f744a05e14e21fc27a69ba Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Tue, 30 Sep 2025 14:28:45 -0500 Subject: [PATCH 09/10] working on api.water --- R/A_nhdplusTools.R | 12 +++++++++++- R/oafeat_tools.R | 40 ++++++++++++++++++++++++++++++++------ vignettes/plot_nhdplus.Rmd | 2 +- 3 files changed, 46 insertions(+), 8 deletions(-) diff --git a/R/A_nhdplusTools.R b/R/A_nhdplusTools.R index 5d4cee7f..770e5b37 100644 --- a/R/A_nhdplusTools.R +++ b/R/A_nhdplusTools.R @@ -101,7 +101,7 @@ assign("arcrest_root", "https://hydro.nationalmap.gov/arcgis/rest/services/", assign("gocnx_ref_base_url", "https://reference.geoconnex.us/", envir = nhdplusTools_env) -assign("usgs_water_root", "https://labs-beta.waterdata.usgs.gov/api/fabric/pygeoapi/", +assign("usgs_water_root", "https://api.water.usgs.gov/fabric/pygeoapi/", envir = nhdplusTools_env) assign("split_flowlines_attributes", @@ -526,6 +526,16 @@ filter_list_kvp <- \(l, key, val, type = NULL, n = NULL) { extract <- `[[` +split_equal_size <- function (x, n) +{ + nr <- try(nrow(x), silent = TRUE) + if (inherits(nr, "try-error") | is.null(nr)) + nr <- try(length(x), silent = TRUE) + if (!inherits(nr, "numeric") & length(nr) != 1) + stop("x can't be interpreted as a data.frame or list") + split(x, rep(1:ceiling(nr/n), each = n, length.out = nr)) +} + #' @title Trim and Cull NULLs #' @description Remove NULL arguments from a list #' @param x a list diff --git a/R/oafeat_tools.R b/R/oafeat_tools.R index f01c0eb2..72e48fdf 100644 --- a/R/oafeat_tools.R +++ b/R/oafeat_tools.R @@ -189,7 +189,9 @@ get_oafeat <- function(base, avail <- filter(avail, id == !!type) base_call <- paste0(base, "collections/", type, "/items") - post_body <- "" + post_body <- list() + + limit <- 1000 if(!is.null(AOI)) { if(is.character(AOI)) { @@ -226,7 +228,8 @@ get_oafeat <- function(base, if(length(ids) == 1) { base_call <- paste0(base_call, paste0("?", id_attribute, "=", ids)) } else { - post_body <- id_filter_cql(ids, id_attribute) + post_body <- list(ids = ids, id_attribute = id_attribute) + limit <- 500 } } @@ -235,7 +238,7 @@ get_oafeat <- function(base, base_call <- paste0(add_sep(base_call), "filter=", filter) } - out <- get_features_paging(base_call, post_body, status = status) + out <- get_features_paging(base_call, post_body, limit = limit, status = status) if(!is.null(t_srs)) out <- st_transform(out, t_srs) @@ -282,10 +285,18 @@ make_request <- function(req, body = "") { out } -get_features_paging <- function(base_call, post_body = "", limit = 1000, status = TRUE) { +get_features_paging <- function(base_call, ids_list = list(), limit = 1000, status = TRUE) { + + if(!identical(ids_list, list())) { + # we will page through ids + ids <- split_equal_size(ids_list$ids, limit) + } base_call <- add_sep(base_call) + + post_body <- "" + offset <- 0 keep_going <- TRUE @@ -296,7 +307,18 @@ get_features_paging <- function(base_call, post_body = "", limit = 1000, status i <- 1 while(keep_going) { - req <- paste0(base_call, "limit=", limit, "&offset=", offset) + + if(!identical(ids_list, list())) { + + req <- base_call + post_body <- id_filter_cql(ids[[i]], ids_list$id_attribute) + req <- paste0(base_call, "limit=", limit) + + } else { + + req <- paste0(base_call, "limit=", limit, "&offset=", offset) + + } out[[i]] <- make_request(req, post_body) @@ -311,7 +333,13 @@ get_features_paging <- function(base_call, post_body = "", limit = 1000, status keep_going <- FALSE } - if(status & keep_going) message("Starting next download from ", offset, ".") + if(status & keep_going) { + if(identical(ids_list, list())) { + message("Starting next download from ", offset, ".") + } else { + message("starting next download from ", i * limit, ".") + } + } i <- i + 1 } diff --git a/vignettes/plot_nhdplus.Rmd b/vignettes/plot_nhdplus.Rmd index a084922c..78277eec 100644 --- a/vignettes/plot_nhdplus.Rmd +++ b/vignettes/plot_nhdplus.Rmd @@ -176,7 +176,7 @@ mapsf::mf_map(bb, type = "base", col = NA, border = NA) maptiles::plot_tiles(tiles, add = TRUE) mapsf::mf_map(bb, type = "base", add = TRUE, col = NA, border = NA) -mapsf::mf_arrow(adjust = bb) +mapsf::mf_arrow(align = bb) mapsf::mf_scale() plot(prep_layer(basin), From 32e152414eb98d1b709850dad423eaeb8c932be6 Mon Sep 17 00:00:00 2001 From: David Blodgett Date: Tue, 30 Sep 2025 14:50:05 -0500 Subject: [PATCH 10/10] version update --- DESCRIPTION | 2 +- NEWS.md | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 1b9b2692..65000add 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nhdplusTools Type: Package Title: NHDPlus Tools -Version: 1.3.2 +Version: 1.4.0 Authors@R: c(person(given = "David", family = "Blodgett", role = c("aut", "cre"), diff --git a/NEWS.md b/NEWS.md index cac5e650..2cd5cea7 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +nhdplusTools 1.4.0 +========== + +This release migrates to a new web service provider for the key function `subset_nhdplus()`. +The change is mostly backward compatible but some minor differences in performance and response +data may be noticed. + +- `subset_nhdplus()` and `get_huc()` migrated to pygeoapi-based web services. + nhdplusTools 1.3.2 ==========