# 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)