diff --git a/DESCRIPTION b/DESCRIPTION index 7769949822225c8cf34cc04a410759ebc2da47cd..ba49b609b9242ec51abce7bdb4c4a5b57629ea1d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: RedModRphree Title: Utilities Leveraging the R Interface to the PHREEQC Geochemical Solver -Version: 0.3.10 +Version: 0.3.11 Authors@R: c(person(given = "Marco", family = "De Lucia", email = "delucia@gfz-potsdam.de", diff --git a/NAMESPACE b/NAMESPACE index be89e038bc4ec58545aef8fe5711591ecfbc8c9c..fbc9121e2766463974f7222658c58a5e19e71810 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -36,6 +36,7 @@ export(MatplotSingle) export(PQCRunAndDump) export(ParseFormula) export(PlotModsInSample) +export(PlotPourbaix) export(Pourbaix) export(RGetPhases) export(RPhreeExt) diff --git a/R/Rphree_Pourbaix.R b/R/Rphree_Pourbaix.R index c34393772064113e13677c8d147202a38d6fd780..9b5f96894027e053cbc80dc45b065df02f42d926 100644 --- a/R/Rphree_Pourbaix.R +++ b/R/Rphree_Pourbaix.R @@ -1,5 +1,6 @@ -### Licence: LGPL version 2.1 -## Time-stamp: "Last modified 2022-07-01 17:07:49 delucia" +## Functions to compute and plot Pourbaix diagrams + +## Time-stamp: "Last modified 2022-07-14 15:55:47 delucia" ##' @title Revert a string ##' @param x character vector. All element of the vector will be @@ -105,18 +106,20 @@ FormSelectedOutput <- function(sim, element) { " -percent_error true", Totals, Species, SI) - return(SelOut) + SelOut } -##' This function uses PHREEQC to calculate and plot a Pourbaix -##' diagram for the specified system at a given T and P. The -##' geochemical system is specified as a single PHREEQC input -##' specifying one "SOLUTION" while the actual values of pe and pH at -##' which the predominance shall be calculated must be specified by -##' the user at input. To try and prevent numerical instabilities, -##' which are quite frequent, the function only actually performs -##' PHREEQC calculations inside the water stability region. +##' This function uses PHREEQC to calculate (and optionally plot) a +##' Pourbaix diagram for the specified system at a given T and P. The +##' geochemical system is defined as a single PHREEQC input specifying +##' one "SOLUTION" while the actual values of pe and pH at which the +##' predominance shall be calculated must be specified by the user as +##' argument in the function call. To try and prevent numerical +##' instabilities, which are quite frequent, the function only +##' actually performs PHREEQC calculations inside the water stability +##' region. Through parameter \code{offset} the user can request a +##' given offset w.r.t. the limiting lines. ##' ##' @title Pourbaix (pe/pH) diagrams using PHREEQC ##' @param element character of length 1. The element for which the @@ -136,18 +139,6 @@ FormSelectedOutput <- function(sim, element) { ##' species are considered ##' @param suppress optional character vector. Species and minerals ##' contained here are not considered in the diagram -##' @param ann.phases Logical, defaults to TRUE. Do we write the -##' formulas of each phase in the (approximate) center of each -##' corresponding diagram region? -##' @param ann.title Optional string. If specified, it will replace -##' main title of the picture -##' @param ann.sub Optional string. If specified, it will replace -##' subtitle of the picture (usually this specifies the -##' concentrations) -##' @param colors Optional character vector containing the colors for -##' each region. If missing, grDevices::terrain.colors is used -##' internally for plotting. If vector is named, the names are -##' matched, if unnamed the colors are used sequentially ##' @param procs how many CPUs should be used for the computations? ##' Defaults to 1. Parallelization is achieved through DoParallel ##' and foreach, and should work on both *nix and Windows using @@ -156,8 +147,8 @@ FormSelectedOutput <- function(sim, element) { ##' limiting lines of water stability ##' @param plot logical, defaults to TRUE. Do you want to plot the ##' diagram? -##' @param ... further graphic parameters passed to the \code{image} -##' for plotting +##' @param ... further parameters passed to \code{PlotPourbaix} for +##' plotting ##' @return invisibly returns a list containing all calculated ##' simulations, input script for PHREEQC simulations, formulas, ##' names, etc @@ -180,8 +171,9 @@ FormSelectedOutput <- function(sim, element) { ##' par(mfrow=c(1,2)) ##' a <- Pourbaix(element="Ca", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, ##' Patm=1, base=base, first="", ann.title="Custom title", db=llnl.dat) -##' b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, -##' Patm=1, base=base, first="", ann.sub="Custom subtitle", db=llnl.dat) +##' b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), +##' Tc= 25, Patm=1, base=base, first="", ann.sub="Custom subtitle", +##' db=llnl.dat, offset=0.2) ##' ##' ## Now we restrict the speciation to a fixed valence state of the element, ##' ## here Fe(2) and Cu(2), only aqueous species: @@ -212,14 +204,14 @@ Pourbaix <- function(element, db, aqonly=FALSE, suppress, - ann.phases=TRUE, - ann.title, - ann.sub, - colors, procs=1, offset=0, plot=TRUE, ...) { + + ## store parameters + store_call <- as.list(match.call(expand.dots = FALSE)) + phreeqc::phrLoadDatabaseString(db) phreeqc::phrSetOutputStringsOn(TRUE) @@ -240,11 +232,14 @@ Pourbaix <- function(element, aa <- phreeqc::phrRunString(c(first, base)) con <- phreeqc::phrGetOutputStrings() + + has_element <- FALSE has_valence <- FALSE if (missing(element)) { utils::capture.output(res <- ReadOutVal(con)[[1]]) selout <- FormSelectedOutput(res) } else { + has_element <- TRUE ## first we distinguish if we want total amount of element or a specific valence if (length(grep("(", element, fixed=TRUE))>0) { ## now we call "valence" the full expression (including the parentheses) @@ -411,36 +406,38 @@ Pourbaix <- function(element, vec <- dat$Aqname } - todisplay <- stats::na.omit(unique(vec)) + displayed <- stats::na.omit(unique(vec)) if (!missing(suppress)) - todisplay <- todisplay[! todisplay %in% suppress] + displayed <- displayed[! displayed %in% suppress] - naqueous <- length(aq <- which(todisplay %in% unique(Aqmaxnam))) + naqueous <- length(aq <- which(displayed %in% unique(Aqmaxnam))) ## fonts 1 for eq phases, font 4 for aqueous - is.aqueous <- rep(1, length(todisplay)) - formulas <- rep("", length(todisplay)) + is.aqueous <- rep(1, length(displayed)) + formulas <- rep("", length(displayed)) is.aqueous[aq] <- 4 if (!aqonly) { - neqphases <- length(eq <- which(todisplay %in% unique(SImaxnam))) + neqphases <- length(eq <- which(displayed %in% unique(SImaxnam))) ## formulas of eqphases - is.pp <- todisplay[is.aqueous==1] + is.pp <- displayed[is.aqueous==1] formulas[is.aqueous==1] <- res$SI$formula[match(is.pp, rownames(res$SI))] - } else eq <- NA + } else { + eq <- NA + } - - numvec <- match(vec, todisplay)## so that vec==todisplay[numvec] - levs <- max(numvec, na.rm=TRUE) + molalities <- paste(rownames(res$tot), res$tot[,1], collapse="; ") + numvec <- match(vec, displayed)## so that vec==displayed[numvec] + levels <- max(numvec, na.rm=TRUE) mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) - centers <- matrix(0, ncol=2, nrow=length(todisplay)) - for (i in seq_along(todisplay)) { + centers <- matrix(0, ncol=2, nrow=length(displayed)) + for (i in seq_along(displayed)) { centers[i, ] <- floor(colMeans(which(mat == i, arr.ind = TRUE))) } ## Fix the names - parsednames <- todisplay + parsednames <- displayed legexpr <- FormulaToExpression(formulas[is.aqueous==1]) Rformulas <- parse(text = paste0("italic(", legexpr,")")) aqlegexpr <- FormulaToExpression(parsednames[is.aqueous==4]) @@ -449,83 +446,155 @@ Pourbaix <- function(element, on.exit() - ## plotting + ret <- list(CALL=store_call, + pe=pe, + pH=pH, + Tc=Tc, + Patm=Patm, + mat=mat, + molalities=molalities, + centers=centers, + displayed=displayed, + parsednames=parsednames, + Rformulas=Rformulas, + is.aqueous=is.aqueous, + aqonly=aqonly, + phreeqc=list(input=c(first, selout), pqc=bigres, db=db, res=res), + #dat=dat, + vec=vec, + aq=aq, + eq=eq, + levels=levels, + suppr=suppr, + has_valence=has_valence, + has_element=has_element, + valence=ifelse(has_valence, valence, ""), + element=ifelse(has_element, element, "")) + if (plot) { - if (missing(colors)) { - region_colors <- grDevices::terrain.colors(levs, alpha = 0.5) - names(region_colors) <- todisplay - } else if (!is.null(names(colors))) { - region_colors <- stats::na.omit(colors[match(todisplay, names(colors))]) - ## we now match the mat to the actual names of colors - numvec <- match(vec, names(region_colors))## so that vec==todisplay[numvec] - levs <- max(numvec, na.rm=TRUE) - mat <- matrix(numvec, ncol=length(pe), nrow=length(pH), byrow=TRUE) - } else { - region_colors <- colors - } - - - par(mar = c(5,5,5,5), las=1, family="serif") - - image(x=pH, y=pe, z=mat, col = region_colors, useRaster=TRUE, axes=FALSE, ...) - - ## Do we want to plot the region names? - if (ann.phases) { - text(pH[centers[,1]], pe[centers[,2]], parsednames, font=is.aqueous, pos=3) - if (!aqonly) { - ii <- which(is.aqueous==1) - ## we need also to check if ii is not empty! - if (length(ii)>0) { - ## we plot these labels - text(pH[centers[ii,1]], pe[centers[ii,2]], Rformulas, font=3, pos=1) + msg("Invoking PlotPourbaix") + PlotPourbaix(ret, ...) + invisible(ret) + } else { + ret + } +} + + + +##' @title Plot a Pourbaix diagram +##' @param pbx list returned by \code{Pourbaix()} +##' @param ann.phases Logical, defaults to TRUE. Do we write the +##' formulas of each phase in the (approximate) center of each +##' corresponding diagram region? +##' @param ann.title Optional string. If specified, it will replace +##' main title of the picture +##' @param ann.sub Optional string. If specified, it will replace +##' subtitle of the picture (usually this specifies the +##' concentrations) +##' @param colors Optional character vector containing the colors for +##' each region. If missing, grDevices::terrain.colors is used +##' internally for plotting. If vector is named, the names are +##' matched, if unnamed the colors are used sequentially +##' @param ... further parameters passed to \code{image} +##' @return nothing +##' @author MDL +##' @export +##' @examples +##' base <- c("SOLUTION 1", +##' "units mol/kgw", +##' "temp 25.0", +##' "pressure 1", +##' "-water 1", +##' "pH 7", +##' "pe 4.0", +##' "Ca 0.025", +##' "Cu 0.001", +##' "Cl 0.05122 charge", +##' "Fe 3.0E-04", +##' "END") +##' +##' b <- Pourbaix(pe=seq(-10,12, length=51), pH=seq(2,12, length=51), +##' Patm=1, Tc=25, base=base, db=llnl.dat, plot=FALSE) +##' +##' mycol <- c("Cu"="grey", "Cl-"="white","Hematite"="light blue", "Nantokite"="red", +##' "Delafossite"="orange", "Ferrite-Cu"="gold", "Magnetite"="light green", +##' "Atacamite"="grey90") +##' +##' par(mfrow=c(1,2)) +##' PlotPourbaix(b, ann.title = "Default colors") +##' PlotPourbaix(b, colors = mycol, ann.title = "Custom colors via named vector") +PlotPourbaix <- function(pbx, + ann.phases=TRUE, + ann.title, + ann.sub, + colors, + ...) { + if (missing(colors)) { + region_colors <- grDevices::terrain.colors(pbx$levels, alpha = 0.5) + names(region_colors) <- pbx$displayed + } else if (!is.null(names(colors))) { + region_colors <- stats::na.omit(colors[match(pbx$displayed, names(colors))]) + ## we now match the mat to the actual names of colors + numvec <- match(pbx$vec, names(region_colors))## so that vec==displayed[numvec] + mat <- matrix(numvec, ncol=length(pbx$pe), nrow=length(pbx$pH), byrow=TRUE) + } else { + region_colors <- colors + } + + par(mar = c(5,5,5,5), las=1, family="serif") + + image(x=pbx$pH, y=pbx$pe, z=pbx$mat, col = region_colors, useRaster=TRUE, xlab="pH", ylab="pe", axes=FALSE, ...) + + ## Do we want to plot the region names? + if (ann.phases) { + text(pbx$pH[pbx$centers[,1]], pbx$pe[pbx$centers[,2]], pbx$parsednames, font=pbx$is.aqueous, pos=3) + if (!pbx$aqonly) { + ii <- which(pbx$is.aqueous==1) + ## we need also to check if ii is not empty! + if (length(ii)>0) { + ## we plot these labels + text(pbx$pH[pbx$centers[ii,1]], pbx$pe[pbx$centers[ii,2]], pbx$Rformulas, font=3, pos=1) } - } } - - if (missing(ann.title)) { - mymain <- paste("Pourbaix diagram for", - ifelse(missing(element), "the specified system", ifelse(has_valence, valence, element)), - "\n T=", Tc, "\u00B0C / P=", Patm, "atm") - } else - mymain <- ann.title - - if (missing(ann.sub)) { - mysub <- paste("Molalities of dissolved elements:", paste(rownames(res$tot),res$tot[,1], collapse="; ")) - } else - mysub <- ann.sub - - title(main = mymain, - sub = mysub) - - axis(1, at=pretty(pH)) - axis(2, at=pretty(pe)) - - ## stability lines of water - abline(0, -1, lty="dashed") - abline(20.8, -1, lty="dashed") - - F <- 96485.3365 - - Eh <- seq(-1,1,0.25) - ppe <- Eh / 2.3 / 8.3145 / { Tc + 273.15 } * F - ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F - ## par(new = TRUE) - ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") - axis(side=4, at = ppe, labels=Eh) - ## axis(3, at=1/{Tgrid+273.15}, labels=(Tgrid)) - - mtext("Eh / V", side=4, line=3, las=0) - box() - } ## plotting + } - invisible(list(mat=mat, - displayed=todisplay, - formula=parsednames, - is.aq=is.aqueous, - phreeqc=list(input=c(first, selout), pqc=bigres, db=db, res=res), - dat=dat, - vec=vec, - aq=aq, - eq=eq)) -} + if (missing(ann.title)) { + mymain <- paste("Pourbaix diagram for", + ifelse(pbx$has_element, "the specified system", ifelse(pbx$has_valence, pbx$valence, pbx$element)), + "\n T=", pbx$Tc, "\u00B0C / P=", pbx$Patm, "atm") + } else { + mymain <- ann.title + } + if (missing(ann.sub)) { + mysub <- paste("Molalities of dissolved elements:", pbx$molalities) + } else { + mysub <- ann.sub + } + + title(main = mymain, sub = mysub) + + axis(1, at=pretty(pbx$pH)) + axis(2, at=pretty(pbx$pe)) + + ## stability lines of water + abline(0, -1, lty="dashed") + abline(20.8, -1, lty="dashed") + + Faraday <- 96485.3365 + + Eh <- seq(-1,1,0.25) + ppe <- Faraday * Eh / 2.3 / 8.3145 / (pbx$Tc + 273.15) + ## Eh <- 2.3*8.3145*{Tc+273.15}* ppe/F + ## par(new = TRUE) + ## plot(0, 0, type = "n", axes = FALSE, bty = "n", xlab = "", ylab = "") + axis(side=4, at = ppe, labels=Eh) + ## axis(3, at=1/{Tgrid+273.15}, labels=(Tgrid)) + + mtext("Eh / V", side=4, line=3, las=0) + box() + + invisible() +} + diff --git a/man/PlotPourbaix.Rd b/man/PlotPourbaix.Rd new file mode 100644 index 0000000000000000000000000000000000000000..5d4f50e236f6d26f86065b7324f4e6aa9e8ca61e --- /dev/null +++ b/man/PlotPourbaix.Rd @@ -0,0 +1,63 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/Rphree_Pourbaix.R +\name{PlotPourbaix} +\alias{PlotPourbaix} +\title{Plot a Pourbaix diagram} +\usage{ +PlotPourbaix(pbx, ann.phases = TRUE, ann.title, ann.sub, colors, ...) +} +\arguments{ +\item{pbx}{list returned by \code{Pourbaix()}} + +\item{ann.phases}{Logical, defaults to TRUE. Do we write the +formulas of each phase in the (approximate) center of each +corresponding diagram region?} + +\item{ann.title}{Optional string. If specified, it will replace +main title of the picture} + +\item{ann.sub}{Optional string. If specified, it will replace +subtitle of the picture (usually this specifies the +concentrations)} + +\item{colors}{Optional character vector containing the colors for +each region. If missing, grDevices::terrain.colors is used +internally for plotting. If vector is named, the names are +matched, if unnamed the colors are used sequentially} + +\item{...}{further parameters passed to \code{image}} +} +\value{ +nothing +} +\description{ +Plot a Pourbaix diagram +} +\examples{ +base <- c("SOLUTION 1", + "units mol/kgw", + "temp 25.0", + "pressure 1", + "-water 1", + "pH 7", + "pe 4.0", + "Ca 0.025", + "Cu 0.001", + "Cl 0.05122 charge", + "Fe 3.0E-04", + "END") + +b <- Pourbaix(pe=seq(-10,12, length=51), pH=seq(2,12, length=51), + Patm=1, Tc=25, base=base, db=llnl.dat, plot=FALSE) + +mycol <- c("Cu"="grey", "Cl-"="white","Hematite"="light blue", "Nantokite"="red", + "Delafossite"="orange", "Ferrite-Cu"="gold", "Magnetite"="light green", + "Atacamite"="grey90") + +par(mfrow=c(1,2)) +PlotPourbaix(b, ann.title = "Default colors") +PlotPourbaix(b, colors = mycol, ann.title = "Custom colors via named vector") +} +\author{ +MDL +} diff --git a/man/Pourbaix.Rd b/man/Pourbaix.Rd index 14d4ba1a6e9e0ecf9ec83135cf16fd89c3a1ea42..7c0219c298ea0e1715a14e2c7365e55d70fe29cf 100644 --- a/man/Pourbaix.Rd +++ b/man/Pourbaix.Rd @@ -15,10 +15,6 @@ Pourbaix( db, aqonly = FALSE, suppress, - ann.phases = TRUE, - ann.title, - ann.sub, - colors, procs = 1, offset = 0, plot = TRUE, @@ -53,22 +49,6 @@ species are considered} \item{suppress}{optional character vector. Species and minerals contained here are not considered in the diagram} -\item{ann.phases}{Logical, defaults to TRUE. Do we write the -formulas of each phase in the (approximate) center of each -corresponding diagram region?} - -\item{ann.title}{Optional string. If specified, it will replace -main title of the picture} - -\item{ann.sub}{Optional string. If specified, it will replace -subtitle of the picture (usually this specifies the -concentrations)} - -\item{colors}{Optional character vector containing the colors for -each region. If missing, grDevices::terrain.colors is used -internally for plotting. If vector is named, the names are -matched, if unnamed the colors are used sequentially} - \item{procs}{how many CPUs should be used for the computations? Defaults to 1. Parallelization is achieved through DoParallel and foreach, and should work on both *nix and Windows using @@ -80,8 +60,8 @@ limiting lines of water stability} \item{plot}{logical, defaults to TRUE. Do you want to plot the diagram?} -\item{...}{further graphic parameters passed to the \code{image} -for plotting} +\item{...}{further parameters passed to \code{PlotPourbaix} for +plotting} } \value{ invisibly returns a list containing all calculated @@ -89,14 +69,16 @@ invisibly returns a list containing all calculated names, etc } \description{ -This function uses PHREEQC to calculate and plot a Pourbaix -diagram for the specified system at a given T and P. The -geochemical system is specified as a single PHREEQC input -specifying one "SOLUTION" while the actual values of pe and pH at -which the predominance shall be calculated must be specified by -the user at input. To try and prevent numerical instabilities, -which are quite frequent, the function only actually performs -PHREEQC calculations inside the water stability region. +This function uses PHREEQC to calculate (and optionally plot) a +Pourbaix diagram for the specified system at a given T and P. The +geochemical system is defined as a single PHREEQC input specifying +one "SOLUTION" while the actual values of pe and pH at which the +predominance shall be calculated must be specified by the user as +argument in the function call. To try and prevent numerical +instabilities, which are quite frequent, the function only +actually performs PHREEQC calculations inside the water stability +region. Through parameter \code{offset} the user can request a +given offset w.r.t. the limiting lines. } \examples{ base <- c("SOLUTION 1", @@ -115,8 +97,9 @@ base <- c("SOLUTION 1", par(mfrow=c(1,2)) a <- Pourbaix(element="Ca", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, Patm=1, base=base, first="", ann.title="Custom title", db=llnl.dat) -b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), Tc= 25, - Patm=1, base=base, first="", ann.sub="Custom subtitle", db=llnl.dat) +b <- Pourbaix(element="Mg", pe=seq(-15,15, length=31), pH=seq(2,12, length=31), + Tc= 25, Patm=1, base=base, first="", ann.sub="Custom subtitle", + db=llnl.dat, offset=0.2) ## Now we restrict the speciation to a fixed valence state of the element, ## here Fe(2) and Cu(2), only aqueous species: