Structure v3

Calibration of forest structure parameters (crown radius allometry and mortality) in TROLL 3.1.8:

Crown radius parameters variation

We gathered all Tallo measurements of crown radius and diameter in a radius of 10-km of Tapajos for trees above 10-cm. We randomly selected 10 trees per class diameter of decimetre a thousand times and computed crown radius intercepts and slopes. We obtained a parameter distribution with a correlation of 0.91, a linear regression of \(b_{CR} = -0.39+0.57 \times a_{CR}\), and a standard variation of 0.08. Dashed lines shows Paracou values, which are lower than the observations.

Code 
tapajos <- tibble(
  site = c("Tapajos"),
  latitude = c(-2.85667),
  longitude = c(-54.9588900),
) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>% 
  st_buffer(10^4)
params <- read_csv("data/species/tallo.csv") %>% 
  filter(longitude <= -39, longitude >= -79, latitude >= -18, latitude <= 10) %>% 
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326) %>% 
  st_intersection(tapajos) %>% 
  filter(!is.na(crown_radius_m), !is.na(stem_diameter_cm)) %>% 
  mutate(dbh_m = stem_diameter_cm * 0.01) %>% 
  filter(dbh_m >= 0.1) %>% 
  mutate(class_dbh_dm = floor(dbh_m * 10)) %>% 
  mutate(rep = list(1:10^3)) %>% 
  unnest(rep) %>% 
  group_by(rep) %>% 
  sample_n(10, replace = T) %>% 
  unique() %>%
  st_drop_geometry() %>% 
  ungroup() %>% 
  nest_by(rep) %>% 
  mutate(coefs = list(coef(lm(log(crown_radius_m) ~ log(dbh_m), data = data)))) %>% 
  mutate(a = coefs[1], b = coefs[2]) %>% 
  select(rep, a, b) %>% 
  filter(a > 0, b > 0)
Code 
ggplot(params, aes(a, b)) +
  geom_vline(xintercept = 2.13, linetype = "dashed", col = "green") +
  geom_hline(yintercept = 0.63, linetype = "dashed", col = "green") +
  geom_point() +
  geom_density2d() +
  theme_bw() +
  xlab(expression(a[CR])) +
  ylab(expression(b[CR])) +
  ggpubr::stat_cor() +
  geom_smooth(col = "blue", method = "lm")

Parameters space

We prepared the parameter space for calibration defining low limit, high limit, and step for a, m, and the error of b. b was computed using the linear equation plus the error. We used \(3 \sigma\) for the error. We obtained a calibration grid of 11,193 simulations. Dashed lines shows Paracou values which are thus lower. The last figure shows the corresponding projected allometries of crown radius.

Code 
pars <- data.frame(parameter = c("a", "error_b", "m"),
                   low = c(1.8, -4*0.08, 0.025),
                   high = c(2.5, 4*0.08, 0.06),
                   by = c(0.05, 0.05, 0.005))
grid <- data.frame(a = seq(pars$low[1], pars$high[1], by = pars$by[1])) %>% 
  mutate(error_b = list(seq(pars$low[2], pars$high[2], by = pars$by[2]))) %>% 
  unnest(error_b) %>% 
  mutate(b = -0.39 + 0.57*a + error_b) %>% 
  mutate(m = list(seq(pars$low[3], pars$high[3], by = pars$by[3]))) %>% 
  unnest(m)
Code 
kable(pars)
parameter low high by
a 1.800 2.50 0.050
error_b -0.320 0.32 0.050
m 0.025 0.06 0.005
Code 
grid %>% 
  ggplot(aes(a, b)) +
  geom_vline(xintercept = 2.13, linetype = "dashed") +
  geom_hline(yintercept = 0.63, linetype = "dashed") +
  geom_point(col = "lightgrey") +
  geom_abline(intercept = -0.37, slope = 0.57, col = "blue") +
  theme_bw() +
  xlab(expression(a[CR])) +
  ylab(expression(b[CR]))

Code 
grid %>% 
  ggplot(aes(a, m)) +
  geom_vline(xintercept = 2.13, linetype = "dashed") +
  geom_hline(yintercept = 0.025, linetype = "dashed") +
  geom_point(col = "lightgrey") +
  theme_bw() +
  xlab(expression(a[CR])) +
  ylab(expression(m))

Code 
grid %>% 
  mutate(dbh = list(seq(1, 150, 1)/100)) %>% 
  unnest(dbh) %>% 
  mutate(cr = exp(a+b*log(dbh))) %>% 
  ggplot(aes(dbh, cr)) +
  geom_line(aes(col = a, group = paste(a, b))) +
  theme_bw() +
  scale_colour_viridis_c(expression(a[CR])) +
  xlab("Diameter (m)") + ylab("Crown radius (m)") +
  geom_line(col = "black", linetype = "dashed",
            data = data.frame(dbh = seq(1, 150, 1)/100) %>% 
              mutate(cr = exp(2.13+0.63*log(dbh))))

Observations

For Tapajos and according to Rice et al. (2004) :

  • 470 live trees per hectare with diameter at breast height (dbh) ≥10
  • Aboveground live biomass was 143.7 ± 5.4 Mg C/ha
  • Stem density -0.058 (±0.003) dbh_cm + 2.912 (±0.076) for dbh < 40
  • Stem density -0.016 (±0.001) dbh_cm + 1.451 (±0.04) for dbh > 40

For Paracou:

  • 612 ± 37 live trees per hectare with diameter at breast height (dbh) ≥10 (from inventory)
  • Aboveground live biomass was 197.5-221.5 Mg C/ha in 2007 -> 209.5 ± 17.0 Mg C/ha (Rutishauser et al. 2010)
  • Stem density is shown below
Code 
abd_par <- read_tsv("outputs/inventories.tsv", show_col_types = FALSE) %>% 
  filter(site == "Paracou", type == "default") %>% 
  filter(dbh > 10) %>% 
  group_by(plot) %>% 
  summarise(abundance = n()/unique(plot_size)) %>% 
  ungroup()
# paste(round(median(abd_par$abundance)), "±", round(sd(abd_par$abundance)))
obs_tap <- data.frame(dbh_class = seq(10, 150, 5)) %>% 
  mutate(abundance = ifelse(dbh_class < 40, 
                            10^(-0.058*dbh_class+2.912), 
                            10^(-0.016*dbh_class+1.451))) %>% 
  mutate(site = "Tapajos")

obs_par <- read_tsv("outputs/inventories.tsv", show_col_types = FALSE) %>% 
  filter(site == "Paracou", type == "default") %>% 
  filter(dbh > 10) %>% 
  mutate(dbh_class = cut(dbh, 
                         breaks = c(seq(10, 155, by = 5)), 
                         labels = c(seq(10, 150, by = 5)))) %>% 
  mutate(dbh_class = as.numeric(as.character(dbh_class))) %>% 
    group_by(plot, dbh_class) %>% 
  summarise(abundance = n()/unique(plot_size)) %>% 
  group_by(dbh_class) %>% 
  summarise(abundance = median(abundance)) %>% 
  mutate(site = "Paracou")
# obs <- bind_rows(obs_tap, obs_par)
obs <- bind_rows(obs_tap, obs_par)
ggplot(obs, (aes(dbh_class, abundance, col = site))) +
  geom_line() +
  theme_bw() +
  scale_y_log10(labels = scales::comma) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  ggtitle(paste("Total:", round(sum(obs$abundance))))

Parameters effect

Code 
calib <- read_tsv("results/calib_str3.tsv") %>% 
  mutate(dbh_class = cut(dbh_cm, 
                         breaks = seq(10, 205, by = 5), 
                         labels = seq(10, 200, by = 5),
                         include.lowest = TRUE)) %>% 
  mutate(dbh_class = as.numeric(as.character(dbh_class))) %>% 
  group_by(site, a, b, m, dbh_class) %>% 
  summarise(abundance = n()/4, agb = sum(agb)/4)
write_tsv(calib, "outputs/calib_structure.tsv")

\(a_{CR}\)

Code 
dat <-  read_tsv("outputs/calib_structure.tsv") %>% 
  group_by(site, a, b, m) %>% 
  summarise(agb = sum(agb)/10^3/2,
            abundance = sum(abundance)) %>% 
  filter(m == 0.025)
cowplot::plot_grid(
 dat %>% 
  ggplot(aes(a, b, col = agb)) +
  geom_point() +
  facet_wrap(~ site, scales = "free_y") +
  theme_bw() +
  scale_color_viridis_c("AGB") +
  xlab(expression(a[CR|m==0.025])) + ylab(expression(b[CR|m==0.025])),
   dat %>% 
  ggplot(aes(a, b, col = abundance)) +
  geom_point() +
  facet_wrap(~ site, scales = "free_y") +
  theme_bw() +
  scale_color_viridis_c("Abundance") +
  xlab(expression(a[CR|m==0.025])) + ylab(expression(b[CR|m==0.025])), 
 nrow = 2)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  group_by(site, a, b, m) %>% 
  summarise(abundance = sum(abundance), 
            agb = sum(agb)/10^3/2) %>% 
  ggplot(aes(abundance, agb, col = a)) +
  geom_point() +
  theme_bw() +
  scale_color_viridis_c(expression(a[CR])) +
  facet_wrap(~ site)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  filter(m == 0.025) %>% 
  mutate(berr = -0.39 + 0.57*a - b) %>% 
  filter(berr > 0, berr < 0.05) %>% 
  ggplot(aes(dbh_class, abundance)) +
  geom_line(aes(col = a, group = paste(a, b, m))) +
  geom_point(data = filter(obs, dbh_class < 150), col = "red") +
  scale_y_log10(labels = scales::comma) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  scale_color_viridis_c(expression(a[CR|m==0.025|berr==0.02])) +
  theme_bw() +
  facet_wrap(~ site, nrow = 2)

\(b_{CR}\)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  group_by(site, a, b, m) %>% 
  summarise(abundance = sum(abundance), 
            agb = sum(agb)/10^3/2) %>% 
  ggplot(aes(abundance, agb, col = b)) +
  geom_point() +
  theme_bw() +
  scale_color_viridis_c(expression(b[CR])) +
  facet_wrap(~ site)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  filter(m == 0.025) %>% 
  mutate(berr = -0.39 + 0.57*a - b) %>% 
  filter(berr > 0, berr < 0.05) %>% 
  ggplot(aes(dbh_class, abundance)) +
  geom_line(aes(col = b, group = paste(a, b, m))) +
  geom_point(data = filter(obs, dbh_class < 150), col = "red") +
  scale_y_log10(labels = scales::comma) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  scale_color_viridis_c(expression(b[CR])) +
  scale_color_viridis_c(expression(b[CR|m==0.025|a==2])) +
  theme_bw() +
  facet_wrap(~ site, nrow = 2)

\(m\)

Code 
dat <-  read_tsv("outputs/calib_structure.tsv") %>% 
  group_by(site, a, b, m) %>% 
  summarise(agb = sum(agb)/10^3/2,
            abundance = sum(abundance)) %>% 
  mutate(berr = -0.39 + 0.57*a - b) %>% 
  filter(berr > 0, berr < 0.05)
cowplot::plot_grid(
 dat %>% 
  ggplot(aes(a, m, col = agb)) +
  geom_point() +
  facet_wrap(~ site, scales = "free_y") +
  theme_bw() +
  scale_color_viridis_c("AGB") +
  xlab(expression(a[CR|berr==0.02])) + ylab(expression(m[berr==0.02])),
   dat %>% 
  ggplot(aes(a, m, col = abundance)) +
  geom_point() +
  facet_wrap(~ site, scales = "free_y") +
  theme_bw() +
  scale_color_viridis_c("Abundance") +
  xlab(expression(a[CR|berr==0.02])) + ylab(expression(m[berr==0.02])),
 nrow = 2)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  group_by(site, a, b, m) %>% 
  summarise(abundance = sum(abundance), 
            agb = sum(agb)/10^3/2) %>% 
  ggplot(aes(abundance, agb, col = m)) +
  geom_point() +
  theme_bw() +
  scale_color_viridis_c(expression(m)) +
  facet_wrap(~ site)

Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  filter(a == 2) %>% 
  mutate(berr = -0.39 + 0.57*a - b) %>% 
  filter(berr > 0, berr < 0.05) %>% 
  ggplot(aes(dbh_class, abundance)) +
  geom_line(aes(col = m, group = paste(a, b, m))) +
  geom_point(data = filter(obs, dbh_class < 150), col = "red") +
  scale_y_log10(labels = scales::comma) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  scale_color_viridis_c(expression(m)) +
  scale_color_viridis_c(expression(m[a==2|berr==0.02])) +
  theme_bw() +
  facet_wrap(~ site, nrow = 2)

Best parameters

Code 
eval <- read_tsv("outputs/calib_structure.tsv") %>% 
  select(site, a, b, m) %>% 
  unique() %>% 
  mutate(dbh_class = list(seq(10, 150, by = 5))) %>% 
  unnest(dbh_class) %>% 
  left_join(read_tsv("outputs/calib_structure.tsv")) %>% 
  mutate_at(c("abundance", "agb"), ~ifelse(is.na(.), 0, .)) %>% 
  rename(simulated = abundance) %>% 
  left_join(rename(obs, observed = abundance)) %>% 
  mutate(observed = ifelse(is.na(observed), 0, observed)) %>% 
  mutate(agb_obs = ifelse(site == "Tapajos", 143.7, 209.5)) %>% 
  mutate(abund_obs = ifelse(site == "Tapajos", 470, 612)) %>% 
  group_by(site, a, b, m) %>% 
  summarise(rmsep_dist = sqrt(mean((observed - simulated)^2)),
            rmsep_agb = sqrt((unique(agb_obs) - sum(agb)/10^3/2)^2),
            rmsep_abund = sqrt((unique(abund_obs) - sum(simulated))^2)
            ) %>% 
  ungroup() %>% 
  mutate(dist_mean_obs = ifelse(site == "Tapajos", 16.2, 30.5)) %>% 
  mutate(agb_obs = ifelse(site == "Tapajos", 143.7, 209.5)) %>% 
  mutate(abund_obs = ifelse(site == "Tapajos", 470, 612)) %>% 
  mutate(rrmsep_dist =  rmsep_dist/dist_mean_obs,
         rrmsep_agb = rmsep_agb/agb_obs,
         rrmsep_abund = rmsep_abund/abund_obs) %>% 
  mutate(all_rmsep_norm = rrmsep_dist + rrmsep_agb + rrmsep_abund) %>% 
  mutate(sim = paste0(site, "_", a, "_", b, "_", m))
best <- eval %>% 
  group_by(site) %>% 
  slice_min(all_rmsep_norm, n = 5)
Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  mutate(sim = paste0(site, "_", a, "_", b, "_", m)) %>% 
  left_join(best) %>% 
  filter(!is.na(all_rmsep_norm)) %>% 
  group_by(site, a, b, m) %>% 
  summarise(abundance = sum(abundance), agb = sum(agb)/10^3/2) %>% 
  gather(variable, value, -site) %>% 
  group_by(site, variable) %>% 
  summarise(median = median(value), minimum = min(value), maximum = max(value)) %>% 
  knitr::kable(caption = "Best structure parameters according to calibration with inventory data.")
Best structure parameters according to calibration with inventory data.
site variable median minimum maximum
Paracou a 1.8500 1.8000 1.9000
Paracou abundance 608.5000 599.7500 622.7500
Paracou agb 210.6042 204.8932 216.6646
Paracou b 0.4445 0.4160 0.4730
Paracou m 0.0550 0.0500 0.0600
Tapajos a 2.5000 2.4000 2.5000
Tapajos abundance 467.2500 453.7500 485.5000
Tapajos agb 152.7362 137.3486 167.9250
Tapajos b 0.8150 0.7080 0.8150
Tapajos m 0.0550 0.0400 0.0600
Code 
read_tsv("outputs/calib_structure.tsv") %>% 
  mutate(sim = paste0(site, "_", a, "_", b, "_", m)) %>% 
  left_join(best %>% slice_min(all_rmsep_norm)) %>% 
  filter(!is.na(all_rmsep_norm)) %>% 
  ggplot(aes(dbh_class, abundance)) +
  geom_line(aes(group = sim)) +
  geom_point(data = obs, col = "red") +
  scale_y_log10(labels = scales::comma) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  theme_bw() +
  facet_wrap(~ site)

Code 
g <- read_tsv("outputs/calib_structure.tsv") %>% 
  mutate(sim = paste0(site, "_", a, "_", b, "_", m)) %>% 
  left_join(best %>% slice_min(all_rmsep_norm)) %>% 
  filter(!is.na(all_rmsep_norm)) %>% 
  ggplot(aes(dbh_class, abundance)) +
  geom_line(aes(group = sim, col = "TROLL")) +
  geom_point(data = obs, aes(col = "field")) +
  xlab("DBH [ cm ]") + ylab(expression(Abundance~"["~ha^{-~1}~"]")) +
  theme_bw() +
  facet_wrap(~ site) +
  scale_colour_manual("", values = c("#3f3f3f", "#0da1f8")) +
  theme(legend.position = "bottom") +
  geom_text(aes(col = "field", label = lab), 
            data = data.frame(dbh_class = 80, abundance = 200,
                              lab = c("470 tress/ha & 143.7 MgC/ha",
                                      "612 tress/ha & 209.5 MgC/ha"),
                              site = c("Tapajos", "Paracou"))) +
  geom_text(aes(col = "TROLL", label = lab), 
            data = data.frame(dbh_class = 80, abundance = 180,
                              lab = c("467 tress/ha & 152.7 MgC/ha",
                                      "608 tress/ha & 210.6 MgC/ha"),
                              site = c("Tapajos", "Paracou")))
ggsave("figures/f01.png", g, dpi = 300, width = 8, height = 5, bg = "white")
g