---
title: "Data expedition 1: spatial data and watershed boundaries"
author: "Kateri Salk, Nicholas Bruns"
date: "1/16/2020"
output: html_document
---
*This is part 1 of 4 in class sessions that will build up to using public hydroinformatic data sets and open source geospatial tools in R to perform spatial analyses of water quality measurements. Today is a warm up, to get all used to R, spatial data, and watersheds.*

#### Objectives by session:

  * display specific watershed boundaries using R and its geospatial features
	  
  * (optional) compare Hydrologic Unit Code nesting in different sized basins
  
  * link water quality data to spatial, "hydrography" datasets. 
  
  * analyze the spatial patterns in water quality measurements that result from fluvial processes

#### We'll use these datasets: 

  * [Watershed Boundary Dataset (WBD), 2.7 GB](http://prd-tnm.s3-website-us-west-2.amazonaws.com/?prefix=StagedProducts/Hydrography/WBD/National/GDB/). 
  
  * [National Hydrography Dataset + (NHD+), 7.3 GB](https://www.epa.gov/waterdata/nhdplus-national-data). select NHDPlusV21_NationalData_Seamless_Geodatabase_Lower48_07.7z  
  
  * [Water Quality Portal (WQP)](https://www.waterqualitydata.us/) via the [dataRetrieval package in R](https://cran.r-project.org/web/packages/dataRetrieval/vignettes/dataRetrieval.html).

#### Additional resources: 
  * [Geocomputation in R](http://prd-tnm.s3-website-us-west-2.amazonaws.com/?prefix=StagedProducts/Hydrography/WBD/National/GDB/). This is a great open source book on all things spatial and R. Everything I know is from here!

  * [HUC visualizer, built by Duke undergrads](http://waterqualityexplorer.rc.duke.edu:3838/explorer/). This is a prototype and may have problems if we all use it together! But.  is a wonderful visualization of the different resolution "quilts" made by the different HUC levels.
  
  

### Today: displaying watershed bounderies
*Homework, setup:* 

1) For reference, we've all already installed R on our machines along with the R packages sf, tidyverse, and tmap by opening R and running this code:

```{r, eval=FALSE}
install.packages(c("sf","tidyverse","tmap"))
```

2) download WBD data, then run this command in R:
```{r message=TRUE, warning=TRUE, paged.print=FALSE}
library(tidyverse)
library(sf) 
wbd_shape_dir <- "~/projects/wqp_hacking/data/WBD_National_GDB/WBD_National_GDB.gdb/" #adjust to wherever you put project data 
st_layers(wbd_shape_dir)
```
You'll need your run of this to command st_layers to produce this output by class time. If you can't, see Nick in his office hours.

*Class activity* 
We will produce a map of the mutiple, nested HUC polygons a point on the landscape sits within. Below, I've done this for my home watershed, and then repeated for Lawrence Kansas. Steps:  

1) Use google maps to get the lat long of a location (drop a pin, right click, and slect "what's here/")
2) Modify and run this command
```{r find_location}
location <- tibble(lat=36.054055, lon=-79.013609) #replace with you location
# location_sf <- st_as_sf(location,coords=c("lon","lat"),crs=4269) #turns a lat long into a "geographically aware", spatial entitity. crs is coordinate reference system
# location_sf <- st_as_sf(location,coords=c("lon","lat"),crs=4269) #
location_sf <- st_as_sf(location,coords=c("lon","lat"),crs=4269,remove = F) #debug update 1:  tmap now needs columns in additioni to geeom, so set remove to F
library(tmap)
tmap_mode("view")
tm_shape(location_sf) + tm_dots()
```

```{r find_locations_huc10, message=FALSE, warning=FALSE}
# huc_10_shapes <- st_read(wbd_shape_dir,layer = "WBDHU10")
# huc_10_shapes <- st_read(wbd_shape_dir,layer = "WBDHU10",drivers = "FileGDB")
huc_10_shapes <- st_read(wbd_shape_dir,layer = "WBDHU10",drivers = "OpenFileGDB")  #debug update 2: specify ESRI drivers. Also, supress warning messages in block (see above, warning=FALSE, message=FALSE)
glimpse(huc_10_shapes)
my_huc10_polygon<- huc_10_shapes[location_sf,] #extracts the polygon that contains your point, this may take a moment
tm_shape(my_huc10_polygon) + tm_borders(col="red") + tm_shape(location_sf) + tm_dots() # map it!
```

Now extract the HUC code to retrieve larger HUCs
```{r extract_the_huc_code}
glimpse(my_huc10_polygon)
my_huc10_code <- my_huc10_polygon %>% 
  pull(HUC10) %>% 
  as.character()  #don't convert to numeric! It will drop the leading 0
my_huc10_code
```
Now truncate this code to the larger, nested polygons
```{r subset_huc_code_for_large_hucs, message=FALSE, warning=FALSE}
my_huc08_code<- substr(my_huc10_code,1,8) 
my_huc08_polygon <- st_read(wbd_shape_dir,layer="WBDHU8",drivers = "OpenFileGDB") %>% 
  filter(HUC8==my_huc08_code) 
tm_shape(my_huc08_polygon) + tm_borders(col="blue") + 
  tm_shape(my_huc10_polygon) + tm_borders(col="red") + 
  tm_shape(location_sf) + tm_dots()
```
Keep going! Extract the polygons for hucs 6,4, and 2, then plot.
```{r extract_all_hucs, message=FALSE, warning=FALSE}
my_huc06_code<- substr(my_huc10_code,1,6) 
my_huc06_polygon <- st_read(wbd_shape_dir,layer="WBDHU6",drivers = "OpenFileGDB") %>%  
  filter(HUC6==my_huc06_code) 

my_huc04_code<- substr(my_huc10_code,1,4) 
my_huc04_polygon <- st_read(wbd_shape_dir,layer="WBDHU4",drivers = "OpenFileGDB") %>%  
  filter(HUC4==my_huc04_code) 

my_huc02_code<- substr(my_huc10_code,1,2) 
my_huc02_polygon <- st_read(wbd_shape_dir,layer="WBDHU2",drivers = "OpenFileGDB") %>%  
  filter(HUC2==my_huc02_code) 
```

Now plot them all together

```{r plot_all_hucs}
tm_shape(my_huc02_polygon) + tm_borders(col="black") + 
tm_shape(my_huc10_polygon) + tm_borders(col="red") + 
tm_shape(my_huc08_polygon) + tm_borders(col="blue") + 
tm_shape(my_huc06_polygon) + tm_borders(col="darkgreen") + 
tm_shape(my_huc04_polygon) + tm_borders(col="purple") + 
tm_shape(my_huc02_polygon) + tm_borders(col="black") + 
  tm_shape(location_sf) + tm_dots()
```

Now repeat with a new location in the Mississippi basin!

```{r everything_together, message=FALSE, warning=FALSE}
location <- tibble(lat=38.962891, lon=-95.245082) #replace with you location

location_sf <- st_as_sf(location,coords=c("lon","lat"),crs=4269,remove=F) #bugfix
huc_10_shapes <- st_read(wbd_shape_dir,layer = "WBDHU10",drivers = "OpenFileGDB")
my_huc10_polygon<- huc_10_shapes[location_sf,] #extracts the polygon that contains your point, this may take a moment
my_huc10_code <- my_huc10_polygon %>% 
  pull(HUC10) %>% 
  as.character()  #don't convert to numeric! It will drop the leading 0

my_huc08_code<- substr(my_huc10_code,1,8) 
my_huc08_polygon <- st_read(wbd_shape_dir,layer="WBDHU8",drivers = "OpenFileGDB") %>%  
  filter(HUC8==my_huc08_code) 

my_huc06_code<- substr(my_huc10_code,1,6) 
my_huc06_polygon <- st_read(wbd_shape_dir,layer="WBDHU6",drivers = "OpenFileGDB") %>%  
  filter(HUC6==my_huc06_code) 

my_huc04_code<- substr(my_huc10_code,1,4) 
my_huc04_polygon <- st_read(wbd_shape_dir,layer="WBDHU4",drivers = "OpenFileGDB") %>%  
  filter(HUC4==my_huc04_code) 

my_huc02_code<- substr(my_huc10_code,1,2) 
my_huc02_polygon <- st_read(wbd_shape_dir,layer="WBDHU2",drivers = "OpenFileGDB") %>%  
  filter(HUC2==my_huc02_code) 
```

```{r}
tm_shape(my_huc02_polygon) + tm_borders(col="black") + 
tm_shape(my_huc04_polygon) + tm_borders(col="purple") + 
tm_shape(my_huc06_polygon) + tm_borders(col="darkgreen") + 
tm_shape(my_huc08_polygon) + tm_borders(col="blue") + 
  tm_shape(my_huc10_polygon) + tm_borders(col="red") + 
  tm_shape(location_sf) + tm_dots()
```