Bakker method

[sklesse]

Hi Andy, here is the code for the bakker function, two files to work with, and two plots.
zoffinal.txt ZOF_d2pith.csv

rwl<-read.rwl("/Users/klesse_adobe/Dropbox/DEEP-tools/Bakker` DBH/zoffinal.txt")

#load file containing metadata
#d2pith must be a data.frame with 4 columns: 1) the same as colnames(rwl), 2) the number of rings estimated to be missing to the pith (PO), 3) the estimated distance to the pith (mm), colnames="d2pith", 4) The DBH (cm) at sampling date, colnames="diam"
#to be consistent with dplR function bai.out() diam should be converted to mm in the d2pith file, change line 102 (*10)

d2pith<-read.csv("/Users/klesse_adobe/zof_d2pith.csv",header=T, sep=",")
d2pith<-d2pith[,-1]



bakker <-function (rwl, d2pith) { 
  if (!is.data.frame(rwl)) 
    stop("'rwl' must be a data.frame")
  
  if (ncol(rwl) != nrow(d2pith)) 
    stop("dimension problem: ", "'ncol(rw)' != 'nrow(d2pith)'")
  
  if (ncol(d2pith) != 4) 
    stop("dimension problem: ", "'ncol(d2pith)' != 4")
  if (!all(d2pith[, 1] == names(rwl))) {
    print(data.frame(rwlNames = names(rwl), seriesID = d2pith[, 1], test = d2pith[, 1] == names(rwl)))
    stop("series ids in 'd2pith' and 'rwl' do not match exactly.")
  }
  if (!is.data.frame(d2pith)) {
    stop("'d2pith' must be a data.frame")
  }
  rwl_d2pith <- d2pith[, "d2pith"]
  rwl_diam <- d2pith[, "diam"]
  rwl_po<-d2pith[,"PO"]
  
  #filling in the missing rings towards the pith is recommended when there are multiple cores per tree. 
  #average multiple BAI time series of a tree, not raw ring width measurements.
  #Otherwise, one can end up with an artifact, i.e. one negative increment when two radii 
  #have different growth rates and the faster growing sample has more rings.
  #The filled rings will be deleted at the end.
    
    rwl_filled<-NULL
    for (i in 1:length(rwl)){
      if(rwl_d2pith[i] !=0){
        temp<-data.frame(rep(rwl_d2pith[i]/rwl_po[i],times=rwl_po[i])) #create a dummy time series with the number of estimated rings missing 
        temp2<-na.omit(rwl[,i])
        select_remove<-c(attr(temp2,"na.action"))
        if(!is.null(select_remove)){
        temp2<-data.frame(temp2,row.names=rownames(rwl)[-select_remove]) 
        }else{
          temp2<- data.frame(temp2,row.names=rownames(rwl))
        }
        rownames(temp)<-c((as.numeric(rownames(temp2))[1]-rwl_po[i]):(as.numeric(rownames(temp2))[1]-1)) 
        colnames(temp)<-colnames(temp2)
        rwl_filled[[i]]<-ts(rbind(temp,temp2),start=as.numeric(rownames(temp))[1]) #paste the dummy time series and ring width measurements
      } else { 
        temp3<-data.frame(rwl[,i],row.names=rownames(rwl))
        rwl_filled[[i]]<-ts(temp3,start=as.numeric(rownames(temp3))[1])
      }
      
      
    }
    names(rwl_filled)<-colnames(rwl)
    rwl_filled<-do.call(cbind,rwl_filled)
    rownames(rwl_filled)<-c(tsp(rwl_filled)[1]:tsp(rwl_filled)[2])

  
 
    
  cum.na <- function(x) { #cumulative sum that changes NAs to zeros. helper function
    x[which(is.na(x))] <- 0
    return(cumsum(x))
  }
  
  IP <- apply(rwl_filled, 2, sum,na.rm=T) #how long is your measured radius including missing distance to pith?
  IH <- rwl_filled
  G <- IH
  for(i in 1:ncol(rwl_filled)){IH[,i] <-cum.na(rwl_filled[,i])}  #cumulative sum
  for(i in 1:ncol(IH)){G[,i]<-(IP[i]-(IH[,i]))/IP[i]} #proportional cumulative sum, reversed
  
  DBHhist<-G
  colnames(DBHhist)<-colnames(rwl)
  for (i in 1:ncol(DBHhist)){
    DBHhist[,i]<-rwl_diam[i]-(G[,i] *rwl_diam[i]) #subtract proportional DBH from measured DBH, the final proportional historic DBH
  }
  
  for(i in 1:ncol(DBHhist)){  #clean zeros, fill with NAs
    temp<-which(is.na(rwl_filled[,i]))
    if(any(which(diff(temp)>1))){
      DBHhist[which(is.na(rwl_filled[,i]))[-which(diff(temp)>1)],i]<-NA
    }else{
      DBHhist[rev(which(is.na(rwl_filled[,i]))),i]<-NA 
    }
  }
  
  if(nrow(rwl_filled)-nrow(rwl)==0){
    DBHhist2<- DBHhist[1:nrow(DBHhist),]
    DBHhist2[is.na(rwl)]<-NA
  }else{
    DBHhist2<- DBHhist[-c(1:(nrow(rwl_filled)-nrow(rwl))),]
    DBHhist2[is.na(rwl)]<-NA
  }
 

  out <- DBHhist
  n.vec <- seq_len(nrow(DBHhist))
  for (i in seq_len(ncol(rwl))) {
    dat <- DBHhist[,i]
    dat2 <- na.omit(dat)

   
    r0 <- dat2/2*10 #if dbh in cm, if dbh in mm, then remove *10; /2 to convert to radius
    bai <- pi * (diff(r0 * r0))
    na <- attributes(dat2)$na.action
    no.na <- n.vec[!n.vec %in% na]
    out[no.na[-1], i] <-  bai
  }
  
  if(nrow(rwl_filled)-nrow(rwl)==0){
    out2<- out[1:nrow(out),]
    out2[is.na(rwl)]<-NA
  }else{
  out2<- out[-c(1:(nrow(rwl_filled)-nrow(rwl))),]
  out2[is.na(rwl)]<-NA #this removes all the filled dummy values from rwl_filled above
  }


  output<-list(DBHhist_raw=DBHhist2,baiBakker_raw=out2)
  #I termed it raw, because at some point back in time I wanted to expand the function
  #to fill the innermost rings of the sample with less rings towards the center with data from the other core(s).
  #Never got to it, though.
  
}
 


######some quick bakker example plots#####

#samples 53 and 54 are from a quite asymmetrically growing tree, the b-sample misses also the innermost 4 cm.
#small differences in the a-sample [,53] that got close to the pith, but note the negative BAI from bai.out in the first few years because the core length + pith offset is larger than 0.5*DBH
silly<-bakker(rwl,d2pith)
pith<-d2pith[,c(1,3)]
diam<-d2pith[,c(1,4)]
diam$diam<-diam$diam*10 #


plot.ts(bai.out(rwl,diam)[,54],main=colnames(rwl)[54],las=1,ylab="BAI mm2",ylim=c(0,6000))
lines(bai.out(rwl)[,54],col=1,lty=2)
abline(h=0,lty=2,col="grey")
lines(bai.in(rwl,pith)[,54],col=2)
lines(bai.in(rwl)[,54],col=2,lty=3)
lines(ts(silly$baiBakker_raw[,54],start=1),col=4,lwd=2)


legend(legend=c("Bakker","outside-in with DBH","inside-out with POE","outside-in no DBH","inside-out no POE"),
       bty="n",x="topleft",lty=c(1,1,1,2,3),col=c(4,1,2,1,2),lwd=c(2,1,1,1,1))


text(0,3000,label=paste("DBH (cm) = ",d2pith[54,4]),adj=0)
text(0,2500,label=paste("2x core length w/o pith offset = ",2/10*(sum(rwl[,54],na.rm=T)),"cm"),adj=0,)
text(0,2000,label=paste("2x pith offset = ",d2pith[54,3]*2/10,"cm"),adj=0)
text(0,1500,label=paste("2x core length w/ pith offset = ",2/10*(sum(rwl[,54],na.rm=T)+d2pith[54,3]),"cm"),adj=0)


plot.ts(bai.out(rwl,diam)[,53],main=colnames(rwl)[53],las=1,ylab="BAI mm2",ylim=c(0,6000))
lines(bai.out(rwl)[,53],col=1,lty=2)
abline(h=0,lty=2,col="grey")
lines(bai.in(rwl,pith)[,53],col=2)
lines(bai.in(rwl)[,53],col=2,lty=3)
lines(ts(silly$baiBakker_raw[,53],start=1),col=4,lwd=2)

legend(legend=c("Bakker","outside-in with DBH","inside-out with POE","outside-in no DBH","inside-out no POE"),
       bty="n",x="topleft",lty=c(1,1,1,2,3),col=c(4,1,2,1,2),lwd=c(2,1,1,1,1))



text(0,3000,label=paste("DBH (cm) = ",d2pith[53,4]),adj=0)
text(0,2500,label=paste("2x core length w/o pith offset = ",2/10*(sum(rwl[,53],na.rm=T)),"cm"),adj=0,)
text(0,2000,label=paste("2x pith offset = ",d2pith[53,3]*2/10,"cm"),adj=0)
text(0,1500,label=paste("2x core length w/ pith offset = ",2/10*(sum(rwl[,53],na.rm=T)+d2pith[53,3]),"cm"),adj=0)