Skip to content

Commit

Permalink
changes to plotting functions to work with spadesCBMpython
Browse files Browse the repository at this point in the history
  • Loading branch information
@camillegiuliano
camillegiuliano committed Nov 14, 2024
1 parent bee1273 commit 3b756f2
Showing 1 changed file with 67 additions and 140 deletions.
207 changes: 67 additions & 140 deletions R/plotting.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,58 +44,26 @@ m3ToBiomIncOnlyPlots <- function(inc) {
#' @importFrom data.table as.data.table merge.data.table
#' @importFrom quickPlot Plot
#' @importFrom terra rast values
spatialPlot <- function(pixelkeep, cbmPools, poolsToPlot, years, masterRaster) {
cbmPools[is.na(cbmPools)] <- 0
colnames(cbmPools)[c(1,3,4)] <- c("simYear", "pixelGroup", "age")
if ("totalCarbon" %in% poolsToPlot) {
totalCarbon <- apply(cbmPools[, SoftwoodMerch:HardwoodBranchSnag], 1, "sum")
cbmPools <- cbind(cbmPools, totalCarbon)
}

if (any(!poolsToPlot %in% colnames(cbmPools))) {
stop("The carbon pool you specified for plotting is not contained in the pool definitions")
}

cbmPools <- as.data.table( cbmPools)
#Build rasters for every year and pool
carbonStacks <- vector(mode = "list", length = length(poolsToPlot))
names(carbonStacks) <- poolsToPlot

for (pool in poolsToPlot) {
carbonStacks[[pool]] <- lapply(years, FUN = function(x, poolsDT = cbmPools,
var = pool,
pixelKeep = pixelkeep) {

poolsDT <- poolsDT[order(pixelGroup)] %>% #order by stand index
.[simYear == x, .(pixelGroup, "var" = get(pool))] #subset by year
#subset pixelKeep
colsToKeep <- c("pixelIndex", paste0("pixelGroup", x))

pixelKeep <- pixelKeep[, colsToKeep, with = FALSE] |>
setnames(c("pixelIndex", "pixelGroup"))
# with=FALSE tells data.table colsToKeep isn't a column name until here it
# is ok...then in 1993 - an extra line gets added from the merge below
# Keep <- poolsDT[pixelKeep, on = c("pixelGroup")]
pixelKeep <- merge(pixelKeep, poolsDT, by = "pixelGroup", all.x = TRUE)
pixelKeep <- pixelKeep[order(pixelIndex)] ## join with pixelKeep; order by rowOrder for raster prep
# pixelKeep <- pixelKeep[poolsDT, on = c("pixelGroup")] %>% #join with pixelKeep
pixels <- values(masterRaster)

plotMaster <- rast(masterRaster)
plotMaster[] <- 0
plotMaster[pixelKeep$pixelIndex] <- pixelKeep$var
# masterRaster[masterRaster == 0] <- NA #Species has zeroes instead of NA. Revisit if masterRaster changes
# masterRaster[!is.na(masterRaster)] <- pixelKeep$var

## name will begin with x if no character assigned
return(plotMaster)
})
}

names(carbonStacks) <- paste0(poolsToPlot)

temp <- unlist(carbonStacks)
quickPlot::Plot(temp, title = paste0(poolsToPlot, " in ", years, " MgC/ha")) ## addTo = "temp",
spatialPlot <- function(cbmPools, years, masterRaster, spatialDT) {
masterRaster <- terra::unwrap(masterRaster)
cbmPools <- as.data.table(cbmPools)
totalCarbon <- apply(cbmPools[, Merch:BranchSnag],
1, "sum")
totalCarbon <- cbind(cbmPools, totalCarbon)
totalCarbon <- totalCarbon[simYear == years,]
t <- spatialDT[, .(pixelIndex, pixelGroup)]
setkey(t, pixelGroup)
setkey(totalCarbon, pixelGroup)
temp <- merge(t, totalCarbon, allow.cartesian=TRUE)
setkey(temp, pixelIndex)
plotM <- terra::rast(masterRaster)
plotM[] <- 0
plotM[temp$pixelIndex] <- temp$totalCarbon
pixSize <- prod(res(masterRaster))/10000
temp[, `:=`(pixTC, totalCarbon * pixSize)]
overallTC <- sum(temp$pixTC)/(nrow(temp) * pixSize)
quickPlot::Plot(plotM, new = TRUE,
title = paste0("Total Carbon in ", years, " in MgC/ha"))
}

#' `carbonOutPlot`
Expand All @@ -110,59 +78,33 @@ spatialPlot <- function(pixelkeep, cbmPools, poolsToPlot, years, masterRaster) {
#' @importFrom ggplot2 aes element_text geom_col geom_line ggplot labs
#' @importFrom ggplot2 scale_fill_discrete scale_y_continuous sec_axis theme xlab
#' @importFrom quickPlot Plot
carbonOutPlot <- function(emissionsProducts, masterRaster) {
# pixelCount <- cbmPools[, .(simYear, pixelGroup, pixelCount)]
# cols <- c("simYear", "pixelGroup")
# productEmissions <- merge(pixelCount, emissionsProducts,
# by = cols)
# totalOutByYr <- productEmissions[, .(
# CO2 = sum(CO2 * (prod(res(masterRaster)) / 10000) * pixelCount),
# CH4 = sum(CH4 * (prod(res(masterRaster)) / 10000) * pixelCount),
# CO = sum(CO * (prod(res(masterRaster)) / 10000) * pixelCount),
# Products = sum(Products * (prod(res(masterRaster)) / 10000) * pixelCount)), by = .(simYear)]
# Units: these are absolute values the total products and emissions, not per ha.

# the CH4 and CO are so small compared to the CO2 that we will add all gases
# and call it emissions
carbonOutPlot <- function(emissionsProducts) {
totalOutByYr <- as.data.table(emissionsProducts)
totalOutByYr[ , Emissions := (CO2 + CH4 + CO)]
cols <- c("CO2", "CH4", "CO")
totalOutByYr[ , (cols) := NULL]
totalOutByYr[, `:=`((cols), NULL)]

## Emissions only
# a <- ggplot(data = totalOutByYr, aes(x = simYear)) +
# + geom_line(aes(y = Emissions), size = 1.5, colour = "#69b3a2")
## Emissions/4 + products
# a + geom_line(aes(y = Products), size = 1.5, colour = "darkred")
absCbyYrProducts <- ggplot(totalOutByYr, aes(x = simYear, y = Products)) +
geom_line(linewidth = 1.5) +
scale_y_continuous(name = "Products in MgC") +
scale_x_continuous(breaks = scales::pretty_breaks()) +
xlab("Simulation Years") + theme_classic() +
theme(axis.title.y = element_text(size = 10),
axis.title.x = element_text(size = 10))

absCbyYrProducts <- ggplot(data = totalOutByYr, aes(x = simYear, y = Products)) +
geom_line(colour = "darkred", size = 1.5) +
geom_line(aes(y = Products), size = 1.5, colour = "#69b3a2") +
scale_y_continuous(
# Features of the first axis
name = "Products in MgC",
# Add a second axis and specify its features
# sec.axis = sec_axis(trans = ~.*coeff, name = "Emissions (CO2+CH4+CO) in MgC")
) +
xlab("Simulation Years") +
theme(axis.title.y = element_text(color = "darkred", size = 13),
axis.title.y.right = element_text(color = "#69b3a2", size = 13, angle = 270))
#ggtitle("Yearly Emissions and Forest Products")
absCbyYrEmissions <- ggplot(data = totalOutByYr, aes(x = simYear, y = Emissions)) +
geom_line(colour = "darkred", size = 1.5) +
geom_line(aes(y = Emissions), size = 1.5, colour = "#69b3a2") +
scale_y_continuous(
# Features of the first axis
name = "Emissions (CO2+CH4+CO) in MgC",
# Add a second axis and specify its features
#sec.axis = sec_axis(trans = ~.*coeff, name = "Emissions (CO2+CH4+CO) in MgC")
) +
xlab("Simulation Years") +
theme(axis.title.y = element_text(color = "darkred", size = 13),
axis.title.y.right = element_text(color = "#69b3a2", size = 13, angle = 270))

quickPlot::Plot(absCbyYrProducts, addTo = "absCbyYrProducts", title = "Yearly Forest Products")
quickPlot::Plot(absCbyYrEmissions, addTo = "absCbyYrEmissions", title = "Yearly Emissions")
geom_line(linewidth = 1.5) +
scale_y_continuous(limits = c(0, NA)) +
scale_x_continuous(breaks = scales::pretty_breaks()) +
labs(x = "Simulation Years", y = expression(paste('Emissions (CO'[2]*'+CH'[4]*'+CO) in MgC'))) +
theme_classic() +
theme(axis.title.y = element_text(size = 10),
axis.title.x = element_text(size = 10))


quickPlot::Plot(absCbyYrProducts, addTo = "absCbyYrProducts",
title = "Yearly Forest Products")
quickPlot::Plot(absCbyYrEmissions, addTo = "absCbyYrEmissions",
title = "Yearly Emissions")
}

#' `NPPplot`
Expand All @@ -178,34 +120,28 @@ carbonOutPlot <- function(emissionsProducts, masterRaster) {
#' @importFrom quickPlot Plot
#' @importFrom raster raster
NPPplot <- function(spatialDT, NPP, masterRaster) {
# Calculate the avgNPP (MgC/ha) by pixel group.
masterRaster <- terra::unwrap(masterRaster)
npp <- as.data.table(copy(NPP))
npp[,avgNPP := mean(NPP), by = c("pixelGroup")]
npp[, `:=`(avgNPP, mean(NPP)), by = c("pixelGroup")]
cols <- c("simYear", "NPP")
avgNPP <- unique(npp[, (cols) := NULL])
# link that to the pixels
avgNPP <- unique(npp[, `:=`((cols), NULL)])
t <- spatialDT[, .(pixelIndex, pixelGroup)]
setkey(t,pixelGroup)
setkey(avgNPP,pixelGroup)
temp <- merge(t, avgNPP, on = "pixelGroup")
setkey(t, pixelGroup)
setkey(avgNPP, pixelGroup)
temp <- merge(t, avgNPP, allow.cartesian=TRUE)
setkey(temp, pixelIndex)
#pixelCount[which(is.na(pixelCount$N)),"N"] <- 0
# temp1 <- temp[which(!is.na(temp$simYear)),.(pixelIndex,NPP)]
# temp1[order(pixelIndex)]
#masterRaster[!masterRaster == 0] <- temp$NPP
plotMaster <- raster(masterRaster)
plotMaster <- terra::rast(masterRaster)
plotMaster[] <- 0
# instead of tC/ha for each pixel,
plotMaster[temp$pixelIndex] <- temp$avgNPP
#pixel size in ha
pixSize <- prod(res(masterRaster))/10000
temp[,pixNPP := avgNPP*pixSize]
overallAvgNpp <- sum(temp$pixNPP)/(nrow(temp)*pixSize)
temp[, `:=`(pixNPP, avgNPP * pixSize)]
overallAvgNpp <- sum(temp$pixNPP)/(nrow(temp) * pixSize)
quickPlot::Plot(plotMaster, new = TRUE,
title = paste0("Pixel-level average NPP MgC/ha/yr.",
"\n Landscape average: ", round(overallAvgNpp,3), " MgC/ha/yr."))
title = paste0("Pixel-level average NPP",
"\n Landscape average: ", round(overallAvgNpp, 3), " MgC/ha/yr."))
}


#' `barPlot`
#'
#' @param cbmPools TODO
Expand All @@ -217,44 +153,35 @@ NPPplot <- function(spatialDT, NPP, masterRaster) {
#' @importFrom data.table as.data.table melt.data.table
#' @importFrom ggplot2 aes geom_col ggplot labs scale_fill_discrete theme_bw
#' @importFrom quickPlot Plot
barPlot <- function(cbmPools, masterRaster) {
#This needs to change to belowground living, aboveground living, soil, snag
#Need to first average values per ha
barPlot <- function(cbmPools) {
cbmPools <- as.data.table(cbmPools)
cbmPools$pixelGroup <- as.character(cbmPools$pixelGroup)
# this ONLY works is the number of modelled pixels does not change per simulation year
pixelNo <- sum(cbmPools$pixelCount / length(unique(cbmPools$simYear))) #Get pixel Sum
pixelNo <- sum(cbmPools$pixelCount/length(unique(cbmPools$simYear)))
cbmPools$simYear <- as.character(cbmPools$simYear)

carbonCompartments <- cbmPools[, .(soil = sum(AboveGroundVeryFastSoil, BelowGroundVeryFastSoil,
AboveGroundFastSoil, BelowGroundFastSoil,
AboveGroundSlowSoil, BelowGroundSlowSoil, MediumSoil),
AGlive = sum(SoftwoodMerch, SoftwoodFoliage, SoftwoodOther,
HardwoodMerch, HardwoodFoliage, HardwoodOther),
BGlive = sum(SoftwoodCoarseRoots, SoftwoodFineRoots,
HardwoodCoarseRoots, HardwoodFineRoots),
snags = sum(SoftwoodStemSnag, SoftwoodBranchSnag,
HardwoodStemSnag, HardwoodBranchSnag),
weight = pixelCount/pixelNo),
AGlive = sum(Merch, Foliage, Other),
BGlive = sum(CoarseRoots,FineRoots),
snags = sum(StemSnag, BranchSnag), weight = pixelCount/pixelNo),
by = .(pixelGroup, simYear)]

outTable <- carbonCompartments[, .(soil = sum(soil * weight),
AGlive = sum(AGlive * weight),
BGlive = sum(BGlive * weight),
snags = sum(snags * weight)),
by = simYear]

outTable <- data.table::melt.data.table(outTable, id.vars = 'simYear',
outTable <- data.table::melt.data.table(outTable, id.vars = "simYear",
measure.vars = c("soil", "AGlive", "BGlive", "snags"),
variable.name = 'pool',
value.name = "carbon")
variable.name = "pool", value.name = "carbon")
outTable$simYear <- as.numeric(outTable$simYear)
outTable$carbon <- as.numeric(outTable$carbon)
barPlots <- ggplot(data = outTable, aes(x = simYear, y = carbon, fill = pool)) +
geom_col(position = "fill") +
scale_fill_discrete(name = "carbon compartment") +
labs(x = "Year", y = "proportion") +
theme_bw()
scale_y_continuous(expand = expansion(mult = c(0, .1))) +
scale_fill_discrete(name = "Carbon Compartment") +
labs(x = "Year", y = "Proportion") + theme_classic() +
guides(fill = guide_legend(title.position= "top", title ="Carbon compartment") ) +
scale_fill_brewer(palette = "Set1", labels = c("Soil", "AGlive", "BGlive", 'snags'))

quickPlot::Plot(barPlots, addTo = 'barPlots', title = "Proportion of C above and below ground compartments.")
quickPlot::Plot(barPlots, addTo = "barPlots", title = "Proportion of C above and below ground compartments.")
}

0 comments on commit 3b756f2

Please sign in to comment.