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)

[AndyBunn]

Great. I'll work on this. Will need a Rd file as well.

[sklesse]

Hi Andy,

I'm pasting the content of the .rd file here. First time for everything, let me know if something is missing or off.

data(bakker.data) would refer to the rwl and d2pith files I attached above. A short description for the required example files could be: "Ring-width series and metadata from a beech site in Zofingen, Switzerland."

Cheers, Stefan

\name{bakker{dplR}} \alias{bakker{dplR}}

\title{ Basal Area Increment (Bakker) } \description{ Convert multiple ring-width series to basal area increment (i.e., ring area) following the proportional method of Bakker (2005). } \usage{ \method{print}{bakker}(rwl, bakker.d2pith) }

\arguments{ \item{rwl}{a data.frame with series as columns and years as rows such as that produced by read.rwl} \item{d2pith}{A mandatory data.frame containing four variables. Variable one (series in the example below) gives the series ID as either characters or factors. These must exactly match colnames(rwl). Variable two (PO) gives the number of rings estimated to be missing to the pith. Variable three (d2pith) gives the estimated distance to the pith (mm). Variable four (diam) gives the diameter at breast height in cm.} } \details{ This converts ring-width series (mm) to ring-area series (mm squared) (aka basal area increments) based on the diameter of the tree, the missing distance to the pith and the missing number of rings to the pith, following the proportional method for reconstructing historical tree diameters by Bakker (2005). It prevents bai.out transformations from producing negative increments when the sum of all ring widths in a series is larger than DBH/2. It prevents bai.in transformations from producing too small values when the sum of all ring widths in a series is smaller than DBH/2.

} \value{ A list containing the following objects: \item{DBHhist_raw}{data.frame with the reconstructed diameter for each series with column names, row names and dimensions of rwl.} \item{baiBakker_raw}{data.frame with the basal area increments for each series with column names, row names and dimensions of rwl.}

} \references{ Bakker, J.D., 2005. A new, proportional method for reconstructing historical tree diameters. Canadian Journal of Forest Research 35, 2515–2520. https://doi.org/10.1139/x05-136

} \author{ Code by Stefan Klesse }

\examples{

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

data(bakker.data) dummy<-bakker(bakker.rwl,bakker.d2pith)

pith<-bakker.d2pith[,c(1,3)] diam<-bakker.d2pith[,c(1,4)] diam$diam<-diam$diam*10 #convert dbh (cm) to dbh (mm)

plot.ts(bai.out(bakker.rwl,diam)[,54],main=colnames(bakker.rwl)[54],las=1,ylab="BAI mm2",ylim=c(0,6000)) lines(bai.out(bakker.rwl)[,54],col=1,lty=2) abline(h=0,lty=2,col="grey") lines(bai.in(bakker.rwl,pith)[,54],col=2) lines(bai.in(bakker.rwl)[,54],col=2,lty=3) lines(ts(dummy$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) = ",bakker.d2pith[54,4]),adj=0) text(0,2500,label=paste("2x core length w/o pith offset = ",2/10(sum(bakker.rwl[,54],na.rm=T)),"cm"),adj=0,) text(0,2000,label=paste("2x pith offset = ",bakker.d2pith[54,3]2/10,"cm"),adj=0) text(0,1500,label=paste("2x core length w/ pith offset = ",2/10*(sum(bakker.rwl[,54],na.rm=T)+bakker.d2pith[54,3]),"cm"),adj=0)

plot.ts(bai.out(bakker.rwl,diam)[,53],main=colnames(bakker.rwl)[53],las=1,ylab="BAI mm2",ylim=c(0,6000)) lines(bai.out(bakker.rwl)[,53],col=1,lty=2) abline(h=0,lty=2,col="grey") lines(bai.in(bakker.rwl,pith)[,53],col=2) lines(bai.in(bakker.rwl)[,53],col=2,lty=3) lines(ts(dummy$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) = ",bakker.d2pith[53,4]),adj=0) text(0,2500,label=paste("2x core length w/o pith offset = ",2/10(sum(bakker.rwl[,53],na.rm=T)),"cm"),adj=0,) text(0,2000,label=paste("2x pith offset = ",bakker.d2pith[53,3]2/10,"cm"),adj=0) text(0,1500,label=paste("2x core length w/ pith offset = ",2/10*(sum(bakker.rwl[,53],na.rm=T)+bakker.d2pith[53,3]),"cm"),adj=0)

}

[AndyBunn]

Ok. Made an initial commit. Still need to fuss with the documentation and data structures a bit to be consistent with other dplR conventions. Note that I'm breaking the data into the rwl file as zof.rwl and the ancillary data into zof.anc. I'll bug you for info about those shortly. But can you give me a site name and species for zof (e.g., gp.rwl is titled "Ponderosa Pine Ring Widths from Gus Pearson Natural Area")? That'll be a start.

[AndyBunn]

@sklesse , You should be able to install the dev version via install.packages("dplR", repos = "https://andybunn.r-universe.dev/") I'll likely muck around with the help files some. I know you mentioned a forthcoming paper. I can add that to the refs section. So species, site name, and approximate location for the data will be a good start.

[sklesse]

Hi Andy, the site name is Zofingen, Switzerland, 47.297 N, 6.963 E, 500 m elevation. It's European beech (Fagus sylvatica). So, "European Beech Ring Widths from Zofingen, Switzerland" should do it.

The data were part of my 2018 paper: Klesse, S., Babst, F., Lienert, S., Spahni, R., Joos, F., Bouriaud, O., Carrer, M., Filippo, A.D., Poulter, B., Trotsiuk, V., Wilson, R., Frank, D.C., 2018. A Combined Tree Ring and Vegetation Model Assessment of European Forest Growth Sensitivity to Interannual Climate Variability. Global Biogeochemical Cycles 32, 1226–1240. https://doi.org/10.1029/2017GB005856

[sklesse]

Regarding the forthcoming paper, I will let you know as soon as it approaches acceptance.

[AndyBunn]

Ok. It's all in there. Take a look and suggest any changes you might want to see. I can do a release of 1.7.7 on CRAN if that's helpful for the paper.

[sklesse]

Somehow R shoots me an error message after downloading your dev version.

library(dplR) Error: package or namespace load failed for ‘dplR’ in dyn.load(file, DLLpath = DLLpath, ...): unable to load shared object '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so': dlopen(/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so, 0x0006): Library not loaded: /opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib Referenced from: <08EE96C1-8CC6-3019-90F1-2E15B59418C9> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so Reason: tried: '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/usr/local/lib/libgfortran.5.dylib' (no such file), '/usr/lib/libgfortran.5.dylib' (no such file, not in dy In addition: Warning message: package ‘dplR’ was built under R version 4.3.3

[AndyBunn]

Via the r-univ?

install.packages("dplR", repos = https://andybunn.r-universe.dev/)

I can see all the binaries.

https://andybunn.r-universe.dev/bin/windows/contrib/4.3/dplR_1.7.7.zip

From: sklesse @.> Date: Friday, May 31, 2024 at 12:56 AM To: AndyBunn/dplR @.> Cc: Andy Bunn @.>, Comment @.> Subject: Re: [AndyBunn/dplR] Bakker method (Issue #21)

Somehow R shoots me an error message after downloading your dev version.

library(dplR) Error: package or namespace load failed for ‘dplR’ in dyn.load(file, DLLpath = DLLpath, ...): unable to load shared object '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so': dlopen(/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so, 0x0006): Library not loaded: /opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib Referenced from: <08EE96C1-8CC6-3019-90F1-2E15B59418C9> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so Reason: tried: '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/usr/local/lib/libgfortran.5.dylib' (no such file), '/usr/lib/libgfortran.5.dylib' (no such file, not in dy In addition: Warning message: package ‘dplR’ was built under R version 4.3.3

— Reply to this email directly, view it on GitHubhttps://github.com/AndyBunn/dplR/issues/21#issuecomment-2141424630, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AC7UCXKVHHRYGFS7J7JIVNDZFAULJAVCNFSM6AAAAABGUJLFMKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNBRGQZDINRTGA. You are receiving this because you commented.Message ID: @.***>

[AndyBunn]

Installs fine from a win machine in my lab via install.packages("dplR", repos = https://andybunn.r-universe.dev/)

So let me know.

From: Andy Bunn @.> Date: Friday, May 31, 2024 at 9:06 AM To: AndyBunn/dplR @.>, AndyBunn/dplR @.> Cc: Comment @.> Subject: Re: [AndyBunn/dplR] Bakker method (Issue #21) Via the r-univ?

install.packages("dplR", repos = https://andybunn.r-universe.dev/)

I can see all the binaries.

https://andybunn.r-universe.dev/bin/windows/contrib/4.3/dplR_1.7.7.zip

From: sklesse @.> Date: Friday, May 31, 2024 at 12:56 AM To: AndyBunn/dplR @.> Cc: Andy Bunn @.>, Comment @.> Subject: Re: [AndyBunn/dplR] Bakker method (Issue #21)

Somehow R shoots me an error message after downloading your dev version.

library(dplR) Error: package or namespace load failed for ‘dplR’ in dyn.load(file, DLLpath = DLLpath, ...): unable to load shared object '/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so': dlopen(/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so, 0x0006): Library not loaded: /opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib Referenced from: <08EE96C1-8CC6-3019-90F1-2E15B59418C9> /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/dplR/libs/dplR.so Reason: tried: '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/System/Volumes/Preboot/Cryptexes/OS/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/opt/gfortran/lib/gcc/aarch64-apple-darwin20.0/12.2.0/libgfortran.5.dylib' (no such file), '/usr/local/lib/libgfortran.5.dylib' (no such file), '/usr/lib/libgfortran.5.dylib' (no such file, not in dy In addition: Warning message: package ‘dplR’ was built under R version 4.3.3

— Reply to this email directly, view it on GitHubhttps://github.com/AndyBunn/dplR/issues/21#issuecomment-2141424630, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AC7UCXKVHHRYGFS7J7JIVNDZFAULJAVCNFSM6AAAAABGUJLFMKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDCNBRGQZDINRTGA. You are receiving this because you commented.Message ID: @.***>

[AndyBunn]

Also @sklesse . I did submit v1.7.7 to CRAN. It was necessary after I transferred the dplR repo to the opendendro GH. So the Bakker func will be live there on CRAN in a few hours. Can sill modify on dev version which will be 1.7.8 here on GH.