1747 lines
67 KiB
R
1747 lines
67 KiB
R
# READ forceps EXPORTS
|
|
|
|
################################################################################
|
|
# load libraries
|
|
|
|
library(ggplot2)
|
|
library(dplyr)
|
|
|
|
################################################################################
|
|
|
|
# 1. LOAD DATA ----
|
|
|
|
# list folders
|
|
dirs <- list.dirs(getwd(), full.names = FALSE, recursive = FALSE)
|
|
dirs <- dirs[grepl(pattern = "PYR", dirs)]
|
|
dirs
|
|
|
|
# list files
|
|
files <- list.files(getwd(), full.names = FALSE, recursive = T, pattern = "PYR")
|
|
files
|
|
|
|
################################################################################
|
|
|
|
# filter files
|
|
|
|
## complete export datamean_stand_export_list
|
|
complete_export <- files[grep("forceps-complete-export", files)]
|
|
|
|
headers_complete_export <- readRDS(file = "headers_complete_export.Rds")
|
|
|
|
## mean stand export data
|
|
mean_stand_export <- files[grep("forceps-mean-stand-export", files)]
|
|
|
|
headers_mean_stand_export <- readRDS(file = "headers_mean_stand_export.Rds")
|
|
|
|
## IBP stand export data
|
|
IBP_export <- files[grep("forceps-IBP-export", files)]
|
|
|
|
headers_IBP_export <- readRDS(file = "headers_IBP_export.Rds")
|
|
|
|
## cut stand export data
|
|
cut_export <- files[grep("forceps-complete-cut-export", files)]
|
|
|
|
headers_cut_export <- readRDS(file = "headers_cut_export.Rds")
|
|
|
|
## scene deadwood export data
|
|
deadwood_export <- files[grep("forceps-scene-deadwood-export", files)]
|
|
|
|
headers_deadwood_export <- readRDS(file = "headers_scene_deadwood_export.Rds")
|
|
|
|
################################################################################
|
|
|
|
# Prepare functions
|
|
## separate filenames to add new columns in DFs
|
|
string_list <- regmatches(x = complete_export,
|
|
gregexpr('[._]',
|
|
complete_export),
|
|
invert = TRUE)
|
|
|
|
# FUN to extract substrings of lists
|
|
string_extractor <- function(lst, n){
|
|
sapply(lst, `[`, n)
|
|
}
|
|
|
|
# 2. OPERATE DATA ----
|
|
|
|
## 2.1. complete_export
|
|
|
|
# read data
|
|
complete_export_list <- lapply(complete_export, read.csv, col.names = headers_complete_export, sep = "\t")
|
|
# names(complete_export_list) <- dirs
|
|
# complete_export_list <- lapply(complete_export_list, subset, patchId != 0) #CHECK TO DO BE AWARE OJO!!!
|
|
|
|
## add more info
|
|
complete_export_list <- Map(cbind, complete_export_list,
|
|
PLOT = string_extractor(string_list, 2))
|
|
complete_export_list <- Map(cbind, complete_export_list,
|
|
PROJECTION = string_extractor(string_list, 4)) # 4th element is projection scenario
|
|
complete_export_list <- Map(cbind, complete_export_list,
|
|
MANAGEMENT = string_extractor(string_list, 3)) # 3th element is management scenario
|
|
|
|
# turn list into DF
|
|
complete_export_list <- bind_rows(complete_export_list)
|
|
|
|
complete_export_list$uniqueID_tree <- paste(
|
|
complete_export_list$PLOT,
|
|
complete_export_list$PROJECTION,
|
|
complete_export_list$MANAGEMENT,
|
|
complete_export_list$id,
|
|
sep = '_')
|
|
|
|
complete_export_list$uniqueID <- paste(
|
|
complete_export_list$PLOT,
|
|
complete_export_list$PROJECTION,
|
|
complete_export_list$MANAGEMENT,
|
|
sep = '_')
|
|
|
|
## 2.2. mean stand_export
|
|
|
|
# read data
|
|
mean_stand_export_list <- lapply(mean_stand_export, read.csv, col.names = headers_mean_stand_export, sep = '\t')
|
|
# names(mean_stand_export_list) <- dirs
|
|
|
|
## add more info
|
|
mean_stand_export_list <- Map(cbind, mean_stand_export_list,
|
|
PLOT = string_extractor(string_list, 2))
|
|
mean_stand_export_list <- Map(cbind, mean_stand_export_list,
|
|
PROJECTION = string_extractor(string_list, 4)) # 4th element is projection scenario
|
|
mean_stand_export_list <- Map(cbind, mean_stand_export_list,
|
|
MANAGEMENT = string_extractor(string_list, 3)) # 3th element is management scenario
|
|
|
|
# turn list into DF
|
|
mean_stand_export_list <- bind_rows(mean_stand_export_list)
|
|
|
|
mean_stand_export_list$uniqueID <- paste(
|
|
mean_stand_export_list$PLOT,
|
|
mean_stand_export_list$PROJECTION,
|
|
mean_stand_export_list$MANAGEMENT,
|
|
sep = "_")
|
|
|
|
## 2.3. IBP stand_export
|
|
|
|
# read data
|
|
IBP_export_list <- lapply(IBP_export, read.table, col.names = headers_IBP_export, sep = "\t")
|
|
|
|
## add more info
|
|
IBP_export_list <- Map(cbind, IBP_export_list,
|
|
PLOT = string_extractor(string_list, 2)) # 4th element is projection scenario
|
|
IBP_export_list <- Map(cbind, IBP_export_list,
|
|
PROJECTION = string_extractor(string_list, 4)) # 4th element is projection scenario
|
|
IBP_export_list <- Map(cbind, IBP_export_list,
|
|
MANAGEMENT = string_extractor(string_list, 3)) # 3th element is management scenario
|
|
|
|
# turn list into DF
|
|
IBP_export_list <- bind_rows(IBP_export_list)
|
|
|
|
IBP_export_list$uniqueID <- paste(
|
|
IBP_export_list$PLOT,
|
|
IBP_export_list$PROJECTION,
|
|
IBP_export_list$MANAGEMENT,
|
|
sep = "_")
|
|
|
|
## 2.4. cut stand_export
|
|
|
|
# read data
|
|
cut_export_list <- lapply(cut_export, read.table, col.names = headers_cut_export, sep = "\t")
|
|
|
|
# empty DF for cut trees
|
|
tp.mat <- data.frame(matrix(0, length(seq(2025,2100, by = 5)), length(headers_cut_export)))
|
|
tp.mat[,1] <- seq(2025,2100, by = 5)
|
|
colnames(tp.mat) <- headers_cut_export
|
|
tp.mat$speciesShortName <- NA
|
|
tp.mat$genotypeName <- NA
|
|
|
|
# loop
|
|
for (df in 1:length(cut_export_list)) {
|
|
|
|
if(nrow(cut_export_list[[df]]) == 0) {cut_export_list[[df]] = tp.mat} else next
|
|
|
|
}
|
|
|
|
## add more info
|
|
cut_export_list <- Map(cbind, cut_export_list,
|
|
PLOT = string_extractor(string_list, 2)) # 2th element is projection scenario
|
|
cut_export_list <- Map(cbind, cut_export_list,
|
|
PROJECTION = string_extractor(string_list, 4)) # 4th element is projection scenario
|
|
cut_export_list <- Map(cbind, cut_export_list,
|
|
MANAGEMENT = string_extractor(string_list, 3)) # 3th element is management scenario
|
|
|
|
# turn list into DF
|
|
cut_export_list <- bind_rows(cut_export_list)
|
|
|
|
cut_export_list$uniqueID <- paste(
|
|
cut_export_list$PLOT,
|
|
cut_export_list$PROJECTION,
|
|
cut_export_list$MANAGEMENT,
|
|
sep = "_")
|
|
|
|
## 2.5. deadwood_expot
|
|
|
|
# read data
|
|
complete_deadwood_list <- lapply(deadwood_export, read.csv, col.names = headers_deadwood_export, sep = "\t")
|
|
# complete_deadwood_list <- lapply(complete_deadwood_list, subset)
|
|
|
|
## add more info
|
|
complete_deadwood_list <- Map(cbind, complete_deadwood_list,
|
|
PLOT = string_extractor(string_list, 2))
|
|
complete_deadwood_list <- Map(cbind, complete_deadwood_list,
|
|
PROJECTION = string_extractor(string_list, 4)) # 4th element is projection scenario
|
|
complete_deadwood_list <- Map(cbind, complete_deadwood_list,
|
|
MANAGEMENT = string_extractor(string_list, 3)) # 3th element is management scenario
|
|
|
|
# turn list into DF
|
|
complete_deadwood_list <- bind_rows(complete_deadwood_list)
|
|
|
|
complete_deadwood_list$uniqueID_tree <- paste(
|
|
complete_deadwood_list$PLOT,
|
|
complete_deadwood_list$PROJECTION,
|
|
complete_deadwood_list$MANAGEMENT,
|
|
complete_deadwood_list$id,
|
|
sep = '_')
|
|
|
|
complete_deadwood_list$uniqueID <- paste(
|
|
complete_deadwood_list$PLOT,
|
|
complete_deadwood_list$PROJECTION,
|
|
complete_deadwood_list$MANAGEMENT,
|
|
sep = '_')
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# 3. PLOT -----
|
|
|
|
# test
|
|
"ggplot(bind_rows(complete_export_list, .id = 'data_frame'),
|
|
aes(col1, col2, colour = data_frame)) +
|
|
geom_line(aes(x = date, y = basalArea_m2, colour = factor(id)),
|
|
show.legend = F) +
|
|
facet_grid(PLOT ~ .)
|
|
|
|
complete_export_list |>
|
|
group_by(speciesShortName, PLOT, date) |>
|
|
summarise(mean_BA = mean(basalArea_m2)) |>
|
|
ggplot(aes(colour = speciesShortName)) +
|
|
geom_line(aes(x = date, y = mean_BA)) +
|
|
facet_grid(. ~ PLOT) + theme_bw()"
|
|
|
|
# ggplot(subset(mean_stand_export_list, MANAGEMENT == "01")) +
|
|
ggplot(mean_stand_export_list) +
|
|
facet_grid(PROJECTION~MANAGEMENT) +
|
|
geom_line(aes(
|
|
x = date,
|
|
y = basalAreaMean,
|
|
colour = PLOT
|
|
))
|
|
|
|
# geom_line(aes(x = date, y = basalAreaMean*20, colour = MANAGEMENT, linetype = PROJECTION)) +
|
|
# scale_y_continuous(
|
|
# name = "nTreesMean",
|
|
# sec.axis = sec_axis(~./20, name="BA"))
|
|
|
|
ggplot(subset(mean_stand_export_list, MANAGEMENT == "01")) +
|
|
facet_wrap(~PLOT) +
|
|
geom_line(aes(x = date, y = basalAreaMean, colour = MANAGEMENT, linetype = PROJECTION))
|
|
# scale_y_continuous(
|
|
# name = "nTreesMean",
|
|
# sec.axis = sec_axis(~./20, name="BA"))
|
|
#
|
|
# mean_stand_export_list |>
|
|
# group_by(speciesShortName, PLOT, date) |>
|
|
# summarise(mean_BA = mean(basalArea_m2)) |>
|
|
# ggplot(aes(colour = speciesShortName)) +
|
|
# geom_line(aes(x = date, y = mean_BA)) +
|
|
# facet_grid(. ~ PLOT) + theme_bw()
|
|
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# 4. Compute Ecosystem services ------------------------------------------------
|
|
|
|
# load the rolling function lag_apply
|
|
source("lag_apply.R")
|
|
|
|
movingWINDOW = 20 # in years
|
|
|
|
# create DF
|
|
services_DF <- data.frame(date = as.numeric(mean_stand_export_list$date),
|
|
PLOT = as.character(mean_stand_export_list$PLOT),
|
|
PROJECTION = as.character(mean_stand_export_list$PROJECTION),
|
|
MANAGEMENT = as.character(mean_stand_export_list$MANAGEMENT),
|
|
uniqueID = as.character(mean_stand_export_list$uniqueID))
|
|
|
|
#-------------------------------------------------------------------------------
|
|
# 4.1. compute ECO-SER ---------------------------------------------------------
|
|
|
|
# 4.1.2. MEAN stand data
|
|
|
|
# Plot tree biomass (totalBiomass_t.ha -> Mg Ha-1 yr-1)
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(mean_stand_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
PlotTreeBiomass = totalBasalArea_m2.ha,
|
|
meanPlotTreeBiomass = lag_apply(totalBasalArea_m2.ha, movingWINDOW, function(x) {mean(x)})
|
|
#, sdPTB = lag_apply(totalBasalArea_m2.ha, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
|
|
# Tree recruitment = N trees (nTrees.ha -> N yr-1)
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(mean_stand_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
N_trees = nTrees.ha,
|
|
meanN_trees = lag_apply(nTrees.ha, movingWINDOW, function(x) {mean(x)})
|
|
#, sdNTREE = lag_apply(nTrees.ha, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
|
|
# some PLOTS ----
|
|
|
|
"ggplot(services_DF, aes(x = date, group = MANAGEMENT)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = PlotTreeBiomass, colour = MANAGEMENT)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = PlotTreeBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date)) +
|
|
facet_grid(PROJECTION ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanPlotTreeBiomass, colour = PLOT)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = uniqueID)) +
|
|
geom_line(aes(y = meanPlotTreeBiomass, colour = uniqueID)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = uniqueID)) +
|
|
facet_grid(PROJECTION ~ MANAGEMENT) +
|
|
geom_line(aes(y = N_trees, colour = uniqueID)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = N_trees, colour = PROJECTION)) +
|
|
theme_bw()
|
|
"
|
|
|
|
ggplot(subset(services_DF, MANAGEMENT == "03" | MANAGEMENT == "01"),
|
|
aes(x = date, group = PROJECTION)) +
|
|
lemon::facet_rep_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = N_trees, colour = PROJECTION), size = 1.5) +
|
|
scale_color_manual(values = c('orange', '#6e2878ff', '#d4cfa6ff')) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = "Ecosystem service",
|
|
title = bquote('Trees' ~ (number/Ha)),
|
|
color = "Projection")
|
|
|
|
|
|
# ----
|
|
|
|
# 4.1.3. Individual tree data
|
|
|
|
# Tree species = N tree species
|
|
temporal <- data.frame(complete_export_list |>
|
|
group_by(date, patchId, uniqueID) |>
|
|
filter(dbh_cm > 10.0) |>
|
|
reframe(SP_N = length(unique(speciesId))) |>
|
|
ungroup() |> group_by(date, uniqueID) |>
|
|
reframe(meanpatchSP_N = mean(SP_N)))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(services_DF, data.frame(
|
|
temporal |> group_by(uniqueID) |>
|
|
reframe(
|
|
date = date,
|
|
N_species = meanpatchSP_N,
|
|
meanN_species = lag_apply(meanpatchSP_N, movingWINDOW, function(x) {mean(x)})
|
|
# , sdSP_N = lag_apply(meanpatchSP_N, movingWINDOW, function(x) {sd(x)})
|
|
)
|
|
))
|
|
|
|
services_DF[is.na(services_DF$N_species), 'N_species'] <- 0
|
|
|
|
'#ggplot(services_DF, aes(x = date, group = uniqueID)) +
|
|
#ggplot(na.omit(subset(services_DF, MANAGEMENT == "01")), aes(x = date)) +
|
|
# geom_line(aes(y = meanN_species,
|
|
# colour = PLOT, linetype = PROJECTION)) +
|
|
# theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = N_species, colour = PROJECTION)) +
|
|
theme_bw()
|
|
'
|
|
|
|
# Sapling growth HEIGHT OR DBH (??)
|
|
|
|
temporal <- complete_export_list |>
|
|
filter(dbh_cm < 10.0) |>
|
|
filter(height_m != 0) |>
|
|
group_by(id, uniqueID, patchId) |>
|
|
reframe(date = date,
|
|
height_m = height_m,
|
|
prev_height_m = lag(height_m)) |> filter(!is.na(prev_height_m)) |>
|
|
mutate(SAP_HEIGTH = height_m - prev_height_m) |> ungroup() |> group_by(date, uniqueID) |>
|
|
reframe(meanPatchSaplingH = mean(SAP_HEIGTH)) |> ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
SaplingHeigthIncrement = meanPatchSaplingH,
|
|
meanSaplingHeigthIncrement = lag_apply(meanPatchSaplingH, movingWINDOW, function(x) {mean(x, na.rm = T)})
|
|
# , sdSaplingHincrement = lag_apply(meanSaplingH, movingWINDOW, function(x) {sd(x)}),
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = MANAGEMENT)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SaplingHeigthIncrement, colour = MANAGEMENT)) +
|
|
theme_bw()
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = SaplingHeigthIncrement, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# Sapling species diversity (Species number) & sapling biomass
|
|
|
|
temporal <- complete_export_list |>
|
|
filter(dbh_cm < 10.0) |>
|
|
filter(height_m != 0) |>
|
|
group_by(id, uniqueID, patchId) |>
|
|
reframe(date = date,
|
|
speciesShortName = speciesShortName,
|
|
vol = volume_m3) |>
|
|
filter(!is.na(vol)) |>
|
|
ungroup() |> group_by(date, uniqueID) |>
|
|
reframe(meanSaplingN = mean(length(unique(speciesShortName))),
|
|
meanSaplingV = mean(vol)) |>
|
|
ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
SaplingSpeciesNumber = meanSaplingN,
|
|
meanSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {mean(x, na.rm = T)}),
|
|
# sdSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {sd(x)}),
|
|
SaplingVolume = meanSaplingV,
|
|
meanSaplingVolume = lag_apply(meanSaplingV, movingWINDOW, function(x) {mean(x, na.rm = T)})
|
|
#, sdSaplingVolume = lag_apply(meanSaplingV, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
services_DF[is.na(services_DF$SaplingSpeciesNumber), 'SaplingSpeciesNumber'] <- 0
|
|
services_DF[is.na(services_DF$SaplingVolume), 'SaplingVolume'] <- 0
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = SaplingVolume, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = SaplingSpeciesNumber, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# GROWTH STABILITY
|
|
temporal <- complete_export_list |>
|
|
filter(dbh_cm > 10.0) |>
|
|
group_by(uniqueID_tree, uniqueID, patchId) |>
|
|
reframe(date = date,
|
|
AGB = aboveGroundBiomass,
|
|
prev_AGB = lag(AGB)) |> filter(!is.na(prev_AGB)) |>
|
|
mutate(AWP = AGB - prev_AGB) |> ungroup() |> group_by(date, uniqueID) |>
|
|
reframe(plotAWP = mean(AWP)) |> ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
meanAWP = lag_apply(plotAWP, movingWINDOW, function(x) {mean(x)}),
|
|
sdAWP = lag_apply(plotAWP, movingWINDOW, function(x) {sd(x)}),
|
|
STABILITY = meanAWP/sdAWP
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = STABILITY, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# ROOT BIOMASS
|
|
temporal <- complete_export_list |> group_by(date, patchId, uniqueID) |>
|
|
reframe(ROOT_BIOMASS = sum(rootBiomass_kg)) |> ungroup() |> group_by(date, uniqueID) |>
|
|
reframe(meanpatchROOT_BIOMASS = mean(ROOT_BIOMASS))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
RootBiomass = meanpatchROOT_BIOMASS,
|
|
meanROOT_BIOMASS = lag_apply(meanpatchROOT_BIOMASS, movingWINDOW, function(x) {mean(x)})
|
|
# , sdROOT_BIOMASS = lag_apply(meanpatchROOT_BIOMASS, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
services_DF[is.na(services_DF$RootBiomass), 'RootBiomass'] <- 0
|
|
|
|
'ggplot(subset(services_DF, meanROOT_BIOMASS < 20000), aes(x = date)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanROOT_BIOMASS, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# Complete fine root biomass RATES
|
|
"doesn't look like a good idea to use this EcoSer"
|
|
|
|
# Fine Root Production (kg yr-1)
|
|
# temporal <- complete_export_list |> group_by(uniqueID, date, patchId) |>
|
|
# reframe(meanpatchFineRootProduction = sum(fineRootBiomass_kg)) |> ungroup() |> group_by(date, uniqueID) |>
|
|
# reframe(date = date,
|
|
# meanpatchFineRootProduction = mean(meanpatchFineRootProduction),
|
|
# prev_finerootbiomass = lag(meanpatchFineRootProduction),
|
|
# FineRootProduction = meanpatchFineRootProduction - prev_finerootbiomass)
|
|
|
|
temporal <- complete_export_list |> group_by(uniqueID, date) |>
|
|
reframe(meanpatchFineRootProduction = sum(fineRootBiomass_kg)) |> ungroup() |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
meanpatchFineRootProduction = meanpatchFineRootProduction,
|
|
prev_finerootbiomass = lag(meanpatchFineRootProduction),
|
|
FineRootProduction = meanpatchFineRootProduction - prev_finerootbiomass) |> filter(!is.na(prev_finerootbiomass)) |>
|
|
ungroup() # |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
FineRootProduction = FineRootProduction,
|
|
meanFineRootProduction = lag_apply(FineRootProduction, movingWINDOW, function(x) {mean(x)})
|
|
# , sdFinRooBioPROD = lag_apply(FineRootProduction, movingWINDOW, function(x) {sd(x)}),
|
|
))
|
|
)
|
|
|
|
services_DF[is.na(services_DF$FineRootProduction), 'FineRootProduction'] <- 0
|
|
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanFineRootProduction, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# Mortality, expressed by Ha⁻¹
|
|
# This is death counts as number of Declinig (slowgrowth) and drougth (sureauDeath)
|
|
|
|
# causeOfDeath == "SlowGrowth Death" | causeOfDeath == "sureauDeath"
|
|
|
|
temporal <- complete_export_list |>
|
|
filter(dbh_cm >= 10.0,
|
|
causeOfDeath == "SlowGrowth Death" | causeOfDeath == "sureauDeath") |>
|
|
group_by(uniqueID, id) |>
|
|
reframe(patchId = patchId,
|
|
date = date,
|
|
status = "DEAD",
|
|
cause = ifelse(causeOfDeath == "SlowGrowth Death", "DECAY", "DR"))
|
|
|
|
# Total mortality per site
|
|
temporal2 <- temporal |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(numberTotalMortality = n()/max(patchId)) |>
|
|
ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal2 <- temporal2[order(temporal2$uniqueID, na.last = NA),]
|
|
|
|
temporal2 <- temporal2 |>
|
|
tidyr::complete(date, uniqueID) |>
|
|
mutate(numberTotalMortality = ifelse(is.na(numberTotalMortality),
|
|
0,
|
|
numberTotalMortality))
|
|
temporal2 <- temporal2 |>
|
|
group_by(uniqueID) |>
|
|
mutate(cumTotalMortality = cumsum(numberTotalMortality))
|
|
|
|
'# bien!
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = numberTotalMortality,
|
|
color = uniqueID),
|
|
show.legend = F)
|
|
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = cumTotalMortality,
|
|
color = uniqueID),
|
|
show.legend = F) + theme_bw()'
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal2 |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
numberTotalMortality = numberTotalMortality,
|
|
rollnumberTotalMortality = lag_apply(numberTotalMortality, movingWINDOW, function(x) {mean(x, na.rm = T)}),
|
|
# sdSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {sd(x)}),
|
|
cumTotalMortality = cumTotalMortality,
|
|
rollcumTotalMortality = lag_apply(cumTotalMortality, movingWINDOW, function(x) {mean(x, na.rm = T)})
|
|
#, sdSaplingVolume = lag_apply(meanSaplingV, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
services_DF[is.na(services_DF$numberTotalMortality), 'numberTotalMortality'] <- 0
|
|
services_DF[is.na(services_DF$cumTotalMortality), 'cumTotalMortality'] <- 0
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = numberTotalMortality, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = cumTotalMortality, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# DECAY mortality per site
|
|
temporal2 <- temporal |>
|
|
filter(cause == "DECAY") |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(numberDecayMortality = n()/max(patchId)) |>
|
|
ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal2 <- temporal2[order(temporal2$uniqueID, na.last = NA),]
|
|
|
|
temporal2 <- temporal2 |>
|
|
tidyr::complete(date, uniqueID) |>
|
|
mutate(numberDecayMortality = ifelse(is.na(numberDecayMortality),
|
|
0,
|
|
numberDecayMortality))
|
|
temporal2 <- temporal2 |>
|
|
group_by(uniqueID) |>
|
|
mutate(cumDecayMortality = cumsum(numberDecayMortality))
|
|
|
|
'# bien!
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = numberDecayMortality,
|
|
color = uniqueID),
|
|
show.legend = F)
|
|
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = cumDecayMortality,
|
|
color = uniqueID),
|
|
show.legend = F) + theme_bw()'
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal2 |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
numberDecayMortality = numberDecayMortality,
|
|
rollnumberDecayMortality = lag_apply(numberDecayMortality, movingWINDOW, function(x) {mean(x, na.rm = T)}),
|
|
# sdSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {sd(x)}),
|
|
cumDecayMortality = cumDecayMortality,
|
|
rollcumDecayMortality = lag_apply(cumDecayMortality, movingWINDOW, function(x) {mean(x, na.rm = T)})
|
|
#, sdSaplingVolume = lag_apply(meanSaplingV, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
services_DF[is.na(services_DF$numberDecayMortality), 'numberDecayMortality'] <- 0
|
|
services_DF[is.na(services_DF$cumDecayMortality), 'cumDecayMortality'] <- 0
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = numberDecayMortality, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = cumDecayMortality, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# DR mortality per site
|
|
temporal2 <- temporal |>
|
|
filter(cause == "DR") |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(numberDRMortality = n()/max(patchId)) |>
|
|
ungroup() |> tidyr::complete(date = 2022:2100, uniqueID, fill = list(val = 0))
|
|
temporal2 <- temporal2[order(temporal2$uniqueID, na.last = NA),]
|
|
|
|
temporal2 <- temporal2 |>
|
|
tidyr::complete(date, uniqueID) |>
|
|
mutate(numberDRMortality = ifelse(is.na(numberDRMortality),
|
|
0,
|
|
numberDRMortality))
|
|
temporal2 <- temporal2 |>
|
|
group_by(uniqueID) |>
|
|
mutate(cumDRMortality = cumsum(numberDRMortality))
|
|
|
|
'# bien!
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = numberDRMortality,
|
|
color = uniqueID),
|
|
show.legend = F)
|
|
|
|
ggplot(temporal2) + geom_line(aes(x = date,
|
|
y = cumDRMortality,
|
|
color = uniqueID),
|
|
show.legend = F) + theme_bw()'
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal2 |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
numberDRMortality = numberDRMortality,
|
|
rollnumberDRMortality = lag_apply(numberDRMortality, movingWINDOW, function(x) {mean(x, na.rm = T)}),
|
|
# sdSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {sd(x)}),
|
|
cumDRMortality = cumDRMortality,
|
|
rollcumDRMortality = lag_apply(cumDRMortality, movingWINDOW, function(x) {mean(x, na.rm = T)})
|
|
#, sdSaplingVolume = lag_apply(meanSaplingV, movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
services_DF[is.na(services_DF$numberDRMortality), 'numberDRMortality'] <- 0
|
|
services_DF[is.na(services_DF$cumDRMortality), 'cumDRMortality'] <- 0
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = numberDRMortality, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = cumDRMortality, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# Drought proxies
|
|
# This is death counts as number of Declinig (slowgrowth) and drougth (sureauDeath)
|
|
|
|
temporal <- complete_export_list |>
|
|
filter(dbh_cm >= 10.0) |>
|
|
group_by(uniqueID, id) |>
|
|
reframe(patchId = patchId,
|
|
date = date,
|
|
DRindex = treeDroughtIndex) |>
|
|
group_by(uniqueID, date) |>
|
|
summarise(meanDR = mean(DRindex))
|
|
|
|
'ggplot(temporal) +
|
|
geom_line(aes(x = date,
|
|
y = meanDR,
|
|
color = uniqueID),
|
|
show.legend = F)'
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
meanDR = meanDR,
|
|
rollmeanDR = lag_apply(meanDR, movingWINDOW, function(x) {mean(x, na.rm = T)}),
|
|
# sdSaplingSpeciesNumber = lag_apply(meanSaplingN, movingWINDOW, function(x) {sd(x)}),
|
|
))
|
|
)
|
|
|
|
# services_DF[is.na(services_DF$meanDR), 'meanDR'] <- 0
|
|
# services_DF[is.na(services_DF$cumTotalMortality), 'cumTotalMortality'] <- 0
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanDR, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = rollmeanDR, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# 4.1.3. IBP tree data -------------------------------------------------------
|
|
|
|
# N of Big Standing Trees (BST (N/Ha), DBH >= 67.5 cm aka VeryVeryLargeTrees)
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
BST_number = numberOfLiveVeryVeryLargeTrees,
|
|
meanBST_N = lag_apply(numberOfLiveVeryVeryLargeTrees, movingWINDOW, function(x) {mean(x)})
|
|
# , sdBST_N = lag_apply(numberOfLiveVeryVeryLargeTrees, movingWINDOW, function(x) {sd(x)}) # DONT USE SD, doen't make any sense
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = BST_number, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# Cumpling index (Thermal regulation potential, aka AlphaClumpingIndex)
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
ClumpingIndex = alphaClumpingIndex,
|
|
meanClumplingInd = lag_apply(alphaClumpingIndex, movingWINDOW, function(x) {mean(x)})
|
|
# , sdClumplingInd = lag_apply(alphaClumpingIndex, movingWINDOW, function(x) {sd(x)}) # DONT USE SD, doen't make any sense
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanClumplingInd, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# N of Standing Dead Trees (DeadBST_N (N/Ha), DBH >= 47.5 cm aka: numberOfStandingDeadVeryLargeTrees + numberOfStandingDeadVeryVeryLargeTrees)
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |>
|
|
group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
DeadBST_N = numberOfStandingDeadVeryLargeTrees +
|
|
numberOfStandingDeadVeryVeryLargeTrees,
|
|
meanDeadBST_N = lag_apply(
|
|
numberOfStandingDeadVeryLargeTrees +
|
|
numberOfStandingDeadVeryVeryLargeTrees,
|
|
movingWINDOW, function(x) {mean(x)})
|
|
# , sdDeadBST_N = lag_apply(numberOfStandingDeadVeryLargeTrees +
|
|
# numberOfStandingDeadVeryVeryLargeTrees, movingWINDOW,
|
|
# function(x) {sd(x)}) # DONT USE SD, doen't make any sense
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanDeadBST_N, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# N of TreM Trees (TREMTrees (N/Ha))
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
TREMTrees_N = numberOfTREMTrees,
|
|
meanTREMTrees_N = lag_apply(numberOfTREMTrees,
|
|
movingWINDOW, function(x) {mean(x)})
|
|
# , sdTREMTrees_N = lag_apply(numberOfTREMTrees, movingWINDOW,
|
|
# function(x) {sd(x)})
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanTREMTrees_N, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# Shannon alpha
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
SP_ShannonIndex = alphaSpeciesShannonIndex,
|
|
meanSP_ShannonIndex = lag_apply(alphaSpeciesShannonIndex,
|
|
movingWINDOW, function(x) {mean(x)})
|
|
# , sdTREMShannon = lag_apply(alphaSpeciesShannonIndex,
|
|
# movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
# ggplot(subset(services_DF, PLOT == "PYR22"), aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = SP_ShannonIndex, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(subset(services_DF, MANAGEMENT == "03" |
|
|
MANAGEMENT == "02" |
|
|
MANAGEMENT == "01"),
|
|
aes(x = date, group = MANAGEMENT)) +
|
|
lemon::facet_rep_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SP_ShannonIndex, colour = MANAGEMENT), size = 1.5) +
|
|
# scale_color_manual(values = c("#f97316", "#1e3a8a")) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = "Ecosystem service",
|
|
title = bquote("Tree species diversity" ~ (H)),
|
|
color = "Management")'
|
|
|
|
# Beta diversity
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(IBP_export_list |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
betaSP_ShannonIndex = betaSpeciesShannonIndex,
|
|
meanbetaSP_ShannonIndex = lag_apply(betaSpeciesShannonIndex,
|
|
movingWINDOW, function(x) {mean(x)})
|
|
# , sdTREMShannon = lag_apply(alphaSpeciesShannonIndex,
|
|
# movingWINDOW, function(x) {sd(x)})
|
|
))
|
|
)
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
# ggplot(subset(services_DF, PLOT == "PYR22"), aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = betaSP_ShannonIndex, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# 4.1.4. CUT tree data ---------------------------------------------------------
|
|
|
|
# 1st, make means and sums per patch and year
|
|
temporal <- cut_export_list |> group_by(uniqueID, patchId, dateOfCut) |>
|
|
reframe(date = dateOfCut,
|
|
aboveGroundBiomass_kg = aboveGroundBiomass_kg,
|
|
id = id,
|
|
patchAGB = sum(aboveGroundBiomass_kg),
|
|
patchVol = sum(volume_m3),
|
|
patchTB = sum(trunkBiomass_kg),
|
|
patchTBratio = patchTB/patchAGB) |> ungroup() |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(date = date,
|
|
AGB = mean(patchAGB),
|
|
Vol = mean(patchVol),
|
|
TB = mean(patchTB),
|
|
TBratio = mean(patchTBratio))
|
|
temporal <- distinct(temporal, uniqueID, AGB, Vol, TB, TBratio, .keep_all = T)
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
temporal <- temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
AGB = AGB,
|
|
Vol = Vol,
|
|
TB = TB,
|
|
TBratio = TBratio,
|
|
cumsumAGB = cumsum(AGB),
|
|
cumsumVol = cumsum(Vol),
|
|
cumsumTB = cumsum(TB))
|
|
cumulative_wood <- temporal
|
|
|
|
'ggplot(temporal, aes(x = date, group = uniqueID)) +
|
|
geom_step(aes(y = cumsumVol, colour = uniqueID), direction = "vh") +
|
|
geom_point(aes(y = cumsumVol, colour = uniqueID)) +
|
|
theme_bw()'
|
|
|
|
temporal <- temporal |> group_by(uniqueID) |>
|
|
tidyr::complete(date = 2022:2100) |>
|
|
tidyr::fill(c(AGB, Vol, TB, TBratio, cumsumAGB, cumsumVol, cumsumTB),
|
|
.direction = "down") |>
|
|
tidyr::fill(c(AGB, Vol, TB, TBratio, cumsumAGB, cumsumVol, cumsumTB),
|
|
.direction = "up")
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
Wood_AboveGroundBiomass = lag_apply(AGB, movingWINDOW, function(x) {mean(x)}),
|
|
# sdpatchAGB = lag_apply(patchAGB, movingWINDOW, function(x) {sd(x)}),
|
|
Wood_Volume = lag_apply(Vol, movingWINDOW, function(x) {mean(x)}),
|
|
# sdpatchVol = lag_apply(patchVol, movingWINDOW, function(x) {sd(x)}),
|
|
Wood_TrunkBiomass = lag_apply(TB, movingWINDOW, function(x) {mean(x)}),
|
|
# sdpatchTB = lag_apply(patchTB, movingWINDOW, function(x) {sd(x)}),
|
|
Wood_TrunkAboveRatio = lag_apply(TBratio, movingWINDOW, function(x) {mean(x)}),
|
|
# sdpatchTBratio = lag_apply(patchTBratio, movingWINDOW, function(x) {sd(x)})
|
|
cumsum_Wood_AboveGroundBiomass = cumsumAGB,
|
|
cumsum_Wood_Volume = cumsumVol,
|
|
cumsum_Wood_TB = cumsumTB
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
# ggplot(subset(services_DF, PLOT == "PYR22"), aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = cumsum_Wood_Volume, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = Wood_AboveGroundBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PLOT)) +
|
|
facet_grid(MANAGEMENT ~ PROJECTION) +
|
|
geom_step(aes(y = cumsum_Wood_Volume, colour = PLOT), direction = "vh") +
|
|
geom_point(aes(y = cumsum_Wood_Volume, colour = PLOT), size = .25) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(.~ MANAGEMENT) +
|
|
geom_line(aes(y = Wood_TrunkBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(MANAGEMENT ~ PLOT) +
|
|
geom_line(aes(y = cumsum_Wood_AboveGroundBiomass, colour = PROJECTION)) +
|
|
geom_line(aes(y = cumsum_Wood_AboveGroundBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = uniqueID)) +
|
|
geom_line(aes(y = meanpatchTBratio, colour = uniqueID)) +
|
|
theme_bw()'
|
|
|
|
# ----
|
|
|
|
# 4.1.5. Deadwood
|
|
|
|
# Litter production. AKA fallenLeafBiomassStock
|
|
temporal <- complete_deadwood_list |> group_by(uniqueID, date) |>
|
|
reframe(meanplotfallenLeafBiomassStock = sum(fallenLeafBiomassStock))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
subset(services_DF, date >= 2023),
|
|
data.frame(temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
FallenLeafBiomass = meanplotfallenLeafBiomassStock,
|
|
meanFallenLeafBiomass = lag_apply(meanplotfallenLeafBiomassStock, movingWINDOW, function(x) {mean(x)})
|
|
# sdFallenLeafBiomass = lag_apply(meanplotfallenLeafBiomassStock, movingWINDOW, function(x) {sd(x)}),
|
|
))
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanFallenLeafBiomass, colour = PROJECTION)) +
|
|
theme_bw()'
|
|
|
|
# SOC. AKA soilOrganicCarbonStock
|
|
|
|
temporal <- complete_deadwood_list |> group_by(uniqueID, date) |>
|
|
reframe(sumSOC = sum(soilOrganicCarbonStock))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(services_DF,
|
|
data.frame(
|
|
temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
currentSOC = sumSOC,
|
|
cumsumSOC = cumsum(sumSOC)) |> ungroup() |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
currentSOC = currentSOC,
|
|
cumsumSOC = cumsumSOC,
|
|
rollingSOC = lag_apply(currentSOC, movingWINDOW, function(x) {mean(x)})
|
|
# , sdSOC = lag_apply(meanSOC, movingWINDOW, function(x) {sd(x)})
|
|
)
|
|
)
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_step(aes(y = cumsumSOC, colour = PROJECTION), direction = "hv") +
|
|
theme_bw()'
|
|
|
|
# Soil Carbon Emissions AKA soilCarbonEmission
|
|
# IDK if this is a flux or a total value
|
|
|
|
# 1st: averaged flux over all patches
|
|
|
|
temporal <- complete_deadwood_list |> group_by(uniqueID, date) |>
|
|
reframe(meanplotSCEmissions = mean(soilCarbonEmission))
|
|
temporal <- temporal[order(temporal$uniqueID, na.last = NA),]
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(
|
|
temporal |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
soilCarbonEmission = meanplotSCEmissions,
|
|
cumsumSCEmission = cumsum(meanplotSCEmissions),
|
|
meanSCEmission = lag_apply(meanplotSCEmissions, movingWINDOW, function(x) {mean(x)})
|
|
# , sdSCEmissions = lag_apply(meanSCEmissions, movingWINDOW, function(x) {sd(x)}),
|
|
)
|
|
)
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(.~ MANAGEMENT) +
|
|
geom_step(aes(y = cumsumSCEmission, colour = PROJECTION), direction = "hv") +
|
|
theme_bw()
|
|
|
|
ggplot(services_DF, aes(x = date, group = MANAGEMENT)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = meanSCEmission, colour = MANAGEMENT), direction = "hv") +
|
|
theme_bw()'
|
|
|
|
# Coarse Woody Debris AKA deadBranchBiomassStock
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(
|
|
complete_deadwood_list |>
|
|
# mutate(pseudoCWD = deadBranchBiomassStock + deadFallenTreeBiomassStock) |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(CWD = mean(deadBranchBiomassStock)) |>
|
|
ungroup() |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
CWD = CWD,
|
|
cumsumCWD = cumsum(CWD),
|
|
meanCWD = lag_apply(CWD, movingWINDOW, function(x) {mean(x)})
|
|
# , sdCWD = lag_apply(CWD, movingWINDOW, function(x) {sd(x)}),
|
|
)
|
|
)
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanCWD, colour = PROJECTION)) +
|
|
theme_bw()
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(.~ MANAGEMENT) +
|
|
geom_step(aes(y = cumsumCWD, colour = PROJECTION), direction = "hv") +
|
|
theme_bw()'
|
|
|
|
# Dead trees biomass AKA deadFallenTreeBiomassStock
|
|
|
|
# ALERT: deadFallenTreeBiomassStock seems very correlated with SOC
|
|
|
|
services_DF <- plyr::join(
|
|
services_DF,
|
|
data.frame(
|
|
complete_deadwood_list |>
|
|
# mutate(pseudoCWD = deadBranchBiomassStock + deadFallenTreeBiomassStock) |>
|
|
group_by(uniqueID, date) |>
|
|
reframe(FallTreeBiomass = mean(deadFallenTreeBiomassStock)) |>
|
|
ungroup() |> group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
FallTreeBiomass = FallTreeBiomass,
|
|
cumsumFallTreeBiomass = cumsum(FallTreeBiomass),
|
|
meanFallTreeBiomass = lag_apply(FallTreeBiomass, movingWINDOW, function(x) {mean(x)})
|
|
# , sdCWD = lag_apply(CWD, movingWINDOW, function(x) {sd(x)}),
|
|
)
|
|
)
|
|
)
|
|
|
|
'ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_line(aes(y = meanFallTreeBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
ggplot(services_DF, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAGEMENT) +
|
|
geom_step(aes(y = cumsumFallTreeBiomass, colour = PROJECTION), direction = "hv") +
|
|
theme_bw()'
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# ------------------------------------------------------------------------------
|
|
# Operate DATA FRAME -----------------------------------------------------------
|
|
|
|
## Add new variable with MANAGEMENT TYPE
|
|
management_list <- read.csv("MANAGEMENT_TYPES.csv", sep = ',')
|
|
management_list$MANAG_TYPE <- as.factor(management_list$MANAG_TYPE)
|
|
management_list$MANAG_TYPE <- ordered(management_list$MANAG_TYPE,
|
|
levels = c('smFD', 'FD+', 'smIFM', 'BAU---', 'BAU--', 'BAU-', 'BAU-b', 'BAU',
|
|
'BAUb', 'BAU+', 'BAU++', 'BAU++a', 'BAU++b', 'BAU++c', 'BAU+++',
|
|
'BAU+++b', 'smCC'),
|
|
labels = c('smFD', 'FD+', 'smIFM', 'BAU---', 'BAU--', 'BAU-', 'BAU-b', 'BAU',
|
|
'BAUb', 'BAU+', 'BAU++', 'BAU++a', 'BAU++b', 'BAU++c', 'BAU+++',
|
|
'BAU+++b', 'smCC')
|
|
)
|
|
str(management_list)
|
|
|
|
# Create matching column
|
|
services_DF$MANAG_INTENSITY <- paste0(services_DF$PLOT, "_", services_DF$MANAGEMENT)
|
|
|
|
services_DF$MANAG_INTENSITY <- management_list$MANAG_TYPE[match(
|
|
unlist(services_DF$MANAG_INTENSITY),
|
|
management_list$UNION)]
|
|
|
|
# Standardize EcoSer (van der Plas 2016)
|
|
|
|
ES_df <- services_DF |>
|
|
group_by(PLOT) |>
|
|
reframe(date = date,
|
|
PROJECTION = factor(PROJECTION, levels = c('BASE', '245', '585'),
|
|
labels = c('Base', '245', '585')),
|
|
MANAGEMENT = MANAGEMENT,
|
|
uniqueID = uniqueID,
|
|
MANAG_INTENSITY = MANAG_INTENSITY,
|
|
PlotTreeBiomass = PlotTreeBiomass,
|
|
N_Trees = N_trees,
|
|
N_TreeSpecies = N_species,
|
|
N_SaplingSpecies = SaplingSpeciesNumber,
|
|
N_SaplingVolume = SaplingVolume,
|
|
RootBiomass = RootBiomass,
|
|
FineRootProduction = FineRootProduction,
|
|
SceneDrougth = meanDR,
|
|
GrowthStability = STABILITY,
|
|
N_BST = BST_number,
|
|
ClumpingIndex = ClumpingIndex,
|
|
N_DeadBST = DeadBST_N,
|
|
N_TREMTrees = TREMTrees_N,
|
|
SpeciesShannonIndex = SP_ShannonIndex,
|
|
BetaSpeciesShannonIndex = betaSP_ShannonIndex,
|
|
WoodAGB = cumsum_Wood_AboveGroundBiomass,
|
|
WoodAGB_V = cumsum_Wood_Volume,
|
|
WoodTrunkBiomass = cumsum_Wood_TB,
|
|
FallenLeafBiomass = FallenLeafBiomass,
|
|
currentSOC = currentSOC,
|
|
CWD = CWD,
|
|
FallTreeBiomass = FallTreeBiomass,
|
|
DecayMortality = cumDecayMortality,
|
|
DRMortality = cumDRMortality,
|
|
TotalMortality = cumTotalMortality
|
|
)
|
|
'write.csv(ES_df, "ECOSER_yearly.csv", row.names = F)'
|
|
|
|
## load the scalling function
|
|
source("scale_values.R")
|
|
|
|
### TODO:
|
|
need to scale values per site (?)
|
|
https://rforhr.com/center.html#review_center
|
|
|
|
SCA_ES_df <- services_DF |>
|
|
group_by(PLOT) |>
|
|
reframe(date = date,
|
|
PROJECTION = factor(PROJECTION, levels = c('BASE', '245', '585'),
|
|
labels = c('Base', '245', '585')),
|
|
MANAGEMENT = MANAGEMENT,
|
|
uniqueID = uniqueID,
|
|
MANAG_INTENSITY = MANAG_INTENSITY,
|
|
SCA_PlotTreeBiomass = scale_values(PlotTreeBiomass),
|
|
SCA_N_Trees = scale_values(N_trees),
|
|
SCA_N_TreeSpecies = scale_values(N_species),
|
|
SCA_N_SaplingSpecies = scale_values(SaplingSpeciesNumber),
|
|
SCA_N_SaplingVolume = scale_values(SaplingVolume),
|
|
SCA_RootBiomass = scale_values(RootBiomass),
|
|
SCA_FineRootProduction = scale_values(FineRootProduction),
|
|
# SCA_SceneDrougth = scale_values(meanDR),
|
|
# SCA_GrowthStability = scale_values(STABILITY),
|
|
SCA_N_BST = scale_values(BST_number),
|
|
SCA_ClumpingIndex = scale_values(ClumpingIndex),
|
|
SCA_N_DeadBST = scale_values(DeadBST_N),
|
|
SCA_N_TREMTrees = scale_values(TREMTrees_N),
|
|
SCA_SpeciesShannonIndex = scale_values(SP_ShannonIndex),
|
|
SCA_BetaSpeciesShannonIndex = scale_values(betaSP_ShannonIndex),
|
|
SCA_WoodAGB = scale_values(cumsum_Wood_AboveGroundBiomass),
|
|
SCA_WoodAGB_V = scale_values(cumsum_Wood_Volume),
|
|
SCA_WoodTrunkBiomass = scale_values(cumsum_Wood_TB),
|
|
SCA_FallenLeafBiomass = scale_values(FallenLeafBiomass),
|
|
SCA_currentSOC = scale_values(currentSOC),
|
|
SCA_CWD = scale_values(CWD),
|
|
SCA_FallTreeBiomass = scale_values(FallTreeBiomass),
|
|
SCA_DecayMortality = scale_values(cumDecayMortality),
|
|
SCA_DRMortality = scale_values(cumDRMortality),
|
|
SCA_TotalMortality = scale_values(cumTotalMortality)
|
|
)
|
|
|
|
## reframe meanDR and STABILITY separately because of NAs
|
|
|
|
SCA_meanDR <- services_DF |>
|
|
filter_at(vars(meanDR), all_vars(!is.na(.))) |>
|
|
group_by(PLOT) |>
|
|
reframe(date = date,
|
|
PROJECTION = factor(PROJECTION, levels = c('BASE', '245', '585'),
|
|
labels = c('Base', '245', '585')),
|
|
MANAGEMENT = MANAGEMENT,
|
|
uniqueID = uniqueID,
|
|
MANAG_INTENSITY = MANAG_INTENSITY,
|
|
SCA_SceneDrougth = scale_values(meanDR)# ,
|
|
# SCA_GrowthStability = scale_values(STABILITY)
|
|
)
|
|
SCA_ES_df <- plyr::join(SCA_ES_df, SCA_meanDR)
|
|
|
|
SCA_STAB <-
|
|
services_DF |>
|
|
filter_at(vars(STABILITY), all_vars(!is.na(.))) |>
|
|
group_by(PLOT) |>
|
|
reframe(date = date,
|
|
PROJECTION = factor(PROJECTION, levels = c('BASE', '245', '585'),
|
|
labels = c('Base', '245', '585')),
|
|
MANAGEMENT = MANAGEMENT,
|
|
uniqueID = uniqueID,
|
|
MANAG_INTENSITY = MANAG_INTENSITY,
|
|
SCA_GrowthStability = scale_values(STABILITY))
|
|
SCA_ES_df <- plyr::join(SCA_ES_df, SCA_STAB)
|
|
|
|
#
|
|
|
|
'
|
|
ggplot(SCA_ES_df, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_line(aes(y = SCA_PlotTreeBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df, aes(x = date, group = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_PlotTreeBiomass, colour = MANAG_INTENSITY)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df, aes(x = date, group = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_SceneDrougth, colour = MANAG_INTENSITY)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df, aes(x = date, group = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_GrowthStability, colour = MANAG_INTENSITY)) +
|
|
theme_bw()
|
|
'
|
|
|
|
'write.csv(SCA_ES_df, "SCALED_ECOSER_yearly.csv", row.names = F)'
|
|
|
|
# Rolling means
|
|
SCA_ES_df_rolling <- SCA_ES_df |>
|
|
group_by(uniqueID) |>
|
|
reframe(date = date,
|
|
PLOT = PLOT,
|
|
PROJECTION = first(PROJECTION),
|
|
MANAG_INTENSITY = first(MANAG_INTENSITY),
|
|
across(starts_with("SCA") & where(~sum(!is.na(.)) > 0),
|
|
lag_apply, n = movingWINDOW,
|
|
callback = function(x) {mean(x)}
|
|
))
|
|
|
|
# Remove NAs created by rolling mean
|
|
SCA_ES_df_rolling <- SCA_ES_df_rolling[!is.na(SCA_ES_df_rolling$SCA_PlotTreeBiomass),]
|
|
summary(SCA_ES_df_rolling)
|
|
|
|
'
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_line(aes(y = SCA_PlotTreeBiomass, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, colour = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION, drop = F) +
|
|
geom_line(aes(y = SCA_PlotTreeBiomass), alpha = .5) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, colour = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION, drop = F) +
|
|
geom_line(aes(y = SCA_N_BST), alpha = .5) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = MANAG_INTENSITY)) +
|
|
facet_grid(PLOT ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_CWD, colour = MANAG_INTENSITY)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PLOT)) +
|
|
facet_grid(MANAG_INTENSITY ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_CWD, colour = PLOT)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PLOT)) +
|
|
facet_grid(MANAG_INTENSITY ~ PROJECTION) +
|
|
geom_line(aes(y = SCA_GrowthStability, colour = PLOT)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_line(aes(y = SCA_GrowthStability, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
facet_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_line(aes(y = SCA_SceneDrougth, colour = PROJECTION)) +
|
|
theme_bw()
|
|
|
|
'
|
|
|
|
'write.csv(SCA_ES_df_rolling, "SCALED_ECOSER_yearly_rolling.csv", row.names = F)'
|
|
|
|
# Average all variables during the period
|
|
SCA_ES_df2 <- SCA_ES_df |>
|
|
group_by(uniqueID) |>
|
|
reframe(PROJECTION = first(PROJECTION),
|
|
MANAG_INTENSITY = first(MANAG_INTENSITY),
|
|
across(starts_with("SCA"), mean))
|
|
summary(SCA_ES_df2)
|
|
|
|
'# write.csv(SCA_ES_df2, "SCALED_ECOSER_mean.csv", row.names = F)'
|
|
|
|
cutoff = SCA_ES_df |>
|
|
select(starts_with("SCA")) |>
|
|
t() |> scale() |>
|
|
dist(diag = T) |> max()/2
|
|
|
|
SCA_ES_df |>
|
|
select(starts_with("SCA")) |>
|
|
t() |> scale() |>
|
|
dist(diag = T) |>
|
|
hclust() |>
|
|
# as.dendrogram() |>
|
|
plot() |> abline(h = cutoff, col = 'red')
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# Multivariate statistics
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# Pure EcoSer dataframe for MFA
|
|
|
|
'example for all averaged years'
|
|
'It is very important to keep the same order in the columns.
|
|
because variables must be ordered by group'
|
|
|
|
# Subset for MFA
|
|
ES_df <- SCA_ES_df[,c(
|
|
'uniqueID', 'PROJECTION', 'MANAG_INTENSITY', # Grouping variables
|
|
'SCA_PlotTreeBiomass', 'SCA_WoodAGB_V', 'SCA_WoodTrunkBiomass', # 3 until here, WP group
|
|
'SCA_RootBiomass', 'SCA_FineRootProduction', 'SCA_FallenLeafBiomass',
|
|
'SCA_currentSOC', # 4 until here, CS group
|
|
'SCA_N_Trees', 'SCA_CWD', 'SCA_DecayMortality', 'SCA_DRMortality', 'SCA_TotalMortality',# 5 until here, MIXED group
|
|
'SCA_N_SaplingSpecies', 'SCA_N_SaplingVolume', 'SCA_N_BST', 'SCA_N_DeadBST',
|
|
'SCA_N_TREMTrees', 'SCA_SpeciesShannonIndex' # 6 untill here, Hab group
|
|
)]
|
|
|
|
'DOWNLOAD SHINY PLOT'
|
|
|
|
ES_df$PROJECTION <- ordered(ES_df$PROJECTION,
|
|
levels = c('Base', '245', '585'),
|
|
labels = c('Base', '245', '585'))
|
|
|
|
# ES_df$MANAGEMENT <- as.factor(string_extractor(string_list2, 3))
|
|
ES_df2 <- as.data.frame(ES_df[, -1])
|
|
rownames(ES_df2) <- ES_df |> pull(uniqueID)
|
|
ES_df2
|
|
|
|
# Remove SCA_ prefix in column names
|
|
colnames(ES_df2) <- gsub("SCA_", "", colnames(ES_df2))
|
|
|
|
# MFA - Multiple Factor Analysis
|
|
# # PLOT #1 ----
|
|
ES_df.mfa <- FactoMineR::MFA(ES_df2,
|
|
group = c(1, 1, 3, 4, 5, 6),
|
|
type = c("n","n","s","s","s","s"),
|
|
ncp = 3,
|
|
name.group = c("Climate",
|
|
"Management",
|
|
"WoodProduction",
|
|
"SoilCarbon",
|
|
"All-purpose",
|
|
"Biodiversity")
|
|
)
|
|
# ------------------------------------------------------------------------------
|
|
# PARA HACER ESTO NO HAY QUE USAR LAS VARIABLES DE MANAGEMENT Y PROJECTION, Y
|
|
# PASARLAS DE FILAS A COLUMNAS (VER 'THE R OPUS' BOOK, CAPITULO MFA). ES DECIR,
|
|
# CADA MANAGEMENT O CADA PROJECTION SERIA UN GRUPO
|
|
|
|
ES_df.mfa$global.pca$var$cor |>
|
|
as_tibble(rownames = "descriptor") |>
|
|
tidyr::separate(descriptor,
|
|
into = c("descriptor", "PLOT", "PROJECTION"),
|
|
sep = "\\.") |>
|
|
ggplot(aes(x = Dim.1, y = Dim.2, color = descriptor)) +
|
|
geom_vline(xintercept = 0, linewidth = 1/10) +
|
|
geom_hline(yintercept = 0, linewidth = 1/10) +
|
|
geom_segment(aes(xend = 0, yend = 0),
|
|
arrow = arrow(length = unit(0.02, units = "in"), ends = "first"),
|
|
linewidth = 1/3, alpha = .5) +
|
|
ggforce::geom_circle(aes(x0 = 0, y0 = 0, r = 1),
|
|
linewidth = 1/10, linetype = 3, color = "grey") +
|
|
coord_equal() +
|
|
paletteer::scale_color_paletteer_d("ggthemes::Classic_20") +
|
|
scale_x_continuous(breaks = c(-1, 0, 1)) +
|
|
scale_y_continuous(breaks = c(-1, 0, 1)) +
|
|
theme_classic() +
|
|
theme(legend.position = "none", strip.text = element_text(size = 6)) +
|
|
|
|
facet_wrap(~descriptor)
|
|
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# Per group
|
|
factoextra::fviz_contrib(ES_df.mfa, "group", axes = 1)
|
|
factoextra::fviz_contrib(ES_df.mfa, "group", axes = 2)
|
|
|
|
print(ES_df.mfa)
|
|
plot(
|
|
ES_df.mfa,
|
|
lab.ind = F,
|
|
partial = c("PYR22_BASE_03", "PYR22_585_03",
|
|
"PYR22_BASE_01", "PYR22_585_01"),
|
|
repel = T
|
|
)
|
|
plot(ES_df.mfa, select = "cos2 0.3", repel = T)
|
|
|
|
# Eigenvectors plot histogram
|
|
factoextra::fviz_screeplot(ES_df.mfa)
|
|
|
|
# Results by group
|
|
groups <- factoextra::get_mfa_var(ES_df.mfa, "group")
|
|
groups$coord
|
|
groups$cos2
|
|
groups$contrib
|
|
|
|
# general plot
|
|
# PLOT #2 ----
|
|
factoextra::fviz_mfa_var(ES_df.mfa, repel = T)
|
|
|
|
# Quantitative variables
|
|
quanti.var <- factoextra::get_mfa_var(ES_df.mfa, "quanti.var")
|
|
quanti.var$contrib
|
|
|
|
# factoextra::fviz_mfa_var(ES_df.mfa, "quanti.var",
|
|
# col.var.sup = "violet", repel = T,
|
|
# geom = c("point", "text"), legend = "bottom")
|
|
|
|
# PLOT #3 A and B ----
|
|
factoextra::fviz_contrib(ES_df.mfa, choice = "quanti.var",
|
|
axes = 1, top = 10)
|
|
|
|
factoextra::fviz_contrib(ES_df.mfa, choice = "quanti.var",
|
|
axes = 2, top = 10)
|
|
|
|
# factoextra::fviz_mfa_var(ES_df.mfa, "quanti.var",
|
|
# col.var = "contrib", gradient.cols = "Blues",
|
|
# col.var.sup = "violet", repel = T,
|
|
# geom = c("point", "text"), legend = "bottom")
|
|
|
|
# Look at cos2 which is the real importance of the variables contribution
|
|
# factoextra::fviz_mfa_ind(ES_df.mfa, col.ind = "cos2",
|
|
# gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
|
|
# repel = T,
|
|
# # geom = c("point", "text"),
|
|
# partial = "all")
|
|
|
|
# contribution on the 1st and 2nd axis
|
|
factoextra::fviz_mfa_axes(ES_df.mfa,
|
|
partial = "all")
|
|
|
|
'factoextra::fviz_mfa_ind(ES_df.mfa,
|
|
habillage = c("MANAG_INTENSITY", "PROJECTION"),
|
|
addEllipses = T, ellipse.type = "confidence",
|
|
repel = T)'
|
|
|
|
# PLOT #4 ----
|
|
factoextra::fviz_mfa_ind(ES_df.mfa,
|
|
habillage = "MANAG_INTENSITY",
|
|
geom = c("point"),
|
|
addEllipses = T, ellipse.type = "confidence",
|
|
repel = T, palette = viridis::plasma(17, direction = -1))
|
|
|
|
# PLOT #5 ----
|
|
factoextra::fviz_mfa_ind(ES_df.mfa,
|
|
habillage = "PROJECTION",
|
|
geom = c("point"),
|
|
addEllipses = T, ellipse.type = "confidence",
|
|
repel = T)
|
|
|
|
factoextra::fviz_mfa_quali_biplot(ES_df.mfa, repel = FALSE, col.var = "#E7B800",
|
|
habillage = "MANAG_INTENSITY", addEllipses = TRUE, ellipse.level = 0.95)
|
|
|
|
# N_Trees hlines ----
|
|
hlines_N_Trees <- SCA_ES_df_rolling |>
|
|
group_by(PLOT, MANAG_INTENSITY, PROJECTION) |>
|
|
summarize(INT = mean(SCA_N_Trees))
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
lemon::facet_rep_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_hline(data = hlines_N_Trees,
|
|
aes(yintercept = hlines_N_Trees$INT, colour = PROJECTION),
|
|
alpha = .5) +
|
|
geom_line(aes(y = SCA_N_Trees, colour = PROJECTION), size = 1.5) +
|
|
scale_color_manual(values = c('#d4cfa6ff', 'orange', '#6e2878ff')) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = "Trees (n / Ha)", color = "Projection")
|
|
|
|
######
|
|
|
|
# Ptreebiomass hlines ----
|
|
hlines_PlotTreeBiomass <- SCA_ES_df_rolling |>
|
|
group_by(PLOT, MANAG_INTENSITY, PROJECTION) |>
|
|
summarize(INT = mean(SCA_PlotTreeBiomass))
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
lemon::facet_rep_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_hline(data = hlines_PlotTreeBiomass,
|
|
aes(yintercept = hlines_PlotTreeBiomass$INT, colour = PROJECTION),
|
|
alpha = .5, linetype = 2) +
|
|
geom_line(aes(y = SCA_PlotTreeBiomass, colour = PROJECTION), size = 1.5) +
|
|
scale_color_manual(values = c('#d4cfa6ff', 'orange', '#6e2878ff')) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = bquote('Basal Area' ~ (m^2/Ha)), color = "Projection")
|
|
######
|
|
|
|
# CWD hlines ----
|
|
hlines_SCA_CWD <- SCA_ES_df_rolling |>
|
|
group_by(PLOT, MANAG_INTENSITY, PROJECTION) |>
|
|
summarize(INT = mean(SCA_CWD))
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
lemon::facet_rep_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_hline(data = hlines_SCA_CWD,
|
|
aes(yintercept = hlines_SCA_CWD$INT, colour = PROJECTION),
|
|
alpha = .5, linetype = 2) +
|
|
geom_line(aes(y = SCA_CWD, colour = PROJECTION), size = 1.5) +
|
|
scale_color_manual(values = c('#d4cfa6ff', 'orange', '#6e2878ff')) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = bquote('Coarse Woody Debris' ~ (kg/Ha)), color = "Projection")
|
|
######
|
|
|
|
# MORTALITY hlines ----
|
|
hlines_DRMortality <- SCA_ES_df_rolling |>
|
|
group_by(PLOT, MANAG_INTENSITY, PROJECTION) |>
|
|
summarize(INT = mean(SCA_DRMortality))
|
|
|
|
ggplot(SCA_ES_df_rolling, aes(x = date, group = PROJECTION)) +
|
|
lemon::facet_rep_grid(PLOT ~ MANAG_INTENSITY) +
|
|
geom_hline(data = hlines_DRMortality,
|
|
aes(yintercept = hlines_DRMortality$INT, colour = PROJECTION),
|
|
alpha = .5) +
|
|
geom_line(aes(y = SCA_DRMortality, colour = PROJECTION), size = 1.5) +
|
|
scale_color_manual(values = c('#d4cfa6ff', 'orange', '#6e2878ff')) +
|
|
theme_bw() +
|
|
theme(axis.text.x = element_text(angle = 90),
|
|
panel.grid = element_blank()) +
|
|
labs(x = "Year", y = "Trees (n / Ha)", color = "Projection")
|
|
|
|
######
|
|
|
|
res <- cor(ES_df2)
|
|
col<- colorRampPalette(c("blue", "white", "red"))(20)
|
|
heatmap(x = res, col = col, symm = TRUE)
|
|
|
|
|
|
######
|
|
'factoextra::fviz_ellipses(ES_df.mfa,
|
|
c("MANAG_INTENSITY", "PROJECTION"),
|
|
repel = T)'
|
|
|
|
|
|
group.mfa.res <- factoextra::get_mfa_var(ES_df.mfa, element = "group")
|
|
FactoMineR::HCPC(group.mfa.res)
|
|
|
|
# HCPC - Hierarchical Clustering on Principal Components
|
|
FactoMineR::HCPC(ES_df.mfa,
|
|
nb.clust = 3, # either 0 or -1 to cut the tree
|
|
min = 3, max = NULL, graph = TRUE)
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# ------------------------------------------------------------------------------
|
|
# Hierarchichal MF
|
|
# Define groups == The plots
|
|
# hierar <- list()
|
|
'esto no es lo que busco'
|
|
# for (site in 1:length(unique(string_extractor(string_list, 2)))) {
|
|
#
|
|
# hierar[[site]] <- which(string_extractor(string_list, 2) %in%
|
|
# unique(string_extractor(string_list, 2))[site])
|
|
# # names(hierar)[[site]] <- unique(string_extractor(string_list, 2))[site]
|
|
#
|
|
# }
|
|
#
|
|
# res.hmfa <- HMFA(ES_df2, H = hierar,
|
|
# group = c(1, 1, 3, 4, 2, 6),
|
|
# type = c("n","n","c","c","c","c"),
|
|
# name.group = c("Climate",
|
|
# "Management",
|
|
# "WoodProduction",
|
|
# "SoilCarbon",
|
|
# "All-purpose",
|
|
# "Biodiversity"))
|
|
#
|
|
#
|
|
# ES_df.hmfa <- FactoMineR::HMFA(ES_df2,
|
|
# group = c(1, 1, 3, 4, 2, 6),
|
|
# type = c("n","n","c","c","c","c"),
|
|
# name.group = c("Climate",
|
|
# "Management",
|
|
# "WoodProduction",
|
|
# "SoilCarbon",
|
|
# "All-purpose",
|
|
# "Biodiversity")
|
|
# )
|
|
# print(ES_df.hmfa)
|
|
|
|
# ------------------------------------------------------------------------------
|
|
# ------------------------------------------------------------------------------
|
|
|
|
# Summed EcoSer
|
|
ES_df3 <- ES_df |> group_by(uniqueID) |>
|
|
reframe(
|
|
PROJECTION = PROJECTION,
|
|
MANAG_INTENSITY = MANAG_INTENSITY,
|
|
WoodProduction = mean(c(SCA_PlotTreeBiomass, SCA_WoodAGB_V,
|
|
SCA_WoodTrunkBiomass)),
|
|
SoilCarbon = mean(c(SCA_RootBiomass, SCA_FineRootProduction,
|
|
SCA_FallenLeafBiomass, SCA_currentSOC)),
|
|
All_purpose = mean(c(SCA_N_Trees, SCA_CWD)),
|
|
Biodiversity = mean(c(SCA_N_SaplingSpecies, SCA_N_SaplingVolume,
|
|
SCA_N_BST, SCA_N_DeadBST,
|
|
SCA_N_TREMTrees, SCA_SpeciesShannonIndex))
|
|
)
|
|
|
|
ES_df.famd <- as.data.frame(ES_df3[, -1])
|
|
rownames(ES_df.famd) <- ES_df3 |> pull(uniqueID)
|
|
|
|
# FAMD approach (that is, as above but without groups)
|
|
'Cannot depict both individuals and variables'
|
|
|
|
res.famd <- FactoMineR::FAMD(ES_df.famd)
|
|
|
|
var <- factoextra::get_famd_var(res.famd)
|
|
var
|
|
|
|
factoextra::fviz_famd_var(res.famd)
|
|
factoextra::fviz_famd_var(res.famd, "quanti.var", repel = T)
|
|
factoextra::fviz_famd_var(res.famd, "quali.var", repel = T)
|
|
|
|
factoextra::fviz_famd_ind(res.famd)
|
|
factoextra::fviz_famd_ind(res.famd, col.ind = 'cos2')
|
|
factoextra::fviz_famd_ind(res.famd,
|
|
habillage = "MANAG_INTENSITY",
|
|
addEllipses = T)
|
|
|
|
# X -1. VARIABLE CONTRIBUTION TO ECOS-SER -------------------------------------
|
|
|
|
Quitar ES relacionados con plantulas y repetir
|
|
|
|
modelo lineal o matriz correlaciones entre grupos y ES
|
|
|
|
|
|
# X. PARETO SOLUTIONS -----
|
|
# (before this point (16-01-2026)) use laptop data -----
|
|
|
|
# X.1. Define objectives or "solutions" for the problem
|
|
# 1) diversity 2) C sink 3) habitat protection 4) Wood production loss
|
|
|
|
"hacer un dendrocorrelograma/PA para ver qué variable (EcoSer) estás más relacionado o es
|
|
representativo de las soluciones propuestas.
|
|
|
|
Según VPLAS en PNAS hay que estandarizar (min-max) los ES y hacer un sumatorio de ellos por
|
|
tipo de ES (Provisioning, Supporting, Regulating or Cultural)"
|
|
|
|
# -----
|
|
|
|
library(GPareto)
|
|
|
|
'test' ("https://cran.r-project.org/web/packages/GPareto/refman/GPareto.html")
|
|
|
|
# Example 1 : Surrogate-based multi-objective Optimization with postprocessing
|
|
|
|
set.seed(25468)
|
|
|
|
d <- 2
|
|
fname <- P2
|
|
|
|
plotParetoGrid(P2) # For comparison
|
|
|
|
# Optimization
|
|
budget <- 25
|
|
lower <- rep(0, d)
|
|
upper <- rep(1, d)
|
|
|
|
omEGO <- easyGParetoptim(fn = fname, budget = budget, lower = lower, upper = upper)
|
|
|
|
# Postprocessing
|
|
plotGPareto(omEGO, add= FALSE, UQ_PF = TRUE, UQ_PS = TRUE, UQ_dens = TRUE)
|
|
|
|
# Example 2 : Surrogate-based multi-objective Optimization including a cheap function
|
|
set.seed(42)
|
|
|
|
d <- 2
|
|
|
|
fname <- P1
|
|
n.grid <- 19
|
|
test.grid <- expand.grid(seq(0, 1, length.out = n.grid), seq(0, 1, length.out = n.grid))
|
|
nappr <- 15
|
|
design.grid <- maximinESE_LHS(lhsDesign(nappr, d, seed = 42)$design)$design
|
|
response.grid <- t(apply(design.grid, 1, fname))
|
|
|
|
mf1 <- km(~., design = design.grid, response = response.grid[,1])
|
|
mf2 <- km(~., design = design.grid, response = response.grid[,2])
|
|
model <- list(mf1, mf2)
|
|
|
|
nsteps <- 1
|
|
lower <- rep(0, d)
|
|
upper <- rep(1, d)
|