-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathnasa.r
129 lines (104 loc) · 5.12 KB
/
nasa.r
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
setwd("/Users/dicook/papers/ClimateChange")
nasa <- readRDS(file.choose()) # nasadata.rds
library("maps")
library("ggplot2")
outlines <- as.data.frame(map("world", plot = FALSE,
xlim = -c(113.8, 56.2), ylim = c(-21.2, 36.2))[c("x","y")])
map <- geom_path(aes(x, y), inherit.aes = FALSE, data = outlines,
alpha = 0.2, legend = FALSE)
# Raw data
qplot(long, lat, data=nasa, geom="tile", fill=surftemp, facets=year~month) + scale_fill_gradient2(midpoint=296.9, mid="white") + map
# Fit models
# Seasonal only
temp_seas_models <- dlply(nasa, c("lat", "long"), function(df) {
# lm(ozone ~ factor(month) - 1, data = df)
lm(surftemp ~ factor(month), data = df)
})
# Including trend
temp_trend_models <- dlply(nasa, c("lat", "long"), function(df) {
lm(surftemp ~ year + factor(month), data = df)
})
locs <- dlply(nasa, c("lat", "long"))
resids <- mdply(cbind(d = locs, m = temp_seas_models), function(d, m) {
d$temp_resid <- d$surftemp - predict(m, newdata = d)
d
})
preds <- mdply(cbind(d = locs, m = temp_models), function(d, m) {
d$temp_pred <- predict(m, newdata = d)
d
})
# De-seasonalized
theme_nothing <- function (base_size = 12) {
structure(list(
panel.grid.minor=theme_blank(),
panel.grid.major=theme_blank(),
panel.background=theme_blank(),
panel.border = theme_blank(),
panel.margin = unit(0, "lines"),
axis.title.x=theme_blank(),
axis.title.y=theme_blank(),
axis.line=theme_blank(),
axis.ticks=theme_blank(),
axis.text.y = theme_blank(),
axis.text.x = theme_blank()#,
#legend.position = "none"
), class = "options")
}
qplot(long, lat, data=resids, geom="tile", fill=temp_resid, facets=year~month) + scale_fill_gradient2(midpoint=0, mid="white", high="red", low="blue") + map + theme_nothing() + coord_map()
ggsave("../images/nasa-colored-map.png", height=10, width=18)
resids <- glyphs(resids, "long", "day", "lat", "temp_resid")
qplot(gx, gy, data = resids, geom = "line", group = gid) + map
# Plot predicted values
qplot(long, lat, data=month_preds, geom="tile", fill=pred, facets=year~month) + scale_fill_gradient2(midpoint=296.9, mid="white") + map
# Plot glyphs
month_grid <- expand.grid(year = 2000, month = 1:12)
month_preds <- ldply(temp_models, function(mod) {
month_grid$pred <- predict(mod, newdata = month_grid)
month_grid
})
year_grid <- expand.grid(year = unique(nasa$year), month = 1)
year_preds <- ldply(temp_models, function(mod) {
year_grid$pred <- predict(mod, newdata = year_grid)
year_grid
})
month_preds <- glyphs(month_preds, "long", "month", "lat", "pred")
qplot(gx, gy, data = month_preds, geom = "path", group = gid) + map
year_preds <- glyphs(year_preds, "long", "year", "lat", "pred")
year_preds <- glyphs(year_preds, "long", "year", "lat", "pred", height=3)
qplot(gx, gy, data = year_preds, geom = "path", group = gid, xlab="", ylab="") + map + coord_map()
ggsave("nasa-deseas-trend.pdf", height=8, width=8)
year_resids <- glyphs(year_preds, "long", "year", "lat", "pred", height=3)
qplot(gx, gy, data = year_preds, geom = "path", group = gid, xlab="", ylab="") + map + coord_map()
ggsave("nasa-deseas-trend.pdf", height=8, width=8)
# Explain icons
nasa.temp<-data.frame(time=nasa$time, temp=nasa$surftemp)
nasa.temp<-cast(nasa.temp, formula=temp~time)
nasa.oneloc <- subset(nasa, x == 22 & y == 21, select=c("id","time","surftemp"))
nasa.twoloc <- subset(nasa, x == 7 & y == 12, select=c("id","time","surftemp"))
nasa.threeloc <- subset(nasa, x == 12 & y == 3, select=c("id","time","surftemp"))
nasa.sub<-rbind(nasa.oneloc, nasa.twoloc, nasa.threeloc)
#nasa.oneloc<-t(as.matrix(nasa.oneloc$surftemp))
#colnames(nasa.oneloc)<-paste("t", 1:72, sep="")
nasa.oneloc$gx<-cos(2*pi*(nasa.oneloc$time/72)) * ((nasa.oneloc$surftemp-293)/(302-293))
nasa.oneloc$gy<-sin(2*pi*(nasa.oneloc$time/72)) * ((nasa.oneloc$surftemp-293)/(302-293))
qplot(time, surftemp, data=nasa.oneloc, geom="path", xlab="Time", ylab="Temperature (Kelvin)")
ggsave("nasa-oneloc-ts.pdf", height=2, width=8)
qplot(gx, gy, data=nasa.oneloc, geom="path") + theme_nothing()
nasa.sub$lon<-rep(c(2,4,6), rep(72,3))
nasa.sub$lat<-rep(1,216)
nasa.sub<-glyphs(nasa.sub, "lon", "time", "lat", "surftemp", width=0.5, height=3, polar=TRUE)
nasa.sub$gid<-as.character(nasa.sub$gid)
nasa.sub$gid[nasa.sub$gid=="2.1"]<-"1.1"
nasa.sub$gid[nasa.sub$gid=="4.1"]<-"2.1"
nasa.sub$gid[nasa.sub$gid=="6.1"]<-"3.1"
nasa.sub$gid<-factor(nasa.sub$gid)
qplot(time, surftemp, data=nasa.sub, geom="path", xlab="Time", ylab="Temperature (Kelvin)", group=id, colour=gid)
ggsave("nasa-glyph-ts.pdf", height=2, width=8)
qplot(gx, gy, data=nasa.sub, geom="polygon", colour=gid, fill=I(NA)) + theme_nothing()
ggsave("nasa-glyph.pdf", height=2, width=6)
nasa.gly<-glyphs(nasa, "long", "time", "lat", "surftemp", height=3, polar=TRUE)
ggplot() + map + geom_polygon(data=nasa.gly, mapping=aes(x=gx, y=gy, group=id), colour=I("black"), fill=I(NA)) + theme_fullframe() + coord_map()
ggsave("nasa-glyph-all.pdf", height=6, width=6)
nasa.gly<-glyphs(nasa, "long", "time", "lat", "surftemp", height=1.5)
ggplot() + map + geom_path(data=nasa.gly, mapping=aes(x=gx, y=gy, group=id), colour=I("black")) + theme_fullframe() + coord_map()
ggsave("nasa-glyph-ts-all.pdf", height=6, width=6)