-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy path07-case-study-climate.Rmd
142 lines (120 loc) · 5.39 KB
/
07-case-study-climate.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
133
134
135
136
137
138
139
---
knit: bookdown::preview_chapter
---
# Climate change
## Carbon dioxide
### About the data
### Temporal trend
***Note to self: Need to update the code to use readr, and get all the stations, acknowledge Halldor***
***HH: we made the wilydata package - the stations and the co2 flask data are included***
```{r load_packages-climate, cache=FALSE, echo = FALSE, message = FALSE, warning = FALSE, results='hide'}
library(gridExtra)
library(dplyr)
library(ggplot2)
library(lubridate)
library(rworldmap)
library(ggmap)
library(wilydata)
# remotes::install_github("heike/wilydata")
```
Get the data:
```{r CO2-data, eval=FALSE}
CO2.ptb<-read.table("https://scrippsco2.ucsd.edu/assets/data/atmospheric/stations/merged_in_situ_and_flask/daily/daily_merge_co2_ptb.csv", sep=",", skip=69)
colnames(CO2.ptb)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
CO2.ptb$lat<-71.3
CO2.ptb$lon<-(-156.6)
CO2.ptb$stn<-"ptb"
```
```{r CO2-0, out.width='99%', fig.align='center', fig.cap='Carbon dioxide levels from each of the measurements: (left) stations ordered from north to south, (right) overlaid.', warning=FALSE, message=FALSE, echo=FALSE, cache=FALSE}
# CO2.ptb<-read.table("http://scrippsco2.ucsd.edu/sites/default/files/data/flask_co2_and_isotopic/daily_co2/fldav_ptb.csv", sep=",", skip=69)
# colnames(CO2.ptb)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
# CO2.ptb$lat<-71.3
# CO2.ptb$lon<-(-156.6)
# CO2.ptb$stn<-"ptb"
#
# CO2.ljo<-read.table("http://scrippsco2.ucsd.edu/sites/default/files/data/flask_co2_and_isotopic/daily_co2/fldav_ljo.csv", sep=",", skip=69)
# colnames(CO2.ljo)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
# CO2.ljo$lat<-32.9
# CO2.ljo$lon<-(-117.3)
# CO2.ljo$stn<-"ljo"
#
# #CO2.mlf<-read.table("http://scrippsco2.ucsd.edu/sites/default/files/data/flask_co2_and_isotopic/daily_co2/fldav_mlf.csv", sep=",", skip=69)
# #colnames(CO2.mlf)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
# #CO2.mlf$lat<-19.5
# #CO2.mlf$lon<-(-155.6)
# #CO2.mlf$stn<-"mlf"
#
# # South pole doesn't have daily flask value (Aug 15 2018, HH)
# # CO2.spo<-read.table("http://scrippsco2.ucsd.edu/sites/default/files/data/flask_co2_and_isotopic/daily_co2/fldav_spo.csv", sep=",", skip=69)
# # colnames(CO2.spo)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
# # CO2.spo$lat<- (-90.0)
# # CO2.spo$lon<-0
# # CO2.spo$stn<-"spo"
#
# CO2.ker<-read.table("http://scrippsco2.ucsd.edu/sites/default/files/data/flask_co2_and_isotopic/daily_co2/fldav_ker.csv", sep=",", skip=69)
# colnames(CO2.ker)<-c("date", "time", "day", "decdate", "n", "flg", "co2")
# CO2.ker$lat<-(-29.2)
# CO2.ker$lon<-(-177.9)
# CO2.ker$stn<-"ker"
#
# CO2.all<-rbind(CO2.ker,CO2.ljo,CO2.ptb) #,CO2.spo)
# CO2.all$date<-as.Date(CO2.all$date)
#
# CO2.all$invlat=-1*CO2.all$lat
# CO2.all$stn=reorder(CO2.all$stn,CO2.all$invlat)
#
# CO2.all.loc <- rbind(CO2.ker[1,],CO2.ljo[1,],CO2.ptb[1,]) #,CO2.spo[1,])
#p1 <- qplot(date, co2, data=subset(CO2.all, flg < 2), colour=stn, geom="line",xlab="Year",ylab="CO2 (ppm)") +
# facet_wrap(~stn, ncol=1) + theme(axis.text.y=element_text(size = 6), legend.position="none")
#p2 <- qplot(date, co2, data=subset(CO2.all, flg < 2), colour=stn, geom="line",xlab="Year",ylab="CO2 (ppm)") +
# theme(axis.text.y=element_text(size = 6), legend.position="none")
#grid.arrange(p1, p2, ncol=2)
```
```{r CO2, out.width='99%', fig.align='center', fig.cap='Carbon dioxide levels from each of the measurements: (left) stations ordered from north to south, (right) overlaid.', warning=FALSE, message=FALSE, echo=FALSE, cache=FALSE}
data("scripps_data_co2")
data("scripps_stations")
# order stations by latitude
scripps_stations$station_id <- tolower(scripps_stations$`Station Code`)
scripps_data_co2 <- scripps_data_co2 %>%
left_join(scripps_stations, by = "station_id")
p1 <- scripps_data_co2 %>%
filter(flag < 2) %>%
mutate(
station_id = reorder(factor(station_id), -Latitude)
) %>%
ggplot(aes(x = date, y = co2, colour=station_id)) +
geom_line() +
xlab("Year") +
ylab("CO2 (ppm)") +
facet_wrap(~station_id, ncol=1) +
theme(axis.text.y=element_text(size = 6), legend.position="none")
p2 <- scripps_data_co2 %>%
filter(flag < 2) %>%
ggplot(aes(x = date, y = co2, colour=station_id)) +
geom_line() +
xlab("Year") +
ylab("CO2 (ppm)") +
theme(axis.text.y=element_text(size = 6), legend.position="none")
grid.arrange(p1, p2, ncol=2)
```
```{r CO2-map, out.width='99%', fig.align='center', fig.cap='Locations of stations.', warning=FALSE, message=FALSE, echo=FALSE, cache=FALSE}
world <- map_data("world")
worldmap <- ggplot(world, aes(x=long, y=lat, group=group)) +
geom_path(color="grey80", size=0.5) + xlab("") + ylab("") +
scale_y_continuous(breaks=(-2:2) * 30) +
scale_x_continuous(breaks=(-4:4) * 45) +
theme_bw() + theme(aspect.ratio=0.6)
worldmap +
# geom_point(data = scripps_stations,
# aes(x = Longitude, y = Latitude, group = 1),
# colour = "red", size = 2, alpha= 0) +
geom_text(data = scripps_data_co2 %>%
select(Latitude, Longitude, station_id)%>% unique(),
aes(x = Longitude, y = Latitude,
label = station_id, group = 1),
colour = "orange", size = 5)
# only show stations for which there is data available
# worldmap +
# geom_point(data=CO2.all.loc, aes(x=lon, y=lat, group=1), colour="red", size=2, alpha=0) +
# geom_text(data=CO2.all.loc, aes(x=lon, y=lat, label=stn, group=1), colour="orange", size=5)
```