## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(rsofun)
library(dplyr)
library(ggplot2)

## ----eval = FALSE-------------------------------------------------------------
# library(rsofun)
# 
# # Output obtained from simulating biomee_p_model_luluc_drivers
# biomee_p_model_luluc_output

## ----eval = FALSE-------------------------------------------------------------
# library(rsofun)
# 
# biomee_p_model_luluc_drivers

## ----eval = TRUE--------------------------------------------------------------
df_drivers <- biomee_p_model_luluc_drivers
df_drivers$params_siml[[1]]$nyeartrend <- 300

## ----eval = TRUE--------------------------------------------------------------
df_drivers$init_lu[[1]]

## ----eval = TRUE--------------------------------------------------------------
df_drivers$luc_forcing[[1]]

## ----eval = TRUE, results = 'hide', warning = FALSE---------------------------
out_sim_1 <- runread_biomee_f(df_drivers)

## ----eval = TRUE, results = 'hide'--------------------------------------------
library(cowplot)
library(purrr)

plot1 <- function(lu_name, variable, out, y_limit=NA, yr_start=1, yr_label_offset=0) {
  tile <- out[[lu_name]][[1]]
  if(lu_name != 'aggregated'){
    tile <- tile$output_annual_tile
  } else {
    tile <- tile$output_annual_cell
  }
          
  if (variable != 'lu_fraction') {
    res <- tile %>%
          ggplot() +
          geom_line(aes(x = year + yr_label_offset,
                        y = get(variable) * lu_fraction)) +
          coord_cartesian(xlim = c(yr_start, NA),
                          ylim = c(0, y_limit)) +
          theme_classic() +
          labs(x = "Year", y = paste(variable, "(", lu_name, ")"))
  } else {
    stopifnot(variable == 'lu_fraction')
    res <- tile %>%
          ggplot() +
          geom_line(aes(x = year + yr_label_offset,
                        y = get(variable))) +
          coord_cartesian(xlim = c(yr_start, NA),
                          ylim = c(0, y_limit)) +
          theme_classic() +
          labs(x = "Year", y = paste(variable, "(", lu_name, ")"))
  }

  return(res)
}

plot_variable <- function(variable, out, yr_start=1, yr_label_offset=0) {
  agg <- out[['aggregated']][[1]]$output_annual_cell
  y_limit <- max(agg[variable] * agg$lu_fraction) * 1.01

  # We remove sitename and aggregated
  tile_names <- names(out)[3:length(names(out))]

  names <- c(tile_names, 'aggregated')

  names |> lmap(\(x) list(plot1(x, variable, out, y_limit, yr_start, yr_label_offset)))
}

plot_variables <- function(variables, out, yr_start=1, yr_label_offset=0) {

  plots <- variables |> lmap(\(x) plot_variable(x, out, yr_start, yr_label_offset))
  cowplot::plot_grid(plotlist = plots, ncol = length(variables))
}

## ----eval = TRUE, results = 'hide', warning = FALSE, fig.width=7, fig.height=7----
plot_variables(c('GPP', 'fastSOM', 'lu_fraction'), out_sim_1)

## ----eval = TRUE, results = 'hide', fig.width=7, fig.height=3-----------------
p1 <- out_sim_1$aggregated[[1]]$output_annual_cell %>%
        ggplot() +
        geom_line(aes(x = year, y = prod_pool_1_C)) +
        theme_classic() +
        labs(x = "Year", y = paste("Prod. pools C (2 years)"))

p2 <- out_sim_1$aggregated[[1]]$output_annual_cell %>%
        ggplot() +
        geom_line(aes(x = year, y = prod_pool_2_C)) +
        theme_classic() +
        labs(x = "Year", y = paste("Prod. pools C (20 years)"))
cowplot::plot_grid(p1, p2, nrow = 1)

## ----eval = TRUE--------------------------------------------------------------
df_drivers <- biomee_p_model_luluc_drivers

# We extend the transient period to 300 years
n_trans <- 300
df_drivers$params_siml[[1]]$nyeartrend <- n_trans

lu_defs <- tibble(
        name      = c('primary', 'secondary', 'urban'),
        fraction  = c(1.0, 0.0, 0.0),
        # We defined which LU can accept vegetation cohorts
        # Technically, this was not necessary since no 'lu_index' (see below) points to the third LU
        # (but it is a good idea to enforce this anyways)
        vegetated = c(TRUE, TRUE, FALSE)
        # Equivalently, presets could have been used instead:
        #preset = c('unmanaged', 'unmanaged', 'urban')
)
n_lu <- length(lu_defs$name)

# Initial transition
# 3 LUs, so each year land-use change is represented by 3x3 transitions.
# The first 3 numbers are the transitions from each LU to the first LU,
# the following three numbers are the transitions from each LU to the second LU, and so on.
initial_transition <- c(0, 0, 0, 0.4, 0, 0, 0, 0, 0)

# Harvest pattern: clear-cut secondary forest every 40 years by an area 
# corresponding to 10% of the tile (0.1). That corresponds to a fourth of the area 
# of secondary forest (0.4/4 = 0.1).
# Clear-cuts are represented by self-transitions of LU 2, of area 0.1.
harvest_pattern <- c(
        rep(rep(0, 9), 39), # null pattern
        c(0, 0, 0, 0, 0.1, 0, 0, 0, 0)
)
# We repeat the harvest pattern enough times
harvest_transitions <- rep(harvest_pattern, 10)

# We define a transition from primary to urban repeated every year
urban_transition <- rep(c(0, 0, 0, 0, 0, 0, 0.0006, 0, 0), n_trans)

# We add all the transition matrix
transitions_matrix <- build_luc_matrix(list(initial_transition, harvest_transitions, urban_transition), n_lu, n_trans)

# We set the LU definitions and transition forcing in the driver:
df_drivers$init_lu[[1]] <- lu_defs
df_drivers$luc_forcing[[1]]  <- transitions_matrix

# Finally, we configure two cohorts with two different species (see 'init_cohort_species'): one for tile 1 and one for tile 2 (see 'lu_index').
# LU 3 is an urban tile: it does not accept any cohort.
df_drivers$init_cohort[[1]] <- tibble(
        init_cohort_species = c(2, 3),      # species
        init_cohort_nindivs = rep(0.05, 2), # initial individual density, individual/m2 ! 1 indiv/m2 = 10.000 indiv/ha
        init_cohort_bl      = rep(0.0,  2), # initial biomass of leaves, kg C/individual
        init_cohort_br      = rep(0.0,  2), # initial biomass of fine roots, kg C/individual
        init_cohort_bsw     = rep(0.05, 2), # initial biomass of sapwood, kg C/individual
        init_cohort_bHW     = rep(0.0,  2), # initial biomass of heartwood, kg C/tree
        init_cohort_seedC   = rep(0.0,  2), # initial biomass of seeds, kg C/individual
        init_cohort_nsc     = rep(0.05, 2), # initial non-structural biomass
        lu_index            = c(1, 2)       # index land use (LU) containing this cohort. 0 (default) means any vegetated tile will contain a copy.
)

## ----eval = TRUE, results = 'hide', warning = FALSE, fig.width=7, fig.height=7----
out_sim_2 <- runread_biomee_f(df_drivers)
plot_variables(c('GPP', 'fastSOM', 'lu_fraction'), out_sim_2)

## ----eval = TRUE, results = 'hide', fig.width=7, fig.height=3-----------------
p1 <- out_sim_2$aggregated[[1]]$output_annual_cell %>%
        ggplot() +
        geom_line(aes(x = year, y = prod_pool_1_C)) +
        theme_classic() +
        labs(x = "Year", y = paste("Prod. pools C (2 years)"))

p2 <- out_sim_2$aggregated[[1]]$output_annual_cell %>%
        ggplot() +
        geom_line(aes(x = year, y = prod_pool_2_C)) +
        theme_classic() +
        labs(x = "Year", y = paste("Prod. pools C (20 years)"))

cowplot::plot_grid(p1, p2, nrow = 1)

## ----eval = TRUE--------------------------------------------------------------
df_drivers <- biomee_p_model_luluc_drivers

# We extend the transient period to 300 years
n_trans <- 300
df_drivers$params_siml[[1]]$nyeartrend <- n_trans

lu_defs <- tibble(
        name      = c('primary', 'crop', 'fallow'),
        fraction  = c(1.0, 0.0, 0.0),
        # Here we use presets to easily configure some parameters for the cropland
        # ('extra_N_input', 'extra_turnover_rate', 'oxidized_litter_fraction')
        # preset 'unmanaged' sets all these parameters to 0
        preset = c('unmanaged', 'cropland', 'unmanaged')
)
n_lu <- length(lu_defs$name)

# Initial transition
# 3 LUs, so each year land-use change is represented by 3x3 transitions.
# The first 3 numbers are the transitions from each LU to the first LU,
# the following three numbers are the transitions from each LU to the second LU, and so on.
initial_transition <- c(0, 0, 0, 0.66, 0, 0, 0.34, 0, 0)

# Fallow pattern:
fallow_pattern <- c(
        c(0, 0, 0, 0, 0, 0.33, 0, 0.33, 0),
        rep(0, 9) # null pattern
)
# We repeat the harvest pattern enough times
fallow_transitions <- rep(harvest_pattern, 150)

# We add all the transition matrix
transitions_matrix <- build_luc_matrix(list(initial_transition, fallow_transitions), n_lu, n_trans)

# We set the LU definitions and transition forcing in the driver:
df_drivers$init_lu[[1]] <- lu_defs
df_drivers$luc_forcing[[1]]  <- transitions_matrix

# We configure three cohorts
df_drivers$init_cohort[[1]] <- tibble(
        init_cohort_species = c(2, 1, 2),      # species
        init_cohort_nindivs = rep(0.05, 3), # initial individual density, individual/m2 ! 1 indiv/m2 = 10.000 indiv/ha
        init_cohort_bl      = rep(0.0,  3), # initial biomass of leaves, kg C/individual
        init_cohort_br      = rep(0.0,  3), # initial biomass of fine roots, kg C/individual
        init_cohort_bsw     = rep(0.05, 3), # initial biomass of sapwood, kg C/individual
        init_cohort_bHW     = rep(0.0,  3), # initial biomass of heartwood, kg C/tree
        init_cohort_seedC   = rep(0.0,  3), # initial biomass of seeds, kg C/individual
        init_cohort_nsc     = rep(0.05, 3), # initial non-structural biomass
        lu_index            = c(1, 2, 3)    # index land use (LU) containing this cohort. 0 (default) means any vegetated tile will contain a copy.
)