pbratios: Processing Pb isotope ratio data

#' Three-isotope plot for PM10 data.
#'
#' The function prints a three-isotope plot on screen and saves the plot as pdf, png and svg.
#' Measured data are compared with Pb isotope ratios measured in PM10 samples from the Northern
#' Hemisphere. Such data are obtained from the dataset \code{pm10nh} from the package
#' \code{pbratios}.
#'
#' @import data.table
#' @import ggplot2
#' @import scales
#' @import MASS
#' @import svglite
#'
#' @param data A data.table generated by corr.mbf function and subsetted to cointain only samples to
#'   be plotted in the three-isotope diagram. Data should contain a column the column
#'   "\code{Pb208206}, \code{Pb207206} with the measured Pb isotope ratios for 208Pb/206Pb and
#'   207Pb/206Pb, respectively.
#'
#' @param factor A factor within the data.table that can be used differentiate data points with
#'   different point shapes. Default value is \code{NULL}.
#'
#' @param save Logical. If \code{TRUE}, the plot is saved in the "output" folder with the file name
#'   "data_gen" in PDF, PNG and SVG formats.
#'
#' @return The function prints a three-isotope plot on screen. Additionally, if \code{save = TRUE},
#'   three files named "data_gen" are saved in the "output" folder with the file name "data_gen" in
#'   PDF, PNG and SVG formats.
#'
#' @examples
#' isoplot_pm10(pm10nya, save = FALSE)
#'
#' @seealso \code{\link{pm10nh}} \code{\link{isoplot_generic}}
#'
#' @export
isoplot_pm10 <- function(data, factor = NULL, save = FALSE) {

  # Loading data for PM10 in Northern Hemisphere
  data(pm10nh, envir = environment())

  # Defining the plot
  x <- "Pb208206"
  y <- "Pb207206"
  area <- "area"

  p <- ggplot() + geom_point(data = pm10nh, aes_string(x = x, y = y, colour = area),
                             shape = 3) +
    # Confidence ellipses
    stat_ellipse(data = pm10nh, geom = "polygon", type = "norm", alpha = 0.2,
                 aes_string(x = x, y = y, fill = area, colour = area)) +
    # Colours
    scale_colour_manual(values = c("green3", "gold", "dodgerblue3", "red3", "darkblue"),
                        guide = FALSE) +
    scale_fill_manual(values = c("green3", "gold", "dodgerblue3", "red3", "darkblue"),
                      guide = FALSE) +
    # Limits of the plot and breaks
    coord_cartesian(xlim = c(1.99, 2.12), ylim = c(0.80, 0.88)) +
    scale_x_continuous(breaks = pretty_breaks(n = 10)) +
    scale_y_continuous(breaks = pretty_breaks(n = 6)) +
    # Text annotations
    annotate("segment", x = 2.039, xend = 2.042, y = 0.840, yend = 0.838,
             colour = "gray", size = 0.5) +
    annotate("text", x = 2.03396, y = 0.8418, label = "USA", size = 3.53) +
    annotate("segment", x = 2.080, xend = 2.086, y = 0.841, yend = 0.836,
             colour = "gray", size = 0.5) +
    annotate("text", x = 2.08594, y = 0.835, label = "Russia", size = 3.53) +
    annotate("segment", x = 2.100, xend = 2.106, y = 0.850, yend = 0.846,
             colour = "gray", size = 0.5) +
    annotate("text", x = 2.10594, y = 0.8444, label = "China", size = 3.53) +
    annotate("segment", x = 2.080, xend = 2.074, y = 0.860, yend = 0.866,
             colour = "gray", size = 0.5) +
    annotate("text", x = 2.07394, y = 0.8678, label = "Canada", size = 3.53) +
    annotate("segment", x = 2.100, xend = 2.094, y = 0.870, yend = 0.876,
             colour = "gray", size = 0.5) +
    annotate("text", x = 2.0936, y = 0.8778, label = "Europe", size = 3.53)

  # Ellipse for the local particulate
  nat_local <- get_ellipse(2.044, 0.823, a = 0.003, b = 0.001, angle = 0)

  p <- p + geom_polygon(data = nat_local, aes(x = xel, y = yel),
                        fill = "sienna", colour = "sienna", alpha = 0.2) +
    annotate("text", x = 2.060, y = 0.823, label = "Local natural\n end-member",
             size = 3.53)

  # Preparing axis labels
  xlab <- bquote(phantom(.) ^ 208 * Pb ~ "/" * phantom(.) ^ 206 * Pb)
  ylab <- bquote(phantom(.) ^ 207 * Pb ~ "/" * phantom(.) ^ 206 * Pb)

  # Adding the data
  p <- p + geom_point(data = data, aes_string(x = x, y = y, shape = factor),
                      colour = "black", fill = "gray80", size = 1.8) +
    scale_shape_manual(values = c(21, 22, 23, 24, 25, 0, 1, 2, 5, 6), drop = FALSE) +
    # Axis labels and theme
    xlab(xlab) + ylab(ylab) +
    theme_linedraw(base_size = 10) +
    theme(
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      legend.title = element_blank(),
      legend.justification = c(0, 1),
      legend.position = c(0.01, 0.99)
    )

  if (save == TRUE) {

    filename <- gsub("\\[.*?\\]", "\\1", deparse(substitute(data)))
    filename <- gsub("\\.", "", filename)

      if (dir.exists("output") == TRUE) {

      } else {
        dir.create("output")
      }

    ggsave(
    filename = paste0("output/", filename, ".pdf"), plot = p, width = 13, height = 8.03,
    units = "cm", device = "pdf", dpi = 600)

    ggsave(
    filename = paste0("output/", filename, ".png"), plot = p, width = 13, height = 8.03,
    units = "cm", device = "png", dpi = 600)

    ggsave(
    filename = paste0("output/", filename, ".svg"), plot = p, width = 13, height = 8.03,
    units = "cm", device = "svg", dpi = 600)

    p

  } else {
    p
  }

}
# End of function

#' Generating data points on an ellipse.
#'
#' The function creates a dataframe with the points required to plot an ellipse.
#'
#' @param x,y Coordinates for the centre of the ellipse. Default are \code{0, 0}.
#'
#' @param a,b Major and minor semi-axis of the ellipse. Default are \code{1, 1}.
#'
#' @param angle Angle between the x-axis and the major semi-axis. Default is \code{pi/3}.
#'
#' @param n Number of data points in the dataset. Default is \code{300}.
#'
#' @return The function returns a dataset of \eqn{x, y} coordinates required to draw an ellipse with
#'   \code{ggplot2}. This function is used internally by the function \code{isoplot.pm10}.
#'
#' @seealso \code{\link{isoplot_pm10}}
#'
get_ellipse <- function (x = 0, y = 0, a = 1, b = 1, angle = pi / 3, n = 300){

  cc <- exp(seq(0, n) * (0 + 2i) * pi / n)

  R <- matrix(c(cos(angle), sin(angle),-sin(angle), cos(angle)), ncol = 2, byrow = T)

  res <- cbind(x = a * Re(cc), y = b * Im(cc)) %*% R

  data.frame(xel = res[, 1] + x, yel = res[, 2] + y)
}
# End of file