Merge 8fdb415 into 93d5325

R/subset_nhdplus.R

@@ -21,6 +21,7 @@
 #' if TRUE only the flowline network and attributes will be returned
 #' @param streamorder integer only streams of order greater than or equal will be downloaded.
 #' Not implemented for local data.
+#' @param out_prj character override the default output CRS of NAD83 lat/lon (EPSG:4269)
 #' @details If \code{\link{stage_national_data}} has been run in the current
 #' session, this function will use the staged national data automatically.
 #'
@@ -116,7 +117,7 @@
 
 subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NULL, bbox = NULL,
                            simplified = TRUE, overwrite = FALSE, return_data = TRUE, status = TRUE,
-                           flowline_only = NULL, streamorder = NULL) {
+                           flowline_only = NULL, streamorder = NULL, out_prj = 4269) {
 
   if(is.null(flowline_only)) {
     if(!is.null(nhdplus_data) && nhdplus_data == "download") {
@@ -126,6 +127,8 @@ subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NUL
     }
   }
 
+  if(!is.null(comids) && length(comids) == 0) stop("comids must be NULL or non-empty")
+
   if (status) message("All intersections performed in latitude/longitude.")
 
   if(any(bbox > 180 | bbox < -180)) stop("invalid bbox entry")
@@ -162,20 +165,27 @@ subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NUL
 
     out_list <- c(get_flowline_subset(nhdplus_data, comids,
                                       output_file, paths$fline_path,
-                                      status))
+                                      status, out_prj))
 
     if(!flowline_only) {
+      tryCatch({
+        out_list <- c(out_list, get_catchment_subset(nhdplus_data, comids,
+                                                     output_file, simplified,
+                                                     paths$catchment_path, status,
+                                                     out_prj))
 
-      out_list <- c(out_list, get_catchment_subset(nhdplus_data, comids,
-                                                   output_file, simplified,
-                                                   paths$catchment_path, status))
+        catch_layer <- get_catchment_layer_name(simplified, nhdplus_data)
 
-      catch_layer <- get_catchment_layer_name(simplified, nhdplus_data)
+        envelope <- sf::st_transform(sf::st_as_sfc(sf::st_bbox(out_list[[catch_layer]])),
+                                     4326)
 
-      envelope <- sf::st_transform(sf::st_as_sfc(sf::st_bbox(out_list[[catch_layer]])),
-                                   4326)
+        intersection_names <- c("NHDArea", "NHDWaterbody")
+      }, error = function(e) {
+        warning(paste("error getting catchment from nhdplus_data\n", e))
+      })
+
+      if(!exists("intersection_names")) intersection_names <- c()
 
-      intersection_names <- c("NHDArea", "NHDWaterbody")
     } else {
       intersection_names <- c()
     }
@@ -209,12 +219,16 @@ subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NUL
       layer <- sf::st_transform(envelope, 4326) %>%
         get_nhdplus_bybox(layer = tolower(layer_name), streamorder = streamorder)
 
-      if(return_data) {
-        out_list[layer_name] <- list(layer)
-      }
+      if(nrow(layer) > 0) {
+        layer <- check_valid(layer, out_prj)
 
-      if(!is.null(output_file)) {
-        sf::write_sf(clean_bbox(layer), output_file, layer_name)
+        if(return_data) {
+          out_list[layer_name] <- list(layer)
+        }
+
+        if(!is.null(output_file)) {
+          sf::write_sf(clean_bbox(layer), output_file, layer_name)
+        }
       }
     }
 
@@ -233,7 +247,8 @@ subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NUL
                                          data_path = nhdplus_data,
                                          envelope = envelope,
                                          output_file = output_file,
-                                         status), intersection_names))
+                                         status = status,
+                                         out_prj = out_prj), intersection_names))
   }
 
   if(return_data) return(out_list)
@@ -242,7 +257,7 @@ subset_nhdplus <- function(comids = NULL, output_file = NULL, nhdplus_data = NUL
 }
 
 intersection_write <- function(layer_name, data_path, envelope,
-                               output_file, status) {
+                               output_file, status, out_prj) {
   out_list <- list()
 
   if (status) message(paste("Reading", layer_name))
@@ -253,11 +268,20 @@ intersection_write <- function(layer_name, data_path, envelope,
   try(intersection_test <- suppressMessages(sf::st_intersects(
     sf::st_transform(layer, 4326), envelope)), silent = TRUE)
 
-  out <- dplyr::filter(layer, lengths(intersection_test) > 0)
+  found <- lengths(intersection_test)
+
+  if(length(found) > 0) {
+    out <- dplyr::filter(layer, found > 0)
+  } else {
+    out <- data.frame()
+  }
 
   if (nrow(out) > 0) {
+
+    out <- check_valid(out, out_prj)
+
     if (status) message(paste("Writing", layer_name))
-    if(is.null(output_file)) {
+    if (is.null(output_file)) {
       return(out)
     } else {
       sf::write_sf(clean_bbox(out), output_file, layer_name)
@@ -415,7 +439,7 @@ get_staged_data <- function(nhdplus_data) {
 #' @title Get subset of flowline data later.
 #' @noRd
 get_flowline_subset <- function(nhdplus_data, comids, output_file,
-                                fline_path, status) {
+                                fline_path, status, out_prj) {
 
   layer_name <- get_flowline_layer_name(nhdplus_data)
 
@@ -443,20 +467,24 @@ get_flowline_subset <- function(nhdplus_data, comids, output_file,
     fline <- dplyr::filter(fline, .data$COMID %in% comids)
   }
 
+  fline <- check_valid(fline, out_prj)
+
   if (status) message(paste("Writing", layer_name))
 
   if(!is.null(output_file)) {
     sf::write_sf(clean_bbox(fline), output_file, layer_name)
   }
   out <- list()
   out[layer_name] <- list(fline)
+
   return(out)
 }
 
 #' @title Get subset of catchment data layer.
 #' @noRd
 get_catchment_subset <- function(nhdplus_data, comids, output_file,
-                                 simplified, catchment_path, status) {
+                                 simplified, catchment_path, status,
+                                 out_prj) {
 
   layer_name <- get_catchment_layer_name(simplified, nhdplus_data)
 
@@ -479,6 +507,8 @@ get_catchment_subset <- function(nhdplus_data, comids, output_file,
 
   }
 
+  catchment <- check_valid(catchment, out_prj)
+
   if (status) message(paste("Writing", layer_name))
 
   if(!is.null(output_file)) {
@@ -493,6 +523,29 @@ clean_bbox <- function(x) {
   if("bbox" %in% names(x) && class(x$bbox[1]) == "list") {
     x$bbox <- sapply(x$bbox, paste, collapse = ",")
   }
+
+  return(x)
+}
+
+check_valid <- function(x, out_prj) {
+
+  x <- sf::st_zm(x)
+
+  if (!all(sf::st_is_valid(x))) {
+
+    message("Found invalid geometry, attempting to fix.")
+
+    try(x <- sf::st_make_valid(x))
+  }
+
+  if (any(grepl("POLYGON", class(sf::st_geometry(x))))) {
+    suppressMessages(suppressWarnings(x <- sf::st_buffer(x, 0)))
+  }
+
+  if (sf::st_crs(x) != sf::st_crs(out_prj)) {
+    x <- sf::st_transform(x, out_prj)
+  }
+
   return(x)
 }
 

tests/testthat/test_get_nldi.R

@@ -53,11 +53,31 @@ test_that("navigation works", {
 
   nav3 <- navigate_nldi(nldi_feature = nldi_nwis,
                         mode = "upstreamMain",
-                        data_source = "flowline",
+                        data_source = "flowlines",
                         distance_km = 10)
 
   expect_is(sf::st_geometry(nav3), "sfc_LINESTRING")
 
+  expect_warning(nav3 <- navigate_nldi(nldi_feature = nldi_nwis,
+                                       mode = "upstreamMain",
+                                       data_source = "flowline",
+                                       distance_km = 10),
+                 "data source specified as flowline or '' is deprecated")
+
+  nav <- navigate_nldi(nldi_feature = nldi_nwis,
+                       mode = "https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-08279500/navigation/UM",
+                       data_source = "dumb",
+                       distance_km = 1)
+
+  expect_equal(nav, dplyr::tibble())
+
+  nav <- navigate_nldi(nldi_feature = nldi_nwis,
+                       mode = "https://labs.waterdata.usgs.gov/api/nldi/linked-data/nwissite/USGS-08279500/navigation/UM",
+                       data_source = "nwissite",
+                       distance_km = 1)
+
+  expect("sf" %in% class(nav), "expected an sf data.frame")
+
   # expect_equal(navigate_nldi(list(featureSource = "wqp",
   #                                 featureID = "TCEQMAIN-16638"),
   #                            mode = "upstreamMain",

tests/testthat/test_subset.R

@@ -18,6 +18,8 @@ test_that("subset errors", {
                               "must provide comids or bounding box"))
 
   nhdplus_path("../NHDPlusV21_National_Seamless.gdb")
+
+  expect_error(subset_nhdplus(comids = integer(0), nhdplus_data = "download"), "comids must be NULL or non-empty")
 })
 
 test_that("subset runs as expected", {
@@ -98,6 +100,12 @@ test_that("subset runs as expected", {
   check_layers <- function() {
     expect_equal(nrow(sf::read_sf(out_file, "CatchmentSP")), 4)
     expect_equal(nrow(sf::read_sf(out_file, "NHDWaterbody")), 1)
+    expect_true(sf::st_crs(sf::read_sf(out_file, "CatchmentSP")) ==
+                  sf::st_crs(4269))
+    expect_true(sf::st_crs(sf::read_sf(out_file, "NHDWaterbody")) ==
+                  sf::st_crs(4269))
+    expect_true(sf::st_crs(sf::read_sf(out_file, "NHDFlowline_Network")) ==
+                  sf::st_crs(4269))
   }
 
   check_layers()
@@ -124,7 +132,8 @@ test_that("subset runs as expected", {
                        output_file = out_file,
                        nhdplus_data = "download",
                        overwrite = FALSE,
-                       status = FALSE, flowline_only = FALSE)
+                       status = FALSE,
+                       flowline_only = FALSE)
 
   check_layers()
 
@@ -171,7 +180,7 @@ test_that("subset by bounding box", {
   sample_data <- system.file("extdata/sample_natseamless.gpkg",
                              package = "nhdplusTools")
 
-  bbox <- st_bbox(c(xmin = -89.4, ymin = 43, xmax = -89.3, ymax = 43.1), crs = st_crs(4326))
+  bbox <- sf::st_bbox(c(xmin = -89.4, ymin = 43, xmax = -89.3, ymax = 43.1), crs = sf::st_crs(4326))
 
   expect_error(subset_nhdplus(bbox = (bbox + c(-200, -200, 200, 200)),
                               nhdplus_data = sample_data,
@@ -198,6 +207,13 @@ test_that("subset by bounding box", {
   check_layers <- function() {
     expect_equal(nrow(sf::read_sf(out_file, "CatchmentSP")), 66)
     expect_equal(nrow(sf::read_sf(out_file, "NHDWaterbody")), 12)
+    expect_true(sf::st_crs(sf::read_sf(out_file, "CatchmentSP")) ==
+                 sf::st_crs(4269))
+    expect_true(sf::st_crs(sf::read_sf(out_file, "NHDWaterbody")) ==
+                 sf::st_crs(4269))
+    expect_true(sf::st_crs(sf::read_sf(out_file, "NHDFlowline_Network")) ==
+                 sf::st_crs(4269))
+
   }
 
   check_layers()
@@ -293,3 +309,21 @@ test_that("by rpu", {
   expect(nrow(subset_rpu(sample_flines, rpu = "07b")), 267)
   expect(nrow(subset_rpu(sample_flines, rpu = "07b", run_make_standalone = TRUE)), 267)
 })
+
+
+test_that("projection check", {
+
+  skip_on_cran()
+
+  out <- tempfile(fileext = ".gpkg")
+
+  unlink(out, recursive = TRUE)
+
+  mr <- nhdplusTools::plot_nhdplus(list(13293970), gpkg = out,
+                                   nhdplus_data = out,
+                                   overwrite = FALSE, actually_plot = FALSE)
+
+
+  expect_true(sf::st_crs(mr$flowline) == sf::st_crs(4269))
+
+})