plot_ecorisk.R

# written by Fabian Stenzel, based on work by Sebastian Ostberg
# 2022-2023 - stenzel@pik-potsdam.de

################# EcoRisk plotting functions ##################

#' Plot distribution of similarity within biomes
#'
#' Function to plot the distribution of similarity within biomes
#'
#' @param data data object with distibution - as returned by
#'             calculateWithInBiomeDiffs for each subcategory of ecorisk.
#'             dim: c(biomes,bins)
#' @param file to write into, if not supplied (default is NULL) write to screen
#' @param biomes_abbrv character vector. abbreviated names of biomes 
#'        (defaults to NULL -> extract dimension names from data)
#' @param scale scaling factor for distribution. defaults to 1
#' @param title character string title for plot, default empty
#' @param legendtitle character string legend title, default empty
#' @param eps write as eps or png (default: FALSE -> png)
#' @param palette color palette to plot EcoRisk with, defaults to the Ostberg
#'        color scheme white-blue-yellow-red
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_biome_internal_distribution <- function(
    data,
    file = NULL,
    biomes_abbrv = NULL,
    scale = 1,
    title = "",
    legendtitle = "",
    eps = FALSE,
    palette = NULL)
  {
  
  if (!is.null(file)) {
    if (eps) {
      file <- strsplit(file, ".", fixed = TRUE)[[1]]
      file <- paste(c(file[seq_len((length(file) - 1))], "eps"), collapse = ".")
      grDevices::ps.options(family = c("Helvetica"), pointsize = 18)
      grDevices::postscript(file,
        horizontal = FALSE, onefile = FALSE, width = 8,
        height = 16, paper = "special"
      )
    } else {
      grDevices::png(file,
        width = 3, height = 6, units = "in", res = 300,
        pointsize = 6, type = "cairo"
      )
    }
  }
  di <- dim(data)
  bins <- di["bin"]
  res <- 1 / bins
  biomes <- di["biome"]
  if (is.null(biomes_abbrv)) biomes_abbrv <- names(data)$biome
  
  if (is.null(palette)) {
    palette <- c("white",RColorBrewer::brewer.pal(9,"YlOrRd"))
  }
  col_index <- floor(seq(res / 2, 1 - res / 2, res) * 10) + 1
  
  graphics::par(mar = c(2, 4, 0, 0), oma = c(0, 0, 0, 0)) # bltr
  graphics::plot(NA,
                 xlim = c(0, 1), ylim = c(0, 20), xlab = "EcoRisk",
                 main = title, axes = FALSE, ylab = ""
  )
  graphics::axis(side = 2, labels = FALSE, at = seq_len(biomes))
  brks <- seq(0, 1, 0.1)
  fields::image.plot(
    legend.only = TRUE, col = palette,
    useRaster = FALSE, breaks = brks, horizontal = TRUE,
    lab.breaks = brks, legend.shrink = 0.925,
    legend.args = list("", side = 3, font = 2, line = 1.5)
  )
  graphics::mtext(biomes_abbrv, side = 2, line = 1, at = seq_len(biomes), las = 2)
  for (b in seq_len(biomes)) {
    graphics::rect(
      xleft = seq(0, 1 - res, res),
      xright = seq(res, 1, res),
      ybottom = b,
      ytop = b + data[b, ] * scale,
      col = palette[col_index]
    )
  }
  if (!is.null(file)) grDevices::dev.off()
}

#' Plot EcoRisk map to file
#'
#' Function to plot a global map of EcoRisk values [0-1] per grid cell to file
#'
#' @param ecorisk_object ecorisk object from which to plot
#' @param plot_dimension which dimension from ecorisk object to plot
#' @param file to write into, if not supplied (default is NULL) write to screen
#' @param focus_biome highlight the biome with this id and desaturate all other
#'                    (default NULL -- no highlight)
#' @param biome_classes to mask the focus_biome from
#' @param title character string title for plot, default empty
#' @param legendtitle character string legend title
#' @param eps write as eps or png
#' @param leg_yes logical. whether to plot legend or not. defaults to TRUE
#' @param leg_scale scaling factor for legend. defaults to 1
#' @param palette color palette to plot EcoRisk with, defaults to the Ostberg
#'        color scheme white-blue-yellow-red
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_ecorisk_map <- function(
    ecorisk_object,
    plot_dimension,
    file = NULL,
    focus_biome = NULL,
    biome_classes = NULL,
    title = "",
    legendtitle = "",
    eps = FALSE,
    title_size = 1,
    leg_yes = TRUE,
    palette = NULL
) {
  data <- ecorisk_object[[plot_dimension]]
  lat <- ecorisk_object$lat
  lon <- ecorisk_object$lon
  if (!is.null(file)) {
    path_write <- dirname(file)
    dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)
    if (eps) {
      file <- strsplit(file, ".", fixed = TRUE)[[1]]
      file <- paste(c(file[seq_len(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"
      )
    }
  }
  brks <- seq(0, 1, 0.1)
  data[data < brks[1]] <- brks[1]
  data[data > brks[length(brks)]] <- brks[length(brks)]
  
  if (is.null(palette)) {
    palette = c("white",RColorBrewer::brewer.pal(9,"YlOrRd"))
  }
  
  if (!is.null(focus_biome)) {
    focus <- data
    focus[!(biome_classes == focus_biome)] <- NA
    palette_low_sat <- grDevices::adjustcolor(palette, alpha.f = 0.25)
    ra_f <- terra::rast(ncols = 720, nrows = 360)
    ra_f[terra::cellFromXY(ra_f, cbind(lon, lat))] <- focus
  }
  
  ra <- terra::rast(ncols = 720, nrows = 360)
  ra[terra::cellFromXY(ra, cbind(lon, lat))] <- data
  range <- range(data)
  extent <- terra::ext(c(-180, 180, -60, 90))
  graphics::par(mar = c(0, 0, 1, 3), oma = c(0, 0, 0, 0), bty = "n")
  
  if (is.null(focus_biome)) {
    terra::plot(ra,
                ext = extent, breaks = brks, col = palette, main = "",
                legend = FALSE, axes = FALSE
    )
  } else {
    terra::plot(ra,
                ext = extent, breaks = brks, col = palette_low_sat,
                main = "", legend = FALSE, axes = FALSE
    )
    terra::plot(ra_f,
                ext = extent, breaks = brks, col = palette, main = "",
                legend = FALSE, axes = FALSE, add = TRUE
    )
  }
  
  title(main = title, line = -2, cex.main = title_size)
  maps::map("world", add = TRUE, res = 0.4, lwd = 0.25, ylim = c(-60, 90))
  
  if (leg_yes) {
    fields::image.plot(
      legend.only = TRUE, col = palette, breaks = brks, zlim = range,
      lab.breaks = brks, legend.shrink = 0.7,
      legend.args = list(legendtitle, side = 3, font = 2, line = 1)
    ) # removed zlim
  }
  if (!is.null(file)) grDevices::dev.off()
}


#' Plot radial EcoRisk plot to screen
#'
#' Function to plot an aggregated radial status of EcoRisk values [0-1]
#' for the different sub-categories to screen
#'
#' @param data EcoRisk data array c([nEcoRiskcomponents],
#'             3[min,mean,max])
#' @param title character string title for plot, default empty
#' @param zoom scaling factor for circle plot. defaults to 1
#' @param type plot type, 'legend1' for variable and color legend,
#'             'legend2' for value legend, or 'regular' (default setting)
#'             for the regular EcoRisk plot
#' @param title_size scaling factor for tile. defaults to 1
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_ecorisk_radial_to_screen <- function(data, # nolint
                                          title = "",
                                          zoom = 1.0,
                                          type = "regular",
                                          title_size = 2,
                                          titleline = -2,
                                          use_quantile = TRUE) {
  ecorisk_dims <- length(data[, 1])
  if (ecorisk_dims == 8) {
    # names <- c(
    #  ecorisk = "ecorisk", deltav = "vegetation\nstructure",
    #  local = "local\nchange", global = "global\nimportance",
    #  balance = "ecosystem\nbalance", cstocks = "carbon\nstocks",
    #  cfluxes = "carbon fluxes", wfluxes = "water fluxes"
    # )

    # take the names of the ecorisk list dimensions, removing "_total"
    names <- gsub("_", "\n", gsub("_total", "", dimnames(data)[[1]]))
    # c(blue-green, yellow, violet, red, blue, orange, green, pink, grey,
    #   purple, green-blue, yellow-orange)

    set <- RColorBrewer::brewer.pal(12, "Set3")
    colz <- set[c(4, 7, 8, 11, 1, 10, 5, 6)]
    #   ecorisk vs         lc        gi        eb        ct       wt         nt
    angles <- matrix(
      c(90, 270, 216, 252, 180, 216, 144, 180, 108, 144, -18, 18, -54, -18, 18, 54),
      byrow = TRUE,
      nrow = length(colz)
    )
  } else {
    stop("Unknown number of dimensions for ecorisk data:", ecorisk_dims)
  }

  graphics::par(oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0))
  graphics::plot(c(-zoom, zoom), c(-zoom, zoom),
    type = "n", axes = FALSE,
    ann = FALSE, asp = 1, main = ""
  )
  graphics::title(main = title, line = titleline, cex.main = title_size)

  if (type == "legend1") {
    circlize::draw.sector(0, 360, rou1 = 1)
    ro <- c(1, 1.1, 0.8, 1.1, 0.8, 1, 1, 1, 1, 1, 1)

    for (i in seq_along(angles[, 1])) {
      circlize::draw.sector(
        start.degree = angles[i, 1] + 90,
        end.degree = angles[i, 2] + 90,
        col = colz[i],
        clock.wise = FALSE,
        rou1 = 0,
        rou2 = ro[i],
        border = "black"
      )
    }

    if (ecorisk_dims == 8) {
      graphics::text(names,
        x = c(1.1, 1.0, 0.2, -0.8, -1.6, -0.25, 0.7, -1.3),
        y = c(-0.15, -0.9, -1.3, -1.3, -0.9, 1.2, 1, 1),
        adj = 0
      )
    } else {
      stop("Unknown number of dimensions for ecorisk data:", ecorisk_dims)
    }

    # line lc
    circlize::draw.sector(
      start.degree = (angles[3, 1] + angles[3, 2]) / 2 + 90,
      end.degree = (angles[3, 1] + angles[3, 2]) / 2 + 90,
      rou1 = 0.7,
      rou2 = 1.1
    )
    # line ecorisk
    circlize::draw.sector(
      start.degree = -9,
      end.degree = -9,
      rou1 = 0.9,
      rou2 = 1.05
    )
    # line eb
    circlize::draw.sector(
      start.degree = (angles[5, 1] + angles[5, 2]) / 2 + 90,
      end.degree = (angles[5, 1] + angles[5, 2]) / 2 + 90, rou1 = 0.7,
      rou2 = 1.1
    )
    circlize::draw.sector(
      start.degree = 180,
      end.degree = 180,
      clock.wise = FALSE,
      rou1 = -1.2,
      rou2 = 1.2,
      border = "black",
      lwd = 2
    )
  } else if (type == "legend2") {
    graphics::text("+", x = 0, y = 0)
    circlize::draw.sector(0, 360, rou1 = 1)
    circlize::draw.sector(0, 360, rou1 = 0.65)
    circlize::draw.sector(0, 360, rou1 = 0.3)
    # sector
    circlize::draw.sector(
      start.degree = -18 + 90,
      end.degree = 18 + 90,
      clock.wise = FALSE,
      rou1 = 0.55,
      border = "black"
    )
    # uncertainty arrow
    circlize::draw.sector(
      start.degree = 90,
      end.degree = 90,
      clock.wise = FALSE,
      rou1 = 0.4,
      rou2 = 0.8,
      border = "black"
    )
    # uncertainty lower
    circlize::draw.sector(
      start.degree = -9 + 90,
      end.degree = 9 + 90,
      clock.wise = FALSE,
      rou1 = 0.8,
      rou2 = 0.8,
      border = "black"
    )
    # uncertainty upper
    circlize::draw.sector(
      start.degree = -9 + 90,
      end.degree = 9 + 90,
      clock.wise = FALSE,
      rou1 = 0.4,
      rou2 = 0.4,
      border = "black"
    )
    # 0.3
    circlize::draw.sector(
      start.degree = 270 - 270,
      end.degree = 270 - 270,
      clock.wise = FALSE,
      rou1 = 0.3,
      rou2 = 1.3,
      border = "black",
      lty = "dashed"
    )
    # 0.65
    circlize::draw.sector(
      start.degree = 280 - 270,
      end.degree = 280 - 270,
      clock.wise = FALSE,
      rou1 = 0.65,
      rou2 = 1.3,
      border = "black",
      lty = "dashed"
    )
    # 1.0
    circlize::draw.sector(
      start.degree = 290 - 270,
      end.degree = 290 - 270,
      clock.wise = FALSE,
      rou1 = 1,
      rou2 = 1.3,
      border = "black",
      lty = "dashed"
    )
    graphics::text(c("0.3", "0.65", "1"),
      x = c(1.4, 1.45, 1.25),
      y = c(0, 0.25, 0.45)
    )

    # plot how the whiskers are calculated
    #   quantile case
    if (use_quantile) {
      graphics::text(c("Q90", "Q50", "Q10"),
        x = c(-0.3, -0.29, -0.26),
        y = c(0.8, 0.48, 0.35),
        cex = 0.7
      )

      # minmeanmax case
    } else {
      graphics::text(c("max", "mean", "min"),
        x = c(-0.3, -0.29, -0.26),
        y = c(0.8, 0.48, 0.35),
        cex = 0.7
      )
    }
  } else if (type == "regular") {
    circlize::draw.sector(180, 360, rou1 = 1, col = "gray80")

    for (i in seq_along(angles[, 1])) {
      mangle <- mean(angles[i, ])
      if (i == 1) mangle <- -98
      dmin <- data[i, 1]
      dmedian <- data[i, 2]
      dmax <- data[i, 3]
      circlize::draw.sector(
        start.degree = angles[i, 1] + 90,
        end.degree = angles[i, 2] + 90, col = colz[i],
        rou1 = dmedian,
        clock.wise = FALSE,
        border = "black"
      )
      # uncertainty arrow
      circlize::draw.sector(
        start.degree = mangle + 90,
        end.degree = mangle + 90,
        clock.wise = FALSE,
        rou1 = dmin,
        rou2 = dmax,
        border = "black"
      )
      circlize::draw.sector(
        start.degree = mangle - 9 + 90,
        end.degree = mangle + 9 + 90,
        clock.wise = FALSE,
        rou1 = dmin,
        rou2 = dmin,
        border = "black"
      )
      # uncertainty upper
      circlize::draw.sector(
        start.degree = mangle - 9 + 90,
        end.degree = mangle + 9 + 90,
        clock.wise = FALSE,
        rou1 = dmax,
        rou2 = dmax,
        border = "black"
      )
    }
    circlize::draw.sector(0, 360, rou1 = 1)
    circlize::draw.sector(0, 360, rou1 = 0.6)
    circlize::draw.sector(0, 360, rou1 = 0.3)
    circlize::draw.sector(
      start.degree = 180,
      end.degree = 180,
      clock.wise = FALSE,
      rou1 = -1.2,
      rou2 = 1.2,
      border = "black",
      lwd = 2
    )
  } else {
    stop(
      "Unknown type ", type,
      ". Please use 'legend1' for variable and color legend,
         'legend2' for value legend, or 'regular' (default setting) ",
      "for the regular ecorisk plot."
    )
  }
}


#' Plot radial EcoRisk plot to file
#'
#' Function to plot an aggregated radial status of EcoRisk values [0-1]
#' for the different sub-categories to file
#'
#' @param data EcoRisk data array c(4/19[biomes],[nEcoRiskcomponents],
#'             3[min,mean,max])
#' @param file to write into
#' @param title character string title for plot, default empty
#' @param type plot type, 'legend1' for variable and color legend,
#'             'legend2' for value legend, or 'regular' (default setting)
#'             for the regular EcoRisk plot
#' @param eps write as eps or png
#' @param leg_yes logical. whether to plot legend or not. defaults to TRUE
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_ecorisk_radial <- function(data,
                                file,
                                title = "",
                                leg_yes = TRUE,
                                eps = FALSE,
                                use_quantile = TRUE) {

  path_write <- dirname(file)
  dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)
  if (length(which(data < 0 | data > 1)) > 0) {
    warning(
      "There are values in data outside the expected EcoRisk range [0..1]." # nolint
    )
  }

  if (eps) {
    file <- strsplit(file, ".", fixed = TRUE)[[1]]
    file <- paste(c(file[seq_len(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"
    )
  }

  # adjust the margins, dependent on whether a legend should be plotted or not
  graphics::par(fig = c(0, 0.7, 0, 1)) # , oma=c(0,0,0,0),mar=c(0,0,0,0))

  # plot main EcoRisk radial
  plot_ecorisk_radial_to_screen(
    data = data, title = title, zoom = 1.0,
    type = "regular"
  )

  if (leg_yes) {
    graphics::par(fig = c(0.7, 1, 0, 0.5), new = TRUE)
    plot_ecorisk_radial_to_screen(
      data = data, title = "", zoom = 1.5,
      type = "legend1"
    )
    graphics::par(fig = c(0.7, 1, 0.5, 1), new = TRUE)
    plot_ecorisk_radial_to_screen(
      data = data, title = "", zoom = 1.5,
      type = "legend2", use_quantile = use_quantile
    )
  }
  grDevices::dev.off()
}


#' Plot timeline of EcoRisk variables to screen
#'
#' Function to plot timeline of EcoRisk variables to screen
#'
#' @param data EcoRisk data array
#'        c(4/19[biomes],8/10[nEcoRiskcomponents],3[min,mean,max],timeslices)
#' @param timerange of the data input
#' @param yrange range for y axis default c(0,1)
#' @param leg_yes plot legend (default TRUE)
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_overtime_to_screen <- function(data,
                                    timerange,
                                    yrange = c(0, 1),
                                    leg_yes = TRUE,
                                    leg_only = FALSE,
                                    varnames = NULL) {
  ecorisk_dims <- dim(data)[1]

  if (is.null(varnames)) {
    if (ecorisk_dims == 10) {
      names <- c(
        ecorisk = "ecorisk", deltav = "vegetation structure",
        local = "local change", global = "global importance",
        balance = "ecosystem balance", cstocks = "carbon stocks",
        cfluxes = "carbon fluxes", wfluxes = "water fluxes",
        nstocks = "nitrogen stocks", nfluxes = "nitrogen fluxes"
      )
      # c(blue-green, yellow, violet, red, blue, orange, green, pink, grey,
      #   purple, green-blue, yellow-orange)
      set <- RColorBrewer::brewer.pal(12, "Set3")
      colz <- set[c(4, 7, 8, 11, 1, 3, 10, 5, 12, 6)]
    } else if (ecorisk_dims == 8) {
      names <- c(
        ecorisk = "ecorisk", deltav = "vegetation structure",
        local = "local change", global = "global importance",
        balance = "ecosystem balance", cstocks = "carbon stocks",
        cfluxes = "carbon fluxes", wfluxes = "water fluxes"
      )
      colz <- c(
        "darkgoldenrod", RColorBrewer::brewer.pal(5, "Greens")[5],
        RColorBrewer::brewer.pal(6, "Set1")[seq(2, 6, by = 2)],
        rev(RColorBrewer::brewer.pal(6, "Oranges")[c(4, 5)]),
        RColorBrewer::brewer.pal(6, "PuBu")[6]
      )
    } else {
      stop("Unknown number of dimensions for ecorisk data: ", ecorisk_dims)
    }
  } else {
    names <- varnames
    colz <- RColorBrewer::brewer.pal(length(names), "Set2")
  }
  years <- timerange[1]:timerange[2]
  if (leg_only) {
    graphics::plot(NA,
      ylim = c(yrange[1], yrange[2]), cex.axis = 1,
      axes = FALSE, xlab = "", ylab = ""
    )

    graphics::legend("center", legend = names, fill = colz, border = colz)
  } else {
    graphics::plot(NA,
      xlim = timerange, ylim = c(yrange[1], yrange[2]),
      cex.axis = 1, xlab = "", ylab = ""
    )
    for (i in seq_len(ecorisk_dims)) {
      if (i == 1) {
        graphics::lines(x = years, y = data[i, 2, ], col = colz[i], lwd = 4)
      } else {
        graphics::lines(x = years, y = data[i, 2, ], col = colz[i], lwd = 2)
      }
    }
    if (leg_yes) graphics::legend("topleft", legend = names, fill = colz)
  }
}

#' Plot timeline of EcoRisk variables as panel to file with 4/16 biomes
#'
#' Function to plot a panel of 4/16 timelines per biome aggregated EcoRisk
#' values [0-1]
#' to file
#'
#' @param data EcoRisk data array c(4/19[biomes],[nEcoRiskcomponents],
#'             3[min,mean,max])
#' @param biome_names names of biomes
#' @param file to write into (if not supplied - default NULL - prints to screen)
#' @param yrange range for y axis (default c(0,1))
#' @param timerange of the data input
#' @param eps write as eps or png
#'
#' @examples
#' \dontrun{
#' plot_biome_internal_distribution(
#'   data = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_ecorisk_over_time_panel <- function(data, 
                                         biome_names,
                                         file = NULL,
                                         yrange = c(0, 1),
                                         timerange,
                                         eps = FALSE,
                                         varnames = NULL) {
  if (!is.null(file)) {
    path_write <- dirname(file)
    dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)
  
    if (length(which(data < 0 | data > 1)) > 0) {
      warning("Values in data outside the expected EcoRisk range [0..1].")
    }
  
    if (eps) {
      file <- strsplit(file, ".", fixed = TRUE)[[1]]
      file <- paste(c(file[seq_len(length(file) - 1)], "eps"), collapse = ".")
      grDevices::ps.options(family = c("Helvetica"), pointsize = 18)
      grDevices::postscript(file,
        horizontal = FALSE, onefile = FALSE, width = 15,
        height = 10, paper = "special"
      )
    } else {
      grDevices::png(file,
        width = 5.25, height = 3.5, units = "in", res = 300,
        pointsize = 6, type = "cairo"
      )
    }
  }
  d <- length(data[, 1, 1, 1])
  graphics::par(oma = c(0, 0, 0, 0), mar = c(3, 2, 0.5, 0))
  if (d == 16 | d == 4) {
    k <- sqrt(d)
    xs <- seq(0, 0.8, length.out = k + 1)
    ys <- seq(0.98, 0, length.out = k + 1)

    for (x in seq_len(k)) {
      for (y in seq_len(k)) {
        if (x == 1 & y == 1) {
          graphics::par(
            fig = c(xs[x], xs[x + 1], ys[y + 1], ys[y]),
            xpd = TRUE
          )
        } else {
          graphics::par(
            fig = c(xs[x], xs[x + 1], ys[y + 1], ys[y]), xpd = TRUE,
            new = TRUE
          )
        }
        plot_overtime_to_screen(
          data = data[(x - 1) * k + y, , , ],
          timerange = timerange, yrange = yrange,
          leg_yes = FALSE, varnames = varnames
        )
        graphics::mtext(
          text = biome_names[(x - 1) * k + y], side = 3,
          line = 0, cex = 1, font = 2
        )
      }
    }
  } else {
    stop(paste("Unknown number of biomes: ", length(data[, 1, 1, 1])))
  }
  # legend

  graphics::par(
    fig = c(0.8, 1, 0.5, 1.0), new = TRUE, oma = c(0, 0, 0, 0),
    mar = c(0, 0, 0, 0)
  )
  graphics::plot(NA, axes = FALSE, ylim = c(0, 1), xlim = c(0, 1))
  if (d == 16) {
    graphics::text(
      x = 0.1,
      y = seq(0.95, 0.05, length.out = length(get_biome_names(1))),
      labels = paste0(get_biome_names(1), " : ", get_biome_names(2)),
      cex = 0.7, adj = 0
    )
  }
  graphics::par(
    fig = c(0.8, 1, 0.0, 0.5), new = TRUE, oma = c(0, 0, 0, 0),
    mar = c(0, 0, 0, 0)
  )

  if (is.null(varnames)) {
    plot_overtime_to_screen(
      data = data[1, , , ], timerange = timerange,
      leg_yes = FALSE, leg_only = TRUE
    )
  } else {
    graphics::plot(NA, axes = FALSE, ylim = c(0, 1), xlim = c(0, 1))
    colz <- RColorBrewer::brewer.pal(length(varnames), "Set2")
    graphics::legend("center", legend = varnames, fill = colz, cex = 1)
  }
  if (!is.null(file)) grDevices::dev.off()
}


#' Plot radial EcoRisk panel to file with 4/16 biomes
#'
#' Function to plot an aggregated radial status of EcoRisk values [0-1]
#' for the different sub-categories to file
#'
#' @param data EcoRisk data array c(4/19[biomes],[nEcoRiskcomponents],
#'             3[min,mean,max])
#' @param biome_names names of biomes
#' @param file to write into (if not supplied - default NULL - prints to screen)
#' @param use_quantile is it quantiles or minmeanmax data? - text for whiskers
#' @param eps write as eps or png
#'
#' @examples
#' \dontrun{
#' plot_ecorisk_radial_panel(
#'   data = ecorisk_disaggregated_full[-c(3,18,19),c(1:5,8,9,13),],
#'   biome_names = biomes,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_ecorisk_radial_panel <- function(data,
                                      biome_names,
                                      file = NULL,
                                      use_quantile = TRUE,
                                      eps = FALSE) {
  if (!is.null(file)) {
    path_write <- dirname(file)
    dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)
  
    if (length(which(data < 0 | data > 1)) > 0) {
      warning("Values in data outside the expected EcoRisk range [0..1].")
    }
  
    if (eps) {
      file <- strsplit(file, ".", fixed = TRUE)[[1]]
      file <- paste(c(file[seq_len(length(file) - 1)], "eps"), collapse = ".")
      grDevices::ps.options(family = c("Helvetica"), pointsize = 18)
      grDevices::postscript(file,
        horizontal = FALSE, onefile = FALSE, width = 15,
        height = 10, paper = "special"
      )
    } else {
      grDevices::png(file,
        width = 5.25, height = 3.5, units = "in", res = 300,
        pointsize = 6, type = "cairo"
      )
    }
  }
  d <- length(data[, 1, 1])
  if (d == 16 | d == 4) {
    k <- sqrt(d)
    xs <- seq(0, 0.6, length.out = k + 1)
    ys <- seq(0.98, 0, length.out = k + 1)
    for (x in seq_len(k)) {
      for (y in seq_len(k)) {
        if (x == 1 & y == 1) {
          graphics::par(
            fig = c(xs[x], xs[x + 1], ys[y + 1], ys[y]),
            xpd = TRUE, oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0)
          )
        } else {
          graphics::par(
            fig = c(xs[x], xs[x + 1], ys[y + 1], ys[y]), xpd = TRUE,
            new = TRUE
          )
        }
        plot_ecorisk_radial_to_screen(
          data = data[(x - 1) * k + y, , ],
          title = "", zoom = 1.0, type = "regular"
        )
        graphics::mtext(
          text = biome_names[(x - 1) * k + y], side = 3,
          line = -0.5, cex = 1, font = 2
        )
      }
    }
  } else {
    stop(paste("Unknown number of biomes: ", length(data[, 1, 1])))
  }

  # legend
  graphics::par(fig = c(0.6, 1, 0.1, 0.6), new = TRUE)
  plot_ecorisk_radial_to_screen(
    data = data[1, , ], title = "",
    zoom = 1.5, type = "legend1"
  )

  graphics::par(fig = c(0.6, 1, 0.5, 1.0), new = TRUE)
  plot_ecorisk_radial_to_screen(
    data = data[1, , ], title = "legend", zoom = 1.5,
    type = "legend2", title_size = 1, use_quantile = use_quantile
  )

  if (!is.null(file)) grDevices::dev.off()
}


#' Plot biomes with mercator projection
#'
#' Function to plot biomes to file (or screen) using mercator projection
#'
#' @param biome_ids biome id as given by classify_biomes
#' @param biome_name_length length of biome names in legend: 1 - abbreviation,
#'        2 - short name, 3 - full biome name
#' @param order_legend legend order: either "plants" to first have forests, then
#'        grasslands, then tundra ..., or "zones" to go from north to south
#'        (default: "plants")
#' @param file to write into (if not supplied - default NULL - prints to screen)
#' @param title character string title for plot, default empty
#' @param title_size size of title in cex units (defaukt: 2)
#' @param leg_yes whether to plot legend (default: True)
#' @param leg_scale size of legend in cex units (default 0.5)
#' @param eps write as eps, replacing png in filename (default: True)
#'
#' @examples
#' \dontrun{
#' plot_biomes_mercator(
#'   biome_ids = biome_classification$biome_id,
#'   file = "./biomes.png"
#' )
#' }
#'
#' @md
#' @export
plot_biomes_mercator <- function(biome_ids,
                        biome_name_length = 1,
                        order_legend = "plants",
                        file = NULL,
                        title = "",
                        title_size = 2,
                        leg_yes = TRUE,
                        leg_scale = 1,
                        eps = FALSE) {
  if (!is.null(file)) {
    path_write <- dirname(file)
    dir.create(file.path(path_write), showWarnings = FALSE, recursive = TRUE)
  
    if (eps) {
      file <- strsplit(file, ".", fixed = TRUE)[[1]]
      file <- paste(c(file[seq_len(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"
      )
    }
  }
  # setting up colors and biome names
  colz <- c(
    # warm
    rev(RColorBrewer::brewer.pal(6, "YlOrBr")),
    rev(RColorBrewer::brewer.pal(9, "YlGn")[c(3, 5, 7, 9)]),
    # cold below forest
    rev(RColorBrewer::brewer.pal(9, "GnBu"))[c(2:4, 6, 8, 9)],
    # "lightblue" # Water
    "white",
    # Rocks & Ice
    "lightgrey",
    # montane Tundra/Grassland
    "pink3"
  )
  
  if (order_legend == "plants") {
    order_legend <- seq_len(19)
  } else if (order_legend == "zones") {
    order_legend <- c(
      1, 2, 9, 10, 11, 3, 4, 5, 6, 12, 13, 14, 7, 8, 15, 16, 17, 18, 19
    )
  } else {
    stop(
      "Unknown value for parameter order_legend (plants or zones) - ",
      "was given as: ", order_legend
    )
  }
  biome_class_cols <- (
    colz[c(1, 2, 7, 8, 9, 10, 13, 12, 3, 4, 5, 14, 15, 16, 19, 11, 6, 17, 18)]
  )
  biome_class_names <- get_biome_names(biome_name_length)
  
  if (!(length(biome_class_names) == length(biome_class_cols))) {
    stop("Size of biome class names and colors do not match -- should be 18.")
  }
  
  # plotting
  brks <- seq(
    min(biome_ids, na.rm = TRUE) - 0.5,
    max(biome_ids, na.rm = TRUE) + 0.5,
    1
  )
  ra <- terra::rast(ncols = 720, nrows = 360)
  range <- range(biome_ids)
  ra[terra::cellFromXY(ra, cbind(lon, lat))] <- biome_ids
  extent <- terra::ext(c(-180, 180, -60, 90))
  graphics::par(mar = c(0, 0, 0, 0), oma = c(0, 0, 0, 0), bty = "n")
  terra::plot(ra,
              ext = extent, breaks = brks, col = biome_class_cols,
              main = "", legend = FALSE, axes = FALSE
  )
  graphics::title(main = title, line = -2, cex.main = title_size)
  if (leg_yes) {
    graphics::legend(
      x = -180, y = 27, legend = biome_class_names[order_legend],
      fill = biome_class_cols[order_legend],
      col = biome_class_cols[order_legend],
      cex = leg_scale, bg = "white", bty = "o"
    )
  }
  maps::map("world", add = TRUE, res = 0.4, lwd = 0.25, ylim = c(-60, 90))
  if (!is.null(file)) grDevices::dev.off()
}


#' Plot radial EcoRisk with 4/16 biomes
#'
#' Function to plot to file (or screen) an aggregated radial status of EcoRisk 
#' values [0-1] for the different sub-categories to file
#'
#' @param data EcoRisk data array c(4[biomes],[nEcoRiskcomponents],