#SOME UTILITIES ##JONES GRIDCELL LOOKUP jones<-function(lat,long) { i<-18-floor(lat/5); j<- 36+ceiling(long/5); jones<-72*(i-1) + j jones } ##jones.inv jones.inv<-function(k) { j<- k%%72 temp<-k-j i<-(temp/72)+1 long<-(j-36)*5-2.5 lat<-(18-i)*5+2.5 jones.inv<-c(lat,long) jones.inv } ts.annavg=function(X,M=6) { if(class(X)=="ts") {start0<-floor(tsp(X)[1]);X<-ts.array(X)} else if( ((class(X)=="matrix")|(class(X)=="data.frame") ) &(ncol(X)==12)) start0<-as.numeric(row.names(X)[1]) else {X start0<-X[1,1]; X<-X[,2:13]} count<-apply(!is.na(X),1,sum);temp<-(count0) x<-c(x, rep(NA,12-m)) years<-length(x)/12 x<-array(x,dim=c(12,years)) annavg<-apply(x,2,mean,na.rm=TRUE) return(annavg) } ##strip #FUNCTION TO LISTFILES AND REMOVE SUFFIXES strip<-function(h,n) { h<-unlist(strsplit(h, "\\.")) h<-array(h,dim=c(n,length(h)/n)) h<-t(h) strip<-h[,1] strip } ##split split<-function(h,char) { y<-unlist(strsplit(h,paste("\\",char,sep=""))) y<-array(y,dim=c(2,length(y)/2)) split<-t(y) split<-split[,1] split } ##smooth and pad filter.combine.pad=function(x,a,M=NA,method="norm") { temp<-!is.na(x) w<-x[temp] N<-length(w) if(is.na(M)) M<-trunc (length(a)/2)+1 if(method=="norm") w<-c (rep(mean(w[1:M]),M),w,rep(mean(w[(N-M+1):N]),M) ) #if(method=="reflect") w<-c ( w[M:1],w,w[M:(N-M+1):N]) if(method=="mannomatic") {w<-c( c (w[1]-(w[(M+1):2]-w[1])),w, c(w[N]- (w[(N-1):(N-M)]-w[N]) ) ) } if(method=="mann2") w<-c (rep(w[1],M),w,rep(w[N],M) ) if(method=="linear"){ fm1=lm (w[1:M]~I(1:M)); fm2= lm (w[N+1-(1:M)]~I(N+1-(1:M))); w<-c (fm1$coef[1]+fm1$coef[2]*( -(M-2):0) ,w,fm2$coef[1]+fm2$coef[2]*( N+M-((M-1):1)) ) } y<-filter(w,a,method="convolution",sides=2) y<-y[!is.na(y)] z<- rep(NA,length(x)) z[temp]<-y if(class(x)=="ts") z<-ts(z,start=tsp(x)[1],freq=tsp(x)[3]) filter.combine.pad<-cbind(x,z) #filter.combine<-y[,2] filter.combine.pad } ##gaussian.filter.weights gaussian.filter.weights<-function(N,year) { #N number of points; year - is number of years sigma<-year/6 i<-( (-(N-1)/2):((N-1)/2) ) /sigma w<- ((2*pi)^-0.5) * exp(-0.5* i^2) gaussian.filter.weights<-w/sum(w) gaussian.filter.weights } #GAUSSIAN FILTER #see http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf #The appropriate Gaussian distribution therefore has standard deviation # s = ? / 6, # where ? is the desired wavelength at which the amplitude of frequency response is 0.5. The # filter weights are obtained by sampling the pdf of the standard normal distribution at # values i/s where i={-(N-1)/s, ..., 0, .. (N-1)/s } # MEKo has values to verify for case of N=21, year=10 #this is verified ##truncated.gauss.weights truncated.gauss.weights<-function(year) { a<-gaussian.filter.weights(2*year+1,year) temp<-(a>0.05*max(a)) a<-a[temp] a<-a/sum(a) truncated.gauss.weights<-a truncated.gauss.weights } maketable=function(X,align=-999,widths=-999){ N=nrow(X);M=ncol(X); out= "d:/temp/temp.txt" if (align[1]== -999) right=rep(TRUE,M) else right=align code=rep("",M);code[right]="

" if(widths[1]== -999) pct= floor(100/M) else pct=widths write.table("",file=out,col.names=FALSE,row.names=FALSE,quote=FALSE) for (j in 1:N){ write.table("",append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE) for(i in 1:M){ write.table( paste("",sep=""),append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE) } write.table("",append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE) } write.table("
",code[i],X[j,i],"
", append=TRUE,file=out,row.names=FALSE,col.names=FALSE,quote=FALSE) } #TO EXTEND SERIES BY PERSISTENCE extend.persist<-function(tree) { extend.persist<-tree for (j in 1:ncol(tree) ) { test<-is.na(tree[,j]) end1<-max ( c(1:nrow(tree)) [!test]) test2<-( c(1:nrow(tree))>end1) & test extend.persist[test2,j]<-tree[end1,j] } extend.persist } #4. FUNCTION TO RETURN VERIFICATION PERIOD STATISTICS #calibration and verificaiton are both vectors of length 2 with start and end period #estimator and observed are time series which may have different start and end periods use0="pairwise.complete.obs" verification.stats<-function(estimator,observed=sparse,calibration=c(1902,1980),verification=c(1854,1901)) { combine<-ts.union(estimator,observed) index.cal<- (calibration[1]-tsp(combine)[1]+1):(calibration[2]-tsp(combine)[1]+1) index.ver<-(verification[1]-tsp(combine)[1]+1):(verification[2]-tsp(combine)[1]+1) xmean<-mean ( combine[index.cal,"observed"],na.rm=TRUE ) mx<-mean( combine[index.ver,"observed"],na.rm=TRUE ) # test.ver<-cov ( combine[index.ver,],use=use0) test<-cov(combine[index.cal,],use=use0) R2.cal<-(test[1,2]*test[2,1])/(test[1,1]*test[2,2]) R2.ver<-(test.ver[1,2]*test.ver[2,1])/(test.ver[1,1]*test.ver[2,2]) RE.cal<-1-sum( (combine[index.cal,"estimator"]- combine[index.cal,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE) RE.ver<-1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (xmean- combine[index.ver,"observed"]) ^2,na.rm=TRUE) RE.risk<- - sum( combine[index.ver,"estimator"]^2,na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE) RE.bias<- 2*length(index.ver) * mean(combine[index.ver,"estimator"],na.rm=TRUE)*mean(combine[index.ver,"observed"],na.rm=TRUE)/sum( (xmean- combine[index.cal,"observed"]) ^2,na.rm=TRUE) RE.covar<- 2* (length(index.ver)-1)* cov(combine[index.ver,],use="pairwise.complete.obs")[1,2]/sum( (xmean- combine[index.ver,"observed"]) ^2,na.rm=TRUE) CE<- 1-sum( (combine[index.ver,"estimator"]- combine[index.ver,"observed"]) ^2,na.rm=TRUE)/sum( (mx- combine[index.ver,"observed"]) ^2,na.rm=TRUE) test<-sign(diff(combine[index.ver,],lag=1)) temp<- ( test[,1]*test[,2] ==1)&!is.na(test[,1])&!is.na(test[,2]) sign.test<-sum(temp)/length(index.ver) test1<-scale(combine[index.ver,],scale=FALSE) test1<-test1[!is.na(test1[,1])&!is.na(test1[,2]),] temp<-sign(test1) temp<-(temp[,1]*temp[,2]==1) test2<-test1[,1]*test1[,2] if(sum(!temp) >0) prod.test<- ( mean(test2[temp])+mean(test2[!temp]))/ sqrt( var(test2[temp])/sum(temp) + var(test2[!temp])/sum(!temp) ) else prod.test<-NA verification.stats<-data.frame(RE.cal,RE.ver,R2.cal,R2.ver,CE)#,sign.test,prod.test,RE.risk,RE.bias,RE.covar) names(verification.stats)<-c( "RE.cal","RE.ver","R2.cal","R2.ver","CE")#,"sign.test","prod.test","RE.risk","RE.bias","RE.covar") verification.stats }