diff --git a/R/modules.R b/R/modules.R index 5beb6f1..f65f92c 100644 --- a/R/modules.R +++ b/R/modules.R @@ -141,6 +141,15 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -154,11 +163,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - } leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data %>% @@ -185,6 +189,16 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -198,11 +212,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - } leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data%>% @@ -228,6 +237,16 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + + } values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -241,12 +260,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - - } leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data%>% @@ -272,6 +285,15 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -285,11 +307,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - } leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data%>% @@ -315,6 +332,15 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -328,11 +354,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - } leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data%>% @@ -358,6 +379,16 @@ nhdplusMod <- function(input, output, session, values){ values$hydro_data <- . + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } + values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) @@ -371,12 +402,6 @@ nhdplusMod <- function(input, output, session, values){ req(class(values$hydro_data)[[1]] == 'sf') - if(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION') { - - values$hydro_data <- convert_sf_geocollection(values$hydro_data) - - } - leaflet::leafletProxy("leaf_map", session) %>% leaflet::addPolygons(data = values$hydro_data%>% sf::st_transform(crs = 4326,proj4string = "+init=epsg:4326"), popup = paste0("

", @@ -399,6 +424,17 @@ nhdplusMod <- function(input, output, session, values){ }) %...>% { values$hydro_data <- . + + if(any(sf::st_geometry_type(values$hydro_data) == 'GEOMETRYCOLLECTION')) { + + values$hydro_data <- values$hydro_data %>% + dplyr::mutate(rowid = dplyr::row_number()) %>% + split(.$rowid) %>% + purrr::map(~convert_sf_geocollection(.x)) %>% + dplyr::bind_rows() + + } + values$out <- list(values$hydro_data) names(values$out) <- paste0(input$location_map, '_',sample(1:10000,size = 1, replace = T)) diff --git a/R/nldi.R b/R/nldi.R index 50a5d62..ee7fe61 100644 --- a/R/nldi.R +++ b/R/nldi.R @@ -182,14 +182,14 @@ server = function(input, output, session){ color = "black", fillOpacity = 0, weight = 3, - opacity = 1) + opacity = 1, popup = paste0("", "DA acres: ", "", scales::comma(as.numeric(round(units::set_units(sf::st_area(values$nldi_data()[[2]]), acres), 1)),1), " Acres")) map_nldi <- leaflet::addPolylines(map_nldi, data = values$nldi_data()[[1]], color = "blue", weight = 3, opacity = 1) - }else if (input$location_map == 'local') { + } else if (input$location_map == 'local') { values$nldi_data <- reactive(get_NLDI_catchments(data_sf,method = 'local')) @@ -199,7 +199,8 @@ server = function(input, output, session){ color = "black", fillOpacity = 0, weight = 3, - opacity = 1) + opacity = 1, popup = paste0("", "DA acres: ", "", scales::comma(as.numeric(round(units::set_units(sf::st_area(values$nldi_data()[[2]]), acres), 1)),1), " Acres", + "
", "", "Total length of Tribs: ", "", round(sum(units::set_units(sf::st_length(values$nldi_data()[[1]]), mi)), 1), " Miles")) map_nldi <- leaflet::addPolylines(map_nldi, data = values$nldi_data()[[1]], diff --git a/R/utils.R b/R/utils.R index d185994..6ecd9b2 100644 --- a/R/utils.R +++ b/R/utils.R @@ -2,7 +2,7 @@ #' Get NLDI #' #' @description This function grabs the upstream tributaries, upstream main stream and basin boundary using -#' the NLDI API. It then combines the NLDI zonal stats to the basin boundary shape, i.e. 'TOT' is the 'total' basin zonal statistic. +#' the NLDI API. #' #' @param point A sf point. #' @noRd @@ -10,45 +10,19 @@ #' get_NLDI <- function(point){ - clat <- point$geometry[[1]][[2]] - clng <- point$geometry[[1]][[1]] - - ids <- paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/position?coords=POINT%28", - clng,"%20", clat, "%29") - - error_ids <- httr::GET(url = ids, - httr::write_disk(path = file.path(tempdir(), - "nld_tmp.json"),overwrite = TRUE)) - - nld <- jsonlite::fromJSON(file.path(tempdir(),"nld_tmp.json")) + comid <- nhdplusTools::discover_nhdplus_id(point) - - nldiURLs <- list(site_data = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/position?coords=POINT%28", - clng,"%20", clat, "%29"), - basin_boundary = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",nld$features$properties$identifier,"/basin"), - UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",nld$features$properties$identifier,"/navigation/UT/flowlines?distance=999"), - UM = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",nld$features$properties$identifier,"/navigation/UM/flowlines?distance=999")) + nldiURLs <- list(basin_boundary = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/basin"), + UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UT/flowlines?distance=999"), + UM = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UM/flowlines?distance=999")) nldi_data <- list() for(n in names(nldiURLs)) { nldi_data[n] <- list(sf::read_sf(nldiURLs[n][[1]])) - print(paste(n, "is of class", class(nldi_data[[n]]), "and has", nrow(nldi_data[[n]]), "features")) } - total_characteristic <- data.frame(COMID = nldi_data$site_data$comid) %>% mutate(ID = dplyr::row_number()) %>% - group_by(ID) %>% - tidyr::nest() %>% - mutate(chars = purrr::map(data, ~nhdplusTools::get_nldi_characteristics(list(featureSource = "comid", featureID = as.character(.$COMID)), - type = 'total'))) %>% tidyr::unnest(c(data, chars)) %>% tidyr::unnest(c(chars)) %>% dplyr::ungroup() %>% - dplyr::select(COMID, characteristic_id, characteristic_value) %>% - tidyr::pivot_wider(names_from = "characteristic_id", values_from = "characteristic_value") %>% - dplyr::mutate(dplyr::across(is.character, as.numeric)) - - - nldi_data[['basin_boundary']] <- nldi_data[['basin_boundary']] %>% - cbind(total_characteristic) nldi_data } @@ -69,54 +43,34 @@ get_NLDI <- function(point){ #' get_NLDI_catchments <- function(point, type = 'local', method = 'all'){ - clat <- point$geometry[[1]][[2]] - clng <- point$geometry[[1]][[1]] - - ids <- paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/position?coords=POINT%28", - clng,"%20", clat, "%29") - - error_ids <- httr::GET(url = ids, - httr::write_disk(path = file.path(tempdir(), - "nld_tmp.json"),overwrite = TRUE)) + comid <- tryCatch(expr = {nhdplusTools::discover_nhdplus_id(point)}, + error = function(e) {NA}) - nld <- jsonlite::fromJSON(file.path(tempdir(),"nld_tmp.json")) + if(is.na(comid)){stop('COMID not found')} if(method == 'all'){ - nldiURLs <- list(UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",nld$features$properties$identifier,"/navigation/UT/flowlines?distance=999")) + nldiURLs <- list(UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UT/flowlines?distance=999")) } else if (method == 'local'){ - nldiURLs <- list(UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",nld$features$properties$identifier,"/navigation/UT/flowlines?distance=0")) + nldiURLs <- list(UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UT/flowlines?distance=0")) } nldi_data <- list() for(n in names(nldiURLs)) { nldi_data[n] <- list(sf::read_sf(nldiURLs[n][[1]])) - print(paste(n, "is of class", class(nldi_data[[n]]), "and has", nrow(nldi_data[[n]]), "features")) } nldi_outlets <- nldi_data$UT$nhdplus_comid - nldi_catch <- nhdplusTools::get_nhdplus(comid = nldi_outlets, - realization = 'catchment') - - local_characteristic <- data.frame(COMID = nldi_outlets) %>% mutate(ID = dplyr::row_number()) %>% - group_by(ID) %>% - tidyr::nest() %>% - mutate(chars = purrr::map(data, ~nhdplusTools::get_nldi_characteristics(list(featureSource = "comid", featureID = as.character(.$COMID)), - type = type))) %>% tidyr::unnest(c(data, chars)) %>% tidyr::unnest(c(chars)) %>% dplyr::ungroup() %>% - dplyr::select(COMID, characteristic_id, characteristic_value) %>% - tidyr::pivot_wider(names_from = "characteristic_id", values_from = "characteristic_value") %>% - dplyr::rename(featureid = 'COMID') %>% - dplyr::mutate(dplyr::across(is.character, as.numeric))%>% - dplyr::mutate(featureid = as.integer(featureid)) - - - nldi_catch <- nldi_catch %>% - dplyr::left_join(local_characteristic, by = c('featureid')) + nldi_catch <- suppressMessages(purrr::map(nldi_outlets,~nhdplusTools::get_nhdplus(comid = ., + realization = 'catchment'))) %>% + dplyr::bind_rows() %>% + sf::st_as_sf(crs = 4326) final_data <- list(nldi_data$UT, nldi_catch) } + #' Base Map #' #' @description A generic leaflet base map used in the shiny apps. @@ -237,7 +191,7 @@ convert_sf_geocollection <- function(x) { data <- x %>% sf::st_drop_geometry() - sf::st_union(sf::st_as_sf(sf::st_collection_extract(x %>% sf::st_union(), c('POLYGON')))) %>% + sf::st_union(sf::st_as_sf(sf::st_collection_extract(x %>% sf::st_zm() %>% sf::st_union(), c('POLYGON')))) %>% sf::st_as_sf() %>% dplyr::bind_cols(data) %>% rename_geometry('geometry') @@ -250,7 +204,7 @@ convert_sf_geocollection <- function(x) { #' @param name character. #' #' @return A sf object with a renamed geometry column. -#' @notes This function was grabbed from [stack overflow](https://gis.stackexchange.com/questions/386584/sf-geometry-column-naming-differences-r) from the legend spacedman. +#' @note This function was grabbed from [stack overflow](https://gis.stackexchange.com/questions/386584/sf-geometry-column-naming-differences-r) from the legend spacedman. rename_geometry <- function(g, name){ current = attr(g, "sf_column") names(g)[names(g)==current] = name diff --git a/man/rename_geometry.Rd b/man/rename_geometry.Rd index 4f328cd..ef10aee 100644 --- a/man/rename_geometry.Rd +++ b/man/rename_geometry.Rd @@ -17,3 +17,6 @@ A sf object with a renamed geometry column. \description{ Rename Geometry Column } +\note{ +This function was grabbed from \href{https://gis.stackexchange.com/questions/386584/sf-geometry-column-naming-differences-r}{stack overflow} from the legend spacedman. +} diff --git a/tests/testthat/test-nldi.R b/tests/testthat/test-nldi.R new file mode 100644 index 0000000..241e48e --- /dev/null +++ b/tests/testthat/test-nldi.R @@ -0,0 +1,55 @@ +test_that("testing get catchment characteristics from nhdplusTools", { + + point <- dplyr::tibble(Long = -115.2312,Lat = 48.83037) + + point <- sf::st_as_sf(point, coords = c("Long", "Lat"), crs = 4326) + + comid <- nhdplusTools::discover_nhdplus_id(point) + + catchment_char <- nhdplusTools::get_catchment_characteristics("CAT_BFI", ids = comid) + + testthat::expect_equal(catchment_char$characteristic_value[[1]], 71) + + + catchment_char <- nhdplusTools::get_catchment_characteristics("TOT_BFI", ids = comid) + + testthat::expect_equal(catchment_char$characteristic_value[[1]], 72.07) + + catchment_char <- nhdplusTools::get_catchment_characteristics("ACC_BFI", ids = comid) + + testthat::expect_equal(catchment_char$characteristic_value[[1]], 72.07) + + +}) + + +test_that("testing upstream main and tribs from nhdplusTools", { + + point <- dplyr::tibble(Long = -115.2312,Lat = 48.83037) + + point <- sf::st_as_sf(point, coords = c("Long", "Lat"), crs = 4326) + + comid <- nhdplusTools::discover_nhdplus_id(point) + + nldiURLs <- list(site_data = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/position?coords=POINT%28", + clng,"%20", clat, "%29"), + basin_boundary = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/basin"), + UT = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UT/flowlines?distance=999"), + UM = paste0("https://labs.waterdata.usgs.gov/api/nldi/linked-data/comid/",comid,"/navigation/UM/flowlines?distance=999")) + + nldi_data <- list() + + + for(n in names(nldiURLs)) { + nldi_data[n] <- list(sf::read_sf(nldiURLs[n][[1]])) + } + + testthat::expect_equal(length(nldi_data), 4) + testthat::expect_equal(length(nldi_data$site_data), 6) + testthat::expect_equal(length(nldi_data$basin_boundary), 1) + testthat::expect_equal(nrow(nldi_data$UT), 71) + testthat::expect_equal(nrow(nldi_data$UM), 27) + + + +})