knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(forester)
library(tidyverse)

How do we calculate the volume of a stand?

We usually only have the following information:

Because the mean height of a stand is very sensitive to different forestry operations, we have adopted a measure of the maximum attained height instead: dominant height.

Firstly, you should now that the dominant height of a stand has had many differing definitions:

Let's create som random data to demonstrate how the two variables compare.

Create 200 dummy data observations following a normal distribution with a mean of 6 metres and a standard deviation of 2.5 metres.

norm_distr <- data.frame(rnorm(200, mean=6, sd=2.5))
norm_distr <- norm_distr %>% rename(height = rnorm.200..mean...6..sd...2.5.)

Plot the data as a histogram

norm_distr %>% ggplot(aes(x=height)) + geom_density()+theme(legend.position="bottom")

We'll calculate the mean and top height, using the mean of the tallest 10% definition. Maybe you have not yet come into contact with stats::quantile.

norm_distr_summary <- norm_distr %>% summarise(mean_stand_height=mean(height), dominant_height=quantile(height,probs=0.9, names=FALSE))
norm_distr_summary

Let's see how they are affected by a low-grading thinning (removing the suppressed trees).

norm_thinned <- norm_distr %>% filter(height > quantile(height, probs=0.3))

norm_thinned_summary <- norm_thinned %>% summarise(mean_stand_height=mean(height), dominant_height=quantile(height,probs=0.9,names=FALSE))
norm_thinned_summary

What do the distributions look like now?

#Create identifying variables
norm_distr <- norm_distr %>% mutate(type="Normal")
norm_thinned <- norm_thinned %>% mutate(type="Thinned")

# Append the tables.
distr_joined <- rbind(norm_distr, norm_thinned)

#Plot the distributions of the stand before and after the low-grading.
distr_joined %>% ggplot(aes(x=height,fill=type))+geom_histogram(alpha=0.5, bins=50)+theme(legend.position="bottom")

Let's calculate the percentual change in our variables:

#mean stand height
((norm_thinned_summary[1,1]/norm_distr_summary[1,1])-1)*100
#dominant stand height
((norm_thinned_summary[1,2]/norm_distr_summary[1,2])-1)*100

That's a large difference in the percentual increase! Dominant height seems like a better variable to track the potential height growth of a stand.

But what does this have to do with the total volume produced in a stand?! Easy tiger, we're getting there!

Height ~ Diameter relationship

Trees, unlike us, have an unlimited amount of organs - each year they grow new appendages in the form of leaves, branches, and a new shoot. To support their height growth, they also have cambial growth - which soon enough dies, and is transformed into Xylem. Trees increase their diameter with age.

Can we model what this relationship looks like?

Here's some data on tree heights (m) and diameters (mm)

#Näslund 1936 - y-1.3~I(dia^2)/I((a + b*dia)^2)
#Spruce a=19.0425458
#Spruce b= 0.2302333
hmod <- function(data){
 1.3 + I(data^2)/I((19.0425458+ 0.2302333*data)^2)
} 
diametersequence <- seq(1,300,15)

diametersequence <- data.frame(diameter=diametersequence)

diametersequence <- diametersequence %>% mutate(height=hmod(diameter))
rm(hmod)
#add randomness
diametersequence <- diametersequence %>% group_by(diameter) %>% mutate(noise= runif(1,min=-0.08*height, max=0.08*height)) %>% mutate(height2=height+noise)

Spruce_height_diameter <- diametersequence %>% select(diameter, height2) %>% rename(height=height2)
library(knitr)
kable(Spruce_height_diameter[1:5,],caption="A number of Spruce heights and diameters from a stand")
Spruce_height_diameter %>% ggplot()+geom_point(aes(x=diameter,y=height))

Näslund (1937) showed that you could satisfactorily model the relationship between the height and diameter of pines above breast height (1.3m) along a model of the form: $$Y \sim 1.3 + \frac{DBH^2}{(\alpha+\beta*DBH^2)}$$

We're going to fit a model of the same form to describe the relation between Spruce height and diameter in a stand.

Spruce_height_diameter <- Spruce_height_diameter %>% mutate(species="Spruce") %>% group_by(species) %>% nest() %>%  mutate(
  model= map(data, ~nls(height ~ 1.3 + I(diameter^2)/I((a+ b*diameter)^2), start=list(a=1,b=1),data=.))) %>% 
  # extract formula from each model, convert to one-sided form, &
  # replace coefficients with fitted values, & store in dataframe
  # as character string 
  rowwise() %>%
  mutate(func = formula(model) %>% 
           as.character() %>% 
           magrittr::extract(3) %>%
           gsub("diameter", ".x", ., fixed = T) %>%
           gsub("a", model$m$getPars()[1], .) %>%
           gsub("b", model$m$getPars()[2], .) %>%
           paste("~", ., collapse = "")) %>%
  ungroup()

print(Spruce_height_diameter)

Let's see what it looks like.

ggplot(Spruce_height_diameter[[2]][[1]])+
  theme(legend.position = "bottom")+
  geom_point(aes(x=diameter,y=height)) + #All the data
    # add formula in each row as a separate geom_function layer
  lapply(seq(1, nrow(Spruce_height_diameter)),
         function(i) geom_function(fun = rlang::as_function(formula(Spruce_height_diameter$func[i])),
                                   aes(colour = Spruce_height_diameter$species[i]))) +

  # change legend name (can also change palette / labels / etc.)
  scale_colour_discrete(name = "Spruce")

Great! Now we can predict the heights of all measured trees in our stand from the representative subset we measured!

We'll create a stand to calculate the volume of. Planted at 1.5 m interval, we have 4400 plants outside Umeå at a latitude of c. 65 °N.

forester::plant_spacing(stems_per_ha = 4400)

We've measured the tree diameters:

#We assume normal diameter distribution
diameters <- data.frame(diameters= rnorm(4000, mean=130, sd=50))

#The stand has been heavily thinned from below (40%)
diameters <- diameters %>% filter(diameters > quantile(diameters, probs=0.4))

From the diameters and our earlier relation between height and diameter from a subset of this stand we calculate the heights of the trees:

Height ~ 1.3 + I(diameter^2)/I((19.1265151455858 + 0.229470870878339 * diameter)^2)

diameters <- diameters %>% mutate(height= 1.3 + diameters^2/(19.1265151455858 + 0.229470870878339 * diameters)^2)

Now that we have the diameters and heights of all stems, we can approximate the volume of the stand with help of the forester::tree.volume function.

diameters <- diameters %>% mutate(volume= tree.volume(dbh.cm=diameters/10, height.m=height,species.latin="Picea abies", latitude = 65))

The total stand standing volume is (dm^3):

diameters %>% summarise("Total volume" = sum(volume))

In order to get m^3 from dm^3, divide by 1'000.

diameters %>% mutate(volume_m3 = volume / 1000) %>% summarise("Total volume"= sum(volume_m3))


Silviculturalist/forester documentation built on April 20, 2024, 2:13 p.m.