Revisting Disturbance Treatments

Objective

Previously when we tried perturbing only the carbon storage pool as we expected to mimic girdling treatments there was little response to the treatment. And little to no difference between the different treatments. Now let’s try a sustained disturbance to the carbon storage pool for 1 month.

library(magrittr)
library(ggplot2)
library(ed4forte)
source(file.path(here::here('C.analysis', 'packages.R')))
INPUT_DIR <- here::here('ED-outputs', 'disturbanceII')
source(file.path(here::here('C.analysis', 'functions.R')))
THEME <- theme_bw()

files <- list.files(INPUT_DIR, '.rds', full.names = TRUE)
objs  <- lapply(files, readRDS)
names(objs) <- gsub(pattern = '.rds', replacement = '', x = basename(files))
NPP_data <- dplyr::bind_rows(mapply(function(obj, scn){
  process_NPP(ed_object = obj, scn_name = scn)
  }, obj = objs, scn = names(objs), SIMPLIFY = FALSE))
AGB_data <- dplyr::bind_rows(mapply(function(obj, scn){
  process_ABG(ed_object = obj, scn_name = scn)
  }, obj = objs, scn = names(objs), SIMPLIFY = FALSE))
nplant_data <- dplyr::bind_rows(mapply(function(obj, scn){
  
  process_nplant <- function(ed_object = obj, scn_name = scn){
    
    as.data.table(ed_object$df_pft[[1]])[  ,list(datetime, NPLANT_PY, pft)] %>%  
      .[ , list(value = mean(NPLANT_PY)), by=c("datetime")] %>% 
      .[ , scn := scn_name]
  }
    process_nplant(ed_object = obj, scn_name = scn)
  
  
}, obj = objs, scn = names(objs), SIMPLIFY = FALSE))
# mort_data <- dplyr::bind_rows(mapply(function(obj, scn){
#   
#   process_nplant <- function(ed_object = obj, scn_name = scn){
#     
#     as.data.table(ed_object$df_cohort[[1]])[  ,list(datetime, NPLANT_PY, pft)] %>%  
#       .[ , list(value = mean(NPLANT_PY)), by=c("datetime")] %>% 
#       .[ , scn := scn_name]
# 
#   }
#     process_nplant(ed_object = obj, scn_name = scn)
#   
#   
# }, obj = objs, scn = names(objs), SIMPLIFY = FALSE))

What about mortality?

Can we tell anything about tree morality based on the plant density? Is this something that is useful to look at?

nplant_data_baseline <- nplant_data[scn == 'baseline', ]
nplant_data <- nplant_data[scn != 'baseline', ]
nplant_data$type <- gsub(pattern = 'harvest-45_|harvest-85_', replacement = '', x = nplant_data$scn)
nplant_data$scn <- gsub(pattern = paste(unique(nplant_data$type), collapse = '|'), 
                        replacement = '', x = nplant_data$scn)
ggplot() + 
  geom_line(data = nplant_data_baseline, aes(datetime, value), color = 'grey') + 
  geom_line(data = nplant_data, aes(datetime, value, color = scn)) + 
  THEME + 
  facet_wrap('type') + 
  labs(title = 'Plant Density', caption = 'baseline = grey', y = "Nplants/m2")

Above Ground Biomass

AGB <- AGB_data[ scn != 'baseline',  sum(value), by=c("datetime","scn")]
names(AGB) <- c('datetime', 'scn', 'value')
AGB$type <- gsub(pattern = 'harvest-45_|harvest-85_', replacement = '', x = AGB$scn)
AGB$scn <- gsub(pattern = paste(unique(AGB$type), collapse = '|'), replacement = '', x = AGB$scn)
ABG_baseline <- AGB_data[ scn == 'baseline',  sum(value), by=c("datetime","scn")]
names(ABG_baseline) <- c('datetime', 'scn', 'value')
ggplot() + 
  geom_line(data = ABG_baseline, aes(datetime, value), color = 'grey') + 
  geom_line(data = AGB, aes(datetime, value, color = scn)) + 
  THEME + 
  facet_wrap('type') + 
  labs(title = 'Above Ground Biomass', caption = 'baseline = grey', y = "kgC/m2")

NPP

NPP <- NPP_data[ scn != 'baseline',  sum(value), by=c("year","scn")]
names(NPP) <- c('year', 'scn', 'value')
NPP$type <- gsub(pattern = 'harvest-45_|harvest-85_', replacement = '', x = NPP$scn)
NPP$scn <- gsub(pattern = paste(unique(NPP$type), collapse = '|'), replacement = '', x = NPP$scn)
NPP_baseline <- NPP_data[ scn == 'baseline',  sum(value), by=c("year","scn")]
names(NPP_baseline) <- c('year', 'scn', 'value')
ggplot() + 
  geom_line(data = NPP_baseline, aes(year, value), color = 'grey') + 
  geom_line(data = NPP, aes(year, value, color = scn)) + 
  THEME + 
  facet_wrap('type') + 
  labs(title = 'Net Primary Production', caption = 'baseline = grey', y = "kgC/m2/yr")

# What is the difference the NPP of the different treatments? 
names(NPP) <- c('year', 'scn', 'treatment_value', 'type')
names(NPP_baseline) <- c('year', 'scn', 'baseline_value')
NPP_difference <- NPP[NPP_baseline[, list(year, baseline_value)], on = 'year']
NPP_difference$diff <- NPP_difference$treatment_value - NPP_difference$baseline_value

ggplot() + 
  geom_line(data = NPP_difference, aes(year, diff, color = scn)) + 
  THEME + 
  facet_wrap('type') + 
  labs(title = 'Net Primary Production', y = "kgC/m2/yr")