################## ## COLLATION FUNCTIONS ######################## #get.crn #get.rwl #get.rcs #read.rwm # for climategate #make.rwl yam,al variation ############## #CHRONOLOGIES ################## get.crn=function(dset) { if(dset=="yamal") { loc="http://www.cru.uea.ac.uk/cru/people/melvin/PhilTrans2008/Column.prn" briffa=read.table(loc,skip=1,fill=TRUE) name0=scan(loc,n=8,what="") name0=outer(name0,c("","count"),function(x,y) paste(x,y,sep=".") ) n=nchar(name0[,1]) name0[,1]=substr(name0[,1],1,n-1) names(briffa)=c("year", c(t(name0) ) ) briffa[briffa== -9999]=NA briffa=briffa[,c("year","Yamal.RCS","Yamal.RCS.count")] briffa=briffa[!is.na(briffa[,2]),] briffa=window(ts(briffa[,2:3],start=briffa[1,1]),start=-202) yamal.crn=briffa[,1]/1000 return(yamal.crn) } if(dset=="urals") { download.file("http://www.climateaudit.info/data/esper/esper.archived.tab","temp.rwl",mode="wb"); load("temp.rwl") urals.crn=esper.archived[,"Pol"] year= c(time(urals.crn))[!is.na(urals.crn)] urals.crn= window(urals.crn,start=min(year),end=max(year) )/1000 rm(esper.archived) return(urals.crn) } if (dset=="hantemirov2002") { loc="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/reconstructions/asia/russia/yamal_2002.txt" hant=read.table(loc,skip=57,fill=TRUE,header=TRUE) # Year Recon Chron Samples #dim(hant) #4064 4 hant=window( ts(hant[,2:4],start=hant[1,1]),start= -202 ) #minimum in measurement data return(hant) } if(dset=="briffa09") { #crn is compound of 13 series: # [1] "YAMALAD" "YAMAL_SF" "YAMAL_ALL" "YAMAL_KHAD" "YAMAL_POR" "YAMAL_JAH" "YAMAL_YAD" #[8] "KHAD" "POR" "JAH" "YAD" "LIVE" "YAMALAD_KHAD" index=c("RCM", "RCS", "STD","CSM") loc="http://www.cru.uea.ac.uk/cru/people/briffa/yamal2009/data/Yamal_CRNS_col.txt" name0=scan(loc,n=30,what="") crn=read.table(loc,skip=30,nrow=1996 +202+1,fill=TRUE,colClasses="numeric") name1=scan(loc,n=22,skip=1996+202+1+30,what="") test=read.table(loc,skip=1996+202+1+30+22,fill=TRUE,colClasses="numeric") # 2199 45 crn=cbind(crn,test[,2:ncol(test)]) # name0=c(name0,name1) #52 year=crn[,1]; crn=crn[,2:ncol(crn)] #104 name1=array(c(unlist(strsplit(name0[grep("STD",name0)],"\\.") )),dim=c(2,13) ) name.briffa=name1[1,] count=crn[,seq(2,104,2)] #52 crn=crn[,seq(1,103,2)];dim(crn) #2199 52 dimnames(crn)[[2]]=name0 temp= crn==-9999;crn[temp]=NA X=ts(crn[,grep("STD",name0)],start=year[1])/1000 dimnames(X)[[2]]=name.briffa return(X) } #briffa09 if(dset=="briffa09.CSM") { ##smoothed version of crn #crn is compound of 13 series: # [1] "YAMALAD" "YAMAL_SF" "YAMAL_ALL" "YAMAL_KHAD" "YAMAL_POR" "YAMAL_JAH" "YAMAL_YAD" #[8] "KHAD" "POR" "JAH" "YAD" "LIVE" "YAMALAD_KHAD" index=c("RCM", "RCS", "STD","CSM") loc="http://www.cru.uea.ac.uk/cru/people/briffa/yamal2009/data/Yamal_CRNS_col.txt" name0=scan(loc,n=30,what="") crn=read.table(loc,skip=30,nrow=1996 +202+1,fill=TRUE,colClasses="numeric") name1=scan(loc,n=22,skip=1996+202+1+30,what="") test=read.table(loc,skip=1996+202+1+30+22,fill=TRUE,colClasses="numeric") # 2199 45 crn=cbind(crn,test[,2:ncol(test)]) # name0=c(name0,name1) #52 year=crn[,1]; crn=crn[,2:ncol(crn)] #104 name1=array(c(unlist(strsplit(name0[grep("CSM",name0)],"\\.") )),dim=c(2,13) ) name.briffa=name1[1,] count=crn[,seq(2,104,2)] #52 crn=crn[,seq(1,103,2)];dim(crn) #2199 52 dimnames(crn)[[2]]=name0 temp= crn==-9999;crn[temp]=NA X=ts(crn[,grep("CSM",name0)],start=year[1])/1000 dimnames(X)[[2]]=name.briffa return(X) } #briffa09.CSM if(dset=="briffa00") { loc<-"http://www.cru.uea.ac.uk/cru/people/briffa/qsr1999/qsrfig1.csv" briffa<-read.csv(loc,fill=TRUE,skip=10,header=TRUE)#dim 1999x30 briffa= ts(briffa[,"f1r"],start=briffa[1,1]) briffa=window(briffa,end=max (time(briffa)[!is.na(briffa)]) ) return(briffa) } #dset briffa00 }#get.crn get.rcs=function(dset="briffa09") { if(dset=="briffa09"){ index=c("RCM", "RCS", "STD","CSM") loc="http://www.cru.uea.ac.uk/cru/people/briffa/yamal2009/data/Yamal_CRNS_col.txt" name0=scan(loc,n=30,what="") crn=read.table(loc,skip=30,nrow=1996 +202+1,fill=TRUE,colClasses="numeric") name1=scan(loc,n=22,skip=1996+202+1+30,what="") test=read.table(loc,skip=1996+202+1+30+22,fill=TRUE,colClasses="numeric") # 2199 45 crn=cbind(crn,test[,2:ncol(test)]) # name0=c(name0,name1) #52 year=crn[,1]; crn=crn[,2:ncol(crn)] #104 name1=array(c(unlist(strsplit(name0[grep("STD",name0)],"\\.") )),dim=c(2,13) ) name.briffa=name1[1,] count=crn[,seq(2,104,2)] #52 crn=crn[,seq(1,103,2)];dim(crn) #2199 52 dimnames(crn)[[2]]=name0 dimnames(count)[[2]]=name0 temp= crn==-9999;crn[temp]=NA test=list() #RCM X=crn[,grep("RCM",name0)] Y=count[,grep("RCM",name0)] dimnames(Y)[[2]]=dimnames(X)[[2]]=name.briffa x=apply(!is.na(X),1,sum) temp= x==0 time0= year[!temp] ####starts in year 2 ???? #### assumed to start in year 1 test[["RCM"]]= X[!temp,] #416 13 test[["RCM.count"]]= Y[!temp,] #416 13 #RCS X=crn[,grep("RCS",name0)] Y=count[,grep("RCS",name0)] dimnames(Y)[[2]]=dimnames(X)[[2]]=name.briffa x=apply(!is.na(X),1,sum) temp= x==0 time0= year[!temp] ####starts in year 2 ???? #### assumed to start in year 1 test[["RCS"]]= X[!temp,] #416 13 test[["RCS.count"]]= Y[!temp,] #416 13 return(test) } #dset } test=get.rcs(dset="briffa09") get.rwl= function(dset,directory="",verbose=NULL) { if(is.null(verbose)) { if (directory=="") loc= paste("http://www.climateaudit.info/data/tree/rwl/",dset,".rwl.tab",sep="") else { loc= paste("http://www.climateaudit.info/data/tree/rwl/",directory,"/",dset,".rwl.tab",sep="") } download.file(loc,"temp.dat",mode="wb") load("temp.dat") return (tree) } if(verbose=="verbose") { if(dset=="yamal") { #w" loc="http://www.cru.uea.ac.uk/cru/people/melvin/PhilTrans2008/YamalADring.raw" download.file(loc,"temp.dat") tree=make.rwl_new("temp.dat") tree$id=factor(tree$id) #252 tree=agef(tree) # save(tree,file="d:/climate/data/yamal/yamal_cru.rwl.tab") tree$rw=tree$rw/10 #to make consistent units return(tree) } #yamal } #verbose } #get.rwl make.info= function(tree) { info=data.frame(id=levels(tree$id)) info$id=as.character(info$id) info$start=c(unlist(tapply(tree$year,tree$id,min))) info$end=c(unlist(tapply(tree$year,tree$id,max))) info$max=info$end-info$start+1 n=nchar(info$id) info$tree=substr(info$id,1,n-1) ;# info$tree=factor(info$tree) info$core=substr(info$id,n,n) return(info) } #last digit 0 can be a check get.hadcru=function(dset= "HadCRUT3") { if (dset=="CRUTEM2") loc="http://www.cru.uea.ac.uk/ftpdata/crutem2.nc" else { loc =paste("http://hadobs.metoffice.com/",casefold(dset),"/data/",dset,".nc",sep="") } download.file(loc,"temp.dat",mode="wb") #18.9 MB v<-open.ncdf("temp.dat") instr <- get.var.ncdf( v, v$var[[1]]) # 1850 2006 if (dset=="CRUTEM2") instr=instr[,36:1,] dim(instr)# [1] 72 36 1899 #this is organized in 72 longitudes from -177.5 to 177.5 and 36 latitudes from -87.5 to 87.5 #VRutem2 is opposite return(instr) } get.gridcell=function(lat,long,Instr=instr) { j= max( (1:36) [ lat> seq(-90,85,5) ]) i= max( (1:72) [ long> seq(-180,175,5) ]) gridcell =ts( Instr[i,j,],start=c(1850,1),freq=12) gridcell=window(gridcell,end=2008.99); #length(gridcell) #1908 return(gridcell) } #X= array(gridcell,dim=c(12,length(gridcell)/12)) #count= apply(!is.na(X)[months,],2,sum) #summer= round(ts( apply(X[months,],2,mean,na.rm=T) ,start=tsp(gridcell)[1]),2) #summer[count14] writeLines(fred,"temp") tree= read.fwf( "temp", widths =c(6,4,rep(c(4,3),10)),colClasses=c("character",rep("numeric",21))) } else { tree= read.fwf( url, widths =c(6,4,rep(c(4,3),10)),colClasses=c("character",rep("numeric",21)))} tree[,1]= factor(tree[,1]) levels(tree[,1])<-sub(" +$", "",levels(tree[,1])); count= tree[,c(1:2,seq(4,22,2))] tree= tree[,c(1:2,seq(3,21,2))] temp=tree[,3:12]==9990; tree[,3:12][temp]=NA id=levels(tree[,1]);K=length(id);K #1022 names(tree)[1:2]=c("id","year") long=NULL for(i in 1:K) { temp=(tree$id== levels(tree$id)[i]) test=tree[temp,] test[1,2]=10*floor(test[1,2]/10) A=outer(test[,2], 0:9, function(x,y) x+y) working=cbind( c(t(A)), c(t(test[,3:12])) ) M=nrow(working);M temp1=!is.na(match(1:M,1:10)) & is.na(working[,2]);sum(temp1) temp2=!is.na(match(1:M,(M-9):M )) & (is.na(working[,2])|(working[,2]== -9999)) ;sum(temp2) working=working[!temp1&!temp2,] working=data.frame( working); working$id=levels(tree[,1])[i] long=rbind(long,working) } names(long)=c("year","rw","id") long$id=factor(long$id) long$age=long$id; levels(long$age)=tapply(long$year,long$id,min) long$age=as.numeric(as.character(long$age)); long$age=long$year-long$age+1 long=long[,c("id","year","age","rw")] return(long) } RCS.chronology<-function (tree,method="nls") { dimnames(tree)[[2]][4]<-"x" #sometimes names are "rw", sometimes "mxd" tree<-tree[!is.na(tree$x),] tree$id<-factor(tree$id) meanx = tapply(tree$x,tree$age,mean) meanx=meanx[!is.na(meanx)]; N=length(meanx);N #189 biwt.meanx = tapply(tree$x,tree$age,biweight.mean) biwt.meanx=biwt.meanx[!is.na(biwt.meanx)]; Nb=length(biwt.meanx);Nb #189 if(method=="nls") { #preferred method fm <- nls(x ~ A+B*exp(-C*age),data = tree, start = list( A=mean(tree$x,na.rm=T)/4,B = mean(tree$x,na.rm=T), C= .01 ), alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10)); B<-coef(fm); fitted<- function(x) B[1]+B[2]*exp(-B[3]*x ) } if(method=="esper") { max0=max(tree$age) p=spline.response(.1*Nb);p # 2.695515e-05 fit<-smooth.spline(1:Nb,biwt.meanx,spar=(1-p)) ; fitted<- approxfun(1:Nb,fit$y) B=fit$y } tree$smooth<-fitted(tree$age) tree$delta<- tree$x/tree$smooth series<- c(unlist( tapply(tree$delta,factor(tree$year),mean,na.rm=TRUE) )) series<-ts(series,start=min(tree$year)) RCS.chronology<-list(series=series,fit=fitted,delta=tree$delta,smooth=tree$smooth,coefficients=B,method=method,biweight.mean=biwt.meanx ,mean=meanx,strata= as.character(levels(tree$id))) RCS.chronology } ##RCS.chronology RCS.chronology<-function (tree,method,parameters=NULL,interval=5,units=1) { dimnames(tree)[[2]][4]<-"x" tree<-tree[!is.na(tree$x),] tree$agegroup<-median(1:interval)+interval*trunc( (tree$age-1)/interval) tree$id<-factor(tree$id) tree$agegroupf<-factor(tree$agegroup) #year5f<-factor(tree$group) year5avg <-tapply(tree$x,tree$agegroupf,mean,na.rm=TRUE)/units agegroup<-as.numeric(names(year5avg)) N<-length(agegroup) z<-cbind(agegroup,year5avg) z<-data.frame(z) if (method=="mod.neg.exp") {fm <- nls( year5avg ~ A+B*exp(-C*agegroup),data = z, start = list(A=mean(z$year5avg)/4, B = mean(z$year5avg), C= .01 ), alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10)); B<-coef(fm); fitted<- B[1]+B[2]*exp(-B[3]*(1:agegroup[N])) } if (method=="nls") {fm <- nls(x ~ A+B*exp(-C*age),data = tree, start = list( A=mean(tree$x,na.rm=T)/4,B = mean(tree$x,na.rm=T), C= .01 ), alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10)); B<-coef(fm); fitted<- B[1]+B[2]*exp(-B[3]*(1:max(tree$age))) } if (method=="linear") {fm<-lm(year5avg~agegroup,data=z); B<-coef(fm); fitted<-B[1]+B[2]*(1:agegroup[N])} if (method=="briffa") { B<-parameters; if (length(B) ==2) fitted<- B[1]+B[2]*(1:agegroup[N]); if (length(B) ==3) fitted<- B[1]+B[2]*exp(-B[3]*(1:agegroup[N])) } if (method=="RCS.spline") { if(!is.null(parameters)) p<-spline.response(parameters[1]*N) else p<-spline.response(.66667*N) fit.model<- smooth.spline( agegroup, year5avg, spar = (1-p)); fitted<-predict(fit.model,(1:agegroup[N]))$y; B<-fitted} fitted<-ts( fitted, start=1, end=agegroup[N]) tree$delta<-NA tree$smooth<-NA for (k in 1:length(levels(tree$id))) { temp<-(tree$id==levels(tree$id)[k]) y<-tree$year[temp] w<-tree$x[temp] n<-length(w); tree$smooth[temp]<-fitted[tree$age[temp]] tree$delta[temp]<-w / tree$smooth[temp] } yearf<-factor(tree$year) series<-tapply(tree$delta,yearf,mean,na.rm=TRUE) series<-ts(series,start=min(tree$year),end=max(tree$year)) RCS.chronology<-list(series,z,tree$delta,tree$smooth,B,method) names(RCS.chronology)<-c("series","hist.curve","delta","smooth","coefficients","method") RCS.chronology } ########## ##8. BIWEIGHT MEAN ################ MAD= function(x) median ( abs( x-median(x))) biweight.mean=function(x) { K=length(x) ; if(K<4) return(NA) else { w=rep(0,K); delta=xnew=median(x); S=MAD(x); if(S==0) return(xnew) else { iter=1; tol=.01;c=6 while (delta