-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlab8_oilSpills.Rmd
133 lines (79 loc) · 2.39 KB
/
lab8_oilSpills.Rmd
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
---
title: 'Intro To Git: Oil Spills'
author: "Ian Ladner"
date: "March 8, 2018"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
```{r loading}
library(tidyverse)
library(sf)
library(spatstat)
library(tmap)
library(leaflet)
library(maptools)
oil <- read_csv("data/oil_spills.csv")
texas <- filter(oil, `Accident State` == "TX" & `Accident Year` < 2017) %>%
group_by(`Accident Year`) %>%
summarise(Loss = sum(`Net Loss (Barrels)`))
colnames(texas)<- c("Year","Loss")
ggplot(texas, aes(x = Year, y = Loss))+geom_col()+theme_classic()
```
```{r}
df_loc <- oil %>%
filter(`Accident State` == "TX" & `Accident Year` == 2016) %>%
select(Latitude, Longitude, `Net Loss (Barrels)`)
colnames(df_loc) <- c("Lat","Long","Loss")
oil_sf <- st_as_sf(df_loc, coords = c("Long","Lat"), crs = 4326)
leaflet(oil_sf) %>%
addTiles() %>%
addMarkers()
```
```{r tmap}
states <- st_read(dsn = "data", layer = "states")
tex_border <- states %>%
filter(STATE_NAME == "Texas") %>%
st_transform(4326)
plot(tex_border)
tm_shape(tex_border)+
tm_polygons()+
tm_shape(oil_sf)+
tm_dots(size = 0.7)
```
Point pattern analysis. R is finicky.
```{r spatialAnalysis}
spill_sp <- as(oil_sf, "Spatial")
spill_ppp <- as(spill_sp, "ppp")
#but what's our bounding window
tx_sp <- as(tex_border, "Spatial")
tx_owin <- as(tx_sp, "owin")
#but outer window and points together
all_ppp <- ppp(spill_ppp$x, spill_ppp$y, window = tx_owin)
```
##A density plot
```{r}
plot(density(all_ppp, sigma = 0.4))
```
```{r quadrantTest}
#nx = number of horz subdivisions, ny = number of vert subdivisions
#null hyp: CSR is demonstrated by data
oil_qt <- quadrat.test(all_ppp, nx = 5, ny = 5)
#What do these numbers mean? Ask Allison.
plot(all_ppp)
plot(oil_qt, add = TRUE, cex = 0.5)
```
G function for nearest neighbor analysis
```{r Gfunction}
r <- seq(0,1,0.1)
oil_gFun <- envelope(all_ppp, fun = Gest, r = r, nsim = 100)
ggplot(oil_gFun, aes(x = r, y = obs))+geom_line(color = "black")+geom_line(aes(x = r, y = theo),color = "red")
```
Nearest neighbor using the L (or K) function which looks at concentration of points on the whole
```{r}
lag2 <- seq(0,3,0.5)
oil_lFun <- envelope(all_ppp, fun = Lest, r = lag2, nsim = 20, global = TRUE)
ggplot(oil_lFun, aes(x = r, y = obs))+geom_line(color = "coral3")+geom_line(aes(x = r, y = theo), color = "darkgreen")
#clustered
```