|
1
|
#! /bin/r
|
|
2
|
|
|
3
|
|
|
4
|
#system("ssh -X turaco")
|
|
5
|
#export LD_LIBRARY_PATH=${LD_LIBRARY_PATH:+$LD_LIBRARY_PATH:}/home/adamw/local/lib
|
|
6
|
#export KMP_DUPLICATE_LIB_OK=TRUE
|
|
7
|
#R
|
|
8
|
### Download and clean up GHCN data for CFR
|
|
9
|
#install.packages("ncdf4",lib="/home/adamw/rlib",configure.args="--with-nc-config=/home/adamw/local/bin/nc-config")
|
|
10
|
|
|
11
|
library(sp);library(rgdal);
|
|
12
|
library(reshape)
|
|
13
|
library(ncdf4)
|
|
14
|
library(geosphere)
|
|
15
|
library(rgeos)
|
|
16
|
library(multicore);library(RSQLite)
|
|
17
|
library(spBayes)
|
|
18
|
library(geoR) #for conventional kriging
|
|
19
|
library(raster)
|
|
20
|
## for review/analysis
|
|
21
|
library(coda)
|
|
22
|
library(lattice)
|
|
23
|
library(MBA)
|
|
24
|
|
|
25
|
### Load functions
|
|
26
|
#source("code/interpolationfunc.r")
|
|
27
|
|
|
28
|
X11.options(type="Xlib")
|
|
29
|
ncores=20 #number of threads to use
|
|
30
|
|
|
31
|
setwd("/home/adamw/acrobates/projects/interp")
|
|
32
|
|
|
33
|
### read in full dataset
|
|
34
|
load("data/ghcn/roi_ghcn.Rdata")
|
|
35
|
load("stroi.Rdata")
|
|
36
|
source("code/GHCN_functions.r")
|
|
37
|
|
|
38
|
## rescale units to C and mm
|
|
39
|
d$value=d$value/10.0
|
|
40
|
|
|
41
|
## add flag for Benoit's training stations to identify training/validation
|
|
42
|
stt=readOGR("data/ghcn/tempdata/ghcn_1507_s.shp","ghcn_1507_s")$STAT_ID
|
|
43
|
d$fit=d$id%in%stt
|
|
44
|
|
|
45
|
### read in and subset to Benoit's dates
|
|
46
|
dates=as.Date(as.character(
|
|
47
|
read.table("data/ghcn/tempdata/dates_interpolation_03052012.txt")$V1),"%Y%m%d")
|
|
48
|
|
|
49
|
##########################################################
|
|
50
|
#### load region of interest
|
|
51
|
roi=readOGR("data/boundaries/statesp020.shp","statesp020")
|
|
52
|
proj4string(roi)=CRS("+proj=longlat +datum=WGS84")
|
|
53
|
roi=roi[roi$STATE=="Oregon",]
|
|
54
|
roi=gUnionCascaded(roi) #dissolve any subparts of roi
|
|
55
|
|
|
56
|
## distance from ROI to buffer for border stations
|
|
57
|
dis=100 #km
|
|
58
|
## buffer roi (transform to azimuthal equidistant with centroid of roi, add 'dis' buffer, then transform back)
|
|
59
|
roic=centroid(roi)
|
|
60
|
|
|
61
|
roib=spTransform(
|
|
62
|
gBuffer(spTransform(roi,
|
|
63
|
CRS(paste("+proj=aeqd +lat_0=",roic[2]," +lon_0=",
|
|
64
|
roic[1]," +ellps=WGS84 +datum=WGS84 +units=m +no_defs ",sep=""))),
|
|
65
|
width=dis*1000),
|
|
66
|
CRS("+proj=longlat +datum=WGS84"))
|
|
67
|
|
|
68
|
### generate knots
|
|
69
|
knots=makegrid(roib,cellsize=0.5)
|
|
70
|
coordinates(knots)=c("x1","x2")
|
|
71
|
proj4string(knots)=proj4string(roi)
|
|
72
|
knots=knots[!is.na(overlay(knots,roi))]
|
|
73
|
|
|
74
|
|
|
75
|
###########################################################3
|
|
76
|
#### interpolation function
|
|
77
|
|
|
78
|
interp<-function(parms,niter=100,istart=60,ithin=1){
|
|
79
|
## load variables from parameter list
|
|
80
|
idate=as.Date(parms$date)
|
|
81
|
ivar=parms$variable
|
|
82
|
type=parms$type
|
|
83
|
model=parms$model
|
|
84
|
|
|
85
|
## subset single date from data
|
|
86
|
td=d[d$date==idate,]
|
|
87
|
td=data.frame(cast(td,id+fit~variable,value="value"))
|
|
88
|
td$elev=stroi$elev[match(td$id,stroi$id)]
|
|
89
|
td[c("lon","lat")]=stroi@data[match(td$id,stroi$id),c("lon","lat")]
|
|
90
|
td$intercept=1
|
|
91
|
|
|
92
|
## drop any points without coordinates
|
|
93
|
td=td[!is.na(td$lat)&!is.na(td$lon),]
|
|
94
|
|
|
95
|
###########################
|
|
96
|
## apply climate correction
|
|
97
|
if(type=="anom"){
|
|
98
|
cppt=brick("data/prism/prism_climate.nc",varname="ppt")
|
|
99
|
cppt=subset(cppt,subset=as.numeric(format(idate,"%m")))
|
|
100
|
cppt[cppt==-99.99]=NA
|
|
101
|
td$cppt=as.numeric(extract(cppt,td[,c("lon","lat")],method="bilinear"))
|
|
102
|
## Tmax
|
|
103
|
ctmax=brick("data/prism/prism_climate.nc",varname="tmax")
|
|
104
|
ctmax=subset(ctmax,subset=as.numeric(format(idate,"%m")))
|
|
105
|
ctmax[ctmax==-99.99]=NA
|
|
106
|
td$ctmax=as.numeric(extract(ctmax,td[,c("lon","lat")],method="bilinear"))
|
|
107
|
## Tmin
|
|
108
|
ctmin=brick("data/prism/prism_climate.nc",varname="tmin")
|
|
109
|
ctmin=subset(ctmin,subset=as.numeric(format(idate,"%m")))
|
|
110
|
ctmin[ctmin==-99.99]=NA
|
|
111
|
td$ctmin=as.numeric(extract(ctmin,td[,c("lon","lat")],method="bilinear"))
|
|
112
|
## Calculate anomalies
|
|
113
|
td$atmax=td$tmax-td$ctmax
|
|
114
|
td$atmin=td$tmin-td$ctmin
|
|
115
|
## Scale for precipitation (add 1 to climate to avoid dividing by 0, will be subtracted off later)
|
|
116
|
## Add 1 to final value to avoid transformation troubles (geoR will only transform if >0
|
|
117
|
td$appt=(td$ppt/(ifelse(td$cppt<0,0,td$cppt)+1))+1
|
|
118
|
}
|
|
119
|
|
|
120
|
### Add region flag to data
|
|
121
|
td$region=ifelse(!is.na(overlay(SpatialPoints(td[,c("lon","lat")]),roi)),"roi","roib")
|
|
122
|
td=td[!is.na(overlay(SpatialPoints(td[,c("lon","lat")]),roib)),]
|
|
123
|
|
|
124
|
#### Extract data for validation
|
|
125
|
fd=td[td$fit&td$region=="roi",] #fitting dataset
|
|
126
|
pd=td[!td$fit&td$region=="roi",] #training dataset
|
|
127
|
od=td[!td$fit&td$region!="roi",] #stations outside roi
|
|
128
|
|
|
129
|
|
|
130
|
### drop missing values
|
|
131
|
fd=fd[!is.na(fd[,as.character(ivar)]),] #&!is.na(fd$tmin)
|
|
132
|
|
|
133
|
### define the formula
|
|
134
|
if(type=="raw"&ivar=="tmax") f1=formula(tmax~intercept)
|
|
135
|
if(type=="anom"&ivar=="tmax"&model=="intercept") f1=formula(atmax~intercept)
|
|
136
|
if(type=="anom"&ivar=="tmax"&model=="full") f1=formula(atmax~lon+lat+elev)
|
|
137
|
|
|
138
|
### Run the model
|
|
139
|
t1=Sys.time()
|
|
140
|
m1=spLM(f1,data=fd,coords=as.matrix(fd[,c("lon","lat")]),
|
|
141
|
starting=list("phi"=2.5,"sigma.sq"=4, "tau.sq"=1),
|
|
142
|
sp.tuning=list("phi"=0.8, "sigma.sq"=0.1, "tau.sq"=0.2),
|
|
143
|
priors=list("phi.Unif"=c(0.1, 4), "sigma.sq.IG"=c(2, 1),
|
|
144
|
"tau.sq.IG"=c(2, 1)),
|
|
145
|
cov.model="exponential",
|
|
146
|
knots=coordinates(knots),
|
|
147
|
n.samples=niter, verbose=TRUE, n.report=1000)
|
|
148
|
t2=Sys.time()
|
|
149
|
|
|
150
|
|
|
151
|
#m1s=mcmc(m1$p.samples)
|
|
152
|
#summary(m1s)
|
|
153
|
#xyplot(m1s,scales=list(y=list(rot=0)),main="Posterior Samples of Model Parameters")
|
|
154
|
|
|
155
|
### Run the model
|
|
156
|
#q=3
|
|
157
|
#nltr=q*(q-1)/2+q
|
|
158
|
|
|
159
|
#m1=spMvLM(list(atmin~1,atmax~1),data=fd,coords=as.matrix(fd[,c("lon","lat")]),
|
|
160
|
# starting=list("phi"=0.6,"sigma.sq"=1, "tau.sq"=1),
|
|
161
|
# sp.tuning=list("phi"=0.01, "sigma.sq"=0.05, "tau.sq"=0.05),
|
|
162
|
# priors=list("phi.Unif"=c(0.1, 3), "sigma.sq.IG"=c(2, 1),
|
|
163
|
# "tau.sq.IG"=c(2, 1)),
|
|
164
|
# cov.model="exponential",
|
|
165
|
# knots=coordinates(knots),
|
|
166
|
# n.samples=100, verbose=TRUE, n.report=100)
|
|
167
|
|
|
168
|
#save to play with later
|
|
169
|
#save(m1,file="m1.Rdata")
|
|
170
|
|
|
171
|
|
|
172
|
#load("m1.Rdata")
|
|
173
|
|
|
174
|
##############################
|
|
175
|
#### make predictions
|
|
176
|
|
|
177
|
### Prediction grid
|
|
178
|
# pgrid=data.frame(coordinates(ctmax))
|
|
179
|
# ## crop to bbox
|
|
180
|
# pgrid=pgrid[pgrid$x>=bbox(roi)[1,1]&pgrid$x<=bbox(roi)[1,2]
|
|
181
|
# &pgrid$y<=bbox(roi)[2,2]&pgrid$y>=bbox(roi)[2,1],]
|
|
182
|
# ## Crop to ROI polygon?
|
|
183
|
# #pgrid=pgrid[!is.na(overlay(SpatialPoints(pgrid),roi)),]
|
|
184
|
# rownames(pgrid)=1:nrow(pgrid)
|
|
185
|
# ncluster=200
|
|
186
|
# pgrid$cluster=cut(as.numeric(rownames(pgrid)),ncluster,labels=1:ncluster)
|
|
187
|
# ## predict for the full grid
|
|
188
|
# m1pg=lapply(1:ncluster,function(i){
|
|
189
|
# print(paste("Starting cluster:",i))
|
|
190
|
# m1pg=spPredict(m1,pred.coords=pgrid[pgrid$cluster==i,c("x","y")],
|
|
191
|
# pred.covars=cbind(intercept=rep(1,nrow(pgrid[pgrid$cluster==i,]))),start=istart,thin=ithin)
|
|
192
|
# ## Generate summaries of y.pred
|
|
193
|
# ty=t(apply(m1pg$y.pred,1,function(i) c(mean=mean(i),sd=sd(i),quantile(i,c(0.025,0.5,0.975)))))
|
|
194
|
# colnames(ty)=paste("y.",c("mean","sd","Q2.5","Q50","Q97.5"),sep="")
|
|
195
|
# ## Generate summaries of w.pred
|
|
196
|
# tw=t(apply(m1pg$w.pred,1,function(i) c(mean=mean(i),sd=sd(i),quantile(i,c(0.025,0.5,0.975)))))
|
|
197
|
# colnames(tw)=paste("w.",c("mean","sd","Q2.5","Q50","Q97.5"),sep="")
|
|
198
|
# return(cbind(ty,tw))
|
|
199
|
# })
|
|
200
|
# m1pgs=do.call(rbind.data.frame,m1pg)
|
|
201
|
|
|
202
|
|
|
203
|
########### Predict only to validation locations
|
|
204
|
m1p=spPredict(m1,pred.coords=pd[,c("lon","lat")],pred.covars=model.matrix(lm(f1,data=pd)),start=istart,thin=ithin)#,
|
|
205
|
|
|
206
|
m1p.y=mcmc(t(m1p$y.pred),start=istart,thin=ithin)
|
|
207
|
m1p.w=mcmc(t(m1$sp.effects.knots),start=istart,thin=ithin)
|
|
208
|
m1p.ys=t(apply(m1p.y,2,quantile,c(0.025,0.5,0.975),na.rm=T))
|
|
209
|
m1p.ws=t(apply(m1p.w,2,quantile,c(0.025,0.5,0.975),na.rm=T))
|
|
210
|
|
|
211
|
### recover original scale if necessary
|
|
212
|
resp=as.character(attr(terms.formula(f1),"variables"))[2]
|
|
213
|
## if modeling anomalies, add the climate back on
|
|
214
|
if(type=="anom") {
|
|
215
|
pd2=cbind.data.frame(pd,aQ2.5=m1p.ys[,"2.5%"],aQ50=m1p.ys[,"50%"],aQ97.5=m1p.ys[,"97.5%"])
|
|
216
|
# FIXME! ctmax below needs to be made general for other variables!
|
|
217
|
pd2$Q2.5=pd2$ctmax+pd2$aQ2.5
|
|
218
|
pd2$Q50=pd2$ctmax+pd2$aQ50
|
|
219
|
pd2$Q97.5=pd2$ctmax+pd2$aQ97.5
|
|
220
|
}
|
|
221
|
## if modeling raw values, don't add the climate back on
|
|
222
|
if(type!="anom") {
|
|
223
|
pd2=cbind.data.frame(pd,Q2.5=m1p.ys[,"2.5%"],Q50=m1p.ys[,"50%"],Q97.5=m1p.ys[,"97.5%"])
|
|
224
|
}
|
|
225
|
|
|
226
|
|
|
227
|
### save model output
|
|
228
|
save(m1,m1p,pd2,file=paste("output/modeloutput_",
|
|
229
|
gsub("-","",idate),"_",ivar,"_",type,"_",model,".Rdata",sep=""))
|
|
230
|
|
|
231
|
|
|
232
|
|
|
233
|
####################################
|
|
234
|
#### Accuracy Assement
|
|
235
|
m1.dic=spDiag(m1,start=istart,thin=ithin);m1.dic
|
|
236
|
m1.accuracy=accuracy(pd2[,as.character(ivar)],pd2$Q50,ppt=ivar=="ppt")$summary
|
|
237
|
|
|
238
|
#### Write summary file
|
|
239
|
m1.summary=cbind.data.frame(date=idate,model=model,var=ivar,type=type,
|
|
240
|
formula=paste(terms(f1),collapse=" "),
|
|
241
|
t(m1.dic$DIC),t(m1.accuracy),n.iter=niter,runtime.hours=as.numeric(difftime(t1,t2,units="hours")))
|
|
242
|
write.csv(m1.summary,file=paste("output/modeloutput_summary_",gsub("-","",idate),"_",ivar,"_",type,".csv",sep=""))
|
|
243
|
|
|
244
|
##################################################
|
|
245
|
### some summary plots
|
|
246
|
pdf(width=11,height=8.5,file=paste("output/modeloutput_",gsub("-","",idate),"_",ivar,"_",type,".pdf",sep=""))
|
|
247
|
|
|
248
|
## Show knots and stations
|
|
249
|
plot(roib);axis(1);axis(2)
|
|
250
|
plot(roi,add=T)
|
|
251
|
#points(knots,pch=3,cex=.5)
|
|
252
|
points(unique(fd[,c('lon','lat')]),cex=.8,pch=16,col="red")
|
|
253
|
points(unique(pd[,c('lon','lat')]),cex=.5,pch=16,col="blue")
|
|
254
|
points(unique(od[,c('lon','lat')]),cex=.5,pch=16,col="black")
|
|
255
|
legend(-118,42,legend=c("Fitting","Validation","Border Stations"),pch=c(16,16,16),col=c("red","blue","black"),bg="white")
|
|
256
|
|
|
257
|
### Chain behavior
|
|
258
|
m1s=mcmc(m1$p.samples)
|
|
259
|
print(xyplot(m1s,scales=list(y=list(rot=0)),main="All Posterior Samples of Model Parameters"))
|
|
260
|
m1s=mcmc(m1$p.samples,start=istart,thin=ithin)
|
|
261
|
print(xyplot(m1s,scales=list(y=list(rot=0)),main="Post-burnin, thinned, Posterior Samples of Model Parameters"))
|
|
262
|
|
|
263
|
|
|
264
|
### Pred vs. obs plots
|
|
265
|
plot(pd2[,as.character(ivar)],pd2$Q50,pch=16,cex=.8,ylim=quantile(c(pd2$Q2.5,pd2$Q97.5),c(0.01,0.995)),xlim=quantile(c(pd2$Q2.5,pd2$Q97.5),c(0.01,0.995)))
|
|
266
|
segments(pd2[,as.character(ivar)],pd2$Q2.5,pd2[,as.character(ivar)],pd2$Q97.5,col="grey")
|
|
267
|
points(pd2[,as.character(ivar)],pd2$Q50,pch=16,cex=.8);abline(0,1,col="red") #xlim=c(0,10)
|
|
268
|
|
|
269
|
### Maps
|
|
270
|
par(mfrow=c(1,3))
|
|
271
|
## observed
|
|
272
|
obs.surf <- mba.surf(cbind(fd[,c("lon","lat")],fd$tmax), no.X=100, no.Y=100, extend=TRUE)$xyz.est
|
|
273
|
image(obs.surf, xaxs = "r", yaxs = "r", main="Observed response",asp=1)
|
|
274
|
points(fd[,c("lon","lat")])
|
|
275
|
plot(roi,add=T)
|
|
276
|
#contour(obs.surf, add=T)
|
|
277
|
plot(roi,add=T,border="darkgreen",lwd=3)
|
|
278
|
## spatial effects
|
|
279
|
w.pred.surf <-
|
|
280
|
mba.surf(cbind(coordinates(knots), m1p.ws[,3]), no.X=100, no.Y=100, extend=TRUE)$xyz.est
|
|
281
|
image(w.pred.surf, xaxs = "r", yaxs = "r", main="Spatial Effects",asp=1)
|
|
282
|
#points(pd[,c("lon","lat")], pch=1, cex=1)
|
|
283
|
points(knots, pch=3, cex=1)
|
|
284
|
#contour(w.pred.surf, add=T)
|
|
285
|
plot(roi,add=T,border="darkgreen",lwd=3)
|
|
286
|
legend(1.5,2.5, legend=c("Obs.", "Knots", "Pred."),
|
|
287
|
pch=c(1,3,19), bg="white")
|
|
288
|
## fitted
|
|
289
|
y.pred.surf <-
|
|
290
|
mba.surf(cbind(pd2[,c("lon","lat")], pd2$Q50), no.X=100, no.Y=100, extend=TRUE)$xyz.est
|
|
291
|
image(y.pred.surf, xaxs = "r", yaxs = "r", main="Predicted response",asp=1)
|
|
292
|
points(pd[,c("lon","lat")], pch=1, cex=1)
|
|
293
|
#points(knots, pch=3, cex=1)
|
|
294
|
#contour(y.pred.surf, add=T)
|
|
295
|
plot(roi,add=T,border="darkgreen",lwd=3)
|
|
296
|
legend(1.5,2.5, legend=c("Obs.", "Knots", "Pred."),
|
|
297
|
pch=c(1,3,19), bg="white")
|
|
298
|
|
|
299
|
dev.off()
|
|
300
|
|
|
301
|
### Return objects
|
|
302
|
return(m1.summary)
|
|
303
|
}
|
|
304
|
|
|
305
|
|
|
306
|
|
|
307
|
###########################################################
|
|
308
|
### define models to run
|
|
309
|
ms=expand.grid(variable="tmax",type=c("anom","raw"),date=dates,model=c("intercept","full"))
|
|
310
|
## drop some
|
|
311
|
ms[!ms$type=="raw"&ms$model=="intercept",]
|
|
312
|
|
|
313
|
|
|
314
|
msl=apply(ms,1,as.list)
|
|
315
|
|
|
316
|
|
|
317
|
### test it out
|
|
318
|
parms=msl[[9]]
|
|
319
|
interp(parms,niter=5000,istart=2000,ithin=3)
|
|
320
|
|
|
321
|
#mresults=lapply(msl[1:3],interp)
|
|
322
|
|
|
323
|
### run it!
|
|
324
|
mresults =mclapply(msl,interp,niter=1000,istart=500,ithin=1,mc.cores=1)
|
|
325
|
|
|
326
|
### drop any errors
|
|
327
|
mresults=mresults[!sapply(mresults,function(x) grepl("Error",x[[1]]))]
|
|
328
|
|
|
329
|
res=do.call(rbind.data.frame,mresults)
|
|
330
|
|
|
331
|
save(mresults,res,file="output/mresults.Rdata")
|
|
332
|
|
|
333
|
q("no")
|
|
334
|
|
|
335
|
|
|
336
|
|
|
337
|
|
|
338
|
|
