Packages

The main packages I’ll be using for this blog are as follows:

• dplyr and ggplot2 for wrangling and plotting
• cowplot: an amazingly useful package for formatting/manipulating/making multi plots
• sf: core package for working with spatial data
• mapview: easy way to make html maps with sf objects
• smoothr: cool package for simplifying or smoothing polygons in R

Let’s load the packages:

# 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)

Option 1: Use existing zipped shapefile

# 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)

Option 2: Use function to download basin for USGS Station location

We’ll load a function that grabs the upstream watershed for a given station ID. Here it’s the SF Feather (or most of the watershed). We plot with mapview at the end and add both layers to the same map with a “+” similar to ggplot2.

First we need to load the function:

# 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)

Customizing the Maps: Using ggplot2

Let’s make a map of all these layers as-is, using ggplot2. The newest version of ggplot natively plots sf objects with geom_sf, so if this doesn’t work, try updating.

gg1b <- ggplot() +
  geom_sf(data=fea1, col="blue3", lwd=0.4) +
  geom_sf(data=fea_huc, fill=NA, color="gray20", lwd=1.2) +
  geom_sf(data=fea_wbodies, fill="cyan3") +
  theme_bw()
gg1b

Customizing a Mapview Map

Great, now let’s add a mapview map, and add some advanced trickery to make it look different. See here for some advanced mapview options.

# make a list of the map layers of interest, first layer is what gets plotted
maplayers <- c( "CartoDB.Positron","OpenTopoMap", "OpenStreetMap", "Esri.WorldImagery","Esri.WorldTopoMap", "Stamen.TonerLite")

mapview(fea1, zcol="streamorde", map.types=maplayers, layer.name=c(paste0("Stream <br> Order"))) + 
  mapview(fea_huc, color="gray", col.region=NA, lwd=2.4,map.types=maplayers) +
  mapview(fea_wbodies, col.region="cyan3", map.types=maplayers)

Crop by Watershed

First let’s crop the data by the watershed boundary…because originally these data were pulled in using a bounding box around the watershed (the X/Y coordinates at the corners of a box that completely encapsulates the polygon). Thankfully the sf package makes it easy to clip data down to the area of interest. From here on out I’ll just show the smaller SF Feather watershed and associated data. One important point: make sure the data are in the same projection before doing intersection/joining type operations.

# check crs:
st_crs(fea1)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

st_crs(fea_huc)

## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"

# make a cropped version:
fea_streamline_clipped <- st_intersection(fea1, fea_huc)
fea_lakes_clipped <- st_intersection(fea_wbodies, fea_huc) # nothing to clip

# plot
clipPlot <- ggplot() +
  geom_sf(data=fea_streamline_clipped, aes(color=streamorde), lwd=0.4) +
  geom_sf(data=fea_huc, fill=NA, color="gray20", lwd=1.2) +
  geom_sf(data=fea_lakes_clipped, fill="cyan3", color="blue3") +
  theme_bw() + scale_color_viridis_c("Stream \n Order")
clipPlot

# NOTE: not uncommon to get an Error:
# "Error in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  : 
#  polygon edge not found"
# SUGGESTION: Just try plotting again, or try dev.off() and then replot, 
# typically resolves.

Smooth the Boundaries!

One final step…we have a nice watershed boundary and a some lake boundaries which may benefit from smoothing or simplification using the smoothr and sf packages. There are many uses/purposes, so if you want to read more, check out the smoothr package page. There are lots of cool functions to play with, I’ll just show a couple here. But know there are a few that are helpful for cleaning/tidying data with holes, jagged edges, etc.

Simplify and Smooth

Note, we’ll be getting warnings here because of smoothing data that is in decimal degrees…easy to switch to a different transformation (i.e., CRS=3310) and rerun. Let’s use the watershed boundary. The first step I"m using is to simplify the shape boundary. Again, there will be some warnings which relate to the projection, no worries, it will still work. We could use many different levels of simplification here. I’m picking a value (dTolerance) that still roughly provides the shape. Feel free to play with different values and see what happens.

library(smoothr)

# simplify
fea_huc_simple <- st_simplify(fea_huc, preserveTopology = TRUE, dTolerance = .007)
#plot(fea_huc_simple$geometry, border="black", lwd=1.5)
#plot(fea_huc$geometry, border="blue", add=T)

# smooth
fea_huc_smooth1 <- smooth(fea_huc_simple, method = "ksmooth")
fea_huc_smooth2 <- smooth(fea_huc_simple, method = "spline")
fea_huc_smooth3 <- smooth(fea_huc_simple, method = "chaikin")

plot(fea_huc$geometry, border=scales::alpha("#FDE725FF", 0.8), lwd=3)
plot(fea_huc_smooth1$geometry, border="#440154FF", lwd=1, add=T)
plot(fea_huc_smooth2$geometry, border="#31688EFF", lwd=1,lty=2, add=T)
plot(fea_huc_smooth3$geometry, border="#35B779FF", lwd=2, lty=5, add=T)
legend(title= paste0("Different Smoothed \n HUC Outlines"),"bottomright", col=c("#FDE725FF", "#440154FF", "#31688EFF", "#35B779FF"), legend = c("original", "k-smooth","spline", "chaikin"), lty = c(1, 1, 2, 4), lwd = c(3, 1, 1, 2), bty = "n", cex=0.8)

Put it all together in a Map

Let’s use all these pieces to make one final map in ggplot2 and mapview to close this loop. I’ll add a trick where we label the lakes by adding a centroid point to the sf object, which does require the purrr package.

library(purrr)

# get COUNTY data for a given state
fea_lakes_clipped <- fea_lakes_clipped %>% 
  mutate(lon=map_dbl(geometry, ~st_centroid(.x)[[1]]), # add cent. for labels
         lat=map_dbl(geometry, ~st_centroid(.x)[[2]])) # add cent. for labels

Ok, here’s the ggplot version:

finalPlot1 <- ggplot() +
  geom_sf(data=fea_streamline_clipped, aes(color=streamorde), lwd=0.4) +
  geom_sf(data=fea_huc_smooth2, fill=NA, color="black", lwd=0.5) +
  geom_sf(data=fea_lakes_clipped, fill="cyan3", color="slateblue", alpha=0.8) +
  geom_label_repel(data=fea_lakes_clipped, aes(x=lon, y=lat, label=gnis_name), nudge_x = .08, nudge_y = -.01, alpha=0.65) +
  theme_bw() + 
  scale_color_viridis_c("Stream \n Order", option = "D", direction = -1,
                        limits = c(1,3), breaks = c(1, 2, 3),
                        guide = guide_colourbar(nbin=3, 
                                                draw.ulim = FALSE, 
                                                draw.llim = FALSE)) + 
  theme(legend.position = c(0.1, 0.77)) +
  labs(title = "South Feather River Basin", 
       subtitle = "(Spline smoothed watershed boundary)")
finalPlot1

finalPlot2 <- ggplot() +
  geom_sf(data=fea_streamline_clipped, aes(color=streamorde), lwd=0.4) +
  geom_sf(data=fea_huc, fill=NA, color="gray30", lwd=0.8, lty=2) +
  geom_sf(data=fea_lakes_clipped, fill="cyan3", color="slateblue", alpha=0.8) +
  geom_label_repel(data=fea_lakes_clipped, aes(x=lon, y=lat, label=gnis_name), nudge_x = .08, nudge_y = -.01, alpha=0.65) +
  theme_bw() + 
  scale_color_viridis_c("Stream \n Order", option = "D", direction = -1,
                        limits = c(1,3), breaks = c(1, 2, 3),
                        guide = guide_colourbar(nbin=3, 
                                                draw.ulim = FALSE, 
                                                draw.llim = FALSE)) + 
  theme(legend.position = c(0.1, 0.8)) +
  labs(title = "South Feather River Basin",
       subtitle = "(HUC watershed boundary)")
finalPlot2

And the mapview:

maplayers <- c( "CartoDB.Positron","OpenTopoMap", "OpenStreetMap", "Esri.WorldImagery","Esri.WorldTopoMap", "Stamen.TonerLite")

mapview(fea_streamline_clipped, zcol="streamorde", map.types=maplayers, layer.name=c(paste0("Stream <br> Order"))) + 
  mapview(fea_huc, color="gray", col.region=NA, lwd=2.4,map.types=maplayers) +
  mapview(fea_wbodies, col.region="cyan3", map.types=maplayers)