About this notebook

This notebook includes R code to reproduce the figures from the publication How do forest fires affect soil greenhouse gas emissions in upland boreal forests? A review accessible through here.

Packages used ggplot2, ggformula, and extrafont. See citation at the end of script.

GitHub repository for the publication using JMP scripts: datainbrief-repo

GitHub repository for this notebook: datainbrief-binder

Figure 2

#Load libraries and pick fonts
library(ggplot2)
library(ggformula)
library("extrafont")

myFont <- "Segoe UI"
windowsFonts(mono  = myFont)

#Load data
data <- read.table(file = "ghg-fluxes.csv",
                   header = TRUE,
                   dec = ".",
                   sep = ",")

#Select relevant variables for fig.1 
myvars <- c("Rs", "Ageclass","Ageclass100", "Severity1",
            "Severity2", "Source", "Region")
newdata <- data[myvars]

#Make sure factors are factors
newdata$Severity<-factor(newdata$Severity1)
newdata$Severity2<-factor(newdata$Severity2)
newdata$Source<-factor(newdata$Source)
newdata$Region<-factor(newdata$Region)

#Check and exclude missing data
colSums(is.na(newdata))
##          Rs    Ageclass Ageclass100   Severity1   Severity2      Source 
##           0           0           0          11           0           0 
##      Region    Severity 
##           0          11
data1 <- na.exclude(newdata)

#Check means that extrapolate 
##We can see that most mean CO2 fluxes after high (H) and low (L) severity
##fires are within a similar range, and H includes a few means that extrapolate that.
dotchart(data1$Rs, group = data1$Severity)

#Omit means that extrapolate (like done the in the paper)
data2 <- subset(data1, Rs < 30)
data2$ID <- seq.int(nrow(data2))
data3 <- data2[ !(data2$ID %in% c(55, 105, 106, 107)), ]
#Create labels for graph
data3$labelsev2 <- factor(data3$Severity2,labels =
                            c("High severity", "Low severity", "Unburned"))
data3$labelsev <- factor(data3$Severity,labels =
                           c("High severity", "Low severity"))
dataH <- subset(data3, labelsev == "High severity")
dataL <- subset(data3, labelsev == "Low severity")

#Plot the graph as in fig. 1
p1 <- ggplot() + 
  geom_jitter(aes(y = Rs,
                  x = Ageclass100,
                  col = labelsev2,
                  shape = labelsev2),
              data = data3,
              size = 1) +
  geom_spline(aes(y = Rs,
                  x = Ageclass100,
                  col = labelsev),
              data = dataH,
              df = 5.5,
              alpha = 1,
              size = 0.5) +
  geom_spline(aes(y = Rs,
                  x = Ageclass100,
                  col = labelsev),
              data = dataL,
              df = 4,
              alpha = 1,
              size = 0.5) +
  labs(x = "Time after fire (years)",
       y = (bquote(''*g ~CO[2]~ m^-2~d^-1*''))) + 
  theme_bw() + 
  theme(legend.position = c(0.155, 0.91),
        legend.text = element_text(size = 10),
        legend.title = element_blank(),
        legend.key.size = unit(0.5,"line"),
        axis.title = element_text(size=12),
        axis.text = element_text(size=10)) +
  theme(text = element_text(size = 12, family="mono")) +
  scale_x_continuous(breaks = seq(0, 100, by = 20)) +
  scale_colour_manual(name = "labelsev2",
                               breaks = c("High severity", "Low severity",
                                          "Unburned"),
                               values = c("black", "gray",
                                          "#085844")) +
  scale_shape_manual(name = "labelsev2",
                               breaks = c("High severity", "Low severity",
                                          "Unburned"),
                               values = c(1, 1, 16)) 
p1

Fit with smoothing spline: mean soil CO2 efflux over time after fire.