---
title: Spatial Data
date: 2018-09-25
subtitle: Working with Spatial Data and the sf package
reading:
- The [Spatial Features Package (sf)](https://r-spatial.github.io/sf/){target='blank'}
tasks:
- Reproject spatial data using `st_transform()`
- Perform spatial operations on spatial data (e.g. intersection and buffering)
- Generate a polygon that includes all land in NY that is within 10km of the Canadian border and calculate the area
- Save your script as a .R or .Rmd in your course repository
---
```{r setup, include=FALSE, purl=F}
source("functions.R")
source("knitr_header.R")
```
# Reading
```{r reading,results='asis',echo=F,purl=F}
md_bullet(rmarkdown::metadata$reading)
```
# Background
Up to this point, we have dealt with data that fits into the tidy format without much effort. Spatial data has many complicating factors that have made handling spatial data in R complicated. Big strides are being made to make spatial data tidy in R.
# Objective
You woke up in the middle of the night after a bad dream terrified of the Canadians. You decide you need to set up military outposts to defend the border. After you tweet your plans, you realize you have no plan. What will you do next?
> 1) Generate a polygon that includes all land in NY that is within 10km of the Canadian border and 2) calculate it's area in km^2. How much land will you need to defend from the Canadians?
# Tasks
```{r tasks,results='asis',echo=F, purl=F}
md_bullet(rmarkdown::metadata$tasks)
```
[ Download starter R script (if desired)](`r output_nocomment`){target="_blank"}
The details below describe one possible approach.
You will need to load the foillowing packages
```{r, message=FALSE}
library(maps)
library(ggplot2)
library(spData)
library(dplyr)
library(sf)
# library(units) #this one is optional, but can help with unit conversions.
```
And datasets:
```{r}
#load 'world' data from spData package
data(world)
# load 'states' boundaries
states <- sf::st_as_sf(maps::map("state", plot = FALSE, fill = TRUE))
plot(world[1])
plot(states)
```
## Projection
```{r}
albers="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs "
```
1. `world` dataset
1. transform to the albers equal area projection defined above as `albers`
2. rename the `geom` column as `geometry` to make it easier to use `ggplot()`
3. filter the world dataset to include only `name_long=="Canada"`
4. buffer canada to 10km (10000m)
2. `states` object
1. transform to the albers equal area projection defined above as `albers`
2. filter the `states` dataset to include only `ID == "new york"`
3. confirm that the filtered object is valid with `st_is_valid()`, if not, use `st_buffer(0)` to apply a zero-width buffer which can solve most validity problems
3. create a 'border' object
1. use `st_intersection()` to intersect the canada buffer with New York (this will be your final polygon)
2. use `st_area()` to calculate the area of this polygon.
Your final result should look like this:
```{r purl=F, echo=F, message=FALSE, warning=FALSE}
canada <- world%>%
st_transform(crs=albers)%>%
rename(geometry=geom)%>%
filter(name_long=="Canada")%>%
st_buffer(10000)
ny=filter(states,ID == "new york")%>%
st_transform(crs=albers)%>%
st_buffer(0)
ny_border <- canada%>%
st_intersection(ny)
ggplot()+
geom_sf(data=ny)+
geom_sf(data=ny_border,fill="red")+
labs(title="New York Land within 10km of Canada")
library(units)
area <-
st_area(ny_border)%>%
set_units(km^2)%>%
round(1)
```
The total area to be defended is `r area` km^2.
Important note: This is a crude dataset meant simply to illustrate the use of intersections and buffers. The two datasets are not adequate for a highly accurate analysis.
Build a leaflet map of the same dataset.
```{r, purl=F, echo=F}
library(leaflet)
library(htmlwidgets)
l=ny_border%>%
st_transform("+proj=longlat +datum=WGS84")%>%
leaflet() %>%
addTiles() %>%
addPolygons(color = "#444444", weight = 1, fillOpacity = 0.5,
fillColor = "red",
popup = paste("Danger Zone"))
f <-"CS05_leaflet.html"
saveWidget(l,file.path(normalizePath(dirname(f)),basename(f)),libdir="externals",selfcontained = T)
```