I have had success accessing unique database IDs and polygons for various lakes in the US. I'm now looking to get polygons for the watersheds associated with these lakes, but have stumbled a bit.

Example- I get a lake polygon from a GPS point using the Hydrolinks package:


Lon<- -97.0508

Lake <-link_to_waterbodies(Lat, Lon, "Name", dataset = c("nhdh", "hydrolakes", "nhdplusv2"), buffer = 50)

The result is Santa Fe Lake in Kansas. (ID = 45827774)

Using the "nhdR" package, I can get the polygon associated with this lake

dt <- nhd_load(c("KS"), c("NHDWaterbody"))
Poly<-dt[dt$Permanent_Identifier == 45827774]

I'm using results from the streamstats API as my "true" watershed (see screen shot from streamstats API for this lake).

The reachcode for the "Poly" object is "11030017001699". And if I try to calculating watershed data from the first 10 digits of the reachcode, the basin doesn't match what I get from the streamstatas API.

enter image description here

I tried using the pour point for this lake (Hydrolakes lists this as -97.05319 37.709), but calculating a watershed based on these points only gets me a small section of the true watershed.


#Pour Points
PLon<- -97.05319
PLat<- 37.709
ws1 <- delineateWatershed(xlocation = PLon, ylocation = PLat, crs = 4326,
                      includeparameters = "true")

So I'm assuming that the accuracy of these pour points is not good enough to use streamstats to calculate entire lake watersheds.

Using tools from nhdplusTools package works for this lake


 #Step 1
  point <- sf::st_sfc(sf::st_point(c(PLon, PLat)), crs = 4326)
  #Get nhdplus_id from point

  #Get Basin 
  nldi <- list(featureSource = "comid", featureID = ID)
  basin <- get_nldi_basin(nldi_feature = nldi)

  leaflet() %>%
    addPolygons(data = basin) %>%
    addCircleMarkers(data = point,
                     #lat = ~lat, lng = ~long,
                     color = "red",
                     radius = 2) %>%
    addTiles() #Add tiles adds defualt map 

Which seems to work well. But I don' always have pour points for lakes, or some of my pour points do not give me a full watershed. Using random GPS points from within a lake polygon gives me unreliable results.

 Lat1<- 37.72797
 Lon1<- -97.82922
 ID = 135039628

Plugging the above coordinates into the "get_nldi_basin" function only gets a small aspect of the entire watershed

And if I used these coords for a small pond:

Lat2<- 39.27623
Lon2<- -94.89211
ID = 159110681

I get a huge watershed

So using the "get_nldi_basin" seems to be limited to my ability to get an accurate pour point. I'd really prefer a method that gives me an accurate watershed boundary for an entire lake using the entire lake polygon, or a match to the permanent NHD ID.

  • 1
    I'm looking into it but hidrolinks is not available for R 3.6. Which one are you using ? I'm not sure about your final goal. Do you want to get a layer of watershed boundary (with ID) given a pair of coords ? – César Arquero Jan 9 at 10:44
  • 1
    I'm using Hydrolinks ‘0.10.0’ running on R 3.6.1. Though you may be able to get the waterbody ID using nhdR with "ID <- nhd_query(Long, Lat, dsn = c("NHDWaterbody"))". I've been avoiding using nhdR, even though they have a function to find upstream networks ('extract_network') as I'm running on a Mac and nhdR requires the use of 7-zip to unpack NHD files – Vint Jan 9 at 12:22
  • 1
    But to make thing easier- I will edit the above question add the ID from Hydrolinks – Vint Jan 9 at 12:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.