# 3DHP Data Access

#### dblodgett@usgs.gov

This vignette shows how to work with a new service based on data from the 3D Hydrography Program (3DHP) released in 2024. 3DHP data uses “mainstem” identifiers as the primary river ID. We’ll work with one for the Baraboo River in Wisconsin.

We can find one to start with at a url like:

`https://reference.geoconnex.us/collections/mainstems/items?filter=name_at_outlet ILIKE '%Baraboo%'`

This url searches the reference mainstem collection for rivers with names like “Baraboo”. Using that url query, we can find that the Baraboo mainstem ID is `https://geoconnex.us/ref/mainstems/359842` and that it flows to the Wisconsin and Mississippi `https://geoconnex.us/ref/mainstems/323742` and `https://geoconnex.us/ref/mainstems/312091` respectively. At the reference mainstem page, we can also find that the outlet NHDPlusV2 COMID is `https://geoconnex.us/nhdplusv2/comid/937070225`.

The code block just below generates starting locations based on information derived from the reference mainstem for the Baraboo and a point location near the Baraboo’s outlet. It queries the 3DHP service for a flowline within a 10m buffer of the point specified just to get us started.

      comid <- "937070225"

      point <- c(-89.441, 43.487) |>
        sf::st_point() |>
        sf::st_sfc(crs = 4326) |> sf::st_sf()

      dm <- c('https://geoconnex.us/ref/mainstems/323742',
              'https://geoconnex.us/ref/mainstems/312091')

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

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

      plot_nhdplus(bbox = sf::st_bbox(flowline), plot_config = list(flowline = list(col = NULL)), zoom = 14)
      plot(pc(flowline), add = TRUE, col = "darkblue", lwd = 2)
      plot(pc(point), pch = "O", add = TRUE)

Now, we can use dataRetrieval to retrieve a basin upstream of this outlet location. We specify the NHDPlusV2 comid here since we know it, but could have used the point location just as well.

With the basin pulled from the `findNLDI()` function, we can pass it to `get_3dhp()` as a Area of Interest for flowlines, waterbodies, and hydrolocations. We can then use the mainstem ids of our two downstream mainstem rivers found above to get the downtream mainstem rivers.

      basin <- dataRetrieval::findNLDI(comid = comid, find = "basin")

      network <- get_3dhp(basin$basin, type = "flowline")
      water <- get_3dhp(basin$basin, type = "waterbody")
      hydrolocation <- get_3dhp(basin$basin, type = "hydrolocation - reach code, external connection")

      down_mains <- get_3dhp(ids = dm, type = "flowline")

      old_par <- par(mar = c(0, 0, 0, 0))
      plot_nhdplus(bbox = sf::st_bbox(basin$basin), flowline_only = TRUE,
                   plot_config = list(flowline = list(col = NULL)), zoom = 10)
      plot(pc(basin$basin), lwd = 2, add = TRUE)
      plot(pc(network), lwd = 0.5, add = TRUE)
      plot(pc(water), lwd = 0.5, border = "skyblue", col = "lightblue", add = TRUE)
      plot(pc(hydrolocation), pch = "o", col = "#80808026", add = TRUE)
      plot(pc(down_mains), lwd = 3, col = "blue", add = TRUE)

      par(old_par)

      old_par <- par(mar = c(0, 0, 0, 0))
      plot(pc(down_mains))
      plot(pc(basin$basin), add = TRUE)

      par(old_par)

Neat. Now we have flowlines, waterbodies, hydrologic locations, a basin boundary, and major rivers downstream to work with. These could be saved to a local file for use in a GIS or used in some further analysis.

Say you want to know how to relate NHD data to 3DHP data. 3DHP hydrolocations include “reachcode” tops and bottoms to hel with that. Here, we lookup the hydrolocation with reachcode `07070004002889`then grab the rest of the reachcode hydrolocations along the mainstem that the reachcode is part of.

      reachcode <- "07070004002889"

      hydrolocation <- get_3dhp(universalreferenceid = reachcode,
                                type = "hydrolocation - reach code, external connection")
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/40/query
    #> universalreferenceid IN ('07070004002889')jsontrue
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/40/query
    #> 5800707,25271902*geojson

      hydrolocation
    #> Simple feature collection with 2 features and 9 fields
    #> Geometry type: POINT
    #> Dimension:     XY
    #> Bounding box:  xmin: -89.44179 ymin: 43.48568 xmax: -89.43904 ymax: 43.48829
    #> Geodetic CRS:  WGS 84
    #>    id3dhp   featuredate                                mainstemid
    #> 1 01HAJW7 1694649600000 https://geoconnex.us/ref/mainstems/359842
    #> 2 01SVVYQ 1694649600000 https://geoconnex.us/ref/mainstems/359842
    #>   universalreferenceid gnisid gnisidlabel featuretype featuretypelabel
    #> 1       07070004002889   <NA>        <NA>          10  Reachcode Start
    #> 2       07070004002889   <NA>        <NA>          11    Reachcode End
    #>   workunitid                   geometry
    #> 1       <NA> POINT (-89.44179 43.48829)
    #> 2       <NA> POINT (-89.43904 43.48568)

      mainstem_points <- get_3dhp(ids = hydrolocation$mainstemid, type = "hydrolocation - reach code, external connection")
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/40/query
    #> mainstemid IN ('https://geoconnex.us/ref/mainstems/359842', 'https://geoconnex.us/ref/mainstems/359842')jsontrue
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/40/query
    #> 5799214,5799215,5799216,5799218,5799219,5799220,5799221,5799222,5799223,5799224,5799225,5799226,5799227,5799228,5799229,5799230,5799231,5799232,5799233,5799234,5799235,5799236,5799237,5799238,5799239,5799240,5799241,5799242,5799243,5799244,5799245,5799246,5799247,5799248,5799249,5799250,5799251,5799258,5799259,5799260,5799261,5799262,5799270,5799271,5799272,5799273,5799274,5799275,5799276,5799277,5799278,5799285,5799286,5799287,5799288,5799289,5799290,5799291,5799292,5799293,5799294,5799295,5799296,5799297,5799298,5799299,5799300,5799301,5799302,5799303,5799304,5799305,5799306,5799307,5799308,5799309,5799869,5799879,5799881,5800191,5800312,5800433,5800442,5800455,5800478,5800484,5800485,5800629,5800633,5800670,5800707,25270409,25270410,25270411,25270413,25270414,25270415,25270416,25270417,25270418,25270419,25270420,25270421,25270422,25270423,25270424,25270425,25270426,25270427,25270428,25270429,25270430,25270431,25270432,25270433,25270434,25270435,25270436,25270437,25270438,25270439,25270440,25270441,25270442,25270443,25270444,25270445,25270446,25270453,25270454,25270455,25270456,25270457,25270465,25270466,25270467,25270468,25270469,25270470,25270471,25270472,25270473,25270480,25270481,25270482,25270483,25270484,25270485,25270486,25270487,25270488,25270489,25270490,25270491,25270492,25270493,25270494,25270495,25270496,25270497,25270498,25270499,25270500,25270501,25270502,25270503,25270504,25271064,25271074,25271076,25271386,25271507,25271628,25271637,25271650,25271673,25271679,25271680,25271824,25271828,25271865,25271902*geojson
      mainstem_lines <- get_3dhp(ids = hydrolocation$mainstemid, type = "flowline")
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/50/query
    #> mainstemid IN ('https://geoconnex.us/ref/mainstems/359842', 'https://geoconnex.us/ref/mainstems/359842')jsontrue
    #> https://hydro.nationalmap.gov/arcgis/rest/services/3DHP_all/MapServer/50/query
    #> 166149,166153,166156,166167,166226,620198,620242,1549503,1549508,1775834,1775872,2002787,2002816,2288867,2288870,2516418,2516420,2516496,2743563,2743647,2743663,2970336,2970404,3197491,3197492,3197514,3197519,3424074,3424099,3651025,3877526,3877599,4104386,4331439,4558446,4558512,4785327,5012065,5012092,5012169,5012174,5012181,5239211,5466418,5466477,5693565,5920117,5920130,5920163,5920221,6374372,6827725,7054423,7054506,7282175,7282198,7282203,7508821,7508844,7508912,7736344,7736366,7884524,7963932,8418465,8647012,8875153,9329690,9329712,9329751,9329758,9329763,9556807,9556840,9784984,9785004,9785058,10013179,10013227,10013285,10013298,10241420,10241494,10468536,10696621,10924194,10924201,10924210,10924307,11152573,11379880,11607428,11607433,11607492,11835528,12063466,12063471,12291441,12291464,12291473,12518287,12518337,12518383,13202390,13430687,13659066,13659070,13659072,13659101,13659163,14115499,14115561,14342674,14342678,14342742,14571352,15026248,15026278,15254035,15481887,15481944,16165429,16165447,16165475,16165491,16165521,16393674,16393676,16621557,16849105,17304165,17304192,17304194,17304201,17304249,17532332,17759883,17987501,18214503,18441913,18670640,18670648,18898970,18898989,19353340,19580366,19580424,19808311,20036860,20036918,20263934,20263936,20492572,20721053,20721069,20948566,21176001,21630720,21630726,21859053,22086419,22314085,22314094,22314130,22542072,22999441,22999515,23226961,23226989,23455619,23684170,23684174,23684225,24240147,24241311,24368049,24595119,24595159,25279723,25508187,25508266,25508280,25734303,25734325,25961872,26189223,26189317,26189337,26416378,26645046,26645051,26873168,26873214,27103181,27103182,27331602,27331618,27331623,27331629,29329570,29556652,29778146,29781477,29784082,30010627*geojson

      old_par <- par(mar = c(0, 0, 0, 0))
      plot(pc(hydrolocation), pch = "o", col = "#808080BF")
      plot(pc(mainstem_lines), lwd = 2, col = "blue", add = TRUE)

      par(old_par)

      old_par <- par(mar = c(0, 0, 0, 0))
      plot_nhdplus(bbox = sf::st_bbox(basin$basin), flowline_only = TRUE,
                   plot_config = list(flowline = list(col = NULL)), zoom = 10)
    #> Spherical geometry (s2) switched off
    #> https://api.water.usgs.gov/geoserver/wmadata/ows
    #> <?xml version="1.0" encoding="UTF-8"?>
    #> <wfs:GetFeature xmlns:wfs="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:gml="http://www.opengis.net/gml" service="WFS" version="1.1.0" outputFormat="application/json" xsi:schemaLocation="http://www.opengis.net/wfs                         http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    #>   <wfs:Query xmlns:feature="https://api.water.usgs.gov/wmadata" typeName="feature:nhdflowline_network" srsName="EPSG:4269">
    #>     <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
    #>       <ogc:And>
    #>         <ogc:BBOX>
    #>           <ogc:PropertyName>the_geom</ogc:PropertyName>
    #>           <gml:Envelope srsName="urn:x-ogc:def:crs:EPSG:4326">
    #>             <gml:lowerCorner>43.378165219 -90.461663253</gml:lowerCorner>
    #>             <gml:upperCorner>43.84757169 -89.434028091</gml:upperCorner>
    #>           </gml:Envelope>
    #>         </ogc:BBOX>
    #>       </ogc:And>
    #>     </ogc:Filter>
    #>   </wfs:Query>
    #> </wfs:GetFeature>
    #> although coordinates are longitude/latitude, st_intersects assumes that they
    #> are planar
    #> Spherical geometry (s2) switched on
    #> Zoom set to: 10
      plot(pc(mainstem_lines), lwd = 2, col = "blue", add = TRUE)
      plot(pc(mainstem_points), pch = "o", col = "#808080BF", add = TRUE)

      par(old_par)

With this functionality, we have what we need to inter operate between older NHD data and newly released 3DHP data. This service is brand new and `nhdplusTools` is just getting going with functionality to support what’s shown here. Please submit issues if you find problems or have questions!!