Saves a subset of the National Seamless database or other nhdplusTools compatible data based on a specified collection of COMIDs. This function uses get_nhdplus for the "download" data source but returns data consistent with local data subsets in a subset file.

subset_nhdplus(
  comids = NULL,
  output_file = NULL,
  nhdplus_data = NULL,
  bbox = NULL,
  simplified = TRUE,
  overwrite = FALSE,
  return_data = TRUE,
  status = TRUE,
  flowline_only = NULL,
  streamorder = NULL,
  out_prj = 4269
)

Arguments

comids

integer vector of COMIDs to include.

output_file

character path to save the output to defaults to the directory of the nhdplus_data.

nhdplus_data

character path to the .gpkg or .gdb containing the national seamless database, a subset of NHDPlusHR, or "download" to use a web service to download NHDPlusV2.1 data. Not required if nhdplus_path has been set or the default has been adopted. See details for more.

bbox

object of class "bbox" as returned by sf::st_bbox in Latitude/Longitude. If no CRS is present, will be assumed to be in WGS84 Latitude Longitude.

simplified

boolean if TRUE (the default) the CatchmentSP layer will be included. Not relevant to the "download" option or NHDPlusHR data.

overwrite

boolean should the output file be overwritten

return_data

boolean if FALSE path to output file is returned silently otherwise data is returned in a list.

status

boolean should the function print status messages

flowline_only

boolean WARNING: experimental if TRUE only the flowline network and attributes will be returned

streamorder

integer only streams of order greater than or equal will be downloaded. Not implemented for local data.

out_prj

character override the default output CRS of NAD83 lat/lon (EPSG:4269)

Value

character path to the saved subset geopackage

Details

This function relies on the National Seamless Geodatabase or Geopackage. It can be downloaded here.

The "download" option of this function should be considered preliminary and subject to revision. It does not include as many layers and may not be available permenently.

Examples

# \donttest{

source(system.file("extdata/sample_data.R", package = "nhdplusTools"))

nhdplus_path(sample_data)

staged_nhdplus <- stage_national_data(output_path = tempdir())
#> Warning: attributes file exists
#> Warning: flowline file exists
#> Warning: catchment already exists.

sample_flines <- readRDS(staged_nhdplus$flowline)

geom_col <- attr(sample_flines, "sf_column")

plot(sample_flines[[geom_col]],
     lwd = 3)

start_point <- sf::st_sfc(sf::st_point(c(-89.362239, 43.090266)),
                          crs = 4326)

plot(start_point, cex = 1.5, lwd = 2, col = "red", add = TRUE)

start_comid <- discover_nhdplus_id(start_point)

comids <- get_UT(sample_flines, start_comid)

plot(dplyr::filter(sample_flines, COMID %in% comids)[[geom_col]],
     add=TRUE, col = "red", lwd = 2)

output_file <- tempfile(fileext = ".gpkg")

subset_nhdplus(comids = comids,
               output_file = output_file,
               nhdplus_data = sample_data,
               overwrite = TRUE,
               status = TRUE)
#> All intersections performed in latitude/longitude.
#> Reading NHDFlowline_Network
#> 168 comids of 168
#> Writing NHDFlowline_Network
#> Reading CatchmentSP
#> 168 comids of 168
#> Writing CatchmentSP
#> Reading NHDArea
#> Writing NHDArea
#> Reading NHDWaterbody
#> Writing NHDWaterbody
#> Reading NHDFlowline_NonNetwork
#> Writing NHDFlowline_NonNetwork
#> Reading Gage
#> Writing Gage
#> Reading Sink
#> No features to write in Sink
#> $NHDFlowline_Network
#> Simple feature collection with 168 features and 136 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -89.58537 ymin: 43.08521 xmax: -89.21254 ymax: 43.30179
#> Geodetic CRS:  NAD83
#> # A tibble: 168 × 137
#>       COMID FDATE               RESOLU…¹ GNIS_ID GNIS_…² LENGT…³ REACH…⁴ FLOWDIR
#>  *    <int> <dttm>              <chr>    <chr>   <chr>     <dbl> <chr>   <chr>  
#>  1 13293750 1999-10-29 00:00:00 Medium   1577073 Yahara…   1.72  070900… With D…
#>  2 13293504 1999-10-29 00:00:00 Medium   1577073 Yahara…   1.41  070900… With D…
#>  3 13294134 1999-10-29 00:00:00 Medium   1577073 Yahara…   0.74  070900… With D…
#>  4 13294128 1999-10-29 00:00:00 Medium   1577073 Yahara…   3.70  070900… With D…
#>  5 13294394 1999-10-29 00:00:00 Medium   1577073 Yahara…   0.077 070900… With D…
#>  6 13293454 1999-10-29 00:00:00 Medium   1577073 Yahara…   0.94  070900… With D…
#>  7 13293430 1999-10-29 00:00:00 Medium   1577073 Yahara…   1.14  070900… With D…
#>  8 13293424 1999-10-29 00:00:00 Medium   1577073 Yahara…   1.27  070900… With D…
#>  9 13294110 1999-10-29 00:00:00 Medium   1577073 Yahara…   7.48  070900… With D…
#> 10 13293398 1999-10-29 00:00:00 Medium   1577073 Yahara…   0.082 070900… With D…
#> # … with 158 more rows, 129 more variables: WBAREACOMI <int>, FTYPE <chr>,
#> #   FCODE <int>, Shape_Length <dbl>, StreamLeve <int>, StreamOrde <int>,
#> #   StreamCalc <int>, FromNode <dbl>, ToNode <dbl>, Hydroseq <dbl>,
#> #   LevelPathI <dbl>, Pathlength <dbl>, TerminalPa <dbl>, ArbolateSu <dbl>,
#> #   Divergence <int>, StartFlag <int>, TerminalFl <int>, DnLevel <int>,
#> #   UpLevelPat <dbl>, UpHydroseq <dbl>, DnLevelPat <dbl>, DnMinorHyd <dbl>,
#> #   DnDrainCou <int>, DnHydroseq <dbl>, FromMeas <dbl>, ToMeas <dbl>, …
#> 
#> $CatchmentSP
#> Simple feature collection with 167 features and 6 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -89.60479 ymin: 43.03507 xmax: -89.20378 ymax: 43.36607
#> Geodetic CRS:  NAD83
#> # A tibble: 167 × 7
#>    GRIDCODE FEATUREID SOURCEFC AreaS…¹ Shape…² Shape…³                      geom
#>  *    <int>     <int> <chr>      <dbl>   <dbl>   <dbl>        <MULTIPOLYGON [°]>
#>  1  1085160  13293454 NHDFlow…  0.328  0.0256  3.64e-5 (((-89.34777 43.20671, -…
#>  2  1085219  13293750 NHDFlow…  6.69   0.148   7.40e-4 (((-89.38126 43.07518, -…
#>  3  1085247  13294134 NHDFlow…  0.261  0.0307  2.88e-5 (((-89.35547 43.18545, -…
#>  4  1085414  13293570 NHDFlow…  0.0514 0.00963 5.69e-6 (((-89.41468 43.15295, -…
#>  5  1085447  13293430 NHDFlow…  1.63   0.0780  1.81e-4 (((-89.35208 43.20989, -…
#>  6  1085448  13293526 NHDFlow…  0.0690 0.0117  7.64e-6 (((-89.34975 43.16686, -…
#>  7  1085463  13293588 NHDFlow…  0.429  0.0404  4.75e-5 (((-89.42608 43.14793, -…
#>  8  1085478  13293614 NHDFlow…  3.15   0.0953  3.48e-4 (((-89.58446 43.13534, -…
#>  9  1085502  13294264 NHDFlow… 24.8    0.348   2.75e-3 (((-89.53449 43.05226, -…
#> 10  1085509  13293384 NHDFlow…  0.291  0.0296  3.23e-5 (((-89.33033 43.27117, -…
#> # … with 157 more rows, and abbreviated variable names ¹​AreaSqKM,
#> #   ²​Shape_Length, ³​Shape_Area
#> 
#> $NHDArea
#> [1] "NHDArea"
#> 
#> $NHDWaterbody
#> [1] "NHDWaterbody"
#> 
#> $NHDFlowline_NonNetwork
#> [1] "NHDFlowline_NonNetwork"
#> 
#> $Gage
#> [1] "Gage"
#> 
#> $Sink
#> NULL
#> 

sf::st_layers(output_file)
#> Driver: GPKG 
#> Available layers:
#>               layer_name geometry_type features fields crs_name
#> 1    NHDFlowline_Network   Line String      168    136    NAD83
#> 2            CatchmentSP Multi Polygon      167      6    NAD83
#> 3                NHDArea       Polygon        1     14    NAD83
#> 4           NHDWaterbody       Polygon       90     21    NAD83
#> 5 NHDFlowline_NonNetwork   Line String       45     12    NAD83
#> 6                   Gage         Point       33     19    NAD83

catchment <- sf::read_sf(output_file, "CatchmentSP")

plot(catchment[[attr(catchment, "sf_column")]], add = TRUE)

waterbody <- sf::read_sf(output_file, "NHDWaterbody")

plot(waterbody[[attr(waterbody, "sf_column")]],
     col = rgb(0, 0, 1, alpha = 0.5), add = TRUE)


# Cleanup temp
sapply(staged_nhdplus, unlink)
#> attributes   flowline  catchment 
#>          0          0          0 
unlink(output_file)

# Download Option:
subset_nhdplus(comids = comids,
               output_file = output_file,
               nhdplus_data = "download",
               overwrite = TRUE,
               status = TRUE, flowline_only = FALSE)
#> All intersections performed in latitude/longitude.
#> Reading NHDFlowline_Network
#> Spherical geometry (s2) switched off
#> Spherical geometry (s2) switched on
#> Writing NHDFlowline_Network
#> Reading CatchmentSP
#> Spherical geometry (s2) switched off
#> Spherical geometry (s2) switched on
#> Writing CatchmentSP
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> Spherical geometry (s2) switched off
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Spherical geometry (s2) switched on
#> $NHDFlowline_Network
#> Simple feature collection with 168 features and 138 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -89.58537 ymin: 43.08522 xmax: -89.21254 ymax: 43.30179
#> Geodetic CRS:  NAD83
#> # A tibble: 168 × 139
#>    id          comid fdate               resol…¹ gnis_id gnis_…² lengt…³ reach…⁴
#>  * <chr>       <int> <dttm>              <chr>   <chr>   <chr>     <dbl> <chr>  
#>  1 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       3.04  070900…
#>  2 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       1.44  070900…
#>  3 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  "15770… "Yahar…   4.84  070900…
#>  4 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       4.40  070900…
#>  5 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  "15770… "Yahar…   0.839 070900…
#>  6 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       4.67  070900…
#>  7 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  "15770… "Yahar…   0.402 070900…
#>  8 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       1.08  070900…
#>  9 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  "15770… "Yahar…   1.12  070900…
#> 10 nhdflowli… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       0.351 070900…
#> # … with 158 more rows, 131 more variables: flowdir <chr>, wbareacomi <int>,
#> #   ftype <chr>, fcode <int>, shape_length <dbl>, streamleve <int>,
#> #   streamorde <int>, streamcalc <int>, fromnode <dbl>, tonode <dbl>,
#> #   hydroseq <dbl>, levelpathi <dbl>, pathlength <dbl>, terminalpa <dbl>,
#> #   arbolatesu <dbl>, divergence <int>, startflag <int>, terminalfl <int>,
#> #   dnlevel <int>, uplevelpat <dbl>, uphydroseq <dbl>, dnlevelpat <dbl>,
#> #   dnminorhyd <dbl>, dndraincou <int>, dnhydroseq <dbl>, frommeas <dbl>, …
#> 
#> $CatchmentSP
#> Simple feature collection with 167 features and 7 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -89.60479 ymin: 43.03507 xmax: -89.20378 ymax: 43.36607
#> Geodetic CRS:  NAD83
#> # A tibble: 167 × 8
#>    id                 gridcode featureid sourcefc    areasqkm shape_le…¹ shape…²
#>  * <chr>                 <int>     <int> <chr>          <dbl>      <dbl>   <dbl>
#>  1 catchmentsp.948294  1085512  13293376 NHDFlowline  47.8        0.465  5.30e-3
#>  2 catchmentsp.948495  1085714  13293378 NHDFlowline   3.30       0.0899 3.66e-4
#>  3 catchmentsp.948432  1085651  13293380 NHDFlowline  19.1        0.247  2.12e-3
#>  4 catchmentsp.948531  1085750  13293382 NHDFlowline  12.1        0.197  1.35e-3
#>  5 catchmentsp.948291  1085509  13293384 NHDFlowline   0.291      0.0296 3.23e-5
#>  6 catchmentsp.948759  1085979  13293386 NHDFlowline   9.94       0.175  1.10e-3
#>  7 catchmentsp.948299  1085517  13293388 NHDFlowline   0.154      0.0187 1.71e-5
#>  8 catchmentsp.948549  1085768  13293390 NHDFlowline   1.41       0.0612 1.57e-4
#>  9 catchmentsp.948386  1085605  13293392 NHDFlowline   0.461      0.0393 5.11e-5
#> 10 catchmentsp.948674  1085893  13293394 NHDFlowline   0.0796     0.0124 8.82e-6
#> # … with 157 more rows, 1 more variable: geometry <MULTIPOLYGON [°]>, and
#> #   abbreviated variable names ¹​shape_length, ²​shape_area
#> 
#> $NHDArea
#> Simple feature collection with 1 feature and 15 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -89.40194 ymin: 43.15016 xmax: -89.36131 ymax: 43.18112
#> Geodetic CRS:  NAD83
#> # A tibble: 1 × 16
#>   id     comid fdate               resol…¹ gnis_id gnis_…² areas…³ eleva…⁴ ftype
#> * <chr>  <int> <dttm>              <chr>   <chr>   <chr>     <dbl>   <int> <chr>
#> 1 nhda… 1.47e7 1999-10-28 23:00:00 Medium  " "     " "        1.36       0 Stre…
#> # … with 7 more variables: fcode <int>, shape_length <dbl>, shape_area <dbl>,
#> #   onoffnet <int>, purpcode <chr>, purpdesc <chr>, geometry <POLYGON [°]>, and
#> #   abbreviated variable names ¹​resolution, ²​gnis_name, ³​areasqkm, ⁴​elevation
#> 
#> $NHDWaterbody
#> Simple feature collection with 90 features and 22 fields
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -89.72879 ymin: 42.9897 xmax: -89.20939 ymax: 43.40395
#> Geodetic CRS:  NAD83
#> # A tibble: 90 × 23
#>    id          comid fdate               resol…¹ gnis_id gnis_…² areas…³ eleva…⁴
#>  * <chr>       <int> <dttm>              <chr>   <chr>   <chr>     <dbl>   <int>
#>  1 nhdwaterb… 1.47e7 1999-10-28 23:00:00 Medium  " "     " "       1.72        0
#>  2 nhdwaterb… 1.36e7 1999-10-28 23:00:00 Medium  " "     " "       0.015       0
#>  3 nhdwaterb… 1.36e7 1999-10-28 23:00:00 Medium  "15769… "Lake …  29.0       236
#>  4 nhdwaterb… 1.36e7 1999-10-28 23:00:00 Medium  " "     " "       0.018       0
#>  5 nhdwaterb… 1.36e7 1999-10-28 23:00:00 Medium  " "     " "       0.016       0
#>  6 nhdwaterb… 1.36e7 1999-10-28 23:00:00 Medium  " "     " "       0.014       0
#>  7 nhdwaterb… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       0.062       0
#>  8 nhdwaterb… 1.33e7 1999-10-28 23:00:00 Medium  "15737… "Schoe…   0.473       0
#>  9 nhdwaterb… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       0.045       0
#> 10 nhdwaterb… 1.33e7 1999-10-28 23:00:00 Medium  " "     " "       0.123       0
#> # … with 80 more rows, 15 more variables: reachcode <chr>, ftype <chr>,
#> #   fcode <int>, shape_length <dbl>, shape_area <dbl>, onoffnet <int>,
#> #   purpcode <chr>, purpdesc <chr>, meandepth <dbl>, lakevolume <dbl>,
#> #   maxdepth <dbl>, meandused <dbl>, meandcode <chr>, lakearea <dbl>,
#> #   geometry <POLYGON [°]>, and abbreviated variable names ¹​resolution,
#> #   ²​gnis_name, ³​areasqkm, ⁴​elevation
#> 
#> $NHDFlowline_NonNetwork
#> Simple feature collection with 45 features and 13 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -89.62331 ymin: 43.04214 xmax: -89.25271 ymax: 43.35874
#> Geodetic CRS:  NAD83
#> # A tibble: 45 × 14
#>    id          comid fdate               resol…¹ gnis_id gnis_…² lengt…³ reach…⁴
#>  * <chr>       <int> <dttm>              <chr>   <chr>   <chr>     <dbl> <chr>  
#>  1 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        1.90 070700…
#>  2 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        1.80 070700…
#>  3 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        1.39 070700…
#>  4 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        2.36 070700…
#>  5 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        1.56 070700…
#>  6 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        2.87 070700…
#>  7 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        2.42 070700…
#>  8 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        2.12 070700…
#>  9 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        2.19 070700…
#> 10 nhdflowli… 1.36e7 1999-10-28 19:00:00 Medium  " "     " "        1.29 070700…
#> # … with 35 more rows, 6 more variables: flowdir <chr>, wbareacomi <int>,
#> #   ftype <chr>, fcode <int>, shape_length <dbl>, geometry <LINESTRING [°]>,
#> #   and abbreviated variable names ¹​resolution, ²​gnis_name, ³​lengthkm,
#> #   ⁴​reachcode
#> 

sf::st_layers(output_file)
#> Driver: GPKG 
#> Available layers:
#>               layer_name geometry_type features fields crs_name
#> 1    NHDFlowline_Network   Line String      168    138    NAD83
#> 2            CatchmentSP Multi Polygon      167      7    NAD83
#> 3                NHDArea       Polygon        1     15    NAD83
#> 4           NHDWaterbody       Polygon       90     22    NAD83
#> 5 NHDFlowline_NonNetwork   Line String       45     13    NAD83

# NHDPlusHR
source(system.file("extdata/nhdplushr_data.R", package = "nhdplusTools"))

up_ids <- get_UT(hr_data$NHDFlowline, 15000500028335)

sub_gpkg <- file.path(work_dir, "sub.gpkg")
sub_nhdhr <- subset_nhdplus(up_ids, output_file = sub_gpkg,
                            nhdplus_data = hr_gpkg, overwrite = TRUE)
#> All intersections performed in latitude/longitude.
#> Reading NHDFlowline
#> 1000 comids of 1427
#> 1427 comids of 1427
#> Writing NHDFlowline
#> Reading NHDPlusCatchment
#> 1000 comids of 1427
#> 1427 comids of 1427
#> Found invalid geometry, attempting to fix.
#> Writing NHDPlusCatchment
#> Reading NHDArea
#> Writing NHDArea
#> Reading NHDWaterbody
#> Writing NHDWaterbody
#> Reading NHDPlusSink
#> No features to write in NHDPlusSink

sf::st_layers(sub_gpkg)
#> Driver: GPKG 
#> Available layers:
#>         layer_name geometry_type features fields crs_name
#> 1      NHDFlowline   Line String     1427     57    NAD83
#> 2 NHDPlusCatchment Multi Polygon     1361      7    NAD83
#> 3          NHDArea       Polygon        8     14    NAD83
#> 4     NHDWaterbody       Polygon      740     15    NAD83
names(sub_nhdhr)
#> [1] "NHDFlowline"      "NHDPlusCatchment" "NHDArea"          "NHDWaterbody"    
#> [5] "NHDPlusSink"     

plot(sf::st_geometry(hr_data$NHDFlowline), lwd = 0.5)
plot(sf::st_geometry(sub_nhdhr$NHDFlowline), lwd = 0.6, col = "red", add = TRUE)


unlink(output_file)
unlink(sub_gpkg)

# }