Mapping with sf

About a year ago, I wrote a series of posts on using the sf and riverdist packages (among others) to do some stream analysis and make stream maps. Since that time, I feel like the capabilities and content for spatial work in R has exploded. There have been many excellent examples, and kudos to the folks tweeting, posting, and sharing all the code! If you haven’t seen some of the really fantastic visualizations and data wrangling that the USGS OWI folks have been putting out (along with the code!), you’re missing out. Same goes for many others who have posted some great stuff. Here’s just a few interesting examples:

• Groundwater levels affected by Mexico earthquake
• Hurricane Harvey
• California Drought
• Beyond basic mapping (USGS)
• Spatial Pipelines with Congressional Districts

Main objectives of the post:

Assume all of these will be using sf:

• How to grab NHD Streamline data using any spatial polygon
• How to grab waterbodies for a given spatial polgyon
• How to quickly add these to a interactive map in mapview or plot with ggplot
• How to modify a polygon (e.g., a species range) with the smoothr package

# load libraries
library(dplyr) # data munging and piping
library(ggplot2) # plotting
library(sf) # spatial everything
library(mapview) # html mapping
library(smoothr) # smoothing spatial objects
library(purrr); # wrangling lists and such
library(ggrepel) # nice labels in ggplot
#library(sbtools)
#library(dataRetrieval)

Get a HUC Watershed Boundary

First, let’s use one of the NHD HUC boundaries that we can use as our background for grabbing other data. I’m just using something I have handy, which is HUC8 boundaries for the American/Yuba/Bear watersheds. Alternatively, thanks to the awesome folks over at USGS OWI (David Watkins, David Blodgett, Laura DeCicco, and many more), there’s already plenty of code we can adapt to use to download all sorts of NHD/Watershed things.

For example, we can use some of the code from very handy function courtesy of Laura DeCicco, in a github repository for an inactive package called hydroMap. To make that work, you’ll need a local USGS gaging station ID, as well as rgdal installed. If you don’t know a NWIS or USGS station ID, you can look them up via the USGS website. For this example, I’ll use 11396200 which is on the SF Feather River.

# read in zipped shapefile with sf
huc8 <- read_sf(unzip("data/h8_AMR_BEA_YUB.zip"), quiet = F) %>%
  st_transform(crs=4326) #%>% 

## Reading layer `h8_AMR_BEA_YUB' from data source `./h8_AMR_BEA_YUB.dbf' using driver `ESRI Shapefile'
## Simple feature collection with 4 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -121.5976 ymin: 38.61455 xmax: -119.9829 ymax: 39.77696
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

# then remove raw files since file is added in memory
file.remove(list.files(pattern = "h8_AMR_BEA_YUB*",recursive = F))

## [1] TRUE TRUE TRUE TRUE

# filter to the Bear only
h8 <- huc8 %>% filter(HU_8_NAME=="Upper Bear") 

plot(h8$geometry, border="orange", lwd=2)

# the original function is here: https://github.com/USGS-R/hydroMap/blob/master/R/plotWSB.R

# function assumes you have sf installed locally
get_basins <- function(sites, filePath = NA){
  library(httr)
  postURL <- "https://cida.usgs.gov/nwc/geoserver/NWC/ows"
  # postURL <- "http://cida-test.er.usgs.gov/nwc/geoserver/NWC/ows"
  filterXML <- paste0('<?xml version="1.0"?>',
                      '<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="shape-zip" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">',
                      '<wfs:Query xmlns:feature="https://owi.usgs.gov/NWC" typeName="feature:epa_basins" srsName="EPSG:4326">')
  
  
  if(length(sites) > 1){
    siteText <- ""
    for(i in sites){
      siteText <- paste0(siteText,'<ogc:PropertyIsEqualTo  matchCase="true">',
                         '<ogc:PropertyName>site_no</ogc:PropertyName>',
                         '<ogc:Literal>',i,'</ogc:Literal>',
                         '</ogc:PropertyIsEqualTo>')
    }
    
    filterXML <- paste0(filterXML,'<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">',
                        '<ogc:Or>',siteText,'</ogc:Or>',
                        '</ogc:Filter>')
    
  } else {
    filterXML <- paste0(filterXML,
                        '<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">',
                        '<ogc:PropertyIsEqualTo matchCase="true">',
                        '<ogc:PropertyName>site_no</ogc:PropertyName>',
                        '<ogc:Literal>',sites,'</ogc:Literal>',
                        '</ogc:PropertyIsEqualTo>',
                        '</ogc:Filter>')
  }
  
  filterXML <- paste0(filterXML,'</wfs:Query>',
                      '</wfs:GetFeature>')
  
  destination = tempfile(pattern = 'basins_shape', fileext='.zip')
  
  file <- POST(postURL, body = filterXML, write_disk(destination, overwrite=T))
  if(is.na(filePath)){
    filePath <- tempdir()
  }
  
  unzip(destination, exdir = filePath)
  basins <- st_read(filePath, layer='epa_basins')
  
  #basins = readOGR(filePath, layer='epa_basins')
  return(basins)
}

# set a station ID or IDs:
ids <- c(11396200) # S Feather
fea_huc <- get_basins(ids) # use function to get boundary

## Reading layer `epa_basins' from data source `/private/var/folders/jv/56h34hjn2l302ypnnsphwnx00000gn/T/Rtmpa0CSzl' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 4 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -121.2132 ymin: 39.54753 xmax: -120.8683 ymax: 39.78771
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

# use mapview to show both
mapview(fea_huc) + mapview(h8)

Get NHD Streamline

I’ve shown this in previous posts, but here’s a slight modification. Now we can simply provide our polygons, and the function will automagically create the bounding box to grab flowlines for whatever level of detail (stream order) we choose. Higher stream order = more detail, with 1 being the finest scale you can get. First we load the function, again courtesy of the USGS folks. Note, the view_polygon should be a sf polygon, which gets transformed to a WGS84 lat/lon projection before downloading.

Here’s the function:

# the function
get_flowlines <- function(streamorder, view_polygon){
  library(httr)
  
  bbox <- st_bbox(st_transform(view_polygon, 4326))
  
  postURL <- "https://cida.usgs.gov/nwc/geoserver/nhdplus/ows"
  
  filterXML <- paste0('<?xml version="1.0"?>',
                      '<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="shape-zip" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">',
                      '<wfs:Query xmlns:feature="https://gov.usgs.cida/nhdplus" typeName="feature:nhdflowline_network" srsName="EPSG:4326">',
                      '<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">',
                      '<ogc:And>',
                      '<ogc:PropertyIsGreaterThan>',
                      '<ogc:PropertyName>streamorde</ogc:PropertyName>',
                      '<ogc:Literal>',streamorder-1,'</ogc:Literal>',
                      '</ogc:PropertyIsGreaterThan>',
                      '<ogc:BBOX>',
                      '<ogc:PropertyName>the_geom</ogc:PropertyName>',
                      '<gml:Envelope>',
                      '<gml:lowerCorner>',bbox[2]," ",bbox[1],'</gml:lowerCorner>',
                      '<gml:upperCorner>',bbox[4]," ",bbox[3],'</gml:upperCorner>',
                      '</gml:Envelope>',
                      '</ogc:BBOX>',
                      '</ogc:And>',
                      '</ogc:Filter>',
                      '</wfs:Query>',
                      '</wfs:GetFeature>')
  
  destination = file.path(tempdir(),"nhdflowline_network.zip")
  file <- POST(postURL, body = filterXML, write_disk(destination, overwrite=T))
  
  filePath <- tempdir()
  print("unzipping...")
  unzip(destination, exdir = filePath)
  
  flowLines <- sf::st_read(filePath, layer = 'nhdflowline_network')
  
  return(flowLines)
}

Then we grab our data…here at two different scales, stream order = 1, and stream order = 3.

# use the function to get streamlines
fea1 <- get_flowlines(1, fea_huc) # stream order 1

## [1] "unzipping..."
## Reading layer `nhdflowline_network' from data source `/private/var/folders/jv/56h34hjn2l302ypnnsphwnx00000gn/T/Rtmpa0CSzl' using driver `ESRI Shapefile'
## Simple feature collection with 397 features and 90 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: -121.2286 ymin: 39.52225 xmax: -120.8529 ymax: 39.8031
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

fea3 <- get_flowlines(3, fea_huc) # stream order 3

## [1] "unzipping..."
## Reading layer `nhdflowline_network' from data source `/private/var/folders/jv/56h34hjn2l302ypnnsphwnx00000gn/T/Rtmpa0CSzl' using driver `ESRI Shapefile'
## Simple feature collection with 84 features and 90 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: -121.2165 ymin: 39.52225 xmax: -120.8667 ymax: 39.79227
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

# Plot the Feather River
plot(fea1$geometry,col="skyblue") # plot stream order 1
plot(fea3$geometry, add=T,col="blue3", lwd=1.4) # stream order 3
# add HUC with some transparency using the scale::alpha() function
plot(fea_huc$geometry, border=scales::alpha("black", 0.5), lwd=2, add=T) 
title(main = "SF Feather Watershed with NHD Streamlines")

# Now same with the Bear River
bear1 <- get_flowlines(1, h8)

## [1] "unzipping..."
## Reading layer `nhdflowline_network' from data source `/private/var/folders/jv/56h34hjn2l302ypnnsphwnx00000gn/T/Rtmpa0CSzl' using driver `ESRI Shapefile'
## Simple feature collection with 2189 features and 90 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: -121.5895 ymin: 38.8825 xmax: -120.6062 ymax: 39.36088
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

bear3 <- get_flowlines(3, h8)

## [1] "unzipping..."
## Reading layer `nhdflowline_network' from data source `/private/var/folders/jv/56h34hjn2l302ypnnsphwnx00000gn/T/Rtmpa0CSzl' using driver `ESRI Shapefile'
## Simple feature collection with 828 features and 90 fields
## geometry type:  LINESTRING
## dimension:      XYZ
## bbox:           xmin: -121.5895 ymin: 38.91882 xmax: -120.6307 ymax: 39.36088
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

plot(bear1$geometry,col="skyblue")
plot(bear3$geometry, add=T,col=scales::alpha("blue3", 0.8), lwd=1.4)
# add HUC with some transparency using the scale::alpha() function
plot(h8$geometry, border=scales::alpha("maroon", 0.8), lwd=2, add=T)
title(main = "Bear River Watershed with NHD Streamlines")

Get Waterbodies

Great, we have a watershed boundary and a streamline. But let’s add some lakes. Same steps, load function, provide polygon/shapefile of interest (in sf format), a threshold for the minimum area (sq km) waterbody you want, and then wait patiently for something to happen.

Here’s the function:

get_waterbodies <- function(view_polygon, fetch_waterbody_areasqkm) {
  # view_polygon, ind_file were in function 
  bbox <- st_bbox(st_transform(view_polygon, 4326))
  
  postURL <- "https://cida.usgs.gov/nwc/geoserver/nhdplus/ows"
  
  filterXML <- paste0('<?xml version="1.0"?>',
                      '<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="http://gov.usgs.cida/nhdplus" typeName="feature:nhdwaterbody" srsName="EPSG:4326">',
                      '<ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">',
                      '<ogc:And>',
                      '<ogc:PropertyIsGreaterThan>',
                      '<ogc:PropertyName>areasqkm</ogc:PropertyName>',
                      '<ogc:Literal>',fetch_waterbody_areasqkm,'</ogc:Literal>',
                      '</ogc:PropertyIsGreaterThan>',
                      '<ogc:BBOX>',
                      '<ogc:PropertyName>the_geom</ogc:PropertyName>',
                      '<gml:Envelope>',
                      '<gml:lowerCorner>',bbox[2]," ",bbox[1],'</gml:lowerCorner>',
                      '<gml:upperCorner>',bbox[4]," ",bbox[3],'</gml:upperCorner>',
                      '</gml:Envelope>',
                      '</ogc:BBOX>',
                      '</ogc:And>',
                      '</ogc:Filter>',
                      '</wfs:Query>',
                      '</wfs:GetFeature>')
  
  out <- httr::POST(postURL, body = filterXML)
  
  sf_waterbodies <- read_sf(rawToChar(out$content)) %>%
    st_transform(st_crs(view_polygon))
  
}

Now let’s use the function!

# GET LAKES with 0.2 sqkm area or greater
bear_wbodies<-get_waterbodies(h8, 0.2)

# take a look at the plot:
plot(bear1$geometry,col="skyblue")
plot(bear3$geometry, add=T,col=scales::alpha("blue3", 0.8), lwd=1.4)
# add HUC with some transparency using the scale::alpha() function
plot(h8$geometry, border=scales::alpha("maroon", 0.8), lwd=2, add=T) 
plot(bear_wbodies$geometry, col="cyan4", border=scales::alpha("blue3",0.8), add=T)
title("Bear River Watershed with Lakes")

This time lets use a mapview map.

# do the same thing for feather
fea_wbodies<-get_waterbodies(fea_huc, 0.1)

# let's make a big mapview map
mapview(fea_huc, border="blue", col.region=NA) + 
  mapview(fea1, color="blue3") +
  mapview(fea_wbodies, col.region="blue3")