MECO and MCOL re-updated from earlier saved state

R/MCOL.R

 # written by Fabian Stenzel
 # 2022 - stenzel@pik-potsdam.de
 
-# todo: optional reading of monthly netcdf CFT outputs
+# todo: add optional reading of monthly netcdf CFT outputs
 
-#' Plot global map of MCOL to file
-#'
-#' Plot global map of MCOL to file with legend colors similar to Haberl et al. 2007
-#'
-#' @param data array containing MCOL percentage value for each gridcell
-#' @param file character string for location/file to save plot to
-#' @param plotyears range of years to plot over time
-#' @param title character string title for plot
-#' @param legendtitle character string legend title
-#' @param zeroThreshold smallest value to be distinguished from 0 in legend,
-#'        both for negative and positive values (default: 0.1)
-#' @param eps write eps file instead of PNG (boolean) - (default: FALSE)
+#' Calculate MCOL based on a PNV run and LU run of LPJmL
 #'
-#' @return none
-#'
-#' @examples
-#' \dontrun{
-#' }
-#' @export
-plotMCOLmap <- function(data, file, title, legendtitle, zeroThreshold = 0.1, eps = FALSE){
-  brks <- c(-400,-200,-100,-50,-zeroThreshold,zeroThreshold,10,20,30,40,50,60,70,80,100)
-  classes <- c("<-200","-200 - -100","-100 - -50",paste0("-50 - -",zeroThreshold),paste0("-",zeroThreshold," - ",zeroThreshold),paste0(zeroThreshold," - 10"),"10 - 20","20 - 30","30 - 40","40 - 50","50 - 60","60 - 70","70 - 80","80 - 100")
-  palette <- c("navy","royalblue3","royalblue1","skyblue1","grey80","yellowgreen","greenyellow","yellow","gold","orange","orangered","orangered4","brown4","black")
-  data[data < brks[1]] <- brks[1]
-  data[data > brks[length(brks)]] <- brks[length(brks)]
-  if (eps) {
-    file = strsplit(file,".",fixed = TRUE)[[1]]
-    file = paste(c(file[1:(length(file) - 1)],"eps"), collapse = ".")
-    ps.options(family = c("Helvetica"), pointsize = 18)
-    postscript(file,horizontal = FALSE, onefile = FALSE, width = 22,  height = 8.5, paper = "special")
-  }else{
-    png(file, width = 7.25, height = 3.5, units = "in", res = 300, pointsize = 6, type = "cairo")
-  }
-  ra <- raster::raster(ncols = 720, nrows = 360)
-  range <- range(data)
-  ra[raster::cellFromXY(ra,cbind(lon,lat))] <-  data
-  extent <- raster::extent(c(-180, 180, -60, 90))
-  par(bty = "n", oma = c(0,0,0,0), mar = c(0,0,0,0), xpd = T)
-  raster::plot(ra, ext = extent, breaks = brks, col = palette, main = "", legend = FALSE, axes = FALSE)
-  title(title, line = -2)
-  maps::map('world', add = TRUE, res = 0.4, lwd = 0.25,ylim = c(-60,90))
-  legend(x = -180, y = 50, fill = palette, border = palette, legend = classes, title = legendtitle)
-  dev.off()
-}
-
-#' Plot absolute MCOL, overtime, maps, and npp into given folder
-#'
-#' Plot to file a comparison over time of global sums of MCOL, NPPpot, NPPeco,
-#' and NPPact, with legend similar to Krausmann et al. 2013
-#'
-#' @param mcolData mcol data list object (returned from calcMCOL) containing
-#'                  mcol, npp_eco_overtime, npp_act_overtime, npp_pot_overtime,
-#'                  npp_bioenergy_overtime, mcol_overtime, npp_harv_overtime,
-#'                  mcol_overtime_perc_piref, mcol_perc, mcol_perc_piref, ynpp
-#' @param file character string for location/file to save plot to
-#' @param firstyr first year of mcol object
-#' @param plotyears range of years to plot over time
-#' @param minVal y-axis minimum value for plot over time (default: 0)
-#' @param maxVal y-axis maximum value for plot over time (default: 100)
-#' @param legendpos position of legend (default: "topleft")
-#' @param mapyear year to plot mcol map for
-#' @param mapyear_buffer +- years around mapyear to average mcol
-#'                       (make sure these years exist in mcol data)
-#' @param highlightyrs year(s) that should be highlighted in overtime plot (default: 2000)
-#' @param ref reference period for mcol ("pi" or "act"), to either use
-#'        mcolData$mcol_overtime_perc_piref or mcolData$mcol_overtime
-#' @param eps write plots as eps, instead of png (default = FALSE)
-#'
-#' @return none
-#'
-#' @examples
-#' \dontrun{
-#' }
-#' @export
-plotMCOLovertime <- function(mcolData, file, firstyr, plotyrs, highlightyrs = 2000, minVal = 0,
-                            maxVal = 100, legendpos = "topleft", ext = FALSE, eps = FALSE, ref = "pi"){
-  lastyr = firstyr + length(mcolData$npp_act_overtime) - 1
-  colz = c("slateblue","gold","green3","darkorange","black","red3","green","brown","yellow","turquoise")
-  if (eps) {
-    file = strsplit(file,".",fixed=TRUE)[[1]]
-    file = paste(c(file[1:(length(file) - 1)],"eps"),collapse=".")
-    ps.options(family = c("Helvetica"), pointsize = 18)
-    postscript(file,horizontal = FALSE, onefile = FALSE, width = 22,  height = 8.5,paper="special")
-  }else{
-    png(file, width=3.5,  height = 3, units = "in", res = 300, pointsize = 6,type="cairo")
-  }
-  par(bty="o",oma=c(0,0,0,0),mar=c(4,5,1,3))
-  plot(NA,ylab="GtC/yr",xlab="Year",xlim=plotyrs,
-       ylim=c(minVal, maxVal),xaxs="i",yaxs="i")
-  grid()
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_pot_overtime,type = "l",col=colz[1])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_act_overtime,type = "l",col=colz[2])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_eco_overtime,type = "l",col=colz[3])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_luc_overtime,type = "l",col=colz[4])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,type = "l",col=colz[5])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$harvest_cft_overtime,type = "l",col=colz[7])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$rharvest_cft_overtime,type = "l",col=colz[8])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$fire_overtime,type = "l",col=colz[9])
-  lines(x=seq(firstyr,lastyr,1),y=mcolData$timber_harvest_overtime,type = "l",col=colz[10])
-  
-  
-  par(bty="n",oma=c(0,0,0,0),mar=c(4,5,1,3), new = T)
-  if (ref == "pi") {
-    plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime_perc_piref,ylab="",xlab="",xlim=plotyrs,
-         ylim=c(0, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
-  } else if (ref == "act") {
-    plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,ylab="",xlab="",xlim=plotyrs,
-         ylim=c(0, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
-  }else stop(paste0("Unknown value for parameter ref: ",ref," - Aborting."))
-
-  axis(side = 4, col = colz[6],col.axis = colz[6])
-  mtext(text = "%", col=colz[6], side = 4,line = 2)
-
-  if (!is.null(highlightyrs)){
-    for (y in highlightyrs){
-      lines(x=c(y,y),y=c(minVal,maxVal),col="grey40")
-    }
-  }
-  legend(legendpos,legend = c("NPPpot (PNV)","NPPact (landuse)","NPPeco","NPPluc","M-COLabs","M-COL [% NPPpi]","harvestc","rharvest","firec","timber_harvest"),col=colz ,lty=1,cex=1)
-  dev.off()
-}
-#' Calculate the ecosystem change metric gamma between 2 simulations/timesteps
-#'
-#' Function to calculate the ecosystem change metric gamma according
-#' to Sykes (1999), Heyder (2011), and Ostberg (2015,2018).
-#' This is a reformulated version in R, not producing 100% similar values
-#' than the C/bash version from Ostberg 2018, but following their methodology.
+#' Function to calculate MCOL based on a PNV run and LU run of LPJmL
 #'
 #' @param inFol_lu folder of landuse scenario run
 #' @param inFol_pnv folder of pnv reference run
@@ -146,7 +22,10 @@ plotMCOLovertime <- function(mcolData, file, firstyr, plotyrs, highlightyrs = 20
 #' @param varnames data.frame with names (outname) and timestep of output files
 #'                 see code of calcMCOL() for an example with default values
 #'
-#' @return list data object containing arrays of ...
+#' @return list data object containing the data to compute MCOL as arrays of 
+#'         ynpp_potential, ynpp, pftnpp_cft, pftnpp_nat, pftnpp_grasslands, 
+#'         pftnpp_bioenergy, harvest_cft, rharvest_cft, fire, timber, cftfrac, 
+#'         fpc, harvest_grasslands, harvest_bioenergy
 #'
 #' @examples
 #' \dontrun{
@@ -407,12 +286,9 @@ readMCOLdata <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T,
 
 } # end of readMCOLdata
 
-#' Calculate the ecosystem change metric gamma between 2 simulations/timesteps
+#' Calculate MCOL
 #'
-#' Function to calculate the ecosystem change metric gamma according
-#' to Sykes (1999), Heyder (2011), and Ostberg (2015,2018).
-#' This is a reformulated version in R, not producing 100% similar values
-#' than the C/bash version from Ostberg 2018, but following their methodology.
+#' Function to calculate MCOL
 #'
 #' @param inFol_lu folder of landuse scenario run
 #' @param inFol_pnv folder of pnv reference run
@@ -443,7 +319,12 @@ readMCOLdata <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T,
 #'        mapping_lpj67420_to_grazing_regions array with a mapping between 67420
 #'        LPJmL cells and the 29 regions
 #'
-#' @return list data object containing arrays of ...
+#' @return list data object containing MCOL and components as arrays: mcol, 
+#'         mcol_overtime, mcol_overtime_perc_piref, mcol_perc, mcol_perc_piref, 
+#'         ynpp_potential, npp_act_overtime, npp_pot_overtime, npp_eco_overtime, 
+#'         harvest_cft_overtime, npp_luc_overtime, rharvest_cft_overtime, 
+#'         fire_overtime, timber_harvest_overtime, harvest_cft, mcol_harvest,
+#'         grassland_scaling_factor_cellwise,  mcol_luc, mcol_luc_piref
 #'
 #' @examples
 #' \dontrun{
@@ -567,7 +448,8 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T,
                 npp_eco_overtime = npp_eco_overtime, #npp_bioenergy_overtime = npp_bioenergy_overtime,
                 harvest_cft_overtime = harvest_cft_overtime, npp_luc_overtime = npp_luc_overtime,
                 rharvest_cft_overtime = rharvest_cft_overtime, fire_overtime = fire_overtime,
-                timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft,
+                timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft, 
+                rharvest_cft = rharvest_cft,
                 grassland_scaling_factor_cellwise = grassland_scaling_factor_cellwise,
                 mcol_harvest = mcol_harvest, mcol_luc = mcol_luc, mcol_luc_piref = mcol_luc_piref))
   }else{
@@ -579,6 +461,7 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T,
                 harvest_cft_overtime = harvest_cft_overtime, npp_luc_overtime = npp_luc_overtime,
                 rharvest_cft_overtime = rharvest_cft_overtime, fire_overtime = fire_overtime,
                 timber_harvest_overtime = timber_harvest_overtime, harvest_cft = harvest_cft,
+                rharvest_cft = rharvest_cft,
                 mcol_harvest = mcol_harvest, mcol_luc = mcol_luc, mcol_luc_piref = mcol_luc_piref))
   }
 }
@@ -588,9 +471,10 @@ calcMCOL <- function(inFol_lu, inFol_pnv, startyr, stopyr, gridbased = T,
 #' Wrapper function to plot absolute mcol, overtime, maps, and npp into given folder
 #'
 #' @param mcolData mcol data list object (returned from calcMCOL) containing
-#'                  mcol, npp_eco_overtime, npp_act_overtime, npp_pot_overtime,
-#'                  npp_bioenergy_overtime, mcol_overtime, npp_harv_overtime,
-#'                  mcol_overtime_perc_piref, mcol_perc, mcol_perc_piref, ynpp
+#'                 mcol, npp_eco_overtime, npp_act_overtime, npp_pot_overtime,
+#'                 npp_bioenergy_overtime, mcol_overtime, npp_harv_overtime,
+#'                 mcol_overtime_perc_piref, mcol_perc, mcol_perc_piref, ynpp
+#'                 all in GtC
 #' @param outFol folder to write into
 #' @param plotyears range of years to plot over time
 #' @param minVal y-axis minimum value for plot over time
@@ -646,3 +530,124 @@ plotMCOL <- function(mcolData, outFol, plotyears, minVal, maxVal, legendpos,
                 legendtitle = "gC/m2",legYes = T)
   
 } # end of plotMCOL
+
+#' Plot global map of MCOL to file
+#'
+#' Plot global map of MCOL to file with legend colors similar to Haberl et al. 2007
+#'
+#' @param data array containing MCOL percentage value for each gridcell
+#' @param file character string for location/file to save plot to
+#' @param plotyears range of years to plot over time
+#' @param title character string title for plot
+#' @param legendtitle character string legend title
+#' @param zeroThreshold smallest value to be distinguished from 0 in legend,
+#'        both for negative and positive values (default: 0.1)
+#' @param eps write eps file instead of PNG (boolean) - (default: FALSE)
+#'
+#' @return none
+#'
+#' @examples
+#' \dontrun{
+#' }
+#' @export
+plotMCOLmap <- function(data, file, title = "", legendtitle = "", zeroThreshold = 0.1, eps = FALSE){
+  brks <- c(-400,-200,-100,-50,-zeroThreshold,zeroThreshold,10,20,30,40,50,60,70,80,100)
+  classes <- c("<-200","-200 - -100","-100 - -50",paste0("-50 - -",zeroThreshold),paste0("-",zeroThreshold," - ",zeroThreshold),paste0(zeroThreshold," - 10"),"10 - 20","20 - 30","30 - 40","40 - 50","50 - 60","60 - 70","70 - 80","80 - 100")
+  palette <- c("navy","royalblue3","royalblue1","skyblue1","grey80","yellowgreen","greenyellow","yellow","gold","orange","orangered","orangered4","brown4","black")
+  data[data < brks[1]] <- brks[1]
+  data[data > brks[length(brks)]] <- brks[length(brks)]
+  if (eps) {
+    file = strsplit(file,".",fixed = TRUE)[[1]]
+    file = paste(c(file[1:(length(file) - 1)],"eps"), collapse = ".")
+    ps.options(family = c("Helvetica"), pointsize = 18)
+    postscript(file,horizontal = FALSE, onefile = FALSE, width = 22,  height = 8.5, paper = "special")
+  }else{
+    png(file, width = 7.25, height = 3.5, units = "in", res = 300, pointsize = 6, type = "cairo")
+  }
+  ra <- raster::raster(ncols = 720, nrows = 360)
+  range <- range(data)
+  ra[raster::cellFromXY(ra,cbind(lon,lat))] <-  data
+  extent <- raster::extent(c(-180, 180, -60, 90))
+  par(bty = "n", oma = c(0,0,0,0), mar = c(0,0,0,0), xpd = T)
+  raster::plot(ra, ext = extent, breaks = brks, col = palette, main = "", legend = FALSE, axes = FALSE)
+  title(title, line = -2)
+  maps::map('world', add = TRUE, res = 0.4, lwd = 0.25,ylim = c(-60,90))
+  legend(x = -180, y = 50, fill = palette, border = palette, legend = classes, title = legendtitle)
+  dev.off()
+}
+
+#' Plot absolute MCOL, overtime, maps, and npp into given folder
+#'
+#' Plot to file a comparison over time of global sums of MCOL, NPPpot, NPPeco,
+#' and NPPact, with legend similar to Krausmann et al. 2013
+#'
+#' @param mcolData mcol data list object (returned from calcMCOL) containing
+#'                  mcol, npp_eco_overtime, npp_act_overtime, npp_pot_overtime,
+#'                  npp_bioenergy_overtime, mcol_overtime, npp_harv_overtime,
+#'                  mcol_overtime_perc_piref, mcol_perc, mcol_perc_piref, ynpp
+#'                 all in GtC
+#' @param file character string for location/file to save plot to
+#' @param firstyr first year of mcol object
+#' @param plotyrs range of years to plot over time
+#' @param highlightyrs year(s) that should be highlighted in overtime plot (default: 2000)
+#' @param minVal y-axis minimum value for plot over time (default: 0)
+#' @param maxVal y-axis maximum value for plot over time (default: 100)
+#' @param legendpos position of legend (default: "topleft")
+#' @param highlightyrs year(s) that should be highlighted in overtime plot (default: 2000)
+#' @param ref reference period for mcol ("pi" or "act"), to either use
+#'        mcolData$mcol_overtime_perc_piref or mcolData$mcol_overtime
+#' @param eps write plots as eps, instead of png (default = FALSE)
+#'
+#' @return none
+#'
+#' @examples
+#' \dontrun{
+#' }
+#' @export
+plotMCOLovertime <- function(mcolData, file, firstyr, plotyrs, highlightyrs = 2000, minVal = 0,
+                            maxVal = 100, legendpos = "topleft", ext = FALSE, eps = FALSE, ref = "pi"){
+  lastyr = firstyr + length(mcolData$npp_act_overtime) - 1
+  colz = c("slateblue","gold","green3","darkorange","black","red3","green","brown","yellow","turquoise")
+  if (eps) {
+    file = strsplit(file,".",fixed=TRUE)[[1]]
+    file = paste(c(file[1:(length(file) - 1)],"eps"),collapse=".")
+    ps.options(family = c("Helvetica"), pointsize = 18)
+    postscript(file,horizontal = FALSE, onefile = FALSE, width = 22,  height = 8.5,paper="special")
+  }else{
+    png(file, width=3.5,  height = 3, units = "in", res = 300, pointsize = 6,type="cairo")
+  }
+  par(bty="o",oma=c(0,0,0,0),mar=c(4,5,1,3))
+  plot(NA,ylab="GtC/yr",xlab="Year",xlim=plotyrs,
+       ylim=c(minVal, maxVal),xaxs="i",yaxs="i")
+  grid()
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_pot_overtime,type = "l",col=colz[1])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_act_overtime,type = "l",col=colz[2])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_eco_overtime,type = "l",col=colz[3])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$npp_luc_overtime,type = "l",col=colz[4])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,type = "l",col=colz[5])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$harvest_cft_overtime,type = "l",col=colz[7])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$rharvest_cft_overtime,type = "l",col=colz[8])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$fire_overtime,type = "l",col=colz[9])
+  lines(x=seq(firstyr,lastyr,1),y=mcolData$timber_harvest_overtime,type = "l",col=colz[10])
+  
+  
+  par(bty="n",oma=c(0,0,0,0),mar=c(4,5,1,3), new = T)
+  if (ref == "pi") {
+    plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime_perc_piref,ylab="",xlab="",xlim=plotyrs,
+         ylim=c(0, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
+  } else if (ref == "act") {
+    plot(x=seq(firstyr,lastyr,1),y=mcolData$mcol_overtime,ylab="",xlab="",xlim=plotyrs,
+         ylim=c(0, 30),type = "l",col=colz[6],xaxs="i", yaxs="i", axes = F)
+  }else stop(paste0("Unknown value for parameter ref: ",ref," - Aborting."))
+
+  axis(side = 4, col = colz[6],col.axis = colz[6])
+  mtext(text = "%", col=colz[6], side = 4,line = 2)
+
+  if (!is.null(highlightyrs)){
+    for (y in highlightyrs){
+      lines(x=c(y,y),y=c(minVal,maxVal),col="grey40")
+    }
+  }
+  legend(legendpos,legend = c("NPPpot (PNV)","NPPact (landuse)","NPPeco","NPPluc","M-COLabs","M-COL [% NPPpi]","harvestc","rharvest","firec","timber_harvest"),col=colz ,lty=1,cex=1)
+  dev.off()
+}