if (!"LongitudeMeasure" %in% colnames(.data) |

!"LatitudeMeasure" %in% colnames(.data) |

!"HorizontalCoordinateReferenceSystemDatumName" %in% colnames(.data)) {

stop("The dataframe does not contain WQP-style latitude and longitude data (column names `HorizontalCoordinateReferenceSystemDatumName`, `LatitudeMeasure`, and `LongitudeMeasure`.")

} else if (!is.null(.data) & inherits(.data, "sf")) {

stop("Your data is already a spatial object.")

}

message("Transforming your data into a spatial object.")

suppressMessages(suppressWarnings({

# Make a reference table for CRS and EPSG codes

# List should include all codes in WQX domain (see HorizontalCoordinateReferenceSystemDatum CSV at https://www.epa.gov/waterdata/storage-and-retrieval-and-water-quality-exchange-domain-services-and-downloads)

epsg_codes <- tidyr::tribble(

~HorizontalCoordinateReferenceSystemDatumName, ~epsg,

"NAD83", 4269,

"WGS84", 4326,

"NAD27", 4267,

"UNKWN", crs, # Unknowns and NAs should go to user supplied crs or default

"Unknown", crs,

"OTHER", crs,

"OLDHI", 4135,

"AMSMA", 4169,

"ASTRO", 4727,

"GUAM", 4675,

"JHNSN", 4725,

"PR", 6139,

"SGEOR", 4138,

"SLAWR", 4136,

"SPAUL", 4137,

"WAKE", 6732,

"WGS72", 6322,

"HARN", 4152

)

# Check the CRS column for NA or "UNKWN" and warn user if any are found

if (any(is.na(.data$HorizontalCoordinateReferenceSystemDatumName)) ||

any(.data$HorizontalCoordinateReferenceSystemDatumName %in% c("UNKWN", "Unknown", "OTHER"))) {

message(paste0("Your WQP dataframe contains observations without a listed coordinate reference system (CRS). For these, we have assigned CRS ", crs, "."))

}

# join our CRS reference table to our original WQP dataframe:

sf <- .data %>%

tibble::rowid_to_column(var = "index") %>%

dplyr::mutate(

lat = as.numeric(LatitudeMeasure),

lon = as.numeric(LongitudeMeasure),

# If `HorizontalCoordinateReferenceSystemDatumName` is NA...

HorizontalCoordinateReferenceSystemDatumName = ifelse(is.na(HorizontalCoordinateReferenceSystemDatumName),

# ... assign it the same crs as the user-supplied crs:

paste0(epsg_codes %>% dplyr::filter(epsg == as.numeric(crs)) %>% .[1, 1]),

# otherwise, preserve the original crs

HorizontalCoordinateReferenceSystemDatumName

)

) %>%

# Add EPSG codes

dplyr::left_join(

x = .,

y = epsg_codes,

by = "HorizontalCoordinateReferenceSystemDatumName"

) %>%

# Group by CRS:

split(f = .$HorizontalCoordinateReferenceSystemDatumName) %>%

# Transform and re-stack:

purrr::map_df(

.x = .,

.f = ~ .x %>%

sf::st_as_sf(

coords = c("lon", "lat"),

crs = unique(.x$epsg)

) %>%

# transform to the selected CRS:

sf::st_transform(sf::st_crs(as.numeric(crs)))

) %>%

dplyr::arrange(index) %>%

dplyr::select(-c(index, epsg))

}))

return(sf)

# function settings that we ensure go back to their original settings

# after the function stops running:

original_s2 <- sf::sf_use_s2() # Store the original s2 setting first

suppressMessages(sf::sf_use_s2(FALSE))

original_timeout <- getOption("timeout")

options(timeout = 30000)

on.exit(options(timeout = original_timeout), add = TRUE)

on.exit(suppressMessages(suppressWarnings(sf::sf_use_s2(original_s2))), add = TRUE)

message("Depending on your data's observation count and its spatial range, the ATTAINS pull may take a while.")

# EPSG we want our ATTAINS data to be in (always 4326 for this function)

our_epsg <- 4326

# If data is already spatial, just make sure it is in the right CRS

# and add an index as the WQP observations' unique identifier...

if (!is.null(.data) & inherits(.data, "sf")) {

if (sf::st_crs(.data)$epsg != our_epsg) {

.data <- .data %>%

sf::st_transform(our_epsg) %>%

dplyr::distinct(geometry, .keep_all = TRUE)

} else {

.data <- .data %>%

dplyr::distinct(geometry, .keep_all = TRUE)

}

} else if (!"LongitudeMeasure" %in% colnames(.data) |

!"LatitudeMeasure" %in% colnames(.data) |

!"HorizontalCoordinateReferenceSystemDatumName" %in% colnames(.data)) {

stop("The dataframe does not contain WQP-style latitude and longitude data (column names `HorizontalCoordinateReferenceSystemDatumName`, `LatitudeMeasure`, and `LongitudeMeasure`.")

} else {

# ... Otherwise transform into a spatial object then do the same thing:

.data <- .data %>%

data.table::data.table(.) %>%

dplyr::distinct(LongitudeMeasure, LatitudeMeasure, .keep_all = TRUE) %>%

# convert dataframe to a spatial object

TADA_MakeSpatial(.data = ., crs = our_epsg)

}

if (is.null(.data) | nrow(.data) == 0) {

stop("There is no data in your `data` object to use as a bounding box for selecting ATTAINS features.")

}

# REST for ATTAINS geospatial data:

baseurls <- c( # ATTAINS catchments:

"https://gispub.epa.gov/arcgis/rest/services/OW/ATTAINS_Assessment/MapServer/3/query?",

# ATTAINS points:

"https://gispub.epa.gov/arcgis/rest/services/OW/ATTAINS_Assessment/MapServer/0/query?",

# ATTAINS lines:

"https://gispub.epa.gov/arcgis/rest/services/OW/ATTAINS_Assessment/MapServer/1/query?",

# ATTAINS polygons:

"https://gispub.epa.gov/arcgis/rest/services/OW/ATTAINS_Assessment/MapServer/2/query?"

)

# function to download ATTAINS features based on a

# specified bbox:

fetch_bbox <- function(baseurls, sf_bbox) {

# starting at feature 1 (i.e., no offset):

offset <- 0

# empty list to store all features in

all_features <- list()

# The ATTAINS API has a limit of 1000 features that can be pulled in at once.

# Therefore, we must split the call into manageable "chunks" using a moving

# window of what features to pull in, then munging all the separate API calls

# together.

repeat {

query <- urltools::param_set(baseurls, key = "geometry", value = sf_bbox) %>%

urltools::param_set(key = "inSR", value = our_epsg) %>%

# Total of 1000 features at a time...

urltools::param_set(key = "resultRecordCount", value = 1000) %>%

# ... starting at the "offset":

urltools::param_set(key = "resultOffset", value = offset) %>%

urltools::param_set(key = "spatialRel", value = "esriSpatialRelIntersects") %>%

urltools::param_set(key = "f", value = "geojson") %>%

urltools::param_set(key = "outFields", value = "*") %>%

urltools::param_set(key = "geometryType", value = "esriGeometryEnvelope") %>%

urltools::param_set(key = "returnGeometry", value = "true") %>%

urltools::param_set(key = "returnTrueCurves", value = "false") %>%

urltools::param_set(key = "returnIdsOnly", value = "false") %>%

urltools::param_set(key = "returnCountOnly", value = "false") %>%

urltools::param_set(key = "returnZ", value = "false") %>%

urltools::param_set(key = "returnM", value = "false") %>%

urltools::param_set(key = "returnDistinctValues", value = "false") %>%

urltools::param_set(key = "returnExtentOnly", value = "false") %>%

urltools::param_set(key = "featureEncoding", value = "esriDefault")

# Fetch features within the offset window and append to list:

features <- suppressMessages(suppressWarnings({

tryCatch(

{

geojsonsf::geojson_sf(url(query))

},

error = function(e) {

NULL

}

)

}))

# Exit loop if no more features or error occurred

if (is.null(features) || nrow(features) == 0) {

break

}

all_features <- c(all_features, list(features))

# once done, change offset by 1000 features:

offset <- offset + 1000

}

all_features <- dplyr::bind_rows(all_features) %>%

# remove duplicate features (precautionary)

dplyr::distinct(.keep_all = TRUE)

}

# function to download ATTAINS features API based on their name

fetch_au <- function(baseurls, assessment_unit_ids) {

# Split the assessment_unit_ids into chunks of 1000

# API cannot handle more than 1000 features

id_chunks <- split(assessment_unit_ids, ceiling(seq_along(assessment_unit_ids) / 1000))

# Query API for a chunk of assessment unit IDs

fetch_chunk <- function(id_chunk) {

where_clause <- paste0("assessmentunitidentifier IN ('", paste(id_chunk, collapse = "','"), "')")

query_params <- list(

where = where_clause,

outFields = "*",

f = "geojson"

)

response <- httr::GET(baseurls, query = query_params)

if (httr::status_code(response) != 200) {

stop("Failed to retrieve data from EPA ATTAINS API.")

}

geojson_data <- httr::content(response, as = "text", encoding = "UTF-8")

sf_object <- sf::st_read(geojson_data, quiet = TRUE)

return(sf_object)

}

# fetch all chunks and combine results

purrr::map_dfr(id_chunks, fetch_chunk)

}

# function used to grab assessment unit "WaterType".

# (sweet spot chunk size wise is 50):

grab_waterbody_type <- function(au_list, chunk_size = 50) {

# Number of chunks needed

num_chunks <- ceiling(length(au_list) / chunk_size)

# split the au_list into chunks

chunks <- split(au_list, ceiling(seq_along(au_list) / chunk_size))

water_types <- vector("list", length = length(chunks))

for (i in 1:length(chunks)) {

dat <- httr::GET(utils::URLencode(paste0("https://attains.epa.gov/attains-public/api/assessmentUnits?assessmentUnitIdentifier=", paste(chunks[[i]], collapse = ",")))) %>%

httr::content(., as = "text", encoding = "UTF-8") %>%

jsonlite::fromJSON(.)

water_types[[i]] <- dat[["items"]] %>%

tidyr::unnest("assessmentUnits") %>%

tidyr::unnest("waterTypes") %>%

dplyr::select(

assessmentUnitIdentifier,

waterTypeCode

)

}

return(dplyr::bind_rows(water_types))

}

# FOR AOIs THAT ARE GREATER THAN 6,000 sqkm, split into "clusters":

if (as.numeric(sf::st_area(sf::st_as_sfc(.data %>% sf::st_bbox(.)))) >= 6e+9) {

# For user-specified AOIs with a large spatial range, create "clusters" of sites

# whose bounding boxes are smaller.

perform_iterative_clustering <- function(points_sf, min_area = 6e+9, max_iterations = 100) {

# fxn to calculate each cluster's bounding box area

bbox_area <- function(df, clust) {

df %>%

dplyr::filter(cluster == clust) %>%

sf::st_bbox() %>%

sf::st_as_sfc() %>%

sf::st_area() %>%

tidyr::as_tibble() %>%

dplyr::mutate(cluster = clust)

}

# cluster maker fxn:

cluster_iteration <- function(points, eps, min_pts, iteration) {

coords <- sf::st_coordinates(points)

fr <- dbscan::frNN(coords, eps = eps)

clusters <- dbscan::dbscan(fr, minPts = min_pts)$cluster

# create unique cluster IDs that include iteration number

cluster_ids <- ifelse(clusters == -1,

paste0("noise_", iteration),

paste0("cluster_", iteration, "_", clusters)

)

points %>%

dplyr::mutate(

cluster = cluster_ids,

iteration = iteration

)

}

# function to check if any clusters still exceed our bbox area threshold

has_large_clusters <- function(points) {

if (nrow(points) == 0) {

return(FALSE)

}

areas <- unique(points$cluster) %>%

purrr::map_dfr(~ bbox_area(df = points, clust = .))

# Returns TRUE/FALSE

any(as.numeric(areas$value) > min_area)

}

# function to split clusters by their area

split_clusters_by_area <- function(points, min_area) {

# calculate areas for all clusters

cluster_areas <- unique(points$cluster) %>%

purrr::map_dfr(~ bbox_area(df = points, clust = .))

# split into "large" and "small" clusters

large_clusters <- cluster_areas %>%

dplyr::filter(as.numeric(value) > min_area)

small_clusters <- cluster_areas %>%

dplyr::filter(as.numeric(value) <= min_area)

# create data frames for both sets

large_points <- points %>%

dplyr::filter(cluster %in% large_clusters$cluster)

small_points <- points %>%

dplyr::filter(cluster %in% small_clusters$cluster)

list(

large = large_points,

small = small_points,

large_areas = large_clusters,

small_areas = small_clusters

)

}

# store the results of the analysis run

all_small_clusters <- list()

current_points <- points_sf %>% dplyr::distinct(geometry)

iteration <- 1

# sequence of eps values to try

eps_sequence <- c(0.25, 0.05, 1, .1)

eps_index <- 1

while (nrow(current_points) > 0 && iteration <= max_iterations) {

# grab current eps value (cycle through sequence)

current_eps <- eps_sequence[eps_index]

eps_index <- (eps_index %% length(eps_sequence)) + 1

# ... perform cluster

clustered_points <- cluster_iteration(

current_points,

eps = current_eps,

# Minimum number of points within a single cluster

# (meaning, a point can be isolated if super far from others)

min_pts = 1,

iteration = iteration

)

# Split results

split_results <- split_clusters_by_area(clustered_points, min_area)

# Store small clusters

if (nrow(split_results$small) > 0) {

all_small_clusters[[paste0("iteration_", iteration)]] <- split_results$small

}

# check if we have any large bbox clusters to keep processing

if (nrow(split_results$large) == 0) {

break

}

# update points for next iteration

current_points <- split_results$large

iteration <- iteration + 1

}

# combine all results

final_clusters <- dplyr::bind_rows(all_small_clusters) %>%

dplyr::arrange(iteration)

# warning if max iterations reached. This should never happen...

if (iteration == max_iterations) {

warning("Maximum iterations reached. Some clusters may still exceed the area threshold.")

}

return(

list(

clusters = final_clusters,

clusters_by_iteration = all_small_clusters,

total_iterations = iteration,

final_eps = current_eps

)

)

}

# grab all unique points

points_sf <- dplyr::distinct(.data, geometry)

# grab initial clusters

init <- perform_iterative_clustering(points_sf = points_sf) %>%

.[["clusters_by_iteration"]] %>%

dplyr::bind_rows()

# if any didn't get "clustered"....

final_cluster_list <- points_sf %>%

dplyr::filter(!geometry %in% init$geometry) %>%

# tack them on to our df and give them a unique "cluster name"

tibble::rowid_to_column(var = "cluster") %>%

dplyr::mutate(cluster = as.character(cluster)) %>%

dplyr::bind_rows(init)

catchment_features <- vector("list", length = length(unique(final_cluster_list$cluster)))

for (i in 1:length(unique(final_cluster_list$cluster))) {

# bounding box of user's WQP data

suppressMessages(suppressWarnings({

bbox <- final_cluster_list %>%

dplyr::filter(cluster == unique(final_cluster_list$cluster)[i]) %>%

sf::st_bbox(.) %>%

toString(.) %>%

urltools::url_encode(.)

}))

catchment_features[[i]] <- fetch_bbox(baseurls = baseurls[1], sf_bbox = bbox)

}

catchment_features <- catchment_features %>%

purrr::keep(~ !is.null(.)) %>%

purrr::keep(~ nrow(.) > 0) %>%

dplyr::bind_rows()

try(

catchment_features <- catchment_features %>%

.[points_sf, ],

silent = TRUE

)

## GRABBING WATER TYPE:

if (length(catchment_features) == 0 || is.null(catchment_features)) {

message("There are no ATTAINS features associated with your WQP observations.")

} else {

# Use ATTAINS API to grab, for each assessment unit, its WaterType.

# Query the API in "chunks" so it doesn't break. sweet spot is ~50:

all_units <- unique(catchment_features$assessmentunitidentifier)

water_types <- grab_waterbody_type(all_units, chunk_size = 50)

try(catchment_features <- dplyr::left_join(catchment_features, water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")))

}

# If only interested in grabbing catchment data, return just the catchments

if (catchments_only == TRUE) {

return(list("ATTAINS_catchments" = catchment_features))

}

points <- fetch_au(

baseurls = baseurls[2],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

)

lines <- fetch_au(

baseurls = baseurls[3],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

)

polygons <- fetch_au(

baseurls = baseurls[4],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

)

try(

points <- points %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

try(

lines <- lines %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

try(

polygons <- polygons %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

final_features <- list(

"ATTAINS_catchments" = catchment_features,

"ATTAINS_points" = points,

"ATTAINS_lines" = lines,

"ATTAINS_polygons" = polygons

)

return(final_features)

# If area is small (< 6e+9 square meters), just use the bbox in one pull:

} else {

# FOR AOIs THAT ARE LESS THAN 6,000 sqkm, grab data in one go:

points_sf <- .data

bbox <- points_sf %>%

sf::st_bbox(.) %>%

# convert bounding box to characters

toString(.) %>%

# encode for use within the API URL

urltools::url_encode(.)

catchment_features <- fetch_bbox(baseurls = baseurls[1], sf_bbox = bbox)

try(

catchment_features <- catchment_features %>%

.[points_sf, ],

silent = TRUE

)

if (length(catchment_features) == 0 || is.null(catchment_features)) {

message("There are no ATTAINS features associated with your WQP observations.")

} else {

## GRABBING WATER TYPE:

# Use ATTAINS API to grab, for each assessment unit, its WaterType.

# Query the API in "chunks" so it doesn't break:

all_units <- unique(catchment_features$assessmentunitidentifier)

water_types <- grab_waterbody_type(all_units, chunk_size = 50)

try(catchment_features <- dplyr::left_join(catchment_features, water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")), silent = TRUE)

}

# If only interested in grabbing catchment data, return just the catchments

if (catchments_only == TRUE) {

return(list("ATTAINS_catchments" = catchment_features))

}

suppressMessages({

suppressWarnings({

# Otherwise, start grabbing the raw ATTAINS features that intersect those

# catchments

points <- NULL

lines <- NULL

polygons <- NULL

# Download associated point, line, and polygon features using catchment bbox

try(

points <- fetch_au(

baseurls = baseurls[2],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

),

silent = TRUE

)

try(

lines <- fetch_au(

baseurls = baseurls[3],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

),

silent = TRUE

)

try(

polygons <- fetch_au(

baseurls = baseurls[4],

assessment_unit_ids = unique(catchment_features$assessmentunitidentifier)

),

silent = TRUE

)

try(points <- points %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

try(lines <- lines %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

try(polygons <- polygons %>% dplyr::left_join(., water_types, by = c("assessmentunitidentifier" = "assessmentUnitIdentifier")),

silent = TRUE

)

final_features <- list(

"ATTAINS_catchments" = catchment_features,

"ATTAINS_points" = points,

"ATTAINS_lines" = lines,

"ATTAINS_polygons" = polygons

)

})

})

return(final_features)

}

# function settings that we ensure go back to their original settings

# after the function stops running:

original_s2 <- sf::sf_use_s2() # Store the original s2 setting first

suppressMessages(sf::sf_use_s2(FALSE))

original_timeout <- getOption("timeout")

options(timeout = 30000)

on.exit(options(timeout = original_timeout), add = TRUE)

on.exit(suppressMessages(suppressWarnings(sf::sf_use_s2(original_s2))), add = TRUE)

suppressMessages(suppressWarnings({

# If data is already spatial, just make sure it is in the right CRS

if (!is.null(.data) & inherits(.data, "sf")) {

if (sf::st_crs(.data)$epsg != 4326) {

geospatial_data <- .data %>%

sf::st_transform(4326)

} else {

geospatial_data <- .data

}

} else {

# ... Otherwise transform into a spatial object then do the same thing:

geospatial_data <- .data %>%

# convert dataframe to a spatial object

TADA_MakeSpatial(.data = ., crs = 4326) %>%

dplyr::mutate(geometry_join = geometry)

}

}))

# Reduce WQP data to unique coordinates

unique_sites <- dplyr::distinct(geospatial_data, geometry)

# If user wants HighRes NHD...

if (resolution %in% c("Hi", "hi")) {

suppressMessages(suppressWarnings({

# Map server for NHDPlus_HR that is used to download features:

nhd_plus_hr_url <- "https://hydro.nationalmap.gov/arcgis/rest/services/NHDPlus_HR/MapServer"

# bounding box of user's WQP data

wqp_bboxes <- unique_sites %>%

sf::st_buffer(1e-07) %>%

dplyr::rowwise() %>%

dplyr::mutate(bbox = purrr::map(geometry, sf::st_bbox)) %>%

sf::st_as_sfc()

# open the nhd_hr - which contains a bunch of layers

nhd_hr <- arcgislayers::arc_open(nhd_plus_hr_url)

# list the layers of the nhdhr object

# select the layer by id from the items list called above (10 is HR catchments)

nhd_hr_catchments <- arcgislayers::get_layer(nhd_hr, 10)

# use bboxes of the sites to return their associated catchments

nhd_catchments_stored <- vector("list", length = length(wqp_bboxes))

for (i in 1:length(wqp_bboxes)) {

try(

nhd_catchments_stored[[i]] <- arcgislayers::arc_select(nhd_hr_catchments,

filter_geom = wqp_bboxes[i],

crs = sf::st_crs(wqp_bboxes[i])

) %>%

sf::st_make_valid(),

silent = TRUE

)

}

nhd_catchments_stored <- nhd_catchments_stored %>%

purrr::keep(~ !is.null(.)) %>%

dplyr::bind_rows() %>%

dplyr::distinct()

try(nhd_catchments_stored <- nhd_catchments_stored %>%

dplyr::select(nhdplusid,

catchmentareasqkm = areasqkm

) %>%

dplyr::mutate(

NHD.nhdplusid = as.character(nhdplusid),

NHD.resolution = "HR",

NHD.catchmentareasqkm = as.numeric(catchmentareasqkm)

) %>%

dplyr::select(NHD.nhdplusid, NHD.resolution, NHD.catchmentareasqkm, geometry), silent = TRUE)

}))

# Empty version of the df will be returned if no associated catchments

# to avoid breaking downstream fxns reliant on catchment info.

if (nrow(nhd_catchments_stored) == 0 && "catchments" %in% features) {

message("No NHD HR features associated with your WQP observations.")

nhd_catchments_stored <- tibble::tibble(

NHD.nhdplusid = character(),

NHD.resolution = character(),

NHD.catchmentareasqkm = numeric()

)

}

if (nrow(nhd_catchments_stored) == 0 && !"catchments" %in% features) {

stop("No NHD HR features associated with your WQP observations.")

}

if (length(features) == 1 && features == "catchments") {

return(nhd_catchments_stored)

}

# Grab flowlines -

if ("flowlines" %in% features && nrow(nhd_catchments_stored) > 0) {

suppressMessages(suppressWarnings({

# use catchments to grab other NHD features

geospatial_aoi <- nhd_catchments_stored %>%

sf::st_as_sfc()

# select the layer by id from the items list (3 is HR flowlines)

nhd_hr_flowlines <- arcgislayers::get_layer(nhd_hr, 3)

# use catchments to return associated flowlines

nhd_flowlines_stored <- vector("list", length = length(geospatial_aoi))

for (i in 1:length(geospatial_aoi)) {

try(

nhd_flowlines_stored[[i]] <- arcgislayers::arc_select(nhd_hr_flowlines,

filter_geom = geospatial_aoi[i],

crs = sf::st_crs(geospatial_aoi[i])

) %>%

sf::st_make_valid(),

silent = TRUE

)

# so all returned meta data binds properly, must transform all columns into characters,

# EXCEPT for the geometry column:

try(

geometry_col <- sf::st_geometry(nhd_flowlines_stored[[i]]),

silent = TRUE

)

try(

nhd_flowlines_stored[[i]] <- nhd_flowlines_stored[[i]] %>%

dplyr::mutate(dplyr::across(dplyr::where(~ !identical(., geometry_col)), ~ as.character(.))),

silent = TRUE

)

}

nhd_flowlines_stored <- nhd_flowlines_stored %>%

purrr::keep(~ !is.null(.)) %>%

purrr::keep(~ !is.character(.)) %>%

dplyr::bind_rows() %>%

dplyr::distinct()

}))

if (length(features) == 1 && features == "flowlines") {

if (length(nhd_flowlines_stored) == 0 || is.null(nhd_flowlines_stored)) {

message("There are no NHD flowlines associated with your WQP observations.")

}

return(nhd_flowlines_stored)

}

if (length(nhd_flowlines_stored) == 0 || is.null(nhd_flowlines_stored)) {

message("There are no NHD flowlines associated with your WQP observations.")

}

}

# Grab waterbodies -

if ("waterbodies" %in% features & nrow(nhd_catchments_stored) > 0) {

suppressMessages(suppressWarnings({

geospatial_aoi <- nhd_catchments_stored %>%

sf::st_as_sfc()

# select the layer by id from the items list called above (9 is HR waterbodies)

nhd_hr_waterbodies <- arcgislayers::get_layer(nhd_hr, 9)

# use catchments to return associated waterbodies

nhd_waterbodies_stored <- vector("list", length = length(geospatial_aoi))

for (i in 1:length(geospatial_aoi)) {

try(

nhd_waterbodies_stored[[i]] <- arcgislayers::arc_select(nhd_hr_waterbodies,

# where = query,

filter_geom = geospatial_aoi[i],

crs = sf::st_crs(geospatial_aoi[i])

) %>%

sf::st_make_valid(),

silent = TRUE

)

# so all returned meta data binds properly, must transform all columns into characters,

# EXCEPT for the geometry column:

try(

geometry_col <- sf::st_geometry(nhd_waterbodies_stored[[i]]),

silent = TRUE

)

try(

nhd_waterbodies_stored[[i]] <- nhd_waterbodies_stored[[i]] %>%

dplyr::mutate(dplyr::across(dplyr::where(~ !identical(., geometry_col)), ~ as.character(.))),

silent = TRUE

)

}

nhd_waterbodies_stored <- nhd_waterbodies_stored %>%

purrr::keep(~ !is.null(.)) %>%

purrr::keep(~ !is.character(.)) %>%

dplyr::bind_rows() %>%

dplyr::distinct()

}))

if (length(features) == 1 && features == "waterbodies") {

if (length(nhd_waterbodies_stored) == 0 || is.null(nhd_waterbodies_stored)) {

message("There are no NHD waterbodies associated with your WQP observations.")

}

return(nhd_waterbodies_stored)

}

if (length(nhd_waterbodies_stored) == 0 || is.null(nhd_waterbodies_stored)) {

message("There are no NHD waterbodies associated with your WQP observations.")

}

}

# Combinations of features selected, and what they return:

if (length(features) == 2 && "catchments" %in% features && "flowlines" %in% features) {

nhd_list <- list(

"NHD_catchments" = nhd_catchments_stored,

"NHD_flowlines" = nhd_flowlines_stored

)

return(nhd_list)

} else if (length(features) == 2 && "catchments" %in% features && "waterbodies" %in% features) {

nhd_list <- list(

"NHD_catchments" = nhd_catchments_stored,

"NHD_waterbodies" = nhd_waterbodies_stored

)

return(nhd_list)

} else if (length(features) == 2 && "flowlines" %in% features && "waterbodies" %in% features) {

nhd_list <- list(

"NHD_flowlines" = nhd_flowlines_stored,

"NHD_waterbodies" = nhd_waterbodies_stored

)

return(nhd_list)

} else if (length(features) == 3 && "catchments" %in% features && "flowlines" %in% features && "waterbodies" %in% features) {

nhd_list <- list(

"NHD_catchments" = nhd_catchments_stored,

"NHD_flowlines" = nhd_flowlines_stored,

"NHD_waterbodies" = nhd_waterbodies_stored

)

} else {

stop("Please select between 'catchments', 'flowlines', 'waterbodies', or any combination for `feature` argument.")

}

# If user wants NHDPlus V2...

} else if (resolution %in% c("Med", "med")) {

suppressMessages(suppressWarnings({

nhd_catchments <- vector("list", length = nrow(unique_sites))

for (i in 1:nrow(unique_sites)) {

# Use {nhdplusTools} to grab associated catchments...

try(

nhd_catchments[[i]] <- nhdplusTools::get_nhdplus(AOI = unique_sites[i, ], realization = "catchment") %>%

sf::st_make_valid() %>%

dplyr::select(

comid = featureid,

catchmentareasqkm = areasqkm

) %>%

dplyr::mutate(

NHD.comid = as.character(comid),

NHD.resolution = "nhdplusV2",

NHD.catchmentareasqkm = as.numeric(catchmentareasqkm)

) %>%

dplyr::select(NHD.comid, NHD.resolution, NHD.catchmentareasqkm, geometry),

silent = TRUE

)

}

nhd_catchments <- nhd_catchments %>%

purrr::keep(~ !is.null(.))

try(nhd_catchments <- dplyr::bind_rows(nhd_catchments) %>%

dplyr::distinct(), silent = TRUE)

# if NHD catchments are not in the correct CRS, transform them

try(if (sf::st_crs(nhd_catchments) != sf::st_crs(geospatial_data)) {

nhd_catchments <- nhd_catchments %>%

sf::st_transform(sf::st_crs(geospatial_data)$epsg)

}, silent = TRUE)

}))

if (nrow(nhd_catchments) == 0 && "catchments" %in% features) {

message("No NHDPlus V2 features associated with your WQP observations.")

nhd_catchments <- tibble::tibble(

NHD.comid = character(),

NHD.resolution = character(),

NHD.catchmentareasqkm = numeric()

)

}

if (nrow(nhd_catchments) == 0 && !"catchments" %in% features) {

stop("No NHDPlus V2 features associated with your WQP observations.")

}

if (length(features) == 1 && features == "catchments") {

return(nhd_catchments)

}

# Grab flowlines -

if ("flowlines" %in% features && nrow(nhd_catchments) > 0) {

suppressMessages(suppressWarnings({

nhd_flowlines <- vector("list", length = nrow(nhd_catchments))

# use catchments to grab other NHD features:

unique_sites <- nhd_catchments

for (i in 1:nrow(unique_sites)) {

# Use {nhdplusTools} to grab associated flowlines...

try(

nhd_flowlines[[i]] <- nhdplusTools::get_nhdplus(AOI = unique_sites[i, ], realization = "flowline") %>%

sf::st_make_valid(),

silent = TRUE

)

try(geometry_col <- sf::st_geometry(nhd_flowlines[[i]]),

silent = TRUE

)