Making FullCAM outputs useful
A quick-reference quide
Charles Coverdale
17 June 2021
1 Introduction
FullCAM is the Australian Government model for estimating carbon sequestration around Australia. It runs as it’s own software program (which must be downloaded and installed on a computer), and could be described as clunky at best.
This code and methodology takes the FullCAM output and transforms it into something useful for estimating carbon sequestration of a given land area.
#Load up a bunch of packages
library(plyr)
library(dplyr)
library(vctrs)
library(tidyr)
library(ggplot2)
library(rgdal)
library(tmap)
library(ggmap)
library(dplyr)
library(sf)
library(ggspatial)
library(rlang)
library(broom)
library(tidyverse)
library(readxl)
library(purrr)
library(zoo)
library(RcppRoll)
library(rmarkdown)2 Importing data from FullCAM
First we need to import data from the FullCAM model. Luckily, the table that is produced when a scenario is run can be exported as a csv. Let’s load this in.
#Read in the exported FullCAM data
FullCAM_output <- read.csv(file="C:/06 R_code/Test_Data_Environmental_Plantings_Method.csv")
knitr::kable(head(FullCAM_output),digits=2, caption ="Raw FullCAM output","simple")| Year | Month | Day | Dec..Year | C.mass.of.trees…tC.ha. | C.mass.of.forest.debris…tC.ha. | CH4.emitted.due.to.fire..tCH4.ha. | N2O.emitted.due.to.fire..tN2O.ha. |
|---|---|---|---|---|---|---|---|
| 2021 | 6 | 30 | 2021.50 | 0.00 | 0.00 | NA | NA |
| 2021 | 7 | 31 | 2021.58 | 0.29 | 0.01 | 0 | 0 |
| 2021 | 8 | 31 | 2021.67 | 0.29 | 0.01 | 0 | 0 |
| 2021 | 9 | 30 | 2021.75 | 0.29 | 0.01 | 0 | 0 |
| 2021 | 10 | 31 | 2021.83 | 0.29 | 0.01 | 0 | 0 |
| 2021 | 11 | 30 | 2021.91 | 0.29 | 0.01 | 0 | 0 |
3 Data cleaning and analysis
The first thing we realise from this table is that FullCAM outputs this data in a terrible format. Let’s do some cleaning, fix up headings, and make the units more meaningful.
#Let's start by removing the 4th column with the clunky date
FullCAM_output_clean = subset(FullCAM_output, select = -c(4))
#Rename the remaining columns to something more useful
FullCAM_output_clean <- FullCAM_output_clean %>%
rename(
C_mass_tree_tC_ha = 4,
C_mass_forest_debris_tC_ha = 5,
CH4_emmited_fire_tCH4_ha = 6,
N2O_emmited_fire_tN2O_ha = 7)
#Now let's convert these measures of carbon into measures of CO2
FullCAM_output_transformed <- FullCAM_output_clean %>%
mutate(
Date = as.Date(with(FullCAM_output_clean, paste(Year, Month, Day,sep="-")), "%Y-%m-%d"),
#Convert tonnes of C into tonnes of CO2 by multiplying by the weight of a carbon dioxide molecule
CO2_stored_in_trees = C_mass_tree_tC_ha*(44/12),
CO2_stored_in_forest_debris = C_mass_forest_debris_tC_ha*(44/12),
#Add data sets together to get the total amount sequestered
Total_CO2_sequestered = CO2_stored_in_trees+CO2_stored_in_forest_debris,
#Find the marginal change in CO2 sequestration each month
CO2_sequestered_month = lead(Total_CO2_sequestered,n=1) - (Total_CO2_sequestered))
FullCAM_output_transformed <- na.omit(FullCAM_output_transformed)
#Now we want to find the marginal change per year, updated on a rolling monthly basis
FullCAM_output_transformed <- FullCAM_output_transformed %>%
mutate(
#Note the allign='right' bit is important here to make sure it starts the rollsum after the first year.
CO2_sequestered_year = RcppRoll::roll_sum(CO2_sequestered_month,12,fill=NA,align='right'))4 Visualising the sequestration rate
Now that we have clean data, let’s see what the sequestration rate of CO2 (rather than carbon) looks like over the crediting period.
plot_FullCAM_annual <- ggplot(data=FullCAM_output_transformed)+
geom_line(aes(x = Date, y = CO2_sequestered_year), col = "#3C33FE",size=1.5) +
#geom_area(aes(x = Date, y = CO2_sequestered_year), fill = "lightgrey",size=1) +
labs(title="Annual sequestration rate using the environmental plantings method",
subtitle = "Analysis from the Hunter Valley, New South Wales",
caption = "Data: FullCAM Model 2021",
x="",
y="Annual CO2 sequestration (tCO2/ha)")+
scale_x_date(date_breaks = "5 years",date_labels="%Y")+
scale_y_continuous(expand = c(0,0)) +
theme(legend.position="none")+
theme_minimal()+
theme(axis.line.x = element_line(color="black", size = 0.3),axis.line.y = element_line(color="black", size = 0.3))+
theme(axis.text=element_text(size=12))+
theme(panel.grid.minor = element_blank())+
theme(plot.title=element_text(face="bold",size=12))+
theme(plot.subtitle=element_text(size=11,margin=margin(0,0,15,0)))+
theme(plot.caption=element_text(size=8))
plot_FullCAM_annual
ggsave("C:/06 R_code/FullCAM_sequestration_plot.png",units="cm",dpi=200, device="png")We see that gum trees grow fasster and faster until around 10 years old, then the rate of growth slows.
5 Cumulative sequestration
Now we visualize the total carbon sequestration over time.
plot_FullCAM_cumulative <- ggplot(data=FullCAM_output_transformed)+
geom_line(aes(x = Date, y = Total_CO2_sequestered), col = "darkblue",size=1.3,fill="darkblue") +
geom_area(aes(x = Date, y = Total_CO2_sequestered), col = "darkblue",size=1.3,fill="#1589FF") +
geom_area(aes(x = Date, y = CO2_stored_in_forest_debris), col = "darkblue", fill="darkblue",size=1.3) +
labs(title="Cumulative CO2 sequestration using the environmental plantings method",
subtitle = "Analysis from the Hunter Valley, New South Wales",
caption = "Data: FullCAM Model 2021",
x="",
y="Total CO2 sequestration (tCO2/ha)")+
scale_x_date(date_breaks = "5 years",date_labels="%Y")+
scale_y_continuous(expand = c(0,0)) +
theme_minimal()+
theme(axis.line.x = element_line(color="black", size = 0.3),axis.line.y = element_line(color="black", size = 0.3))+
theme(axis.text=element_text(size=12))+
theme(panel.grid.minor = element_blank())+
theme(plot.title=element_text(face="bold",size=12))+
theme(plot.subtitle=element_text(size=11,margin=margin(0,0,15,0)))+
theme(plot.caption=element_text(size=8))
plot_FullCAM_cumulative
ggsave("C:/06 R_code/FullCAM_cumulative_plot.png",units="cm",dpi=200, device="png")6 Developing a yearly estimate
These charts are great… but how many tonnes of CO2 could we expect these planting to sequester per hectare per year?
Let’s make a simple table.
# Calculate the input table annually over the crediting period
FullCAM_output_annual <- FullCAM_output_transformed %>%
group_by(Year) %>%
summarise(yearly_total=sum(CO2_sequestered_month))
knitr::kable(FullCAM_output_annual,digits=2,caption = "Annual sequestration per hectare", "simple")| Year | yearly_total |
|---|---|
| 2021 | 0.04 |
| 2022 | 1.61 |
| 2023 | 7.78 |
| 2024 | 14.10 |
| 2025 | 17.75 |
| 2026 | 19.29 |
| 2027 | 19.53 |
| 2028 | 19.11 |
| 2029 | 18.23 |
| 2030 | 17.24 |
| 2031 | 16.19 |
| 2032 | 15.19 |
| 2033 | 14.14 |
| 2034 | 13.19 |
| 2035 | 12.29 |
| 2036 | 11.49 |
| 2037 | 10.69 |
| 2038 | 9.98 |
| 2039 | 9.33 |
| 2040 | 8.75 |
| 2041 | 8.17 |
| 2042 | 7.66 |
| 2043 | 7.19 |
| 2044 | 6.77 |
| 2045 | 6.36 |
| 2046 | 2.83 |
Better. And useful as an input for ACCU models.
Boiling this all down to a simple number - how many tonnes of CO2 can we expect to be sequestered per hectare per year?
#Calculate the mean sequestration per hectare per year
average_per_ha_per_year <- FullCAM_output_annual %>% summarise(mean_ha_yr=mean(yearly_total))
knitr::kable(average_per_ha_per_year,digits=2, caption ="Average sequestration accross all hectares","simple")| mean_ha_yr |
|---|
| 11.34 |
Ace. 11 tCO2/ha is slightly better than the ~10 tCO2/ha that CSIRO forecasts for the average across Australia.
7 Estimating a model
Let’s use this data to develop some more general equations for carbon sequestration. How many x terms do we need to fit a good polynomial? Unsure. Let’s find out by running a bunch of estimates.
#Create a polynomial to estimate the growth rates
plot_FullCAM_annual_model <- ggplot(FullCAM_output_transformed, aes(Date, CO2_sequestered_year))+
geom_smooth(method = "lm", formula = y ~ poly(x, 1), aes(color = "1"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 2), aes(color = "2"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 3), aes(color = "3"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), aes(color = "4"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 5), aes(color = "5"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 6), aes(color = "6"), se = F) +
geom_smooth(method = "lm", formula = y ~ poly(x, 7), aes(color = "7"), se = F) +
labs(title="Annual sequestration rate using the environmental plantings method",
subtitle = "Fitting adequate polynomials to the series",
caption = "Data: FullCAM Model 2021",
x="",
y="Annual CO2 sequestration (tCO2/ha)",
color= "Polynomial")+
scale_x_date(date_breaks = "5 years",date_labels="%Y")+
scale_y_continuous(expand = c(0,0)) +
theme(legend.position="none")+
theme_minimal()+
theme(axis.line.x = element_line(color="black", size = 0.3),axis.line.y = element_line(color="black", size = 0.3))+
theme(axis.text=element_text(size=12))+
theme(panel.grid.minor = element_blank())+
theme(plot.title=element_text(face="bold",size=12))+
theme(plot.subtitle=element_text(size=11,margin=margin(0,0,15,0)))+
theme(plot.caption=element_text(size=8))
plot_FullCAM_annual_model
The polynomial of order 4 seems to be there best here in estimating the function of CO2 sequestration. It’s about the right shape, and has less x terms than the higher order models. Let’s see what it looks like compared to the original series.
plot_FullCAM_annual_model_2 <- ggplot(FullCAM_output_transformed, aes(Date, CO2_sequestered_year))+
geom_line(aes(x = Date, y = CO2_sequestered_year), col = "#3C33FE",size=1) +
geom_smooth(method = "lm", formula = y ~ poly(x, 4), col = "darkgreen", se = F, linetype="dashed") +
labs(title="Annual sequestration rate using the environmental plantings method",
subtitle = "Polynomial fit compared to actual data set",
caption = "Data: FullCAM Model 2021",
x="",
y="Annual CO2 sequestration (tCO2/ha)",
color= "Polynomial")+
scale_x_date(date_breaks = "5 years",date_labels="%Y")+
scale_y_continuous(expand = c(0,0)) +
theme(legend.position="none")+
theme_minimal()+
theme(axis.line.x = element_line(color="black", size = 0.3),axis.line.y = element_line(color="black", size = 0.3))+
theme(axis.text=element_text(size=12))+
theme(panel.grid.minor = element_blank())+
theme(plot.title=element_text(face="bold",size=12))+
theme(plot.subtitle=element_text(size=11,margin=margin(0,0,15,0)))+
theme(plot.caption=element_text(size=8))
plot_FullCAM_annual_model_2
8 Suitability assessment
There’s a bunch of requirements that must be satisfied to be eligible under the ERF method:
The plantings must include tree or mallee species native to the local area.
Mallee species must only be planted where annual rainfall is under 600 mm.
The project area must have been clear of forest cover for at least five years prior to project commencement.
The plantings must have the potential to reach forest cover (20 per cent crown cover consisting of trees that are at least two metres tall).
So how do we get this information? Let’s start with rainfall. We can download climate data from the BOM here: http://www.bom.gov.au/jsp/ncc/climate_averages/rainfall/index.jsp
#It's possible to read this in as a text file
#BOM_rainfall <- read.delim("C:/06 R_code/rainan.txt")
#However, a raster file is far more useful
#BOM_rainfall <- as.data.frame(raster("C:/06 R_code/rainan.txt"))
#head(BOM_rainfall)
#plot(BOM_rainfall)
#sf_join <- sf::st_join(BOM_rainfall, absmapsdata::sa12016, join = st_intersects)
#Let's cut this raster down to the Australia shapefile using the absmapsdata package
#BOM_rainfall_sa12016 <- crop(BOM_rainfall,extent(absmapsdata::sa12016))
#plot(BOM_rainfall_sa12016)
#mask(BOM_rainfall_sa12016,absmapsdata::sa12016)