Support 3DHP_all service.

[dblodgett-usgs]

Feature layers from: https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer should be accessible with a new subset function.

It should support one or all layers and return in more or less the same format as the existing geoserver utility function.

This will end up being useful for accessing data for a variety of applications including indexing.

[dblodgett-usgs]

Have a very basic implementation based loosely on query_usgs_geoserver and patterns that @cheginit has implemented in hyRiver. Needs some evolving, but the basics of requesting all the OBJECTIDS for a query then iterating through the ids 1000 at a time is in place and seems pretty workable. At this point it will query by id3dhp or take a canned "where" query that can contain various SQL filters.

[dblodgett-usgs]

Good progress -- this will be the example.

library(nhdplusTools)
AOI <- sf::st_as_sfc(sf::st_bbox(c(xmin = -89.56684, ymin = 42.99816,
                                   xmax = -89.24681, ymax = 43.17192),
                                 crs = "+proj=longlat +datum=WGS84 +no_defs"))

# get flowlines and hydrolocations
flowlines <- get_3dhp(AOI = AOI, type = "flowline")
hydrolocation <- get_3dhp(AOI = AOI, type = "hydrolocation")
waterbody <- get_3dhp(AOI = AOI, type = "waterbody")

plot(sf::st_geometry(waterbody), col = "lightblue", border = "lightgrey")
plot(sf::st_geometry(flowlines), col = "blue", add = TRUE)
plot(sf::st_geometry(hydrolocation), col = "grey", pch = "+", add = TRUE)

# given universamreferenceid (reachcodes), can query for them but only
# for hydrolocations. This is useful for looking up mainstem ids.

get_3dhp(universalreferenceid = unique(hydrolocation$universalreferenceid),
         type = "hydrolocation")
#> Request failed [504]. Retrying in 1.4 seconds...
#> Request failed [504]. Retrying in 3.6 seconds...
#> No encoding supplied: defaulting to UTF-8.
#> Error in split.default(ids, ceiling(seq_along(ids)/chunk_size)): first argument must be a vector

# given mainstem ids from any source, can query for them in ids.

get_3dhp(ids = unique(flowlines$mainstemid), type = "flowline")
#> Request failed [504]. Retrying in 1 seconds...
#> Request failed [504]. Retrying in 2.2 seconds...
#> Getting features 0 to 1000 of 27874
#> Getting features 1000 to 2000 of 27874
#> Getting features 2000 to 3000 of 27874
#> Getting features 3000 to 4000 of 27874
#> Getting features 4000 to 5000 of 27874
#> Getting features 5000 to 6000 of 27874
#> Getting features 6000 to 7000 of 27874
#> Getting features 7000 to 8000 of 27874
#> Getting features 8000 to 9000 of 27874
#> Getting features 9000 to 10000 of 27874
#> Getting features 10000 to 11000 of 27874
#> Getting features 11000 to 12000 of 27874
#> Getting features 12000 to 13000 of 27874
#> Getting features 13000 to 14000 of 27874
#> Getting features 14000 to 15000 of 27874
#> Getting features 15000 to 16000 of 27874
#> Getting features 16000 to 17000 of 27874
#> Getting features 17000 to 18000 of 27874
#> Getting features 18000 to 19000 of 27874
#> Getting features 19000 to 20000 of 27874
#> Getting features 20000 to 21000 of 27874
#> Getting features 21000 to 22000 of 27874
#> Getting features 22000 to 23000 of 27874
#> Getting features 23000 to 24000 of 27874
#> Getting features 24000 to 25000 of 27874
#> Getting features 25000 to 26000 of 27874
#> Getting features 26000 to 27000 of 27874
#> Getting features 27000 to 27874 of 27874
#> Simple feature collection with 27874 features and 35 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -124.2935 ymin: 25.37482 xmax: -66.77981 ymax: 50.88568
#> Geodetic CRS:  WGS 84
#> # A tibble: 27,874 × 36
#>    id3dhp   featuredate mainstemid                gnisid gnisidlabel featuretype
#>    <chr>          <dbl> <chr>                      <int> <chr>             <int>
#>  1 100SO  1694649600000 https://geoconnex.us/us…      NA <NA>                  5
#>  2 104AU  1694649600000 https://geoconnex.us/us…      NA <NA>                  5
#>  3 104C6  1694649600000 https://geoconnex.us/us…      NA <NA>                  1
#>  4 107GG  1694649600000 https://geoconnex.us/us… 1209828 Minister C…           5
#>  5 109C7  1694649600000 https://geoconnex.us/us…      NA <NA>                  1
#>  6 109I4  1694649600000 https://geoconnex.us/us…      NA <NA>                  1
#>  7 10CDG  1694649600000 https://geoconnex.us/us…      NA <NA>                  1
#>  8 10CE7  1694649600000 https://geoconnex.us/us…      NA <NA>                  1
#>  9 10DVV  1694649600000 https://geoconnex.us/us…      NA <NA>                  5
#> 10 10E2X  1694649600000 https://geoconnex.us/us…      NA <NA>                  3
#> # ℹ 27,864 more rows
#> # ℹ 30 more variables: featuretypelabel <chr>, lengthkm <dbl>,
#> #   waterbodyid3dhp <chr>, flowdirection <int>, flowdirectionlabel <chr>,
#> #   onsurface <int>, onsurfacelabel <chr>, catchmentid3dhp <chr>,
#> #   flowpathid3dhp <chr>, streamlevel <chr>, startflag <chr>,
#> #   terminalflag <chr>, streamorder <chr>, streamcalculator <chr>,
#> #   hydrosequence <chr>, dnhydrosequence <chr>, uphydrosequence <chr>, …

CO <- get_3dhp(ids = "https://geoconnex.us/ref/mainstems/29559",
               type = "flowline")
#> Getting features 0 to 1000 of 4395
#> Getting features 1000 to 2000 of 4395
#> Getting features 2000 to 3000 of 4395
#> Getting features 3000 to 4000 of 4395
#> Getting features 4000 to 4395 of 4395

plot(sf::st_geometry(CO), col = "blue")

# get all the waterbodies along the CO river
CO_wb <- get_3dhp(ids = unique(CO$waterbodyid3dhp), type = "waterbody")

plot(sf::st_geometry(CO_wb[grepl("Powell", CO_wb$gnisidlabel),]),
     col = "blue", border = "NA")

^{Created on 2023-10-19 with reprex v2.0.2}

[dblodgett-usgs]

Some additional reprex

library(dataRetrieval)
library(nhdplusTools)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

pc <- function(x) st_geometry(st_transform(x, 3857))

basin <- findNLDI(location = c(-89.441, 43.487), 
                  find = "basin")

plot_box <- st_bbox(st_buffer(
  st_transform(basin$basin, 5070), 
  dist = units::as_units(-11, "km")))

network <- get_3dhp(st_as_sfc(plot_box), 
                    type = "flowline")
water <- get_3dhp(st_as_sfc(plot_box), 
                  type = "waterbody")

par(mar = c(0, 0, 0, 0))
plot_nhdplus(bbox = plot_box, 
             plot_config = list(flowline = list(col = NULL)), 
             zoom = 13)
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> Zoom set to: 13
plot(pc(network), col = "darkblue", 
     lwd = 1.5, add = TRUE)             
plot(pc(water), col = "lightblue", 
     border = "blue", add = TRUE)

^{Created on 2024-03-24 with reprex v2.0.2}

library(nhdplusTools)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

pc <- function(x) st_geometry(st_transform(x, 3857))

outlet <- st_sfc(st_point(c(-89.441, 43.487)), crs = 4326)

flowline <- get_3dhp(outlet, type = "flowline", buffer = 10)

unique(flowline$mainstemid)
#> [1] "https://geoconnex.us/ref/mainstems/359842"

mainstem <- read_sf(unique(flowline$mainstemid))

downmain <- read_sf(mainstem$downstream_mainstem_id)

dplyr::glimpse(st_drop_geometry(mainstem), width = 55)
#> Rows: 1
#> Columns: 22
#> $ id                           <chr> "359842"
#> $ name_at_outlet_gnis_id       <int> 1561153
#> $ outlet_rf1id                 <int> 18616
#> $ head_nhdpv2_comid            <chr> "https://geoconn…
#> $ head_nhdpv1_comid            <int> 13652654
#> $ fid                          <int> 11347
#> $ outlet_nhdpv2_comid          <chr> "https://geoconn…
#> $ outlet_nhdpv1_comid          <int> 13654326
#> $ uri                          <chr> "https://geoconn…
#> $ head_nhdpv2huc12             <chr> "https://geoconn…
#> $ head_2020huc12               <chr> "070700040102"
#> $ featuretype                  <chr> "['https://www.o…
#> $ outlet_nhdpv2huc12           <chr> "https://geoconn…
#> $ outlet_2020huc12             <chr> "070700040406"
#> $ downstream_mainstem_id       <chr> "https://geoconn…
#> $ lengthkm                     <dbl> 178.5
#> $ superseded                   <lgl> FALSE
#> $ encompassing_mainstem_basins <chr> "['https://geoco…
#> $ outlet_drainagearea_sqkm     <dbl> 1676.1
#> $ new_mainstemid               <chr> ""
#> $ name_at_outlet               <chr> "Baraboo River"
#> $ head_rf1id                   <int> 18627

par(mar = c(0, 0, 0, 0))
plot_nhdplus(bbox = st_bbox(mainstem), 
             plot_config = list(flowline = list(col = NULL)))
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> Zoom set to: 9
plot(pc(flowline), col = "darkblue", lwd = 4, add = TRUE)             
plot(pc(mainstem), col = "blue", add = TRUE)
plot(pc(downmain), col = "blue", lwd = 2, add = TRUE)

^{Created on 2024-03-24 with reprex v2.0.2}

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(hydroloom)
library(nhdplusTools)
#> 
#> Attaching package: 'nhdplusTools'
#> The following object is masked from 'package:hydroloom':
#> 
#>     make_node_topology
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

pc <- function(x) st_geometry(st_transform(x, 3857))

outlet <- st_sfc(st_point(c(-90.387826, 43.602920)), crs = 4326)

flowline <- get_3dhp(outlet, type = "flowline", buffer = 10)

unique(flowline$mainstemid)
#> [1] "https://geoconnex.us/usgs/mainstems/4994552"

mainstem <- read_sf(unique(flowline$mainstemid))

invisible(hydroloom_names(c(id3dhp = "id")))

network <- make_attribute_topology(st_transform(mainstem, 5070),
                                   min_distance = 10)
#> Loading required namespace: future
#> Loading required namespace: future.apply

outlet <- filter(mainstem, 
                 id3dhp == network$id3dhp[network$toid == ""])

next_down <- get_3dhp(outlet, type = "flowline", buffer = 10)
#> st_as_s2(): dropping Z and/or M coordinate

downmain <- filter(next_down, 
                   !mainstemid %in% mainstem$mainstemid)

unique(downmain$mainstemid)
#> [1] "https://geoconnex.us/usgs/mainstems/10208138"

downmain <- read_sf(unique(downmain$mainstemid))

par(mar = c(0, 0, 0, 0))
plot_nhdplus(bbox = st_bbox(c(st_geometry(mainstem), 
                              st_geometry(downmain))), 
             plot_config = list(flowline = list(col = NULL)), 
             zoom = 13)
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> Zoom set to: 13
plot(pc(flowline), col = "darkblue", lwd = 4, add = TRUE)             
plot(pc(mainstem), col = "blue", add = TRUE)
plot(pc(downmain), col = "blue", lwd = 2, add = TRUE)

^{Created on 2024-03-24 with reprex v2.0.2}