###Exercise 1
options(scipen = 50) library(dplyr) library(tidyr) library(ggplot2)
luc<-read_csv(“EDGC/LCC_Chaco_Numbers_Updated.csv”)#import of data set under land_cover name luc_su <- luc %>% # Summarize per ecoregion –> convert right away into ha to match the units of the other data group_by(Ecoregion, Year) %>% summarise(WL_P_ha = sum(WL_P_km2)100, WL_C_ha = sum(WL_C_km2)100, P_C_ha = sum(P_C_km2) *100) %>% ungroup() %>% arrange(Year, Ecoregion) # Print the table head(luc_su,10) # Define the simulation period startYear <- 1986 endYear <- 2050 period <- seq(startYear, endYear, 1) # Define variables for the carbon stocks # Wet Chaco wet_wl <- 104.2 # Carbon in woodlands (Mg C/ha) wet_crop <- 5 # Carbon in crops (Mg C/ha) wet_past_Tree <- 28 # Carbon in pastures with trees (Mg C/ha) wet_past_noTree <- 13.5 # Carbon in pastures without trees (Mg C/ha) # Dry Chaco dry_wl <- 60.8 dry_crop <- 5 dry_past_Tree <- 28 dry_past_noTree <- 6.5 # Very dry Chaco vdry_wl <- 42.6 vdry_crop <- 5 vdry_past_Tree <- 28 vdry_past_noTree <- 6.1 agb_func <- function(lcTS, simPeriod, agb_fromLC, agb_toLC, ecoregion, conversion){ ########################################################################### # lcTS: vector with the estimates of the land cover conversions # simPeriod: vector with a sequence of years for the simulation period # agb_fromLC: AGB stocks in the land cover type that is converted from # agb_toLC: AGB stocks in the land cover type that is converted to # ecoregion: name of the ecoregion, has to be one of “Wet Chaco”, “Dry Chaco”, “Very Dry Chaco” # conversion: name of the conversion, has to be one of “WL_P”, “WL_C”, “P_C” ###########################################################################
#### Define the proportions of AGB that goes into the different stocks p_burn <- 0.55 # Vegetation burned on site p_slash <- 0.35 # Vegetation left on site after clearing p_wood <- 0.08 # Vegetation taken out for forest products p_ret <- 0.02 # Fraction returning to soil as mineral carbon after burning
#### Instantiate the data frame df <- data.frame(matrix(nrow = length(simPeriod), ncol = 0)) df\(year <- period # For this implementation, we assume that no further land-cover change occurred after 2022. Thus, for all years post 2023 we define the conversion to be zero df\)lcc <- c(lcTS, rep(0, length(period) - length(lcTS)))
#### Fill the data frame with values # Losses df\(loss_AGB_tot <- df\)lcc * (agb_fromLC - agb_toLC) df\(loss_burn <- p_burn * df\)loss_AGB_tot df\(loss_slash <- p_slash * df\)loss_AGB_tot df\(loss_wood <- p_wood * df\)loss_AGB_tot df\(loss_ret <- p_ret * df\)loss_AGB_tot # Pools: df\(pool_slash <- 0 df\)pool_wood <- 0 df\(pool_ret <- 0 # Emissions df\)em_slash <- 0 df\(em_burn <- 0 #df\)loss_burn # immediate emissions, we can df\(em_wood <- 0 df\)em_ret <- 0
### Iterate over the years and distribute the pools & emissionss for (i in 1:nrow(df)) { # Calculate the emissions df\(em_burn[i] <- df\)loss_burn[i] # Emissions from burn happen entirely in this year df\(em_slash[i] <- sum(df\)loss_slash[i:ifelse(i > 2, i-1, 1)]) * 0.5 # Emissions from slash happen within 2 years df\(em_wood[i] <- sum(df\)loss_wood[i:ifelse(i > 10, i-9, 1)]) * 0.1 # Emissions from wood products happen within 10 years df\(em_ret[i] <- sum(df\)loss_ret[i:ifelse(i > 100, i-99, 1)]) * 0.01 # Emissions from carbon returning to soil as mineral carbon after burning # Fill the pools df\(pool_slash[i] <- ifelse(i == 1, 0, df\)pool_slash[i-1]) + df\(loss_slash[i] - df\)em_slash[i] df\(pool_wood[i] <- ifelse(i == 1, 0, df\)pool_wood[i-1]) + df\(loss_wood[i] - df\)em_wood[i] df\(pool_ret[i] <- ifelse(i == 1, 0, df\)pool_ret[i-1]) + df\(loss_ret[i] - df\)em_ret[i] } # Round everything to 3 digits df <- df %>% round(2) # Calculate the total emissions df\(em_tot <- df\)em_burn + df\(em_slash + df\)em_wood + df\(em_ret # Add the ecoregion and conversion type to the data frame df\)Ecoregion <- ecoregion df$Conversion <- conversion # Return the data frame return(df) } #### Wet Chaco # Woodland to pasture wet_wl_p <- luc_su %>% filter(Ecoregion == “Wet Chaco”) %>% select(WL_P_ha) %>% pull() wet_wl_p <- agb_func(lcTS = wet_wl_p, simPeriod = period, agb_fromLC = wet_wl, agb_toLC = wet_past_Tree, ecoregion = “Wet Chaco”, conversion = “Woodland to Pasture”) # Woodland to cropland wet_wl_c <- luc_su %>% filter(Ecoregion == “Wet Chaco”) %>% select(WL_C_ha) %>% pull() wet_wl_c <- agb_func(lcTS = wet_wl_c, simPeriod = period, agb_fromLC = wet_wl, agb_toLC = wet_crop, ecoregion = “Wet Chaco”, conversion = “Woodland to Cropland”) # Pasture to cropland wet_p_c <- luc_su %>% filter(Ecoregion == “Wet Chaco”) %>% select(P_C_ha) %>% pull() wet_p_c <- agb_func(lcTS = wet_p_c, simPeriod = period, agb_fromLC = wet_past_Tree, agb_toLC = wet_crop, ecoregion = “Wet Chaco”, conversion = “Pasture to Cropland”) # Rbind the tables wet_all <- rbind(wet_wl_p, wet_wl_c, wet_p_c)
head(wet_all, 25)
dry_wl_p <- luc_su %>% filter(Ecoregion == “Dry Chaco”) %>% select(WL_P_ha) %>% pull() dry_wl_p <- agb_func(lcTS = dry_wl_p, simPeriod = period, agb_fromLC = dry_wl, agb_toLC = dry_past_Tree, ecoregion = “Dry Chaco”, conversion = “Woodland to Pasture”) # Woodland to cropland dry_wl_c <- luc_su %>% filter(Ecoregion == “Dry Chaco”) %>% select(WL_C_ha) %>% pull() dry_wl_c <- agb_func(lcTS = dry_wl_c, simPeriod = period, agb_fromLC = dry_wl, agb_toLC = dry_crop, ecoregion = “Dry Chaco”, conversion = “Woodland to Cropland”) # Pasture to cropland dry_p_c <- luc_su %>% filter(Ecoregion == “Dry Chaco”) %>% select(P_C_ha) %>% pull() dry_p_c <- agb_func(lcTS = dry_p_c, simPeriod = period, agb_fromLC = dry_past_Tree, agb_toLC = dry_crop, ecoregion = “Dry Chaco”, conversion = “Pasture to Cropland”) # Rbind the tables dry_all <- rbind(dry_wl_p, dry_wl_c, dry_p_c)
vdry_wl_p <- luc_su %>% filter(Ecoregion == “Very dry Chaco”) %>% select(WL_P_ha) %>% pull() vdry_wl_p <- agb_func(lcTS = vdry_wl_p, simPeriod = period, agb_fromLC = vdry_wl, agb_toLC = vdry_past_Tree, ecoregion = “Very Dry Chaco”, conversion = “Woodland to Pasture”) # Woodland to cropland vdry_wl_c <- luc_su %>% filter(Ecoregion == “Very dry Chaco”) %>% select(WL_C_ha) %>% pull() vdry_wl_c <- agb_func(lcTS = vdry_wl_c, simPeriod = period, agb_fromLC = vdry_wl, agb_toLC = vdry_crop, ecoregion = “Very Dry Chaco”, conversion = “Woodland to Cropland”) # Pasture to cropland vdry_p_c <- luc_su %>% filter(Ecoregion == “Very dry Chaco”) %>% select(P_C_ha) %>% pull() vdry_p_c <- agb_func(lcTS = vdry_p_c, simPeriod = period, agb_fromLC = vdry_past_Tree, agb_toLC = vdry_crop, ecoregion = “Very Dry Chaco”, conversion = “Pasture to Cropland”) # Rbind the tables vdry_all <- rbind(vdry_wl_p, vdry_wl_c, vdry_p_c) # Rbind all tables, only keep the columns we need em_all <- rbind(wet_all, dry_all, vdry_all) %>% select(year, Ecoregion, Conversion, em_tot)
em_all %>% ggplot() + geom_line(aes(x=year, y=em_tot / 1000, color=Conversion), linewidth=1) + facet_wrap(~Ecoregion) + scale_x_continuous(breaks = seq(1985, 2050, 15)) + scale_y_continuous(name = “Emissions from AGB [thsd. t]”) + theme_light() + theme(axis.title.x = element_blank(), axis.title.y = element_text(size=12, vjust=1.5), axis.text.y = element_text(size=12), axis.text.x = element_text(size=10), legend.position = “bottom”, legend.title = element_blank(), legend.text = element_text(size=12), panel.grid.minor = element_blank()) + geom_vline(xintercept = 2023, linetype=“dashed”, color = “black”)
tot_emission <- em_all %>% group_by(Conversion, Ecoregion) %>% summarise(em_tot = sum(em_tot)) %>% ungroup() tot_emission
This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
summary(cars)
## speed dist
## Min. : 4.0 Min. : 2.00
## 1st Qu.:12.0 1st Qu.: 26.00
## Median :15.0 Median : 36.00
## Mean :15.4 Mean : 42.98
## 3rd Qu.:19.0 3rd Qu.: 56.00
## Max. :25.0 Max. :120.00
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the
code chunk to prevent printing of the R code that generated the
plot.