linking up DBH with core measurements

[ValentineHerr]

@teixeirak,

I couldn't find a straightforward answer in Ian's analysis for calculating DBH so I am going with a simple solution which is:

1. find out what main census year is closest to the last core measurement year (only looking at main census years where the tree was alive (and measured)).
2. assign the dbh to that year, if it exists, in the core measurements or, if it doesn't, extrapolate the dbh to the closest year of core measurement, taking the average tree ring increment of the past few years (e.g. past 6 years if the latest core measurement is 6 year prior to the dbh measurement).
3. calculate the dbh of other years by subtracting the tree ring measurements (cummulatively).

Is that okay?

Now, there is one issue, the only info I have about the cores is the tag number. But sometimes, the tag corresponds to a multi-stem tree. What stem should I take the dbh from in that case? was there a rule to always core the main stem in the coring protocol?

[teixeirak]

From @mcgregorian1 : "The code is in my script "canopy_position_analysis", section 5e, "new method". That code chunk references the trees_all table made from section 4a."

[teixeirak]

I couldn't find a straightforward answer in Ian's analysis for calculating DBH so I am going with a simple solution which is:

1. find out what main census year is closest to the last core measurement year (only looking at main census years where the tree was alive (and measured)).

2. assign the dbh to that year, if it exists, in the core measurements or, if it doesn't, extrapolate the dbh to the closest year of core measurement, taking the average tree ring increment of the past few years (e.g. past 6 years if the latest core measurement is 6 year prior to the dbh measurement).

3. calculate the dbh of other years by subtracting the tree ring measurements (cummulatively).

Is that okay?

This is okay for a first pass, but we should ultimately use the bark allometries, based on our data published here, to account for changes in bark thickness. @mcgregorian1's script does that.

[teixeirak]

Now, there is one issue, the only info I have about the cores is the tag number. But sometimes, the tag corresponds to a multi-stem tree. What stem should I take the dbh from in that case? was there a rule to always core the main stem in the coring protocol?

I think its generally safe to assume main stem, although I want to double check with @gonzalezeb regarding the mortality census cores. It looks like we have cored multiple stems in at least one instance last year:

[mcgregorian1]

I agree my code is a little complicated - I definitely would have written it differently if I was starting it now. However, the main steps I did are here, with the actual code (clarified a bit) below. Let me know if you need something else.

Steps

1. Calculate diameter_nobark for 2008 = DBH.mm.2008-2*bark.depth.mm
2. log-transform both diam_nobark_2008 (x) and bark.depth.mm (y)
3. Fit a linear model, and use model to predict log(bark.depth.mm)
4. Take exponent of bark.depth.mm and make sure predicted values look good.
5. Get mean bark thickness per species in 2008.
   • The equation for calculating old dbh, using 1999 as an example, is
   • dbh1999 = dbh2008 - 2(ring.width2013 - ring.width1999) - 2(bark.depth2008) + 2(bark.depth1999)
   • using the dataset from calculating the regression equations, we can get mean bark thickness per species in 2008.
6. Thus, the only value we're missing is bark depth in 1999.
   • This is ok, because we can calculate from the regression equation per each species (all we need is diam_nobark_1999). Calculate diam_nobark_1999 using
   • diam_nobark_1999 = dbh2008 - 2(bark.depth2008) - 2(sum(ring.width1999:ring.width2008))
7. Calculate bark thickness using regression equation per appropriate sp
   • log(bark.depth.1999) = intercept + log(diam_nobark)*constant
8. Add to solution from #6 to get full dbh1999
   • dbh1999 = diam_nobark_1999 + 2*bark.depth.1999

Code

bark <- read.csv("data/traits/SCBI_bark_depth.csv")
bark <- bark[bark$species %in% unique(trees_all$sp), ]

#1. Calculate diameter_nobark for 2008 = DBH.mm.2008-2*bark.depth.mm
bark$diam_nobark_2008.mm <- bark$DBH.mm.2008 - 2*bark$bark.depth.mm 

#2. log-transform both diam_nobark_2008 (x) and bark.depth.mm (y)
#3. Fit a linear model, and use model to predict log(bark.depth.mm)
bark$predict_barkthick.ln.mm <- NA
bark$predict_barkthick.ln.mm <- 
      ifelse(bark$species == "caco", -1.56+0.416*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "cagl", -0.393+0.268*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "caovl", -2.18+0.651*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "cato", -0.477+0.301*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "fram", 0.418+0.268*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "juni", 0.346+0.279*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "litu", -1.14+0.463*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "qual", -2.09+0.637*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "qupr", -1.31+0.528*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "quru", -0.593+0.292*log(bark$diam_nobark_2008.mm),
      ifelse(bark$species == "quve", 0.245+0.219*log(bark$diam_nobark_2008.mm),
                           bark$predict_barkthick.ln)))))))))))

#4. Take exponent of bark.depth.mm and make sure predicted values look good.
bark$predict_barkthick.mm <- exp(bark$predict_barkthick.ln.mm)
range(bark$predict_barkthick.mm - bark$bark.depth.mm)

#5. Get mean bark thickness per species in 2008.
## The equation for calculating old dbh, using 1999 as an example, is
## dbh1999 = dbh2008 - 2(ring.width2013 - ring.width1999) - 2(bark.depth2008) + 2(bark.depth1999)

## using the dataset from calculating the regression equations, we can get mean bark thickness per species in 2008.

##set up dbh dataframe
dbh <- trees_all[, c(1:4)]
scbi.stem1 <- read.csv(text=getURL("https://raw.githubusercontent.com/SCBI-ForestGEO/SCBI-ForestGEO-Data/master/tree_main_census/data/census-csv-files/scbi.stem1.csv"), stringsAsFactors = FALSE)
scbi.stem1$dbh <- as.numeric(scbi.stem1$dbh)
dbh$dbh2008.mm <- scbi.stem1$dbh[match(dbh$tree, scbi.stem1$tag)]

mean_bark <- aggregate(bark$bark.depth.mm, by=list(bark$species), mean) #mm
colnames(mean_bark) <- c("sp", "mean_bark_2008.mm")

dbh$mean_bark_2008.mm <- ifelse(dbh$sp %in% mean_bark$sp, mean_bark$mean_bark_2008.mm[match(dbh$sp, mean_bark$sp)], mean(bark$bark.depth.mm))
dbh$mean_bark_2008.mm <- round(dbh$mean_bark_2008.mm, 2)

#6. Thus, the only value we're missing is bark depth in 1999.
## This is ok, because we can calculate from the regression equation per each species (all we need is diam_nobark_1999).Calculate diam_nobark_1999 using
## diam_nobark_1999 = dbh2008 - 2*(bark.depth2008) - 2*(sum(ring.width1999:ring.width2008))

##define this column before loop
dbh$diam_nobark_old.mm <- 0

df <- rings #the original file read-in from read.rwl
colnames(df) <- gsub("A", "", colnames(df)) #remove "A"
colnames(df) <- gsub("^0", "", colnames(df)) #remove leading 0

cols <- colnames(df) #define cols for below
colnames(df) <- gsub("^", "x", colnames(df)) #add "x" to make calling colnames below feasible

for (j in seq(along=cols)){
  for (k in seq(along=colnames(df))){
    ring_ind <- cols[[j]]
    ring_col <- colnames(df)[[k]]

    if(j==k){
      #the output of this loop is 3 separate columns for each year's old dbh, hence why it is set to q as a dataframe before being combined below. Pointer_years_simple comes from #4d.
      q <- data.frame(sapply(pointer_years_simple, function(x){
        rw <- df[rownames(df)>=x, ]
        ifelse(dbh$year == x & dbh$tree == ring_ind, 
               dbh$dbh2008.mm - 2*(dbh$mean_bark_2008.mm) - sum(rw[rownames(rw) %in% c(x:2008), ring_col], na.rm=TRUE), 0)
      }))

      q$diam_nobark_old.mm <- q[,1] +q[,2] + q[,3] #add columns together
      # q$dbh_old.mm <- q[,1] +q[,2] + q[,3] + q[,4]
      dbh$diam_nobark_old.mm <- dbh$diam_nobark_old.mm + q$diam_nobark_old.mm #combine with dbh (it's the same order of rows) #mm
    }
  }
}

#7. Calculate bark thickness using regression equation per appropriate sp
## log(bark.depth.1999) = intercept + log(diam_nobark)*constant
## bark.depth.1999 = exp(log(bark.depth.1999))

#the full equation at the bottom is the regression equation for all these species put together. "fagr" is given a bark thickness of 0 because it is negligble
#these equations are the same as above in #3 of this code section
dbh$bark_thick_old.ln.mm <- NA
dbh$bark_thick_old.ln.mm <- 
  ifelse(dbh$sp == "caco", -1.56+0.416*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "cagl", -0.393+0.268*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "caovl", -2.18+0.651*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "cato", -0.477+0.301*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "fram", 0.418+0.268*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "juni", 0.346+0.279*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "litu", -1.14+0.463*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "qual", -2.09+0.637*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "qupr", -1.31+0.528*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "quru", -0.593+0.292*log(dbh$diam_nobark_old.mm),
  ifelse(dbh$sp == "quve", 0.245+0.219*log(dbh$diam_nobark_old.mm),
                      0)))))))))))

dbh$bark_thick_old.mm <- ifelse(dbh$sp == "fagr", 0, exp(dbh$bark_thick_old.ln.mm))

#8. Add to solution from #6 to get full dbh1999
## dbh1999 = diam_nobark_1999 + 2*bark.depth.1999
dbh$dbh_old.mm <- dbh$diam_nobark_old.mm + 2*dbh$bark_thick_old.mm

##NOTE
##The first time I ran this code I was getting NaNs for one tree (140939), because the dbh in 2008 was listed as 16.9. I double-checked this, and that was the second stem, which we obviously didn't core at a size of 1.69 cm (or 2.2 cm in 2013). The dbh is meant to be the first stem. However, there was confusion with the dbh in the field. 

trees_all$dbh_old.mm <- dbh$dbh_old.mm[match(trees_all$tree, dbh$tree) & match(trees_all$year, dbh$year)] #mm
trees_all$dbh_old.cm <- trees_all$dbh_old.mm/10
trees_all$dbh.ln.cm <- log(trees_all$dbh_old.cm)

[teixeirak]

Many thanks, @mcgregorian1 !

[gonzalezeb]

Now, there is one issue, the only info I have about the cores is the tag number. But sometimes, the tag corresponds to a multi-stem tree. What stem should I take the dbh from in that case? was there a rule to always core the main stem in the coring protocol?

I think its generally safe to assume main stem, although I want to double check with @gonzalezeb regarding the mortality census cores. It looks like we have cored multiple stems in at least one instance last year:

For that very reason I asked Alyssa to add the column "cored" in 2018 and 2019 data to indicate if a core was taken or not form a stem. For previous censuses, we would assume that yes, the mail stem was the one cored, but in general, the expectation was/is if the stem is dead, then core it.

[ValentineHerr]

I'll work on calculating dbh like Ian did.

Just for my record, I add plots of the cases where dbh goes <0 when dbh is calculated the simple way (keep removing the tree ring measurement each year). (Vertical dashed line is year 15 after first measurement).

Note: this is after cutting data back to 2000 for trees cored dead and removing first 15 years if an outlier (core measurement > 10) was found then.

[teixeirak]

I started a new issue that may fit better here: #9.

[ValentineHerr]

@mcgregorian1, do you by any chance remember where you got your data/traits/SCBI_bark_depth.csv file from? I would like to source it directly from where it belongs originally but I can't find it.

[ValentineHerr]

Also, @mcgregorian1, did you remove your script mentioned in this comment ?

[teixeirak]

@mcgregorian1, do you by any chance remember where you got your data/traits/SCBI_bark_depth.csv file from? I would like to source it directly from where it belongs originally but I can't find it.

@ValentineHerr, this comes from:

• Dryad data publication accompanying Anderson-Teixeira et al. 2015, Functional Ecology. The former contains data, latter describes allometries.

Ian and I decided it made sense to re-do the allometries.

[ValentineHerr]

Ok thanks @teixeirak , I thouht I could find the data in Github too.

[mcgregorian1]

Also, @mcgregorian1, did you remove your script mentioned in this comment ?

Hi Sorry I changed the name of the script and clarified the other manuscript code files to be more streamlined. The DBH code you're looking for is in climate_sensitivity_analysis.R, section #2d.

[ValentineHerr]

Thanks @mcgregorian1 !

Other question (that you can answer in person here today): About step 5 in your code - Why do you calculate mean bark thickness by species and not use the allometric equation (that take species and dbh as input) to calculate bark thickness in 2008?

[ValentineHerr]

@teixeirak , the allometric equations here do not include 'frni', 'fagr' and 'pist'. What is your preferred approach for these species?

[teixeirak]

FAGR has such thin bark that we assume negligible thickness across species. For FRNI, it probably makes most sense to use the FRAM bark allometry. For PIST, not sure... maybe we can find some bark thickness data, or ask our SMSC student to go measure bark thickness on several pines.

[teixeirak]

For now, just use an equation for any species (except FAGR) for PIST. It should have minimal influence on results.

[teixeirak]

Actually, for PIST, here's a pub with a bit of data: a 28.4 cm PIST in S IL had bark thickness of ~10mm. Let's use an equation for a species with a similar bark thickness. 28.4 cm,

[ValentineHerr]

I found fram to be the closest but maybe we can improve this solution

[teixeirak]

We'll run into this issue for other sites, so I need to think up a general solution.

@gonzalezeb, @cpiponiot, have either of you come across any bark thickness allometries?

[gonzalezeb]

I haven't but by doing a quick look online it seems that there is some literature on bark allometries out there.

Oak allometries for 8 species (Texas)

Various conifers (California)..

Now, if we need bark allometries for many species I think we won't find them.

[teixeirak]

This is now covered in this issue; closing.