title | week | type | subtitle | reading | tasks | |||||
---|---|---|---|---|---|---|---|---|---|---|
Beware the Canadians! |
5 |
Case Study |
Working with Spatial Data and the sf package |
|
|
- The Spatial Features Package (sf){target='blank'}
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.
You woke up in the middle of the night terrified of the Canadians after a bad dream. You decide you need to set up military bases to defend the Canada-NY border. After you tweet your plans, you realize you have no plan. What will you do next?
- 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?
- 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
Download starter R script (if desired){target="_blank"}
You will need to load the following packages
library(spData)
library(sf)
library(tidyverse)
# library(units) #this one is optional, but can help with unit conversions.
#load 'world' data from spData package
data(world)
# load 'states' boundaries from spData package
data(us_states)
# plot(world[1]) #plot if desired
# plot(us_states[1]) #plot if desired
-
world
dataset- transform to the albers equal area projection:
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"
- use
st_set_geometry()
to specify that thegeom
column contains thegeometry
. This will also rename the column fromgeom
togeometry
to make it easier to useggplot()
- filter the world dataset to include only
name_long=="Canada"
- buffer canada to 10km (10000m)
-
us_states
object- transform to the albers equal area projection defined above as
albers
- filter the
us_states
dataset to include onlyNAME == "New York"
- transform to the albers equal area projection defined above as
-
Create a 'border' object
- use
st_intersection()
to intersect the canada buffer with New York (this will be your final polygon) - Plot the border area using
ggplot()
andgeom_sf()
. - use
st_area()
to calculate the area of this polygon. - Convert the units to km^2. You can use
set_units(km^2)
(from theunits
library) or some other method.
- use
-
Do not worry about small waterways, etc. Just use the two datasets listed above.
Your final result should look something like this:
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. Please do not use these data for real military purposes.
Build a leaflet map of the same dataset.
<iframe id="test" style=" height:400px; width:100%;" scrolling="no" frameborder="0" src="CS05_leaflet.html"></iframe>