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:
library(hydrolinks) Lat<-37.70127 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
library(nhdR) 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.
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.
library(streamstats) #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
library(leaflet) library(nhdplusToos) #Step 1 point <- sf::st_sfc(sf::st_point(c(PLon, PLat)), crs = 4326) #Get nhdplus_id from point ID<-discover_nhdplus_id(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.