In tylerbhampton/watershedtools: Delineating Watersheds with R

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  message = FALSE
)

watershedtools

The goal of watershedtools is to delineate watersheds using packages like nhdplusTools and others. Most methods will work best for the continental U.S.

Installation

You can install watershedtools from github using devtools:

devtools::install_github("tylerbhampton/watershedtools")

Test Watershed Delineation with nhdplusTools

library(watershedtools)
library(ggplot2)
library(leaflet)
library(sf)

# Point for Harlan County Lake on the Republican River, Nebraska USA
p=c(-99.189990,40.071521)

wsshape=NHDWBDws(method = "point",point = p,
         WBDstagepath = "D:/clone/USA_WBD_WatershedBoundaryDataset")


leaflet(data = wsshape) %>% leaflet::addPolygons() %>% leaflet::addProviderTiles("Esri.WorldImagery")

Test Watershed Delineation with RQGIS3

Note, you need OSGeo4W installed on your computer to use this package. See the installation tutorial in the OSGeo folder on Github.

Install R Packages

devtools::install_github("jannes-m/RQGIS3")
devtools::install_github("Nowosad/spData")
devtools::install_github("Nowosad/spDataLarge")

Start up RStudio

suppressPackageStartupMessages({
library("raster")
library("dplyr")
library("sf")
library("ggplot2")
library("RQGIS3")
})
set_env("C:/OSGeo4W64")

qgis_session_info()
version

R will identify for you the current version of QGIS running on your computer.

Test using Zion National Park Data

We will use a digital elevation model (DEM) of Zion National Park, Utah, USA, to construct watersheds. The geographic system is in latitude-longitude (WGS84) and elevation is in meters. The most apparent feature on the landscape is the Virgin River, which drains the plateau in the northeast corner of the map and exists to the southwest.

dem=spDataLarge::elevation
plot(dem)

RQGIS functionality

Our first task for watershed deliniation is to fill the sinks in the dem, to make sure all water drains out. We use the find_algorithms function to query the python libraries. We can simply add the string "sinks".

find_algorithms("sinks")

We will use the saga library fillsinks function.

Next, the get_usage and get_args_man functions describe the inputs, parameters, and outputs of the functions. We see that the required input is DEM, a dem, and the output is named RESULT, which gives us our filled dem.

get_usage("saga:fillsinks")
get_args_man("saga:fillsinks")

Finally, we can integrate these with the run_qgis function, defining the algorithm and its parameters. We can supply a file path for the output. Saga packages can only output as .sdat files, not .tif. We will use .tif later.

temppath=tempdir()

demfill = run_qgis(alg = "saga:fillsinks",
                   DEM=dem,
                   RESULT=file.path(temppath, "demfill.sdat"),
                   load_output = FALSE,
                   show_output_paths=FALSE)

run_qgis(alg = "grass7:r.watershed",
         elevation=file.path(temppath, "demfill.sdat"),
         threshold=1000,
         accumulation=file.path(temppath, "acc.tif"),
         drainage=file.path(temppath, "flowdir.tif"),
         load_output = FALSE,
         show_output_paths=FALSE)

flowacc=raster(file.path(temppath, "acc.tif"))
flowacc=calc(flowacc,function(x){ifelse(x>0,log10(x),NA)})
names(flowacc)="acc"
plot(flowacc)

Plotting the log of flow accumulation, we can see several rivers shown in green with high accumulation. For a dem, any cell that accumulates water from outside the boundary is negative.

RQGIS functionality

We will use these calculated raster products to draw watersheds from several points on the landscape. First, we will calculate a list of points on the raster where flow accumulation is greater than our threshold (1000), and then use lapply to snap each of our points to the raster coordinates.

zionpoints=list(
  c(-113.1634, 37.2176),
  c(-113.1,37.22833),
  c(-113.17,37.23083),
  c(-113.11,37.253),
  c(-112.954,37.3525)
)

r.pts <- as.data.frame(rasterToPoints(flowacc))
subset=r.pts[r.pts$acc>3,]

points.s=lapply(1:length(zionpoints),function(i){
  p=zionpoints[[i]]
  row=which.min(abs(subset$x-p[1])+abs(subset$y-p[2]))
  return(subset[row,1:2])
})

Next, for each point, we will use the flow direction raster to create a watershed boundary, and then convert that raster to a vector shapefile.

for(i in 1:length(zionpoints)){
  run_qgis(alg = "grass7:r.water.outlet",
           input=file.path(temppath, "flowdir.tif"),
           coordinates=paste(as.vector(unlist(points.s[[i]])),collapse=","),
           output=file.path(tempdir(), paste0("basin",i,".tif")),
           load_output = FALSE,
           show_output_paths=FALSE)
  run_qgis(alg="grass7:r.to.vect",
           input=file.path(temppath, paste0("basin",i,".tif")),
           type=2L,
           output=file.path(temppath, paste0("basin",i,".shp")),
           load_output = FALSE,
           show_output_paths=FALSE)
}

Finally, we will convert our points to a shape layer, and read in all our basin shapefiles as one joined layer. The boundary of zion national park comes with the spDataLarge package. And to add aesthetic to the map, we will make a hillshade baselayer.

point.d=do.call("rbind",points.s)%>%as.data.frame()
names(point.d)=c("lon","lat")
coordinates(point.d)=c("lon","lat")
crs(point.d)=proj4string(dem)

basins.d=do.call("rbind",lapply(
  1:length(zionpoints),function(i){
    l=st_read(file.path(tempdir(), paste0("basin",i,".shp")),quiet = T)
    l$value=i
    l=l[1,]
    return(l)
  }
))
st_crs(basins.d)=st_crs(zion)

zion=system.file("vector/zion.gpkg", package = "spDataLarge")%>%
  st_read()%>%st_transform(.,proj4string(dem))
hillshade = run_qgis(alg = "grass7:r.relief",
                     input=dem,
                     output=file.path(tempdir(), "hillshade.tif"),
                     load_output = TRUE,
                     show_output_paths = FALSE)


ggplot()+
  geom_tile(data=as.data.frame(rasterToPoints(hillshade)),aes(x=x,y=y,fill=hillshade))+
  scale_fill_gradientn(colors=c("white",1))+
  geom_tile(data=subset(as.data.frame(rasterToPoints(flowacc)),acc>1000 | acc<{-100}),aes(x=x,y=y))+
  geom_sf(data=zion,col="white",fill="transparent")+
  geom_sf(data=basins.d,aes(col=factor(value)),fill="transparent")+
  geom_point(data=as.data.frame(point.d),aes(x=lon,y=lat),size=2,col=2)+
  ggsn::blank()+theme(legend.position = "n")+
  ggtitle("Zion National Park\nNorth Creek Watershed")

Compare RQGIS3 and nhdplustools

zionNHDwatersheds=do.call("rbind",lapply(1:5,function(i){
  NHDWBDws(method = "point",point = zionpoints[[i]],
         WBDstagepath = "B:/MASTERDATA/data_Hydrology/hydrology_shapefiles/USA_WBD_WatershedBoundaryDataset")
}))
zionNHDwatersheds$value=1:5


ggplot()+
  geom_tile(data=as.data.frame(rasterToPoints(hillshade)),aes(x=x,y=y,fill=hillshade))+
  scale_fill_gradientn(colors=c("white",1))+
  geom_tile(data=subset(as.data.frame(rasterToPoints(flowacc)),acc>1000 | acc<{-100}),aes(x=x,y=y))+
  geom_sf(data=zion,col="white",fill="transparent")+
  geom_sf(data=zionNHDwatersheds,aes(col=factor(value)),fill="transparent")+
  geom_point(data=as.data.frame(point.d),aes(x=lon,y=lat),size=2,col=2)+
  ggsn::blank()+theme(legend.position = "n")+
  ggtitle("Zion National Park\nNorth Creek Watershed")