Project

General

Profile

Download (11.5 KB) Statistics
| Branch: | Revision:
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

    
(15-15/15)