forked from AdamWilsonLabEDU/SpatialDataScience
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathCS_05.Rmd
More file actions
155 lines (116 loc) · 4.78 KB
/
CS_05.Rmd
File metadata and controls
155 lines (116 loc) · 4.78 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
---
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)
```
[<i class="fa fa-file-code-o fa-1x" aria-hidden="true"></i> Download starter R script (if desired)](`r output_nocomment`){target="_blank"}
<div class="well">
<button data-toggle="collapse" class="btn btn-primary btn-sm round" data-target="#demo1">Show Hints</button>
<div id="demo1" class="collapse">
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.
</div>
</div>
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.
<div class="extraswell">
<button data-toggle="collapse" class="btn btn-link" data-target="#extras">
Extra time? Try these extra activities...
</button>
<div id="extras" class="collapse">
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)
```
<iframe id="test" style=" height:400px; width:100%;" scrolling="no" frameborder="0" src="CS05_leaflet.html"></iframe>
</div>
</div>