###TREE RING UTILITIES ################################### #1. make.crn - makes time series from *.crn format allowing for some oddball formats. #this is a very old script and could be updated #2. make.rwl_new - makes organized measurement data set from Tucson format (allowing for some oddball formats) #update of make.rwl #3. countrycode - utility used in make.crn, make.rwl to look up country code from identification form #4. RCS.chronology - emulation of RCS chronology with several methods. I usually use method="nls" #Sep 28, 2009: this is a very old program for me and I re-drafted it in improved style today without changing output # older version is still at bottom under RCS.chronology.old #5. jacoby.chronology - emulaiton of "conventional" chronology using modified neg exp. This has been tested carefully. #6. spline chronology - ## spline.response - from LTRR protocols #6. countf #utility to give time series of ring counts #7. agef - utility to calculate ring age from rwl file #8. Biweight Mean; MAD # rbar - from rwl data.frame, returns mean rbar, rbar for tree, rbar time series, eps time series # eps - using Wigley 1984 definition calculates eps from count (n) and rbar # sss - signal subsample strength from Wigley 1984 definition # some older functions seldom used ############# #1. make.crn make.crn<-function(loc,widths0=c(6,4,rep( c(4,3),10)),N= -1,readtype="REG",url0="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/chronologies") { #may need some more quality control on oddballs but will work for STNDRD, SCHW, Long if(url0=="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/updates") url<-file.path(url0,paste(loc,"crn",sep=".")) else url<-file.path(url0,countrycode(loc),paste(loc,"crn",sep=".")) #see if it is a Schweingruber format by checking fifth row schwein<-c("RAW","RES","STD") h<-try(readLines(url,n=5)) if (class(h)=="try-error") make.crn<-NA else { if( substr(h[4],6,6)=="-") readtype<-"LONG" type<-substr(h,83,85) if ( !is.na( match(type[5],schwein) ) ) readtype<-"SCHW"} #read in different ways for STNDRD and SCHW formats #SCHW format has a number of options; only the RAW version is read here #this is controlled by calculating N in Schwein case; otherwise N= -1 in STNDRD if(readtype=="SCHW") { h<-readLines(url); temp<-(substr(h,83,85)=="STD") h=h[temp] h=h[4:length(h)] write.table(h,"temp.dat",row.names=FALSE,quote=FALSE,col.names=FALSE) test=read.fwf("temp.dat",widths=c(6,4,rep(c(4,3),10))) start0=floor(test[1,2]/10)*10 test=test[,seq(3,21,2)];dim(test) test=ts(c(t(test)),start=start0) temp=(test>=9990);test[temp]=NA temp=!is.na(test);r0=range( time(test)[temp]) temp=(time(test)>=r0[1])&(time(test)<=r0[2]) make.crn=ts(test[temp],start=r0[1]) } if(readtype=="LONG") widths0[1:2]<-c(5,5) if(!(readtype=="SCHW")) { h<-read.fwf(url,widths=widths0,skip=3,nrow=N,blank.lines.skip=TRUE,colClasses=c("factor",rep("numeric",21))) #blan.lines.skip added June 2006 start1<-h[1,2] #} temp<-!is.na(h[,1]) h<-h[temp,] #adjusts for blank lines if (h[nrow(h),2]< h[(nrow(h)-1),2]) h<-h[1:(nrow(h)-1),] #this adjusts for some versions which have statistics in last row chron.crn<-h[,seq(3,21,2)] if (start1>2000) print(loc, "ERROR 1 - START AFTER 2000") temp<-(chron.crn==9990) chron.crn[temp]<-NA chron.crn<- c(t(chron.crn)) temp<-!is.na(chron.crn); make.crn<- ts(chron.crn[temp],start=start1) } make.crn } #EXAMPLE #chron.crn<-make.crn("arge048") ################## #2. make.rwl_new make.rwl =function(url,skip0=0,widths0 =c(8,4,rep(6,10))) { # "id","year","rw" trailer.verify<- 999 trailing.marker<- (-9999) h<-try(readLines(url,n=5)) if (class(h)=="try-error") make.rwl<-NA else { tree<-read.fwf(file=url,widths=widths0,fill=TRUE,skip=skip0,blank.lines.skip=TRUE) levels(tree[,1])<-sub(" +$", "",levels(tree[,1])); temp<-is.na(tree[,1]) #typically 0 and next segment is not invoked if(sum(temp)>0) tree<-tree[!temp,] (test3=sapply(tree,class)[2:12]); if (sum(!is.na(match(test3,"factor")))>0) { fred=readLines(url); fred=fred[!(fred=="")]; fred=fred[4:length(fred)]; n=nchar(fred);treeid=substr(fred,1,8);treeid=gsub(" +$","",treeid) fred=substr(fred,9,n) fred=gsub("-9999"," -99 ",fred) write.table(fred,"temp.dat",quote=FALSE,row.names=FALSE,col.names=FALSE) tree=read.table("temp.dat",fill=TRUE) tree[tree== -99]= -9999 tree=data.frame(treeid,tree) } #endif #check for trailing 999 markers; very annoying; also -9999 trailers 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[,1]== levels(tree[,1])[i]) test=tree[temp,] 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) } #endloop long=long[,c(3,1,2)];names(long)=c("id","year","rw") return(long) } #end class-error } #tree<-make.rwl("russ022w") #save(tree,file="d:/climate/data/briffa.2008/TornFin.rwl.tab") ############################ #3. countrycode: utility used in making crn and rwl download.file("http://www.climateaudit.org/scripts/tree/country.code.tab","temp.dat",mode="wb");load("temp.dat") download.file("http://www.climateaudit.org/scripts/tree/country.list.tab","temp.dat",mode="wb");load("temp.dat") countrycode<-function(id) { temp<-match(substr(id,1,4),country.list[,1]) if(is.na(temp)) countrycode<-"northamerica/usa" else countrycode<- country.list[temp,2] countrycode } ############################# #4. RCS.chronology ############################# # this calculates an RCS chronology for a tree ring data set: the script is one of my earlier scripts # it works, but I'd write it more cleanly now. # it permits a couple of different 'methods' - in this case, variant styles of RCS #my preferred method is a nls fit 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) 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) meanx = tapply(tree$x,tree$age,biweight.mean) meanx=meanx[!is.na(meanx)]; N=length(meanx);N #189 p=spline.response(.1*N);p # 2.695515e-05 fit<-smooth.spline(1:N,meanx,spar=(1-p)) ; fitted<- approxfun(1:N,fit$y) B=NA } 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) RCS.chronology } ################################ #5. jacoby.chronology #this fits with straight lines and negative exponentials jacoby.chronology<-function(tree,method.chron="neg.exp") { tree$delta<-NA tree$smooth<-NA dimnames(tree)[[2]][4]<-"x" temp.na<-is.na(tree$x) treef<-factor(tree$id) N<-length(levels(treef)) jac.coef<-array(NA,dim=c(N,3)) for (k in 1:N) { temp<-(tree$id[!temp.na]==levels(treef)[k]) age<-tree$age[!temp.na][temp] w<-tree$x[!temp.na][temp] n<-length(w); Data<-cbind(age,w) Data<-data.frame(Data) fm1<-lm(w~age,data=Data) if (method.chron=="linear") {fit<-coef(fm1)[1]+coef(fm1)[2]*age;jac.coef[k,1:2]<-coef(fm1)} if(method.chron=="neg.exp") { if (fm1$coefficients[2]>0) {fit<-rep(mean(w),length(age));jac.coef[k,1:2]<-c(mean(w),0)} else { fm <- try( nls( w ~ A+B*exp(-C*age),data = Data, start = list( A = fm1$coefficients[1]+fm1$coefficients[2]*200, B= - fm1$coefficients[2]*200, C=0.01 ), alg = "default", trace = TRUE,control=nls.control(maxiter=200, tol=1e-05, minFactor=1e-10))); if (class(fm)=="try-error") {fit<-rep(mean(w),length(age));jac.coef[k,1:2]<-c(mean(w),0)} else {a<-summary(fm)$parameters[,1]; fit<-a[1]+a[2]*exp(-a[3]*age); jac.coef[k,]<-a } } #else } #neg.exp tree$smooth[!temp.na][temp]<-fit tree$delta [!temp.na][temp]<-w/tree$smooth[!temp.na][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)) row.names(jac.coef)=levels(treef) jacoby.chronology<-list(series,tree$delta,tree$smooth,jac.coef) names(jacoby.chronology)<-c("series","delta","smooth","coef") jacoby.chronology } ################### ## spline.chronology ######################## #for methods spline at 67% of tree length and RCS; #does on long-form data.frame tree #returns tree$spline; tree$smooth and chronology (ts) #data.frame must be passed in two columns - year and rw/mxd #call is is tree[,c("id","year","age","rw")] if more than one item ##spline.response http://www.ltrr.arizona.edu/~dmeko/notes_7.pdf # for 100 years, p = 1.5585e-05 #updated from http://www.ltrr.arizona.edu/~dmeko/notes_7.pdf # for 100 years, p = 1.5585e-05 May 2007 spline.response<-function(year) { g<-cos (2*pi/year) spline.response<-1/ ( ( 1-.5)/.5 * (g+2)/ (12*(g-1)^2 ) +1 ) spline.response } ##spline.chronology #tag - if COEFECHA then 32 year spline otherwise 2/3 of tree length spline.chronology<-function(tree,tag="COFECHA") { tree$delta<-NA tree$smooth<-NA dimnames(tree)[[2]][4]<-"x" temp.na<-is.na(tree$x) treef<-factor(tree$id) N<-length(levels(treef)) for (k in 1:N) { temp<-(tree$id[!temp.na]==levels(treef)[k]) y<-tree$year[!temp.na][temp] w<-tree$x[!temp.na][temp] n<-length(w); if (tag=="COFECHA") p<-spline.response(32) if (tag=="briffa") p<-spline.response(.666667*n); z<-smooth.spline(y,w,spar=(1-p)) ; tree$smooth[!temp.na][temp]<-z$y tree$delta [!temp.na][temp]<-w/z$y } yearf<-factor(tree$year) series<-tapply(tree$delta,yearf,mean,na.rm=TRUE) series<-ts(series,start=min(tree$year),end=max(tree$year)) spline.chronology<-list(series,tree$delta,tree$smooth) names(spline.chronology)<-c("series","delta","smooth") spline.chronology } #################### #6. countf - return time series of tree counts countf=function(tree,year=NULL) { if(class(year)=="NULL") year=min(tree$year):max(tree$year) count= ts( c(unlist(tapply(!is.na(tree$rw),factor(tree$year,levels=year),sum))) ,start=min(year)) count[is.na(count)]=0 return(count) } ######################## #7. agef - function to return age of each ring as column in database agef=function(tree) { tree$age=factor(tree$id) levels(tree$age)=tapply(tree$year,tree$id,min) tree$age=as.numeric(as.character(tree$age)) tree$age=tree$year-tree$age+1 tree=tree[,c(1,2,4,3)] return(tree) } ########## ##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 .85 # n*a > (1+(n-1)*a) * .85 # n*a > .85 - a*.85 + n*a* .85 # n*a (1-.85) > .85*(1 - a) # n > (.85/.15) * (1 - a) / a #VARIANCE ADJUSTMENT ACCORDING TO OSBORN ET AL 1999 #page 461: variance reduced according to Osborn et al. (1998) Dendrochronologia #see also NAture 1998 # "Each regional mean thus obtained tended to have greater variance during years when few chronologies were available to contribute to the average; #this effect was corrected for by scaling by the square root of the ffective number of independent samples available in each year # n'= n/(1+(n-1)*rbar) #this is a function to implement rbar-adjustment #rbar0 is a vector of length of the total series, calculated in various ways externally bo.adjust<-function(x,count,rbar){ count.eff<- count/(1+(count-1)*rbar.series);#plot.ts(count.eff) #goes from 2.3 at the beginning to maximum of 5.35 # NN<-max(count,na.rm=T) count.eff.max<- max(count.eff) ;#count.eff.max# 5.35 var.adj<-sqrt ( count.eff/count.eff.max) ; #equation 7 of Osborn et al. #factor goes from about 0.7 at beginning to 1 bo.adjust<-x *var.adj bo.adjust } ##sss: defined Wigley 1984 page 204 and equation 24 second part, discussed on page 206 ##sss= Rbar_{n,N}^2 a_hat in Wigley=rbar below # N is maximum count presumed to be recent ; n is actual sss<-function(rbar,n,N) n*(1+(N-1)*rbar)/(N*(1+(n-1)*rbar)) #equation 24 ##eps: defined Wigley 1984 page 204 and equation 26 ##eps = Rbar_{N}^2 ## Wigley eq 33: eps= SNR/(1+SNR) ## Wigley eq 32 SNR= (N*rbar)/(1-rbar) eps<-function(rbar,N) (N*rbar)/ (1+ (N-1)*rbar) #equation 26 snr<-function(rbar,N) (N*rbar)/ (1 -rbar) #equation 26 #MAKES TEMPLATE IN NXN MATRIX TO ASSIST IN MEAN CORRELATION CALCULATION #this returns a logical matrix with TRUE in left-subdiagonal to simplify mean correlation calculations template<-function(N) { template<-NULL for (k in 1:N) { template<- c(template, rep(NA,k), rep(TRUE,N-k) ) } template<-array(template, dim =c(N,N)) template } #CALCULATES OVERLAPS FOR SPARSE TIME SERIES recon larger than M #this returns a logical matrix of same shape as correlation matrix showing whether the two series have >M overlap overlap<-function (recon,M=20) { test<-!is.na(recon) overlap<-t(test) %*% test overlap[overlapSTART) if (sum(temp)>0) { treetest<-tree[temp,]; treetemp<-ts ( rep(NA,max(treetest[,2])- min(treetest[,2])+1), start=min(treetest[,2]), end=max(treetest[,2]) ) ; test5<-match(treetest[,2],c(tsp(treetemp)[1]:tsp(treetemp)[2])) treetemp[test5]<-treetest[,4] ; recon<-ts.union(recon,treetemp); name<-c(name,j) } } dimnames(recon)[[2]]<- sub(" +$", "", levels(treef)[name])#removes trailing blanks collate<-recon collate } #note: this includes some quality control checking on NGDC data - #the tag verify2 keeps track of few negative values, of which there are a few #the tag start0 and end0 keeps track of start-dates and end-dates. #some "floating series" intentionally have start-dates way into the future #a couple of series had RES as well as STD series and the RES data was manually deleted #these entries are not relevant to the MBH subsets # see script "ngdc.qc" ##make.chron.are make.chron.ars<-function(loc.in) { h<-readLines(loc.in,n=-1) start0<-as.numeric(substr(h[2],13,18)) length0<-as.numeric(substr(h[2],4,8)) h<-h[3:length(h)] fred<-as.numeric (substr(h,22,24)) N<-length(h) test<-c(1:N)[is.na(fred)] #this looks for first alphanumeric header h<-h[1:(min(test)-1)] write.table(h,file="d:/temp/h.txt",row.names=FALSE,quote=FALSE) g<-read.fwf("d:/temp/h.txt",skip=1, widths=rep( c(4,2),13)) g<-g[,seq(1,25,2)] g<-c(t(g)) make.chron.ars<-ts(g,start=start0,end=start0+length0-1) make.chron.ars } ### update details update.details<-function(loc,url0="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/measurements") { if(!(url0=="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/updates")) url0<-file.path(url0,countrycode(loc)) url<-file.path(url0,paste(loc,"rwl",sep=".")) h<-try(readLines(url,n=5)) if (class(h)=="try-error") h<-readLines(file.path("ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/chronologies",countrycode(loc),paste(loc,"crn",sep=".")) ,n=5) dummy<-details[1,] dummy[1,]<-NA dummy$id<-loc dummy$location<-sub(" +$", "",substr(h[1],10,40)) dummy$type<-substr(h[1],62,65) dummy$long<-as.numeric(substr(h[2],52,57)) dummy$long<-dummy$long/100 dummy$lat<-as.numeric(substr(h[2],47,51)) dummy$lat<-dummy$lat/100 dummy$altitude<-as.numeric(substr(h[2],41,44)) dummy$altitude<-dummy$altitude dummy$start<-as.numeric(substr(h[2],67,71)) dummy$start<-dummy$start dummy$end<-as.numeric(substr(h[2],72,76)) dummy$end<-dummy$end dummy$author<-sub(" +$", "",substr(h[3],10,50)) update.details<-dummy update.details } country.out<-function(id) { temp<-match(substr(id,1,4),country.list[,1]) if(is.na(temp)) country.out<-"northamerica" else country.out<- country.list[temp,2] if ((substr(id,1,4)=="cana")|(substr(id,1,4)=="mexi")) country.out<-"northamerica" country.out } ##CHECK MEASUREMENT CODE ##make.rwl #this is a very old script and can probably be improved, #it's been re-edited to post on internet and has not been checked yet against oddball formats #earlier edition was thoroughly checked # a few oddball formats need to be checked out #FUNCTION TO READ RWL FILE AND MAKE LONG-FORM VERSION #returns list # tree, #two verification lists #the parameter long adjusts format when dates are -1000AD are less make.rwl_old<-function(loc,long=FALSE,method="wdcp",skip0=3,url0="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/measurements") { trailer.verify<- 999 trailing.marker<- (-9999) if(!(url0=="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/updates")&!(method=="local")) url0<-file.path(url0,countrycode(loc)) url<-file.path(url0,paste(loc,"rwl",sep=".")) h<-try(readLines(url,n=5)) if (class(h)=="try-error") make.rwl<-NA else { if( substr(h[4],6,6)=="-") long<-TRUE #if (long) widths0<- c(8,4,rep(6,10)) ##else widths0<- c(7,5,rep(6,10)) test<-as.numeric(substr(h,9,12)) skip0<-sum(is.na(test)) # if(!is.numeric(tree[2,4])) tree<-read.fwf(file=url,widths=c(7,6,rep(6,10)),fill=TRUE,skip=3,blank.lines.skip=TRUE) #usually FALSE #if(!is.numeric(tree[1,3])) tree<-read.fwf(file=url,widths=c(7,5,rep(6,10)),fill=TRUE,skip=4,blank.lines.skip=TRUE) #usually FALSE #if(!is.numeric(tree[1,3])) tree<-read.fwf(file=url,widths=c(7,5,rep(6,10)),fill=TRUE,skip=5,blank.lines.skip=TRUE) #usually FALSE tree<-read.fwf(file=url,widths=widths0,fill=TRUE,skip=skip0,blank.lines.skip=TRUE) #slightly different format when long if(class(tree[,2])=="factor") {widths0<- c(8,4,rep(6,10)); tree<-read.fwf(file=url,widths=widths0,fill=TRUE,skip=skip0,blank.lines.skip=TRUE) } tree[,1]<-factor(tree[,1]) levels(tree[,1])<-sub(" +$", "",levels(tree[,1])); #check for oddball input formats temp<-is.na(tree[,1]) #typically 0 and next segment is not invoked if(sum(temp)>0) tree<-tree[!temp,] # if (sum(temp)>1) {h<-readLines(url); # h<-h[!temp]; # write.table(h,file="d:/temp/tree.txt",quote=FALSE,col.names=FALSE,row.names=FALSE); # tree<-read.fwf(file="d:/temp/tree.txt",widths=c(8,4,rep(6,10)),fill=TRUE,skip=3,blank.lines.skip=TRUE) # } #tests for oddball format #class(tree$V4)<-"numeric" #change from integer #tree<-tree[!is.na(tree[,1]),] #remove NA values in ID if any; usually not test3=sapply(tree,class)[2:12]; if (sum(!is.na(match(test3,"factor")))>0) { fred=readLines(url); fred=fred[!(fred=="")]; fred=fred[4:length(fred)]; n=nchar(fred);treeid=substr(fred,1,8);treeid=gsub(" +$","",treeid) fred=substr(fred,9,n) fred=gsub("-9999"," -99 ",fred) write.table(fred,"temp.dat",quote=FALSE,row.names=FALSE,col.names=FALSE) tree=read.table("temp.dat",fill=TRUE) tree[tree== -99]= -9999 tree=data.frame(treeid,tree) } #check for trailing 999 markers; very annoying; also -9999 trailers n<-nrow(tree) temp.marker<-(tree[2:n,1]==tree[1:(n-1),1]) temp.marker<-c(temp.marker,FALSE) temp.marker<-!temp.marker test1<- apply( tree[temp.marker,3:12],1,min,na.rm=TRUE) temp.trailer1<- (test1== -9999) id1<-tree[temp.marker,1][temp.trailer1] test2<- apply( tree[temp.marker,3:12],1,max,na.rm=TRUE) temp.trailer2<- (test2== 999) id2<-tree[temp.marker,1][temp.trailer2] temp.trailer<- ((tree[temp.marker,3:12]==999)|(tree[temp.marker,3:12]== -9999)) &!is.na(tree[temp.marker,3:12]) tree[temp.marker,3:12][temp.trailer]<-NA #usually results in substitutions ##IF MIXED if ( (sum(temp.trailer1)>0) & (sum(temp.trailer2)>0) ) trailer.verify<-0 if ( !(sum(temp.trailer1)>0) & (sum(temp.trailer2)>0) ) trailer.verify<-0 print(paste(loc,trailer.verify)) if (trailer.verify ==0) { temp<- !is.na(match( tree[,1], id2)) tree[temp,3:12]<-10* tree[temp,3:12]} ### m<-length(levels(tree[,1])) tree<-cbind(tree,c(rep(0,n))) #make dummy column line.tree.no<-match(tree[,1],levels(tree[,1])) #this is index for tree year.start<-tapply(tree[,2],tree[,1],min) #NOTE: this doesn't work if ID is 8-digits (nv514) line.start<-year.start[line.tree.no] #this gives v- vector tree[,13]<-line.start #tree<-cbind(tree[gl(n,1,10*n),1],tree[gl(n,1,10*n),2],tree[gl(n,1,10*n),13],stack(tree[1:nrow(tree),3:12])[,1] );##reshapes NGDC to long form tree<-data.frame(tree[gl(n,1,10*n),1],tree[gl(n,1,10*n),2],tree[gl(n,1,10*n),13],stack(tree[1:nrow(tree),3:12])[,1] );##reshapes NGDC to long form temp2<-t(array(rep(1:10,n),dim=c(10,n))); tree<-cbind(tree,c(temp2));##sets tag for column of NGDC entry to add back in long format #reshape into long format with dim = n*10 x 9 #one column keeps indicator for potential reconstruction #id<-as.character(tree[,1])#scrapped tree<-data.frame(tree[,1],tree[,2]+tree[,5]-1,tree[,2]+tree[,5]-tree[,3],tree[,4]); dimnames(tree)[[2]]<-c("id","year","age","rw") #use NA removal to remove trailing values on first and last records, but not missing values with na.marker temp<-!is.na(tree[,4]) tree<-tree[temp,]; #set negative values as NA temp<-(tree[,4]<0) if (sum(temp)>0) { verify<-tree[temp,]; temp1<- (verify[,4] == -999); verify.NA<-verify[temp1,]; verify.negative<-verify[!temp1,]; tree[temp,4]<-NA } else {verify.NA <-NULL; verify.negative<-NULL } #order by tree and by year order.tree<-order(tree[,1],tree[,2]); make.rwl<-tree[order.tree,]; #possibly put into dataframe and edit column names } make.rwl } #tree<-make.rwl("russ022w") check.rwl<-function(loc,long=FALSE,skip0=3,url0="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/measurements") { trailer.verify<- 999 trailing.marker<- (-9999) if(!(url0=="ftp://ftp.ncdc.noaa.gov/pub/data/paleo/treering/updates")) url0<-file.path(url0,countrycode(loc)) url<-file.path(url0,paste(loc,"rwl",sep=".")) h<-try(readLines(url,n=5)) if (class(h)=="try-error") check.rwl<-NA else { if( substr(h[4],6,6)=="-") long<-TRUE #if (long) widths0<- c(8,4,rep(6,10)) ##else widths0<- c(7,5,rep(6,10)) test<-as.numeric(substr(h,9,12)) skip0<-sum(is.na(test)) # if(!is.numeric(tree[2,4])) tree<-read.fwf(file=url,widths=c(7,6,rep(6,10)),fill=TRUE,skip=3,blank.lines.skip=TRUE) #usually FALSE #if(!is.numeric(tree[1,3])) tree<-read.fwf(file=url,widths=c(7,5,rep(6,10)),fill=TRUE,skip=4,blank.lines.skip=TRUE) #usually FALSE #if(!is.numeric(tree[1,3])) tree<-read.fwf(file=url,widths=c(7,5,rep(6,10)),fill=TRUE,skip=5,blank.lines.skip=TRUE) #usually FALSE tree<-read.fwf(file=url,widths=widths0,fill=TRUE,skip=skip0,blank.lines.skip=TRUE) #slightly different format when long if(class(tree[,2])=="factor") {widths0<- c(8,4,rep(6,10)); tree<-read.fwf(file=url,widths=widths0,fill=TRUE,skip=skip0,blank.lines.skip=TRUE) } tree[,1]<-factor(tree[,1]) levels(tree[,1])<-sub(" +$", "",levels(tree[,1])); #check for oddball input formats temp<-is.na(tree[,1]) #typically 0 and next segment is not invoked if(sum(temp)>0) tree<-tree[!temp,] #check for trailing 999 markers; very annoying; also -9999 trailers n<-nrow(tree) temp.marker<-(tree[2:n,1]==tree[1:(n-1),1]) temp.marker<-c(temp.marker,FALSE) temp.marker<-!temp.marker test1<- apply( tree[temp.marker,3:12],1,min,na.rm=TRUE) temp.trailer1<- (test1== -9999) id1<-tree[temp.marker,1][temp.trailer1] test2<- apply( tree[temp.marker,3:12],1,max,na.rm=TRUE) temp.trailer2<- (test2== 999) id2<-tree[temp.marker,1][temp.trailer2] temp.trailer<- ((tree[temp.marker,3:12]==999)|(tree[temp.marker,3:12]== -9999)) &!is.na(tree[temp.marker,3:12]) tree[temp.marker,3:12][temp.trailer]<-NA #usually results in substitutions ##IF MIXED if ( (sum(temp.trailer1)>0) & (sum(temp.trailer2)>0) ) trailer.verify<-0 if ( !(sum(temp.trailer1)>0) & (sum(temp.trailer2)>0) ) trailer.verify<-0 check.rwl<-trailer.verify } check.rwl }