-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgridded.temp.r
134 lines (115 loc) · 4.82 KB
/
gridded.temp.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
130
131
132
133
library(ncdf)
# Download dataset if not present
if (!file.exists("air.mon.anom.nc")) {
url <- "ftp://ftp.cdc.noaa.gov/Datasets/cru/hadcrut2/var/air.mon.anom.nc"
download.file(url, basename(url))
}
fo<-open.ncdf("air.mon.anom.nc")
fo$nvars
lon <- get.var.ncdf(fo, "lon")
lat <- get.var.ncdf(fo, "lat")
time <- get.var.ncdf(fo, "time")
temp <- get.var.ncdf(fo, "air")
temp.usa<-temp[10:24, 8:14,] # lon -130, 65; lat 25, 50
temp.usa.melt<-melt(temp.usa)
colnames(temp.usa.melt)<-c("gridx", "gridy", "time", "temp")
temp.usa.melt$lon<-lon[temp.usa.melt$gridx+9]
temp.usa.melt$lat<-lat[temp.usa.melt$gridy+7]
temp.usa.melt.gly<-glyphs(temp.usa.melt, "lon", "time", "lat", "temp")
qplot(gx, gy, data=temp.usa.melt.gly, group=gid, geom="line")
library(maps)
world <- map_data("world")
world <- getbox(world, xlim = c(-138, -55), ylim = c(22, 52))
world <- ddply(world, "group", function(df) {
if (diff(range(df$long)) < 1e-6) return(NULL)
if (diff(range(df$lat)) < 1e-6) return(NULL)
df
})
colnames(world)[1]<-"lon"
map <- list(
geom_polygon(aes(lon, lat, group = group), inherit.aes = FALSE,
data = world, legend = FALSE, fill = "grey80", colour = "grey90"),
scale_x_continuous(breaks = NA, expand = c(0.02, 0)),
scale_y_continuous(breaks = NA, expand = c(0.02, 0)),
xlab(NULL),
ylab(NULL))
qplot(gx, gy, data=temp.usa.melt.gly, group=gid, geom="line") + map
fo<-open.ncdf("air.2x2.250.mon.anom.land.nc")
fo$nvars
lon <- get.var.ncdf(fo, "lon")
lat <- get.var.ncdf(fo, "lat")
time <- get.var.ncdf(fo, "time")
temp <- get.var.ncdf(fo, "air")
# Convert lon from 1 to 360 to -180 to 180
lon<-(lon+180)%%360 - 180
length(lon)
dim(temp)
library(ggplot2)
# Ok, this now works for transforming the lon, and selecting the USA
temp.usa<-temp[114:150, 20:33,] # lon -130, -65; lat 25, 50
temp.usa.melt<-melt(temp.usa)
colnames(temp.usa.melt)<-c("gridx", "gridy", "time", "temp")
temp.usa.melt$lon<-lon[temp.usa.melt$gridx+113]
temp.usa.melt$lat<-lat[temp.usa.melt$gridy+19]
temp.usa.melt<-subset(temp.usa.melt, !(lon>-70 & lat<37))
x<-ddply(temp.usa.melt, c("lon","lat"), function(df) {
if (nrow(subset(df,is.na(temp))) == nrow(df)) NULL else df
})
temp.usa.melt.gly<-glyphs(temp.usa.melt, "lon", "time", "lat", "temp")
qplot(gx, gy, data=temp.usa.melt.gly, group=gid, geom="line") + map + coord_map() # Testing
p <-ggplot(temp.usa.melt.gly, aes(gx, gy, group = gid)) +
map
p <- p + geom_line(aes(y = lat), colour = "white", size = 1.2)
p <- p + geom_path()
p <- p + ref_boxes # This part is really slow
p <- p + theme_fullframe() #+ coord_map() # but not as slow as this part
ggsave("../images/gistemp-raw.png", width = 8, height = 4)
temp.usa.melt.gly<-glyphs(temp.usa.melt, "lon", "time", "lat", "temp", polar=T)
p <- ggplot(temp.usa.melt.gly, aes(gx, gy, group = gid)) +
map_gistemp
p <- p + add_ref_lines(temp.usa.melt.gly) + add_ref_boxes(temp.usa.melt.gly)
p <- p + geom_path()
# theme_fullframe() + coord_map()
p
ggsave("../images/gistemp-polar-raw.png", width = 8, height = 4)
# Predictions
library(lubridate)
temp.usa.melt$date <- ymd(18800101) + months(temp.usa.melt$time-1)
temp.usa.melt$year <- year(temp.usa.melt$date)
temp.usa.melt$month <- month(temp.usa.melt$date)
# Just use 1950-2010, cleaner data
temp.usa.melt.sub<-subset(temp.usa.melt, year>1950 & year<2011)
temp.usa.melt.sub<-subset(temp.usa.melt.sub, !is.na(temp))
#temp_models %<-cache% dlply(temp.usa.melt.sub, c("lat", "lon"), function(df) {
# lm(temp ~ year + factor(month), data = df)
#}) This uses old data!!!
temp_models <- dlply(temp.usa.melt.sub, c("lat", "lon"), function(df) {
lm(temp ~ year + factor(month), data = df)
})
month_grid <- expand.grid(year = 1950:2010, 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(temp.usa.melt.sub$year), month = 1)
year_preds <- ldply(temp_models, function(mod) {
year_grid$pred <- predict(mod, newdata = year_grid)
year_grid
})
temp.usa.melt.sub.gly<-glyphs(year_preds, "lon", "year", "lat", "pred")
ggplot(temp.usa.melt.sub.gly, aes(gx, gy, group = gid)) +
map_gistemp + geom_line(aes(y = lat), colour = "white", size = 1.5) +
add_ref_lines(temp.usa.melt.sub.gly) + add_ref_boxes(temp.usa.melt.sub.gly) + geom_path()
ggsave("../images/gistemp-pred.png", width = 8, height = 4)
# Coord map seems to make the plotting REALLY SLOW!
temp.usa.melt.sub.gly<-glyphs(year_preds, "lon", "year", "lat", "pred", polar=T)
ggplot(temp.usa.melt.sub.gly, aes(gx, gy, group = gid)) +
map + # Need to get the circles drawn here
geom_path() +
ref_boxes +
theme_fullframe() + coord_map()
ggsave("../images/gistemp-polar-pred.png", width = 8, height = 4)
# Testing/checking
library(plyr)
temp.mean<-ddply(temp.usa.melt, c("lon","lat"), summarise, tm=mean(temp, na.rm=T))
qplot(lon, lat, data=temp.mean, colour=tm) + map + coord_map()