R/util.plot.R

In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry

# CHNOSZ/util.plot.R
# Functions to create and modify plots

thermo.plot.new <- function(xlim,ylim,xlab,ylab,cex=par('cex'),mar=NULL,lwd=par('lwd'),side=c(1,2,3,4),
  mgp=c(1.7,0.3,0),cex.axis=par('cex'),col=par('col'),yline=NULL,axs='i',plot.box=TRUE,
  las=1,xline=NULL, grid = "", col.grid = "gray", ...) {
  # Start a new plot with some customized settings
  thermo <- get("thermo", CHNOSZ)
  # 20120523 store the old par in thermo()$opar
  if(is.null(thermo$opar)) {
    thermo$opar <- par(no.readonly=TRUE)
    assign("thermo", thermo, CHNOSZ)
  }
  # 20090324 mar handling: NULL - a default setting; NA - par's setting
  # 20090413 changed mar of top side from 2 to 2.5
  marval <- c(3, 3.5, 2.5, 1)
  if(identical(mar[1], NA)) marval <- par("mar")
  # 20181007 get mar from the current device (if it exists) and par("mar") is not the default
  if(!is.null(dev.list())) {
    if(!identical(par("mar"), c(5.1, 4.1, 4.1, 2.1))) marval <- par("mar")
  }
  # Assign marval to mar if the latter is NULL or NA
  if(!is.numeric(mar)) mar <- marval
  par(mar=mar,mgp=mgp,tcl=0.3,las=las,xaxs=axs,yaxs=axs,cex=cex,lwd=lwd,col=col,fg=col, ...)
  plot.new()
  plot.window(xlim=xlim,ylim=ylim)
  if(plot.box) box()
  # Labels
  if(is.null(xline)) xline <- mgp[1]
  thermo.axis(xlab,side=1,line=xline,cex=cex.axis,lwd=NULL)
  if(is.null(yline)) yline <- mgp[1]
  thermo.axis(ylab,side=2,line=yline,cex=cex.axis,lwd=NULL)
  # (optional) tick marks
  if(1 %in% side) thermo.axis(NULL,side=1,lwd=lwd, grid = grid, col.grid = col.grid, plot.line = !plot.box)
  if(2 %in% side) thermo.axis(NULL,side=2,lwd=lwd, grid = grid, col.grid = col.grid, plot.line = !plot.box)
  if(3 %in% side) thermo.axis(NULL,side=3,lwd=lwd, plot.line = !plot.box)
  if(4 %in% side) thermo.axis(NULL,side=4,lwd=lwd, plot.line = !plot.box)
}

label.plot <- function(x, xfrac=0.07, yfrac=0.93, paren=FALSE, italic=FALSE, ...) {
  # Make a text label e.g., "(a)" in the corner of a plot
  # xfrac, yfrac: fraction of axis where to put label (default top right)
  # paren: put a parenthesis around the text, and italicize it?
  if(italic) x <- substitute(italic(a), list(a=x))
  if(paren) x <- substitute(group('(',a,')'), list(a=x))
  if(italic | paren) x <- as.expression(x)
  pu <- par('usr')
  text(pu[1]+xfrac*(pu[2]-pu[1]), pu[3]+yfrac*(pu[4]-pu[3]), labels=x, ...)
}

usrfig <- function() {
  # Function to get the figure limits in user coordinates
  # Get plot limits in user coordinates (usr) and as fraction [0,1] of figure region (plt)
  xusr <- par('usr')[1:2]; yusr <- par('usr')[3:4]
  xplt <- par('plt')[1:2]; yplt <- par('plt')[3:4]
  # Linear model to calculate figure limits in user coordinates
  xlm <- lm(xusr ~ xplt); ylm <- lm(yusr ~ yplt)
  xfig <- predict.lm(xlm, data.frame(xplt=c(0, 1)))
  yfig <- predict.lm(ylm, data.frame(yplt=c(0, 1)))
  return(list(x=xfig, y=yfig))
}

label.figure <- function(x, xfrac=0.05, yfrac=0.95, paren=FALSE, italic=FALSE, ...) {
  # Function to add labels outside of the plot region  20151020
  f <- usrfig()
  # Similar to label.plot(), except we have to set xpd=TRUE here
  opar <- par(xpd=NA)
  if(italic) x <- substitute(italic(a), list(a=x))
  if(paren) x <- substitute(group('(',a,')'), list(a=x))
  if(italic | paren) x <- as.expression(x)
  text(f$x[1]+xfrac*(f$x[2]-f$x[1]), f$y[1]+yfrac*(f$y[2]-f$y[1]), labels=x, ...)
  par(opar)
}

water.lines <- function(eout, which=c('oxidation','reduction'),
  lty=2, lwd=1, col=par('fg'), plot.it=TRUE) {

  # Draw water stability limits for Eh-pH, logfO2-pH, logfO2-T or Eh-T diagrams
  # (i.e. redox variable is on the y axis)

  # Get axes, T, P, and xpoints from output of affinity() or equilibrate()
  if(missing(eout)) stop("'eout' (the output of affinity(), equilibrate(), or diagram()) is missing")
  # Number of variables used in affinity()
  nvar1 <- length(eout$vars)
  # If these were on a transect, the actual number of variables is less
  dim <- dim(eout$loga.equil[[1]]) # for output from equilibrate()
  if(is.null(dim)) dim <- dim(eout$values[[1]]) # for output from affinity()
  nvar2 <- length(dim)
  # We only work on diagrams with 1 or 2 variables
  if(!nvar1 %in% c(1,2) | !nvar2 %in% c(1,2)) return(NA)

  # If needed, swap axes so redox variable is on y-axis
  # Also do this for 1-D diagrams 20200710
  if(is.na(eout$vars[2])) eout$vars[2] <- "nothing"
  swapped <- FALSE
  if(eout$vars[2] %in% c("T", "P", "nothing")) {
    eout$vars <- rev(eout$vars)
    eout$vals <- rev(eout$vals)
    swapped <- TRUE
  }
  xaxis <- eout$vars[1]
  yaxis <- eout$vars[2]
  xpoints <- eout$vals[[1]]
  # Make xaxis "nothing" if it is not pH, T, or P 20201110
  # (so that horizontal water lines can be drawn for any non-redox variable on the x-axis)
  if(!identical(xaxis, "pH") & !identical(xaxis, "T") & !identical(xaxis, "P")) xaxis <- "nothing"

  # T and P are constants unless they are plotted on one of the axes
  T <- eout$T
  if(eout$vars[1]=="T") T <- envert(xpoints, "K")
  P <- eout$P
  if(eout$vars[1]=="P") P <- envert(xpoints, "bar")
  # logaH2O is 0 unless given in eout$basis
  iH2O <- match("H2O", rownames(eout$basis))
  if(is.na(iH2O)) logaH2O <- 0 else logaH2O <- as.numeric(eout$basis$logact[iH2O])
  # pH is 7 unless given in eout$basis or plotted on one of the axes
  iHplus <- match("H+", rownames(eout$basis))
  if(eout$vars[1]=="pH") pH <- xpoints
  else if(!is.na(iHplus)) {
    minuspH <- eout$basis$logact[iHplus]
    # Special treatment for non-numeric value (happens when a buffer is used, even for another basis species)
    if(can.be.numeric(minuspH)) pH <- -as.numeric(minuspH) else pH <- NA
  }
  else pH <- 7

  # O2state is gas unless given in eout$basis
  iO2 <- match("O2", rownames(eout$basis))
  if(is.na(iO2)) O2state <- "gas" else O2state <- eout$basis$state[iO2]
  # H2state is gas unles given in eout$basis
  iH2 <- match("H2", rownames(eout$basis))
  if(is.na(iH2)) H2state <- "gas" else H2state <- eout$basis$state[iH2]

  # Where the calculated values will go
  y.oxidation <- y.reduction <- NULL
  if(xaxis %in% c("pH", "T", "P", "nothing") & yaxis %in% c("Eh", "pe", "O2", "H2")) {
    # Eh/pe/logfO2/logaO2/logfH2/logaH2 vs pH/T/P
    if('reduction' %in% which) {
      logfH2 <- logaH2O # usually 0
      if(yaxis=="H2") {
        logK <- suppressMessages(subcrt(c("H2", "H2"), c(-1, 1), c("gas", H2state), T=T, P=P, convert=FALSE))$out$logK
        # This is logfH2 if H2state=="gas", or logaH2 if H2state=="aq"
        logfH2 <- logfH2 + logK
        y.reduction <- rep(logfH2, length.out=length(xpoints))
      } else {
        logK <- suppressMessages(subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", O2state, "gas"), T=T, P=P, convert=FALSE))$out$logK 
        # This is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
        logfO2 <- 2 * (logK - logfH2 + logaH2O)
        if(yaxis=="O2") y.reduction <- rep(logfO2, length.out=length(xpoints))
        else if(yaxis=="Eh") y.reduction <- convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O)
        else if(yaxis=="pe") y.reduction <- convert(convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O), "pe", T=T)
      }
    }
    if('oxidation' %in% which) {
      logfO2 <- logaH2O # usually 0
      if(yaxis=="H2") {
        logK <- suppressMessages(subcrt(c("H2O", "O2", "H2"), c(-1, 0.5, 1), c("liq", "gas", H2state), T=T, P=P, convert=FALSE))$out$logK 
        # This is logfH2 if H2state=="gas", or logaH2 if H2state=="aq"
        logfH2 <- logK - 0.5*logfO2 + logaH2O
        y.oxidation <- rep(logfH2, length.out=length(xpoints))
      } else {
        logK <- suppressMessages(subcrt(c("O2", "O2"), c(-1, 1), c("gas", O2state), T=T, P=P, convert=FALSE))$out$logK 
        # This is logfO2 if O2state=="gas", or logaO2 if O2state=="aq"
        logfO2 <- logfO2 + logK
        if(yaxis=="O2") y.oxidation <- rep(logfO2, length.out=length(xpoints))
        else if(yaxis=="Eh") y.oxidation <- convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O)
        else if(yaxis=="pe") y.oxidation <- convert(convert(logfO2, 'E0', T=T, P=P, pH=pH, logaH2O=logaH2O), "pe", T=T)
      }
    }
  } else return(NA)

  # Now plot the lines
  if(plot.it) {
    if(swapped) {
      if(nvar1 == 1 | nvar2 == 2) {
        # Add vertical lines on 1-D diagram 20200710
        abline(v = y.oxidation[1], lty=lty, lwd=lwd, col=col)
        abline(v = y.reduction[1], lty=lty, lwd=lwd, col=col)
      } else {
        # xpoints above is really the ypoints
        lines(y.oxidation, xpoints, lty=lty, lwd=lwd, col=col)
        lines(y.reduction, xpoints, lty=lty, lwd=lwd, col=col)
      }
    } else {
      lines(xpoints, y.oxidation, lty=lty, lwd=lwd, col=col)
      lines(xpoints, y.reduction, lty=lty, lwd=lwd, col=col)
    }
  }
  # Return the values
  return(invisible(list(xpoints=xpoints, y.oxidation=y.oxidation, y.reduction=y.reduction, swapped=swapped)))
}

mtitle <- function(main, line=0, spacing=1, ...) {
  # Make a possibly multi-line plot title 
  # Useful for including expressions on multiple lines 
  # 'line' is the margin line of the last (bottom) line of the title
  len <- length(main)
  for(i in 1:len) mtext(main[i], line = line + (len - i)*spacing, ...)
}

# Get colors for range of ZC values 20170206
ZC.col <- function(z) {
  # Scale values to [1, 1000]
  z <- z * 999/diff(range(z))
  z <- round(z - min(z)) + 1
  # Diverging (blue - light grey - red) palette
  # dcol <- colorspace::diverge_hcl(1000, c = 100, l = c(50, 90), power = 1)
  # Use precomputed values
  file <- system.file("extdata/cpetc/bluered.txt", package = "CHNOSZ")
  dcol <- read.table(file, as.is=TRUE)[[1]]
  # Reverse the palette so red is at lower ZC (more reduced)
  rev(dcol)[z]
}

# Function to add axes and axis labels to plots,
#   with some default style settings (rotation of numeric labels)
# With the default arguments (no labels specified), it plots only the axis lines and tick marks
#   (used by diagram() for overplotting the axis on diagrams filled with colors).
thermo.axis <- function(lab=NULL,side=1:4,line=1.5,cex=par('cex'),lwd=par('lwd'),col=par('col'), grid = "", col.grid="gray", plot.line=FALSE) {
  if(!is.null(lwd)) {
    for(thisside in side) {

      ## Get the positions of major tick marks
      at <- axis(thisside,labels=FALSE,tick=FALSE) 
      # Get nicer divisions for axes that span exactly 15 units 20200719
      if(thisside %in% c(1,3)) lim <- par("usr")[1:2]
      if(thisside %in% c(2,4)) lim <- par("usr")[3:4]
      if(abs(diff(lim)) == 15) at <- seq(lim[1], lim[2], length.out = 6)
      if(abs(diff(lim)) == 1.5) at <- seq(lim[1], lim[2], length.out = 4)
      # Make grid lines
      if(grid %in% c("major", "both") & thisside==1) abline(v = at, col=col.grid)
      if(grid %in% c("major", "both") & thisside==2) abline(h = at, col=col.grid)
      ## Plot major tick marks and numeric labels
      do.label <- TRUE
      if(missing(side) | (missing(cex) & thisside %in% c(3,4))) do.label <- FALSE
      # col and col.ticks: plot the tick marks but no line (we make it with box() in thermo.plot.new()) 20190416
      # mat: don't plot ticks at the plot limits 20190416
      if(thisside %in% c(1, 3)) pat <- par("usr")[1:2]
      if(thisside %in% c(2, 4)) pat <- par("usr")[3:4]
      mat <- setdiff(at, pat)
      if(plot.line) axis(thisside, at=mat, labels=FALSE, tick=TRUE, lwd=lwd, col.axis=col, col=col)
      else axis(thisside, at=mat, labels=FALSE, tick=TRUE, lwd=lwd, col.axis=col, col = NA, col.ticks = col)
      # Plot only the labels at all major tick points (including plot limits) 20190417
      if(do.label) axis(thisside, at=at, tick=FALSE, col=col)

      ## Plot minor tick marks
      # The distance between major tick marks
      da <- abs(diff(at[1:2]))
      # Distance between minor tick marks
      di <- da / 4
      if(!da %% 3) di <- da / 3
      else if(da %% 2 | !(da %% 10)) di <- da / 5
      # Number of minor tick marks
      if(thisside %in% c(1,3)) {
        ii <- c(1,2) 
        myasp <- par('xaxp')
      } else {
        ii <- c(3,4)
        myasp <- par('yaxp')
      }
      myusr <- par('usr')[ii]
      daxis <- abs(diff(myusr))
      nt <- daxis / di + 1
      ## If nt isn't an integer, it probably
      ## means the axis limits don't correspond
      ## to major tick marks (expect problems)
      ##at <- seq(myusr[1],myusr[2],length.out=nt)
      # Start from (bottom/left) of axis?
      bl <- 1
      #if(myasp[2]==myusr[2]) bl <- 2
      # Is forward direction (top/right)?
      tr <- 1
      if(xor(myusr[2] < myusr[1] , bl==2)) tr <- -1
      #at <- myusr[bl] + tr * di * seq(0:(nt-1))
      # Well all of that doesn't work in a lot of cases,
      # where none of the axis limits correspond to
      # major tick marks. perhaps the following will work
      at <- myusr[1] + tr * di * (0:(nt-1))
      # Apply an offset
      axt <- axTicks(thisside)[1]
      daxt <- (axt - myusr[1])/di
      daxt <- (daxt-round(daxt))*di
      at <- at + daxt
      ## Get the positions of major tick marks and make grid lines
      if(grid %in% c("minor", "both") & thisside==1) abline(v = at, col=col.grid, lty = 3)
      if(grid %in% c("minor", "both") & thisside==2) abline(h = at, col=col.grid, lty = 3)
      tcl <- par('tcl') * 0.5
      at <- setdiff(at, pat)
      if(plot.line) axis(thisside,labels=FALSE,tick=TRUE,lwd=lwd,col.axis=col,at=at,tcl=tcl, col = col)
      else axis(thisside,labels=FALSE,tick=TRUE,lwd=lwd,col.axis=col,at=at,tcl=tcl, col = NA, col.ticks = col)
    }
  }
  # Rotate labels on side axes
  for(thisside in side) {
    if(thisside %in% c(2,4)) las <- 0 else las <- 1
    if(!is.null(lab)) mtext(lab,side=thisside,line=line,cex=cex,las=las)
  }
}

# Function to add transparency to given color 20220223
add.alpha <- function(col, alpha) {
  x <- col2rgb(col)
  newcol <- rgb(x[1], x[2], x[3], maxColorValue = 255)
  newcol <- paste0(newcol, alpha)
  newcol
}