plot_biocol.R

# written by Fabian Stenzel
# 2022-2023 - stenzel@pik-potsdam.de

################# BioCol plot functions  ###################

#' Plot absolute BioCol, overtime, maps, and npp into given folder
#'
#' Wrapper function to plot absolute biocol, overtime, maps, and npp into given
#' folder
#'
#' @param biocol_data biocol data list object (returned from calc_biocol)
#'        containing biocol_overtime, biocol_overtime_abs, 
#'        biocol_overtime_abs_frac_piref, biocol_overtime_frac_piref, 
#'        biocol_overtime_frac, biocol_overtime_abs_frac, npp_harv_overtime, 
#'        npp_luc_overtime, npp_act_overtime, npp_pot_overtime, npp_eco_overtime, 
#'        harvest_grasslands_overtime, harvest_bioenergy_overtime, 
#'        harvest_cft_overtime, rharvest_cft_overtime, fire_overtime, 
#'        timber_harvest_overtime, wood_harvest_overtime, biocol, biocol_frac, 
#'        npp, biocol_frac_piref, npp_potential, npp_ref, harvest_cft, 
#'        rharvest_cft, biocol_harvest, biocol_luc - all in GtC
#' @param path_write folder to write outputs into
#' @param plotyears range of years to plot over time
#' @param min_val y-axis minimum value for plot over time
#' @param max_val y-axis maximum value for plot over time
#' @param legendpos position of legend
#' @param start_year first year of biocol_data object
#' @param details show all harvest components or not
#' @param mapyear year to plot biocol map for
#' @param mapyear_buffer +- years around mapyear to average biocol
#'        (make sure these years exist in biocol_data) - default: 5
#' @param highlightyear year(s) that should be highlighted in overtime plot
#' @param eps write plots as eps, instead of png (default = FALSE)
#'
#' @examples
#' \dontrun{
#' plot_biocol(
#'    biocol_data = biocol,
#'    path_write = "~/BioCol_plots/",
#'    plotyears = c(1980,2014),
#'    min_val = 0,
#'    max_val = 90,
#'    legendpos = "left",
#'    start_year = 1980,
#'    mapyear = 2000,
#'    highlightyear = 2000,
#'    eps = FALSE
#' )
#' }
#'
#' @md
#' @export
plot_biocol <- function(
    biocol_data,
    path_write,
    plotyears,
    min_val,
    max_val,
    legendpos,
    details = FALSE,
    start_year,
    mapyear,
    mapyear_buffer = 5,
    highlightyear,
    eps = FALSE) {
  mapindex <- mapyear - start_year
  message("Plotting BioCol figures")
  dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)

  plot_global(
    data = rowMeans(
      biocol_data$biocol[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "BioCol_absolute_", mapyear, ".png"),
    type = "exp",
    title = "",
    # paste0("BioCol_abs in ", mapyear),
    pow2min = 0,
    pow2max = 12,
    legendtitle = "GtC",
    leg_yes = TRUE,
    only_pos = FALSE,
    eps = eps
  )

  plot_global(
    data = rowMeans(
      biocol_data$biocol_luc[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "BioCol_luc_", mapyear, ".png"),
    type = "exp",
    title = "",
    # paste0("BioCol_luc in ", mapyear),
    pow2min = 0,
    pow2max = 12,
    legendtitle = "GtC",
    leg_yes = TRUE,
    only_pos = FALSE,
    eps = eps
  )

  plot_global(
    data = rowMeans(
      biocol_data$biocol_harvest[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "BioCol_harv_", mapyear, ".png"),
    type = "exp",
    title = "",
    # paste0("BioCol_harv in ", mapyear),
    pow2min = 0,
    pow2max = 12,
    legendtitle = "GtC",
    leg_yes = TRUE,
    only_pos = FALSE,
    eps = eps
  )

  plot_biocol_ts(
    biocol_data = biocol_data,
    file = paste0(
      path_write, "BioCol_overtime_LPJmL_", plotyears[1], "-", plotyears[2], ".png" # nolint
    ),
    first_year = start_year,
    plot_years = plotyears,
    min_val = min_val,
    ref = "pi",
    legendpos = legendpos,
    details = details,
    max_val = max_val,
    eps = eps,
    highlight_years = highlightyear
  )

  plot_global(
    data = rowMeans(
      biocol_data$biocol_frac[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "BioCol_frac_LPJmL_", mapyear, ".png"),
    legendtitle = "frac of NPPpot",
    type = "lin",
    min = -1,
    max = 1,
    col_pos = "Reds",
    col_neg = "Blues",
    leg_yes = TRUE,
    eps = FALSE,
    n_legend_ticks = 11
  )

  plot_global(
    data = rowMeans(
      biocol_data$biocol_frac_piref[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "BioCol_frac_piref_LPJmL_", mapyear, ".png"),
    title = "",
    legendtitle = "frac of NPPref",
    type = "lin",
    min = -1,
    max = 1,
    col_pos = "Reds",
    col_neg = "Blues",
    leg_yes = TRUE,
    eps = FALSE,
    n_legend_ticks = 11
  )

  plot_global(
    data = rowMeans(
      biocol_data$npp[, (mapindex - mapyear_buffer):(mapindex + mapyear_buffer)] # nolint
    ),
    file = paste0(path_write, "NPP_LPJmL_", mapyear, ".png"),
    type = "lin",
    only_pos = TRUE,
    title = "",
    legendtitle = "gC/m2",
    leg_yes = TRUE,
    min = 0,
    max = 1800
  )
}


#' Plot global map of BioCol to file
#'
#' Plot global map of BioCol to file with legend colors similar to
#' Haberl et al. 2007
#'
#' @param data array containing BioCol 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 zero_threshold smallest value to be distinguished from 0 in legend,
#'        both for negative and positive values (default: 0.1)
#' @param haberl_legend use color palette from Haberl et al.? (default: FALSE)
#' @param eps write eps file instead of PNG (boolean) - (default: FALSE)
#'
#' @examples
#' \dontrun{
#' plot_biocol_map(
#'    data = biocol$biocol_frac[,"2000"]*100,
#'    file = "./BioCol_map_yr2000.png",
#' )
#' }
#'
#' @md
#' @export
plot_biocol_map <- function(
    data,
    file,
    title = "",
    legendtitle = "",
    zero_threshold = 0.001,
    eps = FALSE,
    haberl_legend = FALSE) {
  path_write <- dirname(file)
  dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)

  if (haberl_legend) {
    brks <- c(
      -400, -200, -100, -50, -zero_threshold,
      zero_threshold, 10, 20, 30, 40, 50, 60, 70, 80, 100
    )
    classes <- c(
      "<-200", "-200 - -100", "-100 - -50",
      paste0("-50 - -", zero_threshold),
      paste0("-", zero_threshold, " - ", zero_threshold),
      paste0(zero_threshold, " - 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"
    )
  } else {
    brks <- c(
      -400, seq(-100, -10, 10), -zero_threshold,
      zero_threshold, seq(10, 100, 10), 400
    ) / 100
    classes <- c(
      "<-1", "-1 - -0.9", "-0.9 - -0.8", "-0.8 - -0.7",
      "-0.7 - -0.6", "-0.6 - -0.5", "-0.5 - -0.4", "-0.4 - -0.3",
      "-0.3 - -0.2", "-0.2 - -0.1", paste("-0.1 - -", zero_threshold),
      paste("-", zero_threshold, " - ", zero_threshold),
      paste(zero_threshold, " - 0.1"), "0.1 - 0.2", "0.2 - 0.3",
      "0.3 - 0.4", "0.4 - 0.5", "0.5 - 0.6", "0.6 - 0.7",
      "0.7 - 0.8", "0.8 - 0.9", "0.9 - 1", ">1"
    )
    palette <- grDevices::colorRampPalette(rev(
      RColorBrewer::brewer.pal(11, "RdBu")
    ))(length(brks) - 1)
  }

  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 = ".")
    grDevices::ps.options(family = c("Helvetica"), pointsize = 18)
    grDevices::postscript(file,
      horizontal = FALSE, onefile = FALSE, width = 22,
      height = 8.5, paper = "special"
    )
  } else {
    grDevices::png(file,
      width = 7.25, height = 3.5, units = "in", res = 300,
      pointsize = 6, type = "cairo"
    )
  }
  ra <- terra::rast(ncols = 720, nrows = 360)
  range <- range(data)
  ra[terra::cellFromXY(ra, cbind(lon, lat))] <- data
  extent <- terra::ext(-180, 180, -60, 90)
  graphics::par(bty = "n", oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), xpd = TRUE)
  terra::plot(ra,
    ext = extent, breaks = brks, col = palette, main = "",
    legend = FALSE, axes = FALSE
  )
  graphics::title(title, line = -2)
  maps::map("world", add = TRUE, res = 0.4, lwd = 0.25, ylim = c(-60, 90))
  graphics::legend(
    x = -180, y = 50, fill = palette, border = palette,
    legend = classes, title = legendtitle
  )
  grDevices::dev.off()
}

#' Plot absolute BioCol, overtime, maps, and npp into given folder
#'
#' Plot to file a comparison over time of global sums of BioCol, NPPpot, NPPeco,
#' and NPPact, with legend similar to Krausmann et al. 2013
#'
#' @param biocol_data biocol data list object (returned from calc_biocol)
#' containing biocol, npp_eco_overtime, npp_act_overtime, npp_pot_overtime,
#' npp_bioenergy_overtime, biocol_overtime, npp_harv_overtime,
#' biocol_overtime_perc_piref, biocol_perc, biocol_perc_piref, npp
#' all in GtC
#' @param file character string for location/file to save plot to
#' @param first_year first year of biocol object
#' @param plot_years range of years to plot over time
#' @param highlight_years year(s) that should be highlighted in overtime plot
#' (default: 2000)
#' @param min_val y-axis minimum value for plot over time (default: 0)
#' @param max_val y-axis maximum value for plot over time (default: 100)
#' @param legendpos position of legend (default: "topleft")
#' @param highlight_years year(s) that should be highlighted in overtime plot
#' (default: 2000)
#' @param details show all harvest components or not, (default: FALSE)
#' @param ref reference period for biocol ("pi" or "act"), to either use
#'        biocol_data$biocol_overtime_perc_piref or biocol_data$biocol_overtime
#' @param eps write plots as eps, instead of png (default = FALSE)
#'
#' @examples
#' \dontrun{
#' plot_biocol_ts(
#'   biocol_data = biocol_data,
#'   file = "./BioCol_overtime_LPJmL_1550-2015.png"),
#'   first_year = 1550,
#'   plot_years = c(1550,2015),
#'   min_val = 0,
#'   max_val = 80,
#'   ref = "pi",
#'   legendpos = "topleft",
#'   details = TRUE,
#'   max_val = max_val,
#'   highlight_years = c(1900,2000)
#'   )
#' }
#'
#' @md
#' @export
plot_biocol_ts <- function(
    biocol_data,
    file,
    first_year,
    plot_years,
    highlight_years = 2000,
    details = FALSE,
    min_val = 0,
    max_val = 100,
    max_val_right = 0.45,
    legendpos = "topleft",
    ext = FALSE,
    eps = FALSE,
    ref = "pi") {
  path_write <- dirname(file)
  dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)

  last_year <- first_year + length(biocol_data$npp_act_overtime) - 1
  colz <- c(
    "slateblue", "gold", "green3", "grey60", "red3",
    "black", "grey40", "darkorange", "#ff8c009d", "blue",
    "yellow", "turquoise", "darkgreen"
  )

  if (eps) {
    file <- strsplit(file, ".", fixed = TRUE)[[1]]
    file <- paste(c(file[1:(length(file) - 1)], "eps"), collapse = ".")
    grDevices::ps.options(family = c("Helvetica"), pointsize = 18)
    grDevices::postscript(file,
      horizontal = FALSE, onefile = FALSE, width = 22,
      height = 8.5, paper = "special"
    )
  } else {
    grDevices::png(file,
      width = 3.5, height = 3, units = "in", res = 300,
      pointsize = 6, type = "cairo"
    )
  }

  graphics::par(bty = "o", oma = c(0, 0, 0, 0), mar = c(4, 5, 1, 3))
  graphics::plot(NA,
    ylab = "GtC/yr", xlab = "Year", xlim = plot_years,
    ylim = c(min_val, max_val), xaxs = "i", yaxs = "i"
  )
  graphics::grid()
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$npp_pot_overtime,
    type = "l",
    col = colz[1]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$npp_act_overtime,
    type = "l",
    col = colz[2]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$npp_eco_overtime,
    type = "l",
    col = colz[3]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$npp_luc_overtime,
    type = "l",
    col = colz[4]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$biocol_overtime_abs,
    type = "l",
    col = colz[6]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$biocol_overtime,
    type = "l",
    col = colz[7]
  )
  graphics::lines(
    x = seq(first_year, last_year, 1),
    y = biocol_data$npp_harv_overtime,
    type = "l",
    col = colz[5]
  )
  if (details) {
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$rharvest_cft_overtime,
      type = "l",
      col = colz[10]
    )
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$fire_overtime,
      type = "l", col = colz[11]
    )
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$timber_harvest_overtime,
      type = "l",
      col = colz[12]
    )
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$wood_harvest_overtime,
      type = "l",
      col = colz[13]
    )
  }
  graphics::par(bty = "n", oma = c(0, 0, 0, 0), mar = c(4, 5, 1, 3), new = TRUE)
  if (ref == "pi") {
    graphics::plot(
      x = seq(first_year, last_year, 1),
      y = biocol_data$biocol_overtime_abs_frac_piref,
      ylab = "",
      xlab = "",
      xlim = plot_years,
      ylim = c(0, max_val_right),
      type = "l",
      col = colz[8],
      xaxs = "i",
      yaxs = "i",
      axes = FALSE
    )
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$biocol_overtime_frac_piref,
      ylab = "",
      xlab = "",
      xlim = plot_years,
      ylim = c(0, 0.4),
      col = colz[9],
      xaxs = "i",
      yaxs = "i",
      axes = FALSE
    )
  } else if (ref == "act") {
    graphics::plot(
      x = seq(first_year, last_year, 1),
      y = biocol_data$biocol_overtime_abs_frac,
      ylab = "",
      xlab = "",
      xlim = plot_years,
      ylim = c(0, max_val_right),
      type = "l",
      col = colz[8],
      xaxs = "i",
      yaxs = "i",
      axes = FALSE
    )
    graphics::lines(
      x = seq(first_year, last_year, 1),
      y = biocol_data$biocol_overtime_frac,
      ylab = "",
      xlab = "",
      xlim = plot_years,
      ylim = c(0, 0.4),
      col = colz[9],
      xaxs = "i",
      yaxs = "i",
      axes = FALSE
    )
  } else {
    stop(paste0("Unknown value for parameter ref: ", ref, " - Aborting."))
  }

  graphics::axis(side = 4, col = colz[8], col.axis = colz[8])
  graphics::mtext(
    text = "fraction of NPPref",
    col = colz[8],
    side = 4,
    line = 2
  )

  if (!is.null(highlight_years)) {
    for (y in highlight_years) {
      graphics::lines(x = c(y, y), y = c(min_val, max_val), col = "grey40")
    }
  }
  if (details) {
    graphics::legend(
      legendpos,
      legend = c(
        "NPPpot (PNV)", "NPPact (landuse)", "NPPeco", "NPPluc", "NPPharv",
        "HANPP abs sum", "HANPP sum", "BioCol abs sum [frac NPPref]", "BioCol sum [frac NPPref]", # nolint:line_length_linter
        "rharvest", "firec", "timber_harvest", "wood_harvest"
      ), col = colz, lty = 1, cex = 1
    )
  } else {
    graphics::legend(
      legendpos,
      legend = c(
        "NPPpot (PNV)", "NPPact (landuse)", "NPPeco", "NPPluc", "NPPharv",
        "HANPP abs sum", "HANPP sum", "BioCol abs sum [frac NPPref]", "BioCol sum [frac NPPref]" # nolint:line_length_linter
      ), col = colz[seq_len(9)], lty = 1, cex = 1
    )
  }
  grDevices::dev.off()
}