In emf-creaf/medfateland: Mediterranean Landscape Simulation
library(kableExtra)
library(sf)
library(knitr)
library(medfate)
library(ggplot2)
library(cowplot)
library(tidyverse)
evalDIR = "~/OneDrive/mcaceres_work/model_development/medfate_evaluation/RegionalLevelEvaluation"
Introduction
ver <- "2.9.3" # packageVersion("medfate")
files <- c(paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN23_v", ver, ".rds"),
paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN23_v", ver, ".rds"),
paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN34_v", ver, ".rds"),
paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN34_v", ver, ".rds"),
paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN24_v", ver, ".rds"),
paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN24_v", ver, ".rds"))
n <- length(files)
dfdbh_vec <- vector("list", n)
dfh_vec <- vector("list", n)
dfbai_sp_vec <- vector("list", n)
dfbai_plot_vec <- vector("list", n)
dfbadead_sp_vec <- vector("list", n)
dfbadead_plot_vec <- vector("list", n)
dfbarecr_sp_vec <- vector("list", n)
dfbarecr_plot_vec <- vector("list", n)
dfbatot_sp_vec <- vector("list", n)
dfbatot_plot_vec <- vector("list", n)
dfNdead_sp_vec <- vector("list", n)
dfNdead_plot_vec <- vector("list", n)
dfNrecr_sp_vec <- vector("list", n)
dfNrecr_plot_vec <- vector("list", n)
dfNtot_sp_vec <- vector("list", n)
dfNtot_plot_vec <- vector("list", n)
dfshcov_plot_vec <- vector("list", n)
dfshheight_plot_vec <- vector("list", n)
dflai_allom_plot_vec <- vector("list", n)
dflai_state_plot_vec <- vector("list", n)
IDs_23 <- NULL
IDs_34 <- NULL
IDs_24 <- NULL
for(i in 1:n) {
post <- readRDS(files[i])
if(i==1) IDs_23 <- unique(post$dfbai_plot$id)
else if(i==2) IDs_23 <- IDs_23[IDs_23 %in% post$dfbai_plot$id]
else if(i==3) IDs_34 <- unique(post$dfbai_plot$id)
else if(i==4) IDs_34 <- IDs_34[IDs_34 %in% post$dfbai_plot$id]
else if(i==5) IDs_24 <- unique(post$dfbai_plot$id)
else if(i==6) IDs_24 <- IDs_24[IDs_24 %in% post$dfbai_plot$id]
dfdbh_vec[[i]] <- post$dfdbh
dfh_vec[[i]] <- post$dfh
dfbai_sp_vec[[i]] <- post$dfbai_sp
dfbai_plot_vec[[i]] <- post$dfbai_plot
dfbadead_sp_vec[[i]] <- post$dfbadead_sp
dfbadead_plot_vec[[i]] <- post$dfbadead_plot
dfbarecr_sp_vec[[i]] <- post$dfbarecr_sp
dfbarecr_plot_vec[[i]] <- post$dfbarecr_plot
dfbatot_sp_vec[[i]] <- post$dfbatot_sp
dfbatot_plot_vec[[i]] <- post$dfbatot_plot
dfNdead_sp_vec[[i]] <- post$dfNdead_sp
dfNdead_plot_vec[[i]] <- post$dfNdead_plot
dfNrecr_sp_vec[[i]] <- post$dfNrecr_sp
dfNrecr_plot_vec[[i]] <- post$dfNrecr_plot
dfNtot_sp_vec[[i]] <- post$dfNtot_sp
dfNtot_plot_vec[[i]] <- post$dfNtot_plot
dfshcov_plot_vec[[i]] <- post$dfshcov_plot
dfshheight_plot_vec[[i]] <- post$dfshheight_plot
dflai_allom_plot_vec[[i]] <- post$dflai_allom_plot
dflai_state_plot_vec[[i]] <- post$dflai_state_plot
}
for(i in 1:n) {
if(i %in% 1:2) IDs <- IDs_23
else if(i %in% 3:4) IDs <- IDs_34
else if(i %in% 5:6) IDs <- IDs_24
dfdbh_vec[[i]]$version <- as.character(dfdbh_vec[[i]]$version)
dfh_vec[[i]]$version <- as.character(dfh_vec[[i]]$version)
dfbai_sp_vec[[i]]$version <- as.character(dfbai_sp_vec[[i]]$version)
dfbai_plot_vec[[i]]$version <- as.character(dfbai_plot_vec[[i]]$version)
dfbadead_sp_vec[[i]]$version <- as.character(dfbadead_sp_vec[[i]]$version)
dfbadead_plot_vec[[i]]$version <- as.character(dfbadead_plot_vec[[i]]$version)
dfbarecr_sp_vec[[i]]$version <- as.character(dfbarecr_sp_vec[[i]]$version)
dfbarecr_plot_vec[[i]]$version <- as.character(dfbarecr_plot_vec[[i]]$version)
dfbatot_sp_vec[[i]]$version <- as.character(dfbatot_sp_vec[[i]]$version)
dfbatot_plot_vec[[i]]$version <- as.character(dfbatot_plot_vec[[i]]$version)
dfNdead_sp_vec[[i]]$version <- as.character(dfNdead_sp_vec[[i]]$version)
dfNdead_plot_vec[[i]]$version <- as.character(dfNdead_plot_vec[[i]]$version)
dfNrecr_sp_vec[[i]]$version <- as.character(dfNrecr_sp_vec[[i]]$version)
dfNrecr_plot_vec[[i]]$version <- as.character(dfNrecr_plot_vec[[i]]$version)
dfNtot_sp_vec[[i]]$version <- as.character(dfNtot_sp_vec[[i]]$version)
dfNtot_plot_vec[[i]]$version <- as.character(dfNtot_plot_vec[[i]]$version)
dfshcov_plot_vec[[i]]$version <- as.character(dfshcov_plot_vec[[i]]$version)
dfshheight_plot_vec[[i]]$version <- as.character(dfshheight_plot_vec[[i]]$version)
dflai_allom_plot_vec[[i]]$version <- as.character(dflai_allom_plot_vec[[i]]$version)
dflai_state_plot_vec[[i]]$version <- as.character(dflai_state_plot_vec[[i]]$version)
dfdbh_vec[[i]] <- dfdbh_vec[[i]] |> filter(id %in% IDs)
dfh_vec[[i]] <- dfh_vec[[i]] |> filter(id %in% IDs)
dfbai_sp_vec[[i]] <- dfbai_sp_vec[[i]] |> filter(id %in% IDs)
dfbai_plot_vec[[i]] <- dfbai_plot_vec[[i]] |> filter(id %in% IDs)
dfbadead_sp_vec[[i]] <- dfbadead_sp_vec[[i]] |> filter(id %in% IDs)
dfbadead_plot_vec[[i]] <- dfbadead_plot_vec[[i]] |> filter(id %in% IDs)
dfbarecr_sp_vec[[i]] <- dfbarecr_sp_vec[[i]] |> filter(id %in% IDs)
dfbarecr_plot_vec[[i]] <- dfbarecr_plot_vec[[i]] |> filter(id %in% IDs)
dfbatot_sp_vec[[i]] <- dfbatot_sp_vec[[i]] |> filter(id %in% IDs)
dfbatot_plot_vec[[i]] <- dfbatot_plot_vec[[i]] |> filter(id %in% IDs)
dfNdead_sp_vec[[i]] <- dfNdead_sp_vec[[i]] |> filter(id %in% IDs)
dfNdead_plot_vec[[i]] <- dfNdead_plot_vec[[i]] |> filter(id %in% IDs)
dfNrecr_sp_vec[[i]] <- dfNrecr_sp_vec[[i]] |> filter(id %in% IDs)
dfNrecr_plot_vec[[i]] <- dfNrecr_plot_vec[[i]] |> filter(id %in% IDs)
dfNtot_sp_vec[[i]] <- dfNtot_sp_vec[[i]] |> filter(id %in% IDs)
dfNtot_plot_vec[[i]] <- dfNtot_plot_vec[[i]] |> filter(id %in% IDs)
dfshcov_plot_vec[[i]] <- dfshcov_plot_vec[[i]] |> filter(id %in% IDs)
dfshheight_plot_vec[[i]] <- dfshheight_plot_vec[[i]] |> filter(id %in% IDs)
dflai_allom_plot_vec[[i]] <- dflai_allom_plot_vec[[i]] |> filter(id %in% IDs)
dflai_state_plot_vec[[i]] <- dflai_state_plot_vec[[i]] |> filter(id %in% IDs)
}
dfdbh <- bind_rows(dfdbh_vec)
dfh <- bind_rows(dfh_vec)
dfbai_sp <- bind_rows(dfbai_sp_vec) |> sf::st_as_sf()
dfbai_plot<- bind_rows(dfbai_plot_vec) |> sf::st_as_sf()
dfbadead_sp <- bind_rows(dfbadead_sp_vec) |> sf::st_as_sf()
dfbadead_plot<- bind_rows(dfbadead_plot_vec) |> sf::st_as_sf()
dfbarecr_sp <- bind_rows(dfbarecr_sp_vec) |> sf::st_as_sf()
dfbarecr_plot<- bind_rows(dfbarecr_plot_vec) |> sf::st_as_sf()
dfbatot_sp <- bind_rows(dfbatot_sp_vec) |> sf::st_as_sf()
dfbatot_plot<- bind_rows(dfbatot_plot_vec) |> sf::st_as_sf()
dfNdead_sp <- bind_rows(dfNdead_sp_vec) |> sf::st_as_sf()
dfNdead_plot<- bind_rows(dfNdead_plot_vec) |> sf::st_as_sf()
dfNrecr_sp <- bind_rows(dfNrecr_sp_vec) |> sf::st_as_sf()
dfNrecr_plot<- bind_rows(dfNrecr_plot_vec) |> sf::st_as_sf()
dfNtot_sp <- bind_rows(dfNtot_sp_vec) |> sf::st_as_sf()
dfNtot_plot<- bind_rows(dfNtot_plot_vec) |> sf::st_as_sf()
dfshcov_plot<- bind_rows(dfshcov_plot_vec) |> sf::st_as_sf()
dfshheight_plot<- bind_rows(dfshheight_plot_vec) |> sf::st_as_sf()
dflai_allom_plot<- bind_rows(dflai_allom_plot_vec) |> sf::st_as_sf()
dflai_state_plot<- bind_rows(dflai_state_plot_vec) |> sf::st_as_sf()
load("~/OneDrive/EMF_datasets/PoliticalBoundaries/Products/Catalunya/catalonia.rdata")
cat_sfc = sf::st_as_sfc(cat.contour)
evaluation_stats<-function(df) {
df$version <- as.character(df$version)
df$period <- factor(df$period, levels = c("IFN23", "IFN34", "IFN24"))
df_stats <- sf::st_drop_geometry(df) |>
group_by(version, period, transpirationMode) |>
summarize(n = n(),
Obs = mean(obs, na.rm=TRUE),
Pred = mean(pred, na.rm=TRUE),
Bias = mean(error, na.rm=TRUE),
Biasrel = 100*Bias/abs(Obs),
# MAE = mean(abs(error), na.rm=TRUE),
# MAErel = 100*MAE/abs(Obs),
RMSE = sqrt(mean(error^2, na.rm=TRUE)),
RMSErel = 100*RMSE/abs(Obs),
R2 = cor(obs,pred, use = "complete.obs")^2, .groups="drop")
return(as.data.frame(df_stats))
}
Goal
The aim of this article is to provide an assessment of the performance of fordyn (using either the basic and advanced sub-model) for the prediction of forest dynamics in Catalonia (NE of Spain). To this aim, we simulate forest dynamics between surveys of the Spanish National Forest Inventory and compare the model predictions of forest development against inventory data for a set of permanent plots. The evaluation focuses first on the growth (in diameter and height) of surviving trees, then turning the attention to the basal area of dead trees and overall changes in terms of basal area and density. Then, we evaluate changes in stand leaf area index and, finally, changes in shrub cover and mean shrub height.
plot_scatter_tree<-function(df_bas, df_adv, points = TRUE,
xylim = c(0,1.5), errorlim = c(-1,1),
initiallim = c(0,50), var = "diameter"){
if(var=="diameter") {
df_bas$ini <- df_bas$DBH_ini
df_adv$ini <- df_adv$DBH_ini
} else {
df_bas$ini <- df_bas$H_ini
df_adv$ini <- df_adv$H_ini
}
df<- bind_rows(df_bas, df_adv)
g1<-ggplot(df, aes(x = pred, y= obs))
if(points) {
g1 <- g1 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
}
g1 <- g1 + geom_abline(intercept = 0, slope = 1, col ="black")+
xlab(paste0("Predicted ", var, " increment (cm/yr)"))+
ylab(paste0("Observed ", var, " increment (cm/yr)"))+
xlim(xylim)+ ylim(xylim)+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = "none")
g2<-ggplot(df, aes(x = ini, y= error))
if(points) {
g2 <- g2 + geom_point(aes(col=transpirationMode), size=0.1, alpha = 0.05)
}
g2 <- g2 + geom_abline(intercept = 0, slope = 0, col ="black")+
geom_smooth(aes(col=transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
xlab(paste0("Initial ", var, " (cm)"))+
ylab(paste0("Error ", var, " increment (cm/yr)"))+
ylim(errorlim)+xlim(initiallim)+
scale_fill_discrete("")+
scale_color_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = c(0.80,0.85))
plot_grid(g1, g2,nrow=1)
}
plot_scatter_bai<-function(dfbai_bas, dfbai_adv, quantity = "basal area increment", xylim = c(0,2), errorlim = c(-1,1)) {
df<- bind_rows(dfbai_bas, dfbai_adv)
g1<-ggplot(df, aes(x = pred, y= obs))+
geom_point(aes(col=transpirationMode), size=0.3, alpha = 0.3)+
geom_abline(intercept = 0, slope = 1, col ="black")+
xlab(paste0("Predicted ", quantity," (m2/ha/yr)"))+
ylab(paste0("Observed ", quantity," (m2/ha/yr)"))+
xlim(xylim)+ ylim(xylim)+
theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none")
g2<-ggplot(df, aes(x = BA_ini, y= error))+
geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3, data = dfbai_bas)+
geom_abline(intercept = 0, slope = 0, col ="black")+
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+
xlab("Initial basal area (m2/ha)")+
ylab(paste0("Error ", quantity," (m2/ha/yr)"))+
xlim(c(0,30))+ylim(errorlim)+
scale_fill_discrete("")+
scale_color_discrete("")+
theme_bw()+ theme(panel.grid = element_blank(), legend.position = c(0.8,0.85))
plot_grid(g1, g2, nrow=1)
}
plot_cov_clim_tree<-function(df_bas, df_adv, points = TRUE, xylim = c(0,1), errorlim = c(-1,1), var = "diameter"){
df_bas$PPET[df_bas$PPET>2] = 2
df_adv$PPET[df_adv$PPET>2] = 2
df <- bind_rows(df_bas, df_adv)
g1<-ggplot(df, aes(x = PPET, y=pred))
if(points) {
g1 <- g1 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
g1 <- g1 + geom_point(aes(x = PPET, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas)
}
g1 <- g1 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_smooth(aes(x = PPET, y = obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+
ylab(paste0("Annual ", var, " increment (cm/yr)"))+ ylim(xylim)+
xlab("Moisture index (P/PET)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(legend.position = "none",
legend.background = element_blank(),
panel.grid = element_blank())
g2<-ggplot(df, aes(x = PPET, y=error))
if(points) {
g2 <- g2 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
}
g2 <- g2 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_abline(intercept = 0, slope = 0, col ="black")+
ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+
xlab("Moisture index (P/PET)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = c(0.8,0.85))
g3<-ggplot(df, aes(x = MAT, y=pred))
if(points) {
g3 <- g3 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
g3 <- g3 + geom_point(aes(x = MAT, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas)
}
g3 <- g3 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_smooth(aes(x = MAT, y= obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+
ylab(paste0("Annual ", var, " increment (cm/yr)"))+ylim(xylim)+
xlab("Mean annual temperature (ºC)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(legend.position = "none",
legend.background = element_blank(),
panel.grid = element_blank())
g4<-ggplot(df, aes(x = MAT, y=error))
if(points) {
g4 <- g4 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
}
g4 <- g4 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_abline(intercept = 0, slope = 0, col ="black")+
ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+
xlab("Mean annual temperature (ºC)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = "none")
g5<-ggplot(df, aes(x = light_ini, y=pred))
if(points) {
g5 <- g5 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
g5 <- g5 + geom_point(aes(x = light_ini, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas)
}
g5 <- g5 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_smooth(aes(x = light_ini, y= obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+
ylab(paste0("Annual ", var," increment (cm/yr)"))+ylim(xylim)+
xlab("Available PAR (%)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(legend.position = "none",
legend.background = element_blank(),
panel.grid = element_blank())
g6<-ggplot(df, aes(x = light_ini, y=error))
if(points) {
g6 <- g6 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05)
}
g6 <- g6 +
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+
geom_abline(intercept = 0, slope = 0, col ="black")+
ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+
xlab("Available PAR (%)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = "none")
plot_grid(g1, g2, g3, g4, g5,g6, nrow=3)
}
plot_cov_clim_bai<-function(df_bas, df_adv, quantity = "Basal area increment", ylim = c(0,2), errorlim = c(-2,2)){
df_bas$PPET[df_bas$PPET>2] = 2
df_adv$PPET[df_adv$PPET>2] = 2
df <- bind_rows(df_bas, df_adv)
g1<-ggplot(df, aes(x = PPET, y=pred))+
geom_point(aes(x = PPET, y= obs), size=0.3, col="gray", alpha = 0.3, data = df_bas)+
geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+
geom_smooth(aes(x = PPET, y = obs), col = "black", method="loess", formula = y ~ x, data = df_bas)+
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+
ylab(paste0(quantity," (m2/ha/yr)"))+ ylim(ylim)+
xlab("Moisture index (P/PET)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(legend.position = "none",
legend.background = element_blank(),
panel.grid = element_blank())
g2<-ggplot(df, aes(x = PPET, y=error))+
geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+
geom_smooth(aes(col = transpirationMode, fill= transpirationMode), method="loess", formula = y ~ x)+
geom_abline(intercept = 0, slope = 0, col ="black")+
ylab(paste0(quantity," error (m2/ha/yr)"))+ylim(errorlim)+
xlab("Moisture index (P/PET)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = c(0.8,0.85))
g3<-ggplot(df, aes(x = MAT, y=pred))+
geom_point(aes(x = MAT, y= obs), size=0.3, col="gray", alpha = 0.3, data = df_bas)+
geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+
geom_smooth(aes(x = MAT, y = obs), col = "black", method="loess", formula = y ~ x, data = df_bas)+
geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+
ylab(paste0(quantity," (m2/ha/yr)"))+ ylim(ylim)+
xlab("Mean annual temperature (ºC)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(legend.position = "none",
legend.background = element_blank(),
panel.grid = element_blank())
g4<-ggplot(df, aes(x = MAT, y=error))+
geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+
geom_smooth(aes(col = transpirationMode, fill= transpirationMode), method="loess", formula = y ~ x)+
geom_abline(intercept = 0, slope = 0, col ="black")+
ylab(paste0(quantity," error (m2/ha/yr)"))+ ylim(errorlim)+
xlab("Mean annual temperature (ºC)")+
scale_color_discrete("")+ scale_fill_discrete("")+
theme_bw()+
theme(panel.grid = element_blank(), legend.position = "none")
plot_grid(g1, g2, g3, g4, nrow=2)
}
bai_error_map <-function(data_bai, upper = 1, lower = -1){
data_bai = data_bai[!is.na(data_bai$error),]
data_bai$error[data_bai$error > upper] = upper
data_bai$error[data_bai$error < lower] = lower
ggplot(data_bai)+
geom_sf(aes(col=error), size=0.5)+
geom_sf(data = cat_sfc)+
scale_color_fermenter("m2/ha/yr",palette = 1, breaks = c(lower,lower*0.5,lower*0.25,lower*0.1, 0, upper*0.1,upper*0.25,upper*0.5,upper),
type = 'div', direction = 2, na.value = 'transparent')+
# scale_color_gradient2("m2/ha",low="blue", mid="yellow", high="red")+
theme_bw()+
theme(panel.grid = element_blank())
# theme(panel.grid = element_blank(), legend.position = c(0.85,0.25))
}
Simulation procedure
We selected permanent plots between the second (IFN2) and the fourth (IFN4) without signs of management (i.e. the presence of stumps) and whose basal area did not decrease more than 10% between inventory surveys (to avoid the effect of disturbances).
Soil physical properties were drawn from SoilGrids (Hengl 2016), complemented by rock fragment content estimates derived from surface stoniness measurements in forest plots.
fordyn simulations were conducted for different periods:
- IFN2 - IFN3 (~ 10 years)
- IFN3 - IFN4 (~ 15 years)
- IFN2 - IFN4 (~ 25 years)
The actual simulated period varied depending on the sampling years of the target plot. Daily weather data were obtained via interpolation on plot's coordinates using package meteoland.
Default species-specific parameters were modified using the results of the meta-modelling exercise and the growth calibration exercise. These two exercises do not provide values for all the main species included here, so it is expected that evaluation results are worse for those species not included in those exercises.
Simulations were done for both the basic (i.e. Granier) and advanced (i.e. Sperry) transpiration/photosynthesis sub-models. On a server with 20 parallel threads, computational times for the longest simulation period (25 years) are around 4 hours (i.e. 2.5 min/plot) for the basic sub-model, versus around 6 days (i.e. 1.5 hr/plot) for the advanced submodel.
In the following sections, we provide the bias, root mean squared error (in absolute and relative terms) and explained variance (R-squared) of growth and mortality predictions at the tree-level and stand-level obtained by simulations with both sub-models. Scatter plots are provided for the IFN2-IFN4 simulation to represent the relationship between predicted and observed values, as well as the factors that may influence the direction and magnitude of prediction error (i.e. initial values, environmental conditions, ...).
Detailed results of growth evaluation by tree species are provided in the last section.
Growth of surviving trees
Comparison of diameter/height growth of trees (DBH >= 7.5) that survived between surveys.
Annual diameter increment
Overall predictive capacity to predict diameter increase (cm/yr):
dfdbh |>
evaluation_stats() |>
kbl() |>
kable_styling()
Predictive capacity plots (IFN2-IFN4):
dfdbh_bas <- dfdbh |> filter(transpirationMode == "Granier",
period == "IFN24")
dfdbh_adv <- dfdbh |> filter(transpirationMode == "Sperry",
period == "IFN24")
plot_scatter_tree(dfdbh_bas, dfdbh_adv, points = TRUE)
Relationship between diameter increase and climatic variables (MAT, P/PET and available PAR; IFN2 - IFN4):
plot_cov_clim_tree(dfdbh_bas, dfdbh_adv, points = TRUE)
Annual height increment
Overall predictive capacity to predict height increase (cm/yr):
dfh %>%
evaluation_stats() %>%
kbl() %>%
kable_styling()
Predictive capacity plots (IFN2-IFN4):
dfh_bas <- dfh |> filter(transpirationMode == "Granier",
period == "IFN24")
dfh_adv <- dfh |> filter(transpirationMode == "Sperry",
period == "IFN24")
plot_scatter_tree(dfh_bas, dfh_adv, points = TRUE, xylim = c(0,20), errorlim = c(-20,20), initiallim = c(0, 2000),
var = "height")
Relationship between height increase and climatic variables (MAT, P/PET and available PAR; IFN2 - IFN4):
plot_cov_clim_tree(dfh_bas, dfh_adv, points = TRUE, xylim = c(0,20), errorlim = c(-20,20), var = "height")
Stand-level basal area increment
Comparison of basal area increment of surviving trees does not take into account changes in density. In other words, densities from the first inventory are used to calculate stand-level basal area of surviving trees. Hence, the comparison is meant to evaluate the effect of diameter increment of surviving trees in terms of annual stand basal area increments (m2/ha/yr) for the period evaluated.
Predictive capacity table:
dfbai_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfbai_bas <- dfbai_plot |> filter(transpirationMode == "Granier",
period == "IFN24")
dfbai_adv <- dfbai_plot |> filter(transpirationMode == "Sperry",
period == "IFN24")
plot_scatter_bai(dfbai_bas, dfbai_adv)
Relationship between basal area increase and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfbai_bas, dfbai_adv, errorlim = c(-1,1), ylim = c(0,1))
Spatial error distribution (IFN2 - IFN4):
p1<-bai_error_map(dfbai_bas)+labs(title="Basic")
p2<-bai_error_map(dfbai_adv)+labs(title="Advanced")
plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"),
get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Mortality
Basal area reduction
Annual reduction of basal area (m2/ha/yr) due to trees (DBH >= 7.5) that died during the evaluation period against model's mortality prediction. In both cases, basal area is calculated using the initial diameter of the trees, so that density reductions are the only prediction that is actually evaluated (and not the possible growth of those trees during the simulation).
Overall predictive capacity:
dfbadead_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfdead_bas <- dfbadead_plot |> filter(transpirationMode == "Granier",
period == "IFN24")
dfdead_adv <- dfbadead_plot |> filter(transpirationMode == "Sperry",
period == "IFN24")
plot_scatter_bai(dfdead_bas, dfdead_adv,
quantity = "dead basal area", xylim = c(0,0.5), errorlim = c(-0.5,0.5))
Relationship between dead basal area and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfdead_bas, dfdead_adv,
quantity = "Dead basal area", ylim = c(0,0.5), errorlim = c(-0.5,0.5))
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dfdead_bas, upper = 0.5, lower = -0.5)+labs(title="Basic")
p2<-bai_error_map(dfdead_bas, upper = 0.5, lower = -0.5)+labs(title="Advanced")
plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"),
get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Density reduction
Annual reduction of density (ind/ha/yr) due to trees (DBH >= 7.5) that died during the evaluation period against model's mortality prediction. This is very similar to evaluating the reduction in basal area
Overall predictive capacity:
dfNdead_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Ingrowth
Basal area increase
Annual increase of basal area (m2/ha/yr) due to ingrowth of trees with diameters between 7.5 cm and 12.5 cm during the evaluated period.
Predictive capacity:
dfbarecr_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfrecr_bas <- dfbarecr_plot |> filter(transpirationMode == "Granier",
period == "IFN24")
dfrecr_adv <- dfbarecr_plot |> filter(transpirationMode == "Sperry",
period == "IFN24")
plot_scatter_bai(dfrecr_bas, dfrecr_adv,
quantity = "ingrowth basal area", xylim = c(0,0.15), errorlim = c(-0.15,0.15))
Relationship between ingrowth basal area and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfrecr_bas, dfrecr_adv,
quantity = "Ingrowth basal area", ylim = c(0,0.15), errorlim = c(-0.15,0.15))
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dfrecr_bas, upper = 0.5, lower = -0.5)+labs(title="Basic")
p2<-bai_error_map(dfrecr_bas, upper = 0.5, lower = -0.5)+labs(title="Advanced")
plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"),
get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Density increase
Annual increase of density (ind/ha/yr) due to ingrowth of trees with diameters between 7.5 cm and 12.5 cm during the evaluated period.
Predictive capacity:
dfNrecr_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Overall stand-level change
Basal area changes
This includes annual changes in basal area (m2/ha/yr) due to growth of surviving trees, mortality reductions and ingrowth derived from sapling growth. In the observed data, basal area changes include also the incorporation of trees into large diameter classes that results from the variable-radius sampling design. Since it takes into account all processes together, this evaluation is the most rellevant of all.
Overall predictive capacity:
dfbatot_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Predictive capacity plot (IFN2 - IFN4):
dftot_bas <- dfbatot_plot |> filter(transpirationMode == "Granier",
period == "IFN24")
dftot_adv <- dfbatot_plot |> filter(transpirationMode == "Sperry",
period == "IFN24")
p1<-plot_scatter_bai(dftot_bas, dftot_adv,
quantity = "basal area change", xylim = c(-0.5,2), errorlim = c(-2,2))
p1
Relationship between overall basal area change and climatic variables (MAT and P/PET; IFN2 - IFN4):
p2<-plot_cov_clim_bai(dftot_bas, dftot_adv,
quantity = "Basal area change", ylim = c(-0.5,2), errorlim = c(-2,2))
p2
p0<-cowplot::plot_grid(p1,p2, nrow=2, rel_heights = c(0.33,0.66))
ggsave(paste0(evalDIR,"/plots/BAtot_eval.png"),p0, width = 8, height=12)
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dftot_bas)+labs(title="Basic sub-model")
p2<-bai_error_map(dftot_adv)+labs(title="Advanced sub-model")
pm<-plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"),
get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
ggsave(paste0(evalDIR,"/plots/BAtot_error_map.png"),pm, width = 13, height=6, bg="white")
pm
Density changes
This includes annual changes in density (ind/ha/yr) due to growth of surviving trees, mortality reductions and ingrowth derived from sapling growth. In the observed data, changes include also the incorporation of trees into large diameter classes that results from the variable-radius sampling design.
Overall predictive capacity:
dfNtot_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Changes in leaf area index of trees > 7.5 cm
Overall predictive capacity using allometric equations:
dflai_allom_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Overall predictive capacity using state variables:
dflai_state_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Shrub cover and mean height
Shrub percent cover changes
Overall predictive capacity:
dfshcov_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Shrub mean height changes
Overall predictive capacity:
dfshheight_plot %>% evaluation_stats() %>%
kbl() %>%
kable_styling()
Detailed evaluation by tree species
In the following we provide detailed evaluation results for the most important species.
sp100 = c("Abies alba", "Fagus sylvatica", "Pinus halepensis", "Pinus nigra",
"Pinus sylvestris", "Pinus uncinata", "Quercus faginea", "Quercus ilex",
"Quercus pubescens", "Quercus suber")
for(i in 1:length(sp100)) {
sp_name = sp100[i]
res <- knitr::knit_child('_RLEvaluation_child.Rmd', quiet = TRUE)
cat(res, sep = '\n')
}
emf-creaf/medfateland documentation built on April 17, 2025, 5:43 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(kableExtra) library(sf) library(knitr) library(medfate) library(ggplot2) library(cowplot) library(tidyverse) evalDIR = "~/OneDrive/mcaceres_work/model_development/medfate_evaluation/RegionalLevelEvaluation"
Introduction
ver <- "2.9.3" # packageVersion("medfate") files <- c(paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN23_v", ver, ".rds"), paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN23_v", ver, ".rds"), paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN34_v", ver, ".rds"), paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN34_v", ver, ".rds"), paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_granier_IFN24_v", ver, ".rds"), paste0(evalDIR, "/data/output/Catalunya/", ver, "/IFNeval_results_sperry_IFN24_v", ver, ".rds")) n <- length(files) dfdbh_vec <- vector("list", n) dfh_vec <- vector("list", n) dfbai_sp_vec <- vector("list", n) dfbai_plot_vec <- vector("list", n) dfbadead_sp_vec <- vector("list", n) dfbadead_plot_vec <- vector("list", n) dfbarecr_sp_vec <- vector("list", n) dfbarecr_plot_vec <- vector("list", n) dfbatot_sp_vec <- vector("list", n) dfbatot_plot_vec <- vector("list", n) dfNdead_sp_vec <- vector("list", n) dfNdead_plot_vec <- vector("list", n) dfNrecr_sp_vec <- vector("list", n) dfNrecr_plot_vec <- vector("list", n) dfNtot_sp_vec <- vector("list", n) dfNtot_plot_vec <- vector("list", n) dfshcov_plot_vec <- vector("list", n) dfshheight_plot_vec <- vector("list", n) dflai_allom_plot_vec <- vector("list", n) dflai_state_plot_vec <- vector("list", n) IDs_23 <- NULL IDs_34 <- NULL IDs_24 <- NULL for(i in 1:n) { post <- readRDS(files[i]) if(i==1) IDs_23 <- unique(post$dfbai_plot$id) else if(i==2) IDs_23 <- IDs_23[IDs_23 %in% post$dfbai_plot$id] else if(i==3) IDs_34 <- unique(post$dfbai_plot$id) else if(i==4) IDs_34 <- IDs_34[IDs_34 %in% post$dfbai_plot$id] else if(i==5) IDs_24 <- unique(post$dfbai_plot$id) else if(i==6) IDs_24 <- IDs_24[IDs_24 %in% post$dfbai_plot$id] dfdbh_vec[[i]] <- post$dfdbh dfh_vec[[i]] <- post$dfh dfbai_sp_vec[[i]] <- post$dfbai_sp dfbai_plot_vec[[i]] <- post$dfbai_plot dfbadead_sp_vec[[i]] <- post$dfbadead_sp dfbadead_plot_vec[[i]] <- post$dfbadead_plot dfbarecr_sp_vec[[i]] <- post$dfbarecr_sp dfbarecr_plot_vec[[i]] <- post$dfbarecr_plot dfbatot_sp_vec[[i]] <- post$dfbatot_sp dfbatot_plot_vec[[i]] <- post$dfbatot_plot dfNdead_sp_vec[[i]] <- post$dfNdead_sp dfNdead_plot_vec[[i]] <- post$dfNdead_plot dfNrecr_sp_vec[[i]] <- post$dfNrecr_sp dfNrecr_plot_vec[[i]] <- post$dfNrecr_plot dfNtot_sp_vec[[i]] <- post$dfNtot_sp dfNtot_plot_vec[[i]] <- post$dfNtot_plot dfshcov_plot_vec[[i]] <- post$dfshcov_plot dfshheight_plot_vec[[i]] <- post$dfshheight_plot dflai_allom_plot_vec[[i]] <- post$dflai_allom_plot dflai_state_plot_vec[[i]] <- post$dflai_state_plot } for(i in 1:n) { if(i %in% 1:2) IDs <- IDs_23 else if(i %in% 3:4) IDs <- IDs_34 else if(i %in% 5:6) IDs <- IDs_24 dfdbh_vec[[i]]$version <- as.character(dfdbh_vec[[i]]$version) dfh_vec[[i]]$version <- as.character(dfh_vec[[i]]$version) dfbai_sp_vec[[i]]$version <- as.character(dfbai_sp_vec[[i]]$version) dfbai_plot_vec[[i]]$version <- as.character(dfbai_plot_vec[[i]]$version) dfbadead_sp_vec[[i]]$version <- as.character(dfbadead_sp_vec[[i]]$version) dfbadead_plot_vec[[i]]$version <- as.character(dfbadead_plot_vec[[i]]$version) dfbarecr_sp_vec[[i]]$version <- as.character(dfbarecr_sp_vec[[i]]$version) dfbarecr_plot_vec[[i]]$version <- as.character(dfbarecr_plot_vec[[i]]$version) dfbatot_sp_vec[[i]]$version <- as.character(dfbatot_sp_vec[[i]]$version) dfbatot_plot_vec[[i]]$version <- as.character(dfbatot_plot_vec[[i]]$version) dfNdead_sp_vec[[i]]$version <- as.character(dfNdead_sp_vec[[i]]$version) dfNdead_plot_vec[[i]]$version <- as.character(dfNdead_plot_vec[[i]]$version) dfNrecr_sp_vec[[i]]$version <- as.character(dfNrecr_sp_vec[[i]]$version) dfNrecr_plot_vec[[i]]$version <- as.character(dfNrecr_plot_vec[[i]]$version) dfNtot_sp_vec[[i]]$version <- as.character(dfNtot_sp_vec[[i]]$version) dfNtot_plot_vec[[i]]$version <- as.character(dfNtot_plot_vec[[i]]$version) dfshcov_plot_vec[[i]]$version <- as.character(dfshcov_plot_vec[[i]]$version) dfshheight_plot_vec[[i]]$version <- as.character(dfshheight_plot_vec[[i]]$version) dflai_allom_plot_vec[[i]]$version <- as.character(dflai_allom_plot_vec[[i]]$version) dflai_state_plot_vec[[i]]$version <- as.character(dflai_state_plot_vec[[i]]$version) dfdbh_vec[[i]] <- dfdbh_vec[[i]] |> filter(id %in% IDs) dfh_vec[[i]] <- dfh_vec[[i]] |> filter(id %in% IDs) dfbai_sp_vec[[i]] <- dfbai_sp_vec[[i]] |> filter(id %in% IDs) dfbai_plot_vec[[i]] <- dfbai_plot_vec[[i]] |> filter(id %in% IDs) dfbadead_sp_vec[[i]] <- dfbadead_sp_vec[[i]] |> filter(id %in% IDs) dfbadead_plot_vec[[i]] <- dfbadead_plot_vec[[i]] |> filter(id %in% IDs) dfbarecr_sp_vec[[i]] <- dfbarecr_sp_vec[[i]] |> filter(id %in% IDs) dfbarecr_plot_vec[[i]] <- dfbarecr_plot_vec[[i]] |> filter(id %in% IDs) dfbatot_sp_vec[[i]] <- dfbatot_sp_vec[[i]] |> filter(id %in% IDs) dfbatot_plot_vec[[i]] <- dfbatot_plot_vec[[i]] |> filter(id %in% IDs) dfNdead_sp_vec[[i]] <- dfNdead_sp_vec[[i]] |> filter(id %in% IDs) dfNdead_plot_vec[[i]] <- dfNdead_plot_vec[[i]] |> filter(id %in% IDs) dfNrecr_sp_vec[[i]] <- dfNrecr_sp_vec[[i]] |> filter(id %in% IDs) dfNrecr_plot_vec[[i]] <- dfNrecr_plot_vec[[i]] |> filter(id %in% IDs) dfNtot_sp_vec[[i]] <- dfNtot_sp_vec[[i]] |> filter(id %in% IDs) dfNtot_plot_vec[[i]] <- dfNtot_plot_vec[[i]] |> filter(id %in% IDs) dfshcov_plot_vec[[i]] <- dfshcov_plot_vec[[i]] |> filter(id %in% IDs) dfshheight_plot_vec[[i]] <- dfshheight_plot_vec[[i]] |> filter(id %in% IDs) dflai_allom_plot_vec[[i]] <- dflai_allom_plot_vec[[i]] |> filter(id %in% IDs) dflai_state_plot_vec[[i]] <- dflai_state_plot_vec[[i]] |> filter(id %in% IDs) } dfdbh <- bind_rows(dfdbh_vec) dfh <- bind_rows(dfh_vec) dfbai_sp <- bind_rows(dfbai_sp_vec) |> sf::st_as_sf() dfbai_plot<- bind_rows(dfbai_plot_vec) |> sf::st_as_sf() dfbadead_sp <- bind_rows(dfbadead_sp_vec) |> sf::st_as_sf() dfbadead_plot<- bind_rows(dfbadead_plot_vec) |> sf::st_as_sf() dfbarecr_sp <- bind_rows(dfbarecr_sp_vec) |> sf::st_as_sf() dfbarecr_plot<- bind_rows(dfbarecr_plot_vec) |> sf::st_as_sf() dfbatot_sp <- bind_rows(dfbatot_sp_vec) |> sf::st_as_sf() dfbatot_plot<- bind_rows(dfbatot_plot_vec) |> sf::st_as_sf() dfNdead_sp <- bind_rows(dfNdead_sp_vec) |> sf::st_as_sf() dfNdead_plot<- bind_rows(dfNdead_plot_vec) |> sf::st_as_sf() dfNrecr_sp <- bind_rows(dfNrecr_sp_vec) |> sf::st_as_sf() dfNrecr_plot<- bind_rows(dfNrecr_plot_vec) |> sf::st_as_sf() dfNtot_sp <- bind_rows(dfNtot_sp_vec) |> sf::st_as_sf() dfNtot_plot<- bind_rows(dfNtot_plot_vec) |> sf::st_as_sf() dfshcov_plot<- bind_rows(dfshcov_plot_vec) |> sf::st_as_sf() dfshheight_plot<- bind_rows(dfshheight_plot_vec) |> sf::st_as_sf() dflai_allom_plot<- bind_rows(dflai_allom_plot_vec) |> sf::st_as_sf() dflai_state_plot<- bind_rows(dflai_state_plot_vec) |> sf::st_as_sf() load("~/OneDrive/EMF_datasets/PoliticalBoundaries/Products/Catalunya/catalonia.rdata") cat_sfc = sf::st_as_sfc(cat.contour)
evaluation_stats<-function(df) { df$version <- as.character(df$version) df$period <- factor(df$period, levels = c("IFN23", "IFN34", "IFN24")) df_stats <- sf::st_drop_geometry(df) |> group_by(version, period, transpirationMode) |> summarize(n = n(), Obs = mean(obs, na.rm=TRUE), Pred = mean(pred, na.rm=TRUE), Bias = mean(error, na.rm=TRUE), Biasrel = 100*Bias/abs(Obs), # MAE = mean(abs(error), na.rm=TRUE), # MAErel = 100*MAE/abs(Obs), RMSE = sqrt(mean(error^2, na.rm=TRUE)), RMSErel = 100*RMSE/abs(Obs), R2 = cor(obs,pred, use = "complete.obs")^2, .groups="drop") return(as.data.frame(df_stats)) }
Goal
The aim of this article is to provide an assessment of the performance of fordyn (using either the basic and advanced sub-model) for the prediction of forest dynamics in Catalonia (NE of Spain). To this aim, we simulate forest dynamics between surveys of the Spanish National Forest Inventory and compare the model predictions of forest development against inventory data for a set of permanent plots. The evaluation focuses first on the growth (in diameter and height) of surviving trees, then turning the attention to the basal area of dead trees and overall changes in terms of basal area and density. Then, we evaluate changes in stand leaf area index and, finally, changes in shrub cover and mean shrub height.
plot_scatter_tree<-function(df_bas, df_adv, points = TRUE, xylim = c(0,1.5), errorlim = c(-1,1), initiallim = c(0,50), var = "diameter"){ if(var=="diameter") { df_bas$ini <- df_bas$DBH_ini df_adv$ini <- df_adv$DBH_ini } else { df_bas$ini <- df_bas$H_ini df_adv$ini <- df_adv$H_ini } df<- bind_rows(df_bas, df_adv) g1<-ggplot(df, aes(x = pred, y= obs)) if(points) { g1 <- g1 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) } g1 <- g1 + geom_abline(intercept = 0, slope = 1, col ="black")+ xlab(paste0("Predicted ", var, " increment (cm/yr)"))+ ylab(paste0("Observed ", var, " increment (cm/yr)"))+ xlim(xylim)+ ylim(xylim)+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none") g2<-ggplot(df, aes(x = ini, y= error)) if(points) { g2 <- g2 + geom_point(aes(col=transpirationMode), size=0.1, alpha = 0.05) } g2 <- g2 + geom_abline(intercept = 0, slope = 0, col ="black")+ geom_smooth(aes(col=transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ xlab(paste0("Initial ", var, " (cm)"))+ ylab(paste0("Error ", var, " increment (cm/yr)"))+ ylim(errorlim)+xlim(initiallim)+ scale_fill_discrete("")+ scale_color_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = c(0.80,0.85)) plot_grid(g1, g2,nrow=1) }
plot_scatter_bai<-function(dfbai_bas, dfbai_adv, quantity = "basal area increment", xylim = c(0,2), errorlim = c(-1,1)) { df<- bind_rows(dfbai_bas, dfbai_adv) g1<-ggplot(df, aes(x = pred, y= obs))+ geom_point(aes(col=transpirationMode), size=0.3, alpha = 0.3)+ geom_abline(intercept = 0, slope = 1, col ="black")+ xlab(paste0("Predicted ", quantity," (m2/ha/yr)"))+ ylab(paste0("Observed ", quantity," (m2/ha/yr)"))+ xlim(xylim)+ ylim(xylim)+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none") g2<-ggplot(df, aes(x = BA_ini, y= error))+ geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3, data = dfbai_bas)+ geom_abline(intercept = 0, slope = 0, col ="black")+ geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+ xlab("Initial basal area (m2/ha)")+ ylab(paste0("Error ", quantity," (m2/ha/yr)"))+ xlim(c(0,30))+ylim(errorlim)+ scale_fill_discrete("")+ scale_color_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = c(0.8,0.85)) plot_grid(g1, g2, nrow=1) }
plot_cov_clim_tree<-function(df_bas, df_adv, points = TRUE, xylim = c(0,1), errorlim = c(-1,1), var = "diameter"){ df_bas$PPET[df_bas$PPET>2] = 2 df_adv$PPET[df_adv$PPET>2] = 2 df <- bind_rows(df_bas, df_adv) g1<-ggplot(df, aes(x = PPET, y=pred)) if(points) { g1 <- g1 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) g1 <- g1 + geom_point(aes(x = PPET, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas) } g1 <- g1 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_smooth(aes(x = PPET, y = obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+ ylab(paste0("Annual ", var, " increment (cm/yr)"))+ ylim(xylim)+ xlab("Moisture index (P/PET)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(legend.position = "none", legend.background = element_blank(), panel.grid = element_blank()) g2<-ggplot(df, aes(x = PPET, y=error)) if(points) { g2 <- g2 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) } g2 <- g2 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_abline(intercept = 0, slope = 0, col ="black")+ ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+ xlab("Moisture index (P/PET)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = c(0.8,0.85)) g3<-ggplot(df, aes(x = MAT, y=pred)) if(points) { g3 <- g3 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) g3 <- g3 + geom_point(aes(x = MAT, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas) } g3 <- g3 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_smooth(aes(x = MAT, y= obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+ ylab(paste0("Annual ", var, " increment (cm/yr)"))+ylim(xylim)+ xlab("Mean annual temperature (ºC)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(legend.position = "none", legend.background = element_blank(), panel.grid = element_blank()) g4<-ggplot(df, aes(x = MAT, y=error)) if(points) { g4 <- g4 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) } g4 <- g4 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_abline(intercept = 0, slope = 0, col ="black")+ ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+ xlab("Mean annual temperature (ºC)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none") g5<-ggplot(df, aes(x = light_ini, y=pred)) if(points) { g5 <- g5 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) g5 <- g5 + geom_point(aes(x = light_ini, y= obs), size=0.1, col="gray", alpha = 0.05, data = df_bas) } g5 <- g5 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_smooth(aes(x = light_ini, y= obs), col = "black", method="gam", formula = y ~ s(x, bs = "cs"), data = df_bas)+ ylab(paste0("Annual ", var," increment (cm/yr)"))+ylim(xylim)+ xlab("Available PAR (%)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(legend.position = "none", legend.background = element_blank(), panel.grid = element_blank()) g6<-ggplot(df, aes(x = light_ini, y=error)) if(points) { g6 <- g6 + geom_point(aes(col = transpirationMode), size=0.1, alpha = 0.05) } g6 <- g6 + geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="gam", formula = y ~ s(x, bs = "cs"))+ geom_abline(intercept = 0, slope = 0, col ="black")+ ylab(paste0("Annual ", var, " increment error (cm)"))+ylim(errorlim)+ xlab("Available PAR (%)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none") plot_grid(g1, g2, g3, g4, g5,g6, nrow=3) }
plot_cov_clim_bai<-function(df_bas, df_adv, quantity = "Basal area increment", ylim = c(0,2), errorlim = c(-2,2)){ df_bas$PPET[df_bas$PPET>2] = 2 df_adv$PPET[df_adv$PPET>2] = 2 df <- bind_rows(df_bas, df_adv) g1<-ggplot(df, aes(x = PPET, y=pred))+ geom_point(aes(x = PPET, y= obs), size=0.3, col="gray", alpha = 0.3, data = df_bas)+ geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+ geom_smooth(aes(x = PPET, y = obs), col = "black", method="loess", formula = y ~ x, data = df_bas)+ geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+ ylab(paste0(quantity," (m2/ha/yr)"))+ ylim(ylim)+ xlab("Moisture index (P/PET)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(legend.position = "none", legend.background = element_blank(), panel.grid = element_blank()) g2<-ggplot(df, aes(x = PPET, y=error))+ geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+ geom_smooth(aes(col = transpirationMode, fill= transpirationMode), method="loess", formula = y ~ x)+ geom_abline(intercept = 0, slope = 0, col ="black")+ ylab(paste0(quantity," error (m2/ha/yr)"))+ylim(errorlim)+ xlab("Moisture index (P/PET)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = c(0.8,0.85)) g3<-ggplot(df, aes(x = MAT, y=pred))+ geom_point(aes(x = MAT, y= obs), size=0.3, col="gray", alpha = 0.3, data = df_bas)+ geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+ geom_smooth(aes(x = MAT, y = obs), col = "black", method="loess", formula = y ~ x, data = df_bas)+ geom_smooth(aes(col = transpirationMode, fill = transpirationMode), method="loess", formula = y ~ x)+ ylab(paste0(quantity," (m2/ha/yr)"))+ ylim(ylim)+ xlab("Mean annual temperature (ºC)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(legend.position = "none", legend.background = element_blank(), panel.grid = element_blank()) g4<-ggplot(df, aes(x = MAT, y=error))+ geom_point(aes(col = transpirationMode), size=0.3, alpha = 0.3)+ geom_smooth(aes(col = transpirationMode, fill= transpirationMode), method="loess", formula = y ~ x)+ geom_abline(intercept = 0, slope = 0, col ="black")+ ylab(paste0(quantity," error (m2/ha/yr)"))+ ylim(errorlim)+ xlab("Mean annual temperature (ºC)")+ scale_color_discrete("")+ scale_fill_discrete("")+ theme_bw()+ theme(panel.grid = element_blank(), legend.position = "none") plot_grid(g1, g2, g3, g4, nrow=2) }
bai_error_map <-function(data_bai, upper = 1, lower = -1){ data_bai = data_bai[!is.na(data_bai$error),] data_bai$error[data_bai$error > upper] = upper data_bai$error[data_bai$error < lower] = lower ggplot(data_bai)+ geom_sf(aes(col=error), size=0.5)+ geom_sf(data = cat_sfc)+ scale_color_fermenter("m2/ha/yr",palette = 1, breaks = c(lower,lower*0.5,lower*0.25,lower*0.1, 0, upper*0.1,upper*0.25,upper*0.5,upper), type = 'div', direction = 2, na.value = 'transparent')+ # scale_color_gradient2("m2/ha",low="blue", mid="yellow", high="red")+ theme_bw()+ theme(panel.grid = element_blank()) # theme(panel.grid = element_blank(), legend.position = c(0.85,0.25)) }
Simulation procedure
We selected permanent plots between the second (IFN2) and the fourth (IFN4) without signs of management (i.e. the presence of stumps) and whose basal area did not decrease more than 10% between inventory surveys (to avoid the effect of disturbances).
Soil physical properties were drawn from SoilGrids (Hengl 2016), complemented by rock fragment content estimates derived from surface stoniness measurements in forest plots.
fordyn simulations were conducted for different periods:
- IFN2 - IFN3 (~ 10 years)
- IFN3 - IFN4 (~ 15 years)
- IFN2 - IFN4 (~ 25 years)
The actual simulated period varied depending on the sampling years of the target plot. Daily weather data were obtained via interpolation on plot's coordinates using package meteoland.
Default species-specific parameters were modified using the results of the meta-modelling exercise and the growth calibration exercise. These two exercises do not provide values for all the main species included here, so it is expected that evaluation results are worse for those species not included in those exercises.
Simulations were done for both the basic (i.e. Granier) and advanced (i.e. Sperry) transpiration/photosynthesis sub-models. On a server with 20 parallel threads, computational times for the longest simulation period (25 years) are around 4 hours (i.e. 2.5 min/plot) for the basic sub-model, versus around 6 days (i.e. 1.5 hr/plot) for the advanced submodel.
In the following sections, we provide the bias, root mean squared error (in absolute and relative terms) and explained variance (R-squared) of growth and mortality predictions at the tree-level and stand-level obtained by simulations with both sub-models. Scatter plots are provided for the IFN2-IFN4 simulation to represent the relationship between predicted and observed values, as well as the factors that may influence the direction and magnitude of prediction error (i.e. initial values, environmental conditions, ...).
Detailed results of growth evaluation by tree species are provided in the last section.
Growth of surviving trees
Comparison of diameter/height growth of trees (DBH >= 7.5) that survived between surveys.
Annual diameter increment
Overall predictive capacity to predict diameter increase (cm/yr):
dfdbh |> evaluation_stats() |> kbl() |> kable_styling()
Predictive capacity plots (IFN2-IFN4):
dfdbh_bas <- dfdbh |> filter(transpirationMode == "Granier", period == "IFN24") dfdbh_adv <- dfdbh |> filter(transpirationMode == "Sperry", period == "IFN24") plot_scatter_tree(dfdbh_bas, dfdbh_adv, points = TRUE)
Relationship between diameter increase and climatic variables (MAT, P/PET and available PAR; IFN2 - IFN4):
plot_cov_clim_tree(dfdbh_bas, dfdbh_adv, points = TRUE)
Annual height increment
Overall predictive capacity to predict height increase (cm/yr):
dfh %>% evaluation_stats() %>% kbl() %>% kable_styling()
Predictive capacity plots (IFN2-IFN4):
dfh_bas <- dfh |> filter(transpirationMode == "Granier", period == "IFN24") dfh_adv <- dfh |> filter(transpirationMode == "Sperry", period == "IFN24") plot_scatter_tree(dfh_bas, dfh_adv, points = TRUE, xylim = c(0,20), errorlim = c(-20,20), initiallim = c(0, 2000), var = "height")
Relationship between height increase and climatic variables (MAT, P/PET and available PAR; IFN2 - IFN4):
plot_cov_clim_tree(dfh_bas, dfh_adv, points = TRUE, xylim = c(0,20), errorlim = c(-20,20), var = "height")
Stand-level basal area increment
Comparison of basal area increment of surviving trees does not take into account changes in density. In other words, densities from the first inventory are used to calculate stand-level basal area of surviving trees. Hence, the comparison is meant to evaluate the effect of diameter increment of surviving trees in terms of annual stand basal area increments (m2/ha/yr) for the period evaluated.
Predictive capacity table:
dfbai_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfbai_bas <- dfbai_plot |> filter(transpirationMode == "Granier", period == "IFN24") dfbai_adv <- dfbai_plot |> filter(transpirationMode == "Sperry", period == "IFN24") plot_scatter_bai(dfbai_bas, dfbai_adv)
Relationship between basal area increase and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfbai_bas, dfbai_adv, errorlim = c(-1,1), ylim = c(0,1))
Spatial error distribution (IFN2 - IFN4):
p1<-bai_error_map(dfbai_bas)+labs(title="Basic") p2<-bai_error_map(dfbai_adv)+labs(title="Advanced") plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"), get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Mortality
Basal area reduction
Annual reduction of basal area (m2/ha/yr) due to trees (DBH >= 7.5) that died during the evaluation period against model's mortality prediction. In both cases, basal area is calculated using the initial diameter of the trees, so that density reductions are the only prediction that is actually evaluated (and not the possible growth of those trees during the simulation).
Overall predictive capacity:
dfbadead_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfdead_bas <- dfbadead_plot |> filter(transpirationMode == "Granier", period == "IFN24") dfdead_adv <- dfbadead_plot |> filter(transpirationMode == "Sperry", period == "IFN24") plot_scatter_bai(dfdead_bas, dfdead_adv, quantity = "dead basal area", xylim = c(0,0.5), errorlim = c(-0.5,0.5))
Relationship between dead basal area and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfdead_bas, dfdead_adv, quantity = "Dead basal area", ylim = c(0,0.5), errorlim = c(-0.5,0.5))
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dfdead_bas, upper = 0.5, lower = -0.5)+labs(title="Basic") p2<-bai_error_map(dfdead_bas, upper = 0.5, lower = -0.5)+labs(title="Advanced") plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"), get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Density reduction
Annual reduction of density (ind/ha/yr) due to trees (DBH >= 7.5) that died during the evaluation period against model's mortality prediction. This is very similar to evaluating the reduction in basal area
Overall predictive capacity:
dfNdead_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Ingrowth
Basal area increase
Annual increase of basal area (m2/ha/yr) due to ingrowth of trees with diameters between 7.5 cm and 12.5 cm during the evaluated period.
Predictive capacity:
dfbarecr_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Predictive capacity plots (IFN2 - IFN4):
dfrecr_bas <- dfbarecr_plot |> filter(transpirationMode == "Granier", period == "IFN24") dfrecr_adv <- dfbarecr_plot |> filter(transpirationMode == "Sperry", period == "IFN24") plot_scatter_bai(dfrecr_bas, dfrecr_adv, quantity = "ingrowth basal area", xylim = c(0,0.15), errorlim = c(-0.15,0.15))
Relationship between ingrowth basal area and climatic variables (MAT and P/PET; IFN2 - IFN4):
plot_cov_clim_bai(dfrecr_bas, dfrecr_adv, quantity = "Ingrowth basal area", ylim = c(0,0.15), errorlim = c(-0.15,0.15))
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dfrecr_bas, upper = 0.5, lower = -0.5)+labs(title="Basic") p2<-bai_error_map(dfrecr_bas, upper = 0.5, lower = -0.5)+labs(title="Advanced") plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"), get_legend(p1),nrow=1, rel_widths = c(1,1,0.25))
Density increase
Annual increase of density (ind/ha/yr) due to ingrowth of trees with diameters between 7.5 cm and 12.5 cm during the evaluated period.
Predictive capacity:
dfNrecr_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Overall stand-level change
Basal area changes
This includes annual changes in basal area (m2/ha/yr) due to growth of surviving trees, mortality reductions and ingrowth derived from sapling growth. In the observed data, basal area changes include also the incorporation of trees into large diameter classes that results from the variable-radius sampling design. Since it takes into account all processes together, this evaluation is the most rellevant of all.
Overall predictive capacity:
dfbatot_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Predictive capacity plot (IFN2 - IFN4):
dftot_bas <- dfbatot_plot |> filter(transpirationMode == "Granier", period == "IFN24") dftot_adv <- dfbatot_plot |> filter(transpirationMode == "Sperry", period == "IFN24") p1<-plot_scatter_bai(dftot_bas, dftot_adv, quantity = "basal area change", xylim = c(-0.5,2), errorlim = c(-2,2)) p1
Relationship between overall basal area change and climatic variables (MAT and P/PET; IFN2 - IFN4):
p2<-plot_cov_clim_bai(dftot_bas, dftot_adv, quantity = "Basal area change", ylim = c(-0.5,2), errorlim = c(-2,2)) p2
p0<-cowplot::plot_grid(p1,p2, nrow=2, rel_heights = c(0.33,0.66)) ggsave(paste0(evalDIR,"/plots/BAtot_eval.png"),p0, width = 8, height=12)
Spatial distribution of errors (IFN2 - IFN4):
p1<-bai_error_map(dftot_bas)+labs(title="Basic sub-model") p2<-bai_error_map(dftot_adv)+labs(title="Advanced sub-model") pm<-plot_grid(p1+theme(legend.position = "none"),p2+theme(legend.position = "none"), get_legend(p1),nrow=1, rel_widths = c(1,1,0.25)) ggsave(paste0(evalDIR,"/plots/BAtot_error_map.png"),pm, width = 13, height=6, bg="white") pm
Density changes
This includes annual changes in density (ind/ha/yr) due to growth of surviving trees, mortality reductions and ingrowth derived from sapling growth. In the observed data, changes include also the incorporation of trees into large diameter classes that results from the variable-radius sampling design.
Overall predictive capacity:
dfNtot_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Changes in leaf area index of trees > 7.5 cm
Overall predictive capacity using allometric equations:
dflai_allom_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Overall predictive capacity using state variables:
dflai_state_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Shrub cover and mean height
Shrub percent cover changes
Overall predictive capacity:
dfshcov_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Shrub mean height changes
Overall predictive capacity:
dfshheight_plot %>% evaluation_stats() %>% kbl() %>% kable_styling()
Detailed evaluation by tree species
In the following we provide detailed evaluation results for the most important species.
sp100 = c("Abies alba", "Fagus sylvatica", "Pinus halepensis", "Pinus nigra", "Pinus sylvestris", "Pinus uncinata", "Quercus faginea", "Quercus ilex", "Quercus pubescens", "Quercus suber") for(i in 1:length(sp100)) { sp_name = sp100[i] res <- knitr::knit_child('_RLEvaluation_child.Rmd', quiet = TRUE) cat(res, sep = '\n') }
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.