[CHNOSZ-commits] r746 - in pkg/CHNOSZ: . R inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Thu Sep 29 04:40:26 CEST 2022
Author: jedick
Date: 2022-09-29 04:40:26 +0200 (Thu, 29 Sep 2022)
New Revision: 746
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/util.data.R
pkg/CHNOSZ/R/util.expression.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/logB_to_OBIGT.Rd
Log:
Minor formatting changes
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-09-20 06:33:14 UTC (rev 745)
+++ pkg/CHNOSZ/DESCRIPTION 2022-09-29 02:40:26 UTC (rev 746)
@@ -1,6 +1,6 @@
-Date: 2022-09-20
+Date: 2022-09-29
Package: CHNOSZ
-Version: 1.9.9-37
+Version: 1.9.9-38
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2022-09-20 06:33:14 UTC (rev 745)
+++ pkg/CHNOSZ/R/diagram.R 2022-09-29 02:40:26 UTC (rev 746)
@@ -74,12 +74,12 @@
n.balance <- bal$n.balance
balance <- bal$balance
} else n.balance <- rep(1, length(eout$values))
- # in case all coefficients are negative (for bimetal() examples) 20200713
- # e.g. H+ for minerals with basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
+ # In case all coefficients are negative (for bimetal() examples) 20200713
+ # e.g. H+ for minerals with basis(c("Fe+2", "Cu+", "hydrogen sulfide", "oxygen", "H2O", "H+"))
if(all(n.balance < 0)) n.balance <- -n.balance
}
} else if(type %in% rownames(eout$basis)) {
- # to calculate the loga of basis species at equilibrium
+ # To calculate the loga of basis species at equilibrium
if(!missing(groups)) stop("can't plot equilibrium activities of basis species for grouped species")
if(isTRUE(alpha) | is.character(alpha)) stop("equilibrium activities of basis species not available with alpha = TRUE")
plot.loga.basis <- TRUE
@@ -86,40 +86,40 @@
} else if(type=="loga.equil" & !"loga.equil" %in% names(eout)) stop("'eout' is not the output from equil()")
else if(!type %in% c("loga.equil", "loga.balance")) stop(type, " is not a valid diagram type")
- ## consider a different number of species if we're grouping them together
+ ## Consider a different number of species if we're grouping them together
ngroups <- length(groups)
- ## keep the values we plot in one place so they can be modified, plotted and eventually returned
- # unless something happens below, we'll plot the loga.equil from equilibrate()
+ ## Keep the values we plot in one place so they can be modified, plotted and eventually returned
+ # Unless something happens below, we'll plot the loga.equil from equilibrate()
plotvals <- eout$loga.equil
plotvar <- "loga.equil"
- ## handle loga.balance here (i.e. solubility calculations)
+ ## Handle loga.balance here (i.e. solubility calculations)
if(type=="loga.balance") {
plotvals <- list(eout$loga.balance)
plotvar <- "loga.balance"
}
- ## modify default arguments with add = TRUE
+ ## Modify default arguments with add = TRUE
if(add) {
- # change missing fill.NA to transparent
+ # Change missing fill.NA to transparent
if(missing(fill.NA)) fill.NA <- "transparent"
- # change missing limit.water to FALSE 20200819
+ # Change missing limit.water to FALSE 20200819
if(missing(limit.water)) limit.water <- FALSE
}
- ## number of dimensions (T, P or chemical potentials that are varied)
+ ## Number of dimensions (T, P or chemical potentials that are varied)
# length(eout$vars) - the number of variables = the maximum number of dimensions
# length(dim(eout$values[[1]])) - nd=1 if it was a transect along multiple variables
nd <- min(length(eout$vars), length(dim(eout$values[[1]])))
- ## deal with output from affinity()
+ ## Deal with output from affinity()
if(eout.is.aout) {
- # plot property from affinity(), divided by balancing coefficients
+ # Plot property from affinity(), divided by balancing coefficients
plotvals <- lapply(1:length(eout$values), function(i) {
- # we divide by the balancing coefficients if we're working with affinities
- # this is not normalizing the formulas! it's balancing the reactions...
- # normalizing the formulas is done below
+ # We divide by the balancing coefficients if we're working with affinities
+ # This is not normalizing the formulas! it's balancing the reactions...
+ # normalizing the formulas is done below
eout$values[[i]] / n.balance[i]
})
names(plotvals) <- names(eout$values)
@@ -128,7 +128,7 @@
plotvar <- "rank_affinity"
message(paste("diagram: plotting average affinity ranking for", length(plotvals), "groups"))
} else if(plotvar=="A") {
- # we change 'A' to 'A/(2.303RT)' so the axis label is made correctly
+ # We change 'A' to 'A/(2.303RT)' so the axis label is made correctly
# 20171027 use parentheses to avoid ambiguity about order of operations
plotvar <- "A/(2.303RT)"
if(nd==2 & type=="auto") message("diagram: using maximum affinity method for 2-D diagram")
@@ -137,10 +137,10 @@
} else message(paste("diagram: plotting", plotvar, " / n.balance"))
}
- ## use molality instead of activity if the affinity calculation include ionic strength 20171101
+ ## Use molality instead of activity if the affinity calculation include ionic strength 20171101
molality <- "IS" %in% names(eout)
- ## when can normalize and as.residue be used
+ ## When can normalize and as.residue be used
if(normalize | as.residue) {
if(normalize & as.residue) stop("'normalize' and 'as.residue' can not both be TRUE")
if(!eout.is.aout) stop("'normalize' or 'as.residue' can be TRUE only if 'eout' is the output from affinity()")
@@ -149,37 +149,37 @@
else message("diagram: using 'as.residue' in calculation of predominant species")
}
- ## sum affinities or activities of species together in groups 20090524
- # using lapply/Reduce 20120927
+ ## Sum affinities or activities of species together in groups 20090524
+ # Use lapply/Reduce 20120927
if(!missing(groups)) {
- # loop over the groups
+ # Loop over the groups
plotvals <- lapply(groups, function(ispecies) {
- # remove the logarithms
+ # Remove the logarithms ...
if(eout.is.aout) act <- lapply(plotvals[ispecies], function(x) 10^x)
# and, for activity, multiply by n.balance 20170207
else act <- lapply(seq_along(ispecies), function(i) eout$n.balance[ispecies[i]] * 10^plotvals[[ispecies[i]]])
- # sum the activities
+ # then, sum the activities
return(Reduce("+", act))
})
- # restore the logarithms
+ # Restore the logarithms
plotvals <- lapply(plotvals, function(x) log10(x))
- # we also combine the balancing coefficients for calculations using affinity
+ # Combine the balancing coefficients for calculations using affinity
if(eout.is.aout) n.balance <- sapply(groups, function(ispecies) sum(n.balance[ispecies]))
}
- ## calculate the equilibrium logarithm of activity of a basis species
+ ## Calculate the equilibrium logarithm of activity of a basis species
## (such that affinities of formation reactions are zero)
if(plot.loga.basis) {
ibasis <- match(type, rownames(eout$basis))
- # the logarithm of activity used in the affinity calculation
+ # The logarithm of activity used in the affinity calculation
is.loga.basis <- can.be.numeric(eout$basis$logact[ibasis])
if(!is.loga.basis) stop(paste("the logarithm of activity for basis species", type, "is not numeric - was a buffer selected?"))
loga.basis <- as.numeric(eout$basis$logact[ibasis])
- # the reaction coefficients for this basis species
+ # The reaction coefficients for this basis species
nu.basis <- eout$species[, ibasis]
- # the logarithm of activity where affinity = 0
+ # The logarithm of activity where affinity = 0
plotvals <- lapply(1:length(eout$values), function(x) {
- # eout$values is a strange name for affinity ... should be named something like eout$affinity ...
+ # NB. eout$values is a strange name for affinity ... should be named something like eout$affinity ...
loga.basis - eout$values[[x]]/nu.basis[x]
})
plotvar <- type
@@ -186,31 +186,31 @@
}
## alpha: plot fractional degree of formation
- # scale the activities to sum=1 ... 20091017
- # allow scaling by balancing component 20171008
+ # Scale the activities to sum=1 ... 20091017
+ # Allow scaling by balancing component 20171008
if(isTRUE(alpha) | is.character(alpha)) {
- # remove the logarithms
+ # Remove the logarithms
act <- lapply(plotvals, function(x) 10^x)
if(identical(alpha, "balance")) for(i in 1:length(act)) act[[i]] <- act[[i]] * eout$n.balance[i]
- # sum the activities
+ # Sum the activities
sumact <- Reduce("+", act)
- # divide activities by the total
+ # Divide activities by the total
alpha <- lapply(act, function(x) x/sumact)
plotvals <- alpha
plotvar <- "alpha"
}
- ## identify predominant species
+ ## Identify predominant species
predominant <- NA
H2O.predominant <- NULL
if(plotvar %in% c("loga.equil", "alpha", "A/(2.303RT)", "rank_affinity") & type!="saturation") {
pv <- plotvals
- # some additional steps for affinity values, but not for equilibrated activities
+ # Some additional steps for affinity values, but not for equilibrated activities
if(eout.is.aout) {
for(i in 1:length(pv)) {
# TODO: The equilibrium.Rmd vignette shows predominance diagrams using
- # 'normalize' and 'as.residue' settings that are consistent with equilibrate(),
- # but what is the derivation of these equations?
+ # 'normalize' and 'as.residue' settings that are consistent with equilibrate(),
+ # but what is the derivation of these equations?
if(normalize) pv[[i]] <- (pv[[i]] + eout$species$logact[i] / n.balance[i]) - log10(n.balance[i])
else if(as.residue) pv[[i]] <- pv[[i]] + eout$species$logact[i] / n.balance[i]
}
@@ -224,23 +224,23 @@
predominant <- which.pmax(pv)
}
- # show water stability region
+ # Show water stability region
if((is.null(limit.water) | isTRUE(limit.water)) & nd==2) {
wl <- water.lines(eout, plot.it=FALSE)
- # proceed if water.lines produced calculations for this plot
+ # Proceed if water.lines produced calculations for this plot
if(!identical(wl, NA)) {
H2O.predominant <- predominant
- # for each x-point, find the y-values that are outside the water stability limits
+ # For each x-point, find the y-values that are outside the water stability limits
for(i in seq_along(wl$xpoints)) {
ymin <- min(c(wl$y.oxidation[i], wl$y.reduction[i]))
ymax <- max(c(wl$y.oxidation[i], wl$y.reduction[i]))
if(!wl$swapped) {
- # the actual calculation
+ # The actual calculation
iNA <- eout$vals[[2]] < ymin | eout$vals[[2]] > ymax
# assign NA to the predominance matrix
H2O.predominant[i, iNA] <- NA
} else {
- # as above, but x- and y-axes are swapped
+ # As above, but x- and y-axes are swapped
iNA <- eout$vals[[1]] < ymin | eout$vals[[1]] > ymax
H2O.predominant[iNA, i] <- NA
}
@@ -249,12 +249,12 @@
}
}
- ## create some names for lines/fields if they are missing
+ ## Create some names for lines/fields if they are missing
is.pname <- FALSE
onames <- names
if(identical(names, FALSE) | identical(names, NA)) names <- ""
else if(!is.character(names)) {
- # properties of basis species or reactions?
+ # Properties of basis species or reactions?
if(eout$property %in% c("G.basis", "logact.basis")) names <- rownames(eout$basis)
else {
if(!missing(groups)) {
@@ -262,24 +262,24 @@
else names <- names(groups)
}
else names <- as.character(eout$species$name)
- # remove non-unique organism or protein names
+ # Remove non-unique organism or protein names
if(all(grepl("_", names))) {
is.pname <- TRUE
- # everything before the underscore (the protein)
+ # Everything before the underscore (the protein)
pname <- gsub("_.*$", "", names)
- # everything after the underscore (the organism)
+ # Everything after the underscore (the organism)
oname <- gsub("^.*_", "", names)
- # if the pname or oname are all the same, use the other one as identifying name
+ # If the pname or oname are all the same, use the other one as identifying name
if(length(unique(pname))==1) names <- oname
if(length(unique(oname))==1) names <- pname
}
- # append state to distinguish ambiguous species names
+ # Append state to distinguish ambiguous species names
isdup <- names %in% names[duplicated(names)]
if(any(isdup)) names[isdup] <- paste(names[isdup],
" (", eout$species$state[isdup], ")", sep="")
}
}
- # numeric values indicate a subset 20181007
+ # Numeric values indicate a subset 20181007
if(all(is.numeric(onames))) {
if(isTRUE(all(onames > 0))) names[-onames] <- ""
else if(isTRUE(all(onames < 0))) names[-onames] <- ""
@@ -286,10 +286,10 @@
else stop("numeric 'names' should be all positive or all negative")
}
- ## apply formatting to chemical formulas 20170204
+ ## Apply formatting to chemical formulas 20170204
if(all(grepl("_", names))) is.pname <- TRUE
if(format.names & !is.pname) {
- # check if names are a deparsed expression (used in mix()) 20200718
+ # Check if names are a deparsed expression (used in mix()) 20200718
parsed <- FALSE
if(any(grepl("paste\\(", names))) {
exprnames <- parse(text = names)
@@ -297,15 +297,15 @@
parsed <- TRUE
} else {
exprnames <- as.expression(names)
- # get formatted chemical formulas
+ # Get formatted chemical formulas
for(i in seq_along(exprnames)) {
- # don't try to format the names if they have "+" followed by a character
- # (created by mix()ing diagrams with format.names = FALSE);
- # expr.species() can't handle it 20200722
+ # Don't try to format the names if they have "+" followed by a character
+ # (created by mix()ing diagrams with format.names = FALSE);
+ # expr.species() can't handle it 20200722
if(!grepl("\\+[a-zA-Z]", names[i])) exprnames[[i]] <- expr.species(exprnames[[i]])
}
}
- # apply bold or italic
+ # Apply bold or italic
bold <- rep(bold, length.out = length(exprnames))
italic <- rep(italic, length.out = length(exprnames))
for(i in seq_along(exprnames)) {
@@ -312,20 +312,20 @@
if(bold[i]) exprnames[[i]] <- substitute(bold(a), list(a=exprnames[[i]]))
if(italic[i]) exprnames[[i]] <- substitute(italic(a), list(a=exprnames[[i]]))
}
- # only use the expression if it's different from the unformatted names
+ # Only use the expression if it's different from the unformatted names
if(parsed | !identical(as.character(exprnames), names)) names <- exprnames
}
- ## where we'll put extra output for predominance diagrams (namesx, namesy)
+ ## Where we'll put extra output for predominance diagrams (namesx, namesy)
out2D <- list()
- ### now on to the plotting ###
+ ### Now we're getting to the plotting ###
if(plot.it) {
- ### general plot parameters ###
+ ### General plot parameters ###
- ## handle line type/width/color arguments
+ ## Handle line type/width/color arguments
if(is.null(lty)) {
if(type=="loga.balance" | nd==2) lty <- 1
else lty <- 1:ngroups
@@ -337,7 +337,7 @@
if(nd==0) {
### 0-D diagram - bar graph of properties of species or reactions
- # plot setup
+ # Plot setup
if(missing(ylab)) ylab <- axis.label(plotvar, units="", molality=molality)
barplot(unlist(plotvals), names.arg=names, ylab=ylab, cex.names=cex.names, col=col, ...)
if(!is.null(main)) title(main=main)
@@ -347,19 +347,19 @@
### 1-D diagram - lines for properties or chemical activities
xvalues <- eout$vals[[1]]
if(missing(xlim)) xlim <- range(xvalues) # TODO: this is backward if the vals are not increasing
- # initialize the plot
+ # Initialize the plot
if(!add) {
if(missing(xlab)) xlab <- axis.label(eout$vars[1], basis=eout$basis, molality=molality)
if(missing(ylab)) {
ylab <- axis.label(plotvar, units="", molality=molality)
if(plotvar == "rank_affinity") ylab <- "Average affinity ranking"
- # use ppb, ppm, ppt (or log ppb etc.) for converted values of solubility 20190526
+ # Use ppb, ppm, ppt (or log ppb etc.) for converted values of solubility 20190526
if(grepl("solubility.", eout$fun, fixed=TRUE)) {
ylab <- strsplit(eout$fun, ".", fixed=TRUE)[[1]][2]
ylab <- gsub("log", "log ", ylab)
}
}
- # to get range for y-axis, use only those points that are in the xrange
+ # To get range for y-axis, use only those points that are in the xrange
if(is.null(ylim)) {
isx <- xvalues >= min(xlim) & xvalues <= max(xlim)
xfun <- function(x) x[isx]
@@ -369,12 +369,12 @@
if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, cex=cex, mar=mar, yline=yline, side=side, ...)
else plot(0, 0, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, ...)
}
- # draw the lines
+ # Draw the lines
spline.n <- 256 # the number of values at which to calculate splines
if(is.null(spline.method) | length(xvalues) > spline.n) {
for(i in 1:length(plotvals)) lines(xvalues, plotvals[[i]], col=col[i], lty=lty[i], lwd=lwd[i])
} else {
- # plot splines instead of lines connecting the points 20171116
+ # Plot splines instead of lines connecting the points 20171116
spline.x <- seq(xlim[1], xlim[2], length.out=spline.n)
for(i in 1:length(plotvals)) {
spline.y <- splinefun(xvalues, plotvals[[i]], method=spline.method)(spline.x)
@@ -385,35 +385,35 @@
# 20120521: use legend.x=NA to label lines rather than make legend
if(is.na(legend.x)) {
maxvals <- do.call(pmax, pv)
- # label placement and rotation
+ # Label placement and rotation
dx <- rep(dx, length.out = length(plotvals))
dy <- rep(dy, length.out = length(plotvals))
srt <- rep(srt, length.out = length(plotvals))
- # don't assign to adj becuase that messes up the missing test below
+ # Don't assign to adj becuase that messes up the missing test below
alladj <- rep(adj, length.out = length(plotvals))
for(i in 1:length(plotvals)) {
# y-values for this line
myvals <- as.numeric(plotvals[[i]])
- # don't take values that lie close to or above the top of plot
+ # Don't take values that lie close to or above the top of plot
myvals[myvals > ylim[1] + 0.95*diff(ylim)] <- ylim[1]
- # if we're adding to a plot, don't take values that are above the top of this plot
+ # If we're adding to a plot, don't take values that are above the top of this plot
if(add) {
this.ylim <- par("usr")[3:4]
myvals[myvals > this.ylim[1] + 0.95*diff(this.ylim)] <- this.ylim[1]
}
- # the starting x-adjustment
+ # The starting x-adjustment
thisadj <- alladj[i]
- # if this line has any of the overall maximum values, use only those values
+ # If this line has any of the overall maximum values, use only those values
# (useful for labeling straight-line affinity comparisons 20170221)
is.max <- myvals==maxvals
if(any(is.max) & plotvar != "alpha") {
- # put labels on the median x-position
+ # Put labels on the median x-position
imax <- median(which(is.max))
} else {
- # put labels on the maximum of the line
+ # Put labels on the maximum of the line
# (useful for labeling alpha plots)
imax <- which.max(myvals)
- # try to avoid the sides of the plot; take care of reversed x-axis
+ # Avoid the sides of the plot; take care of reversed x-axis
if(missing(adj)) {
if(sign(diff(xlim)) > 0) {
if(xvalues[imax] > xlim[1] + 0.8*diff(xlim)) thisadj <- 1
@@ -424,13 +424,13 @@
}
}
}
- # also include y-offset (dy) and y-adjustment (labels bottom-aligned with the line)
- # .. and srt (string rotation) 20171127
+ # Also include y-offset (dy) and y-adjustment (labels bottom-aligned with the line)
+ # ... and srt (string rotation) 20171127
text(xvalues[imax] + dx[i], plotvals[[i]][imax] + dy[i], labels=names[i], adj=c(thisadj, 0), cex=cex.names, srt=srt[i], font=font, family=family)
}
} else legend(x=legend.x, lty=lty, legend=names, col=col, cex=cex.names, lwd=lwd, ...)
}
- # add a title
+ # Add a title
if(!is.null(main)) title(main=main)
} else if(nd==2) {
@@ -437,10 +437,10 @@
### 2-D diagram - fields indicating species predominance, or contours for other properties
- ### functions for constructing predominance area diagrams
- ## color fill function
+ ### Functions for constructing predominance area diagrams
+ ## Color fill function
fill.color <- function(xs, ys, out, fill, nspecies) {
- # handle min/max reversal
+ # Handle min/max reversal
if(xs[1] > xs[length(xs)]) {
tc <- out
t <- numeric()
@@ -461,20 +461,19 @@
out <- tc
ys <- t
}
- # the z values
+ # The z values
zs <- out
for(i in 1:nrow(zs)) zs[i,] <- out[nrow(zs)+1-i,]
zs <- t(zs)
breaks <- c(-1, 0, 1:nspecies) + 0.5
- # use fill.NA for NA values
+ # Use fill.NA for NA values
zs[is.na(zs)] <- 0
image(x=xs, y=ys, z=zs, col=c(fill.NA, fill), add=TRUE, breaks=breaks, useRaster=TRUE)
}
- ## curve plot function
- # 20091116 replaced plot.curve with plot.line; different
- # name, same functionality, *much* faster
+ ## Curve plot function
+ # 20091116 replaced plot.curve with plot.line; different name, same functionality, *much* faster
plot.line <- function(out, xlim, ylim, dotted, col, lwd, xrange) {
- # plot boundary lines between predominance fields
+ # Plot boundary lines between predominance fields
vline <- function(out, ix) {
ny <- nrow(out)
xs <- rep(ix, ny*2+1)
@@ -481,10 +480,10 @@
ys <- c(rep(ny:1, each=2), 0)
y1 <- out[, ix]
y2 <- out[, ix+1]
- # no line segment inside a stability field
+ # No line segment inside a stability field
iy <- which(y1==y2)
ys[iy*2] <- NA
- # no line segment at a dotted position
+ # No line segment at a dotted position
iyd <- rowSums(sapply(dotted, function(y) ys%%y==0)) > 0
ys[iyd] <- NA
return(list(xs=xs, ys=ys))
@@ -495,10 +494,10 @@
xs <- c(0, rep(1:nx, each=2))
x1 <- out[iy, ]
x2 <- out[iy+1, ]
- # no line segment inside a stability field
+ # No line segment inside a stability field
ix <- which(x1==x2)
xs[ix*2] <- NA
- # no line segment at a dotted position
+ # No line segment at a dotted position
ixd <- rowSums(sapply(dotted, function(x) xs%%x==0)) > 0
xs[ixd] <- NA
return(list(xs=xs, ys=ys))
@@ -515,7 +514,7 @@
}
rx <- (xlim[2] - xlim[1]) / (ncol(out) - 1)
ry <- (ylim[2] - ylim[1]) / (nrow(out) - 1)
- # vertical lines
+ # Vertical lines
xs <- ys <- NA
for(ix in 1:(ncol(out)-1)) {
vl <- vline(out,ix)
@@ -527,7 +526,7 @@
ys <- clipfun(ys, ylim)
if(!is.null(xrange)) xs <- clipfun(xs, xrange)
lines(xs, ys, col=col, lwd=lwd)
- # horizontal lines
+ # Horizontal lines
xs <- ys <-NA
for(iy in 1:(nrow(out)-1)) {
hl <- hline(out, iy)
@@ -540,12 +539,12 @@
if(!is.null(xrange)) xs <- clipfun(xs, xrange)
lines(xs, ys, col=col, lwd=lwd)
}
- ## new line plotting function 20170122
+ ## New line plotting function 20170122
contour.lines <- function(predominant, xlim, ylim, lty, col, lwd) {
- # the x and y values
+ # The x and y values
xs <- seq(xlim[1], xlim[2], length.out=dim(predominant)[1])
ys <- seq(ylim[1], ylim[2], length.out=dim(predominant)[2])
- # reverse any axis that has decreasing values
+ # Reverse any axis that has decreasing values
if(diff(xlim) < 0) {
predominant <- predominant[nrow(predominant):1, ]
xs <- rev(xs)
@@ -554,7 +553,7 @@
predominant <- predominant[, ncol(predominant):1]
ys <- rev(ys)
}
- # the categories (species/groups/etc) on the plot
+ # The categories (species/groups/etc) on the plot
zvals <- na.omit(unique(as.vector(predominant)))
if(is.null(lty.aq) & is.null(lty.cr)) {
@@ -561,20 +560,20 @@
# DEFAULT method: loop over species
for(i in 1:(length(zvals)-1)) {
- # get the "z" values
+ # Get the "z" values
# Use + 0 trick to convert T/F to 1/0 20220524
z <- (predominant == zvals[i]) + 0
z[is.na(z)] <- 0
- # use contourLines() instead of contour() in order to get line coordinates 20181029
+ # Use contourLines() instead of contour() in order to get line coordinates 20181029
cLines <- contourLines(xs, ys, z, levels = 0.5)
if(length(cLines) > 0) {
- # loop in case contourLines returns multiple lines
+ # Loop in case contourLines returns multiple lines
for(k in 1:length(cLines)) {
- # draw the lines
+ # Draw the lines
lines(cLines[[k]][2:3], lty = lty[zvals[i]], col = col[zvals[i]], lwd = lwd[zvals[i]])
}
}
- # mask species to prevent double-plotting contour lines
+ # Mask species to prevent double-plotting contour lines
predominant[z == zvals[i]] <- NA
}
@@ -585,24 +584,24 @@
for(i in 1:(length(zvals)-1)) {
for(j in (i+1):length(zvals)) {
z <- predominant
- # draw contours only for this pair
+ # Draw contours only for this pair
z[!z %in% c(zvals[i], zvals[j])] <- NA
- # give them neighboring values (so we get one contour line)
+ # Give them neighboring values (so we get one contour line)
z[z==zvals[i]] <- 0
z[z==zvals[j]] <- 1
- # use contourLines() instead of contour() in order to get line coordinates 20181029
+ # Use contourLines() instead of contour() in order to get line coordinates 20181029
cLines <- contourLines(xs, ys, z, levels=0.5)
if(length(cLines) > 0) {
- # loop in case contourLines returns multiple lines
+ # Loop in case contourLines returns multiple lines
for(k in 1:length(cLines)) {
- # draw the lines
+ # Draw the lines
mylty <- lty[zvals[i]]
if(!is.null(lty.cr)) {
- # use lty.cr for cr-cr boundaries 20190530
+ # Use lty.cr for cr-cr boundaries 20190530
if(all(grepl("cr", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.cr
}
if(!is.null(lty.aq)) {
- # use lty.aq for aq-aq boundaries 20190531
+ # Use lty.aq for aq-aq boundaries 20190531
if(all(grepl("aq", eout$species$state[c(zvals[i], zvals[j])]))) mylty <- lty.aq
}
lines(cLines[[k]][2:3], lty = mylty, col = col[zvals[i]], lwd = lwd[zvals[i]])
@@ -614,12 +613,12 @@
}
}
- ## label plot function
+ ## To add labels
plot.names <- function(out, xs, ys, xlim, ylim, names, srt, min.area) {
- # calculate coordinates for field labels
- # revisions: 20091116 for speed, 20190223 work with user-specified xlim and ylim
+ # Calculate coordinates for field labels
+ # Revisions: 20091116 for speed, 20190223 work with user-specified xlim and ylim
namesx <- namesy <- rep(NA, length(names))
- # even if 'names' is NULL, we run the loop in order to generate namesx and namesy for the output 20190225
+ # Even if 'names' is NULL, we run the loop in order to generate namesx and namesy for the output 20190225
area.plot <- length(xs) * length(ys)
for(i in seq_along(groups)) {
this <- which(out==i, arr.ind=TRUE)
@@ -626,13 +625,13 @@
if(length(this)==0) next
xsth <- xs[this[, 2]]
ysth <- rev(ys)[this[, 1]]
- # use only values within the plot range
+ # Use only values within the plot range
rx <- range(xlim)
ry <- range(ylim)
xsth <- xsth[xsth >= rx[1] & xsth <= rx[2]]
ysth <- ysth[ysth >= ry[1] & ysth <= ry[2]]
if(length(xsth)==0 | length(ysth)==0) next
- # skip plotting names if the fields are too small 20200720
+ # Skip plotting names if the fields are too small 20200720
area <- max(length(xsth), length(ysth))
frac.area <- area / area.plot
if(!frac.area >= min.area) next
@@ -639,7 +638,7 @@
namesx[i] <- mean(xsth)
namesy[i] <- mean(ysth)
}
- # fields that really exist on the plot
+ # Fields that really exist on the plot
if(!is.null(names)) {
cex <- rep(cex.names, length.out = length(names))
col <- rep(col.names, length.out = length(names))
@@ -653,13 +652,13 @@
text(namesx[i] + dx[i], namesy[i] + dy[i], labels=names[i], cex=cex[i], col=col[i], font=font[i], family=family[i], srt = srt[i])
}
}
- return(list(namesx=namesx, namesy=namesy))
+ return(list(namesx = namesx, namesy = namesy))
}
- ### done with predominance diagram functions
- ### now on to the diagram itself
+ ### Done with supporting functions
+ ### Now we really get to make the diagram itself
- # colors to fill predominance fields
+ # Colors to fill predominance fields
if(is.null(fill)) {
if(length(unique(eout$species$state)) == 1) fill <- "transparent"
else fill <- ifelse(grepl("cr,cr", eout$species$state), "#DEB88788", ifelse(grepl("cr", eout$species$state), "#FAEBD788", "#F0F8FF88"))
@@ -667,8 +666,8 @@
else if(isTRUE(fill[1]=="rainbow")) fill <- rainbow(ngroups)
else if(isTRUE(fill[1] %in% c("heat", "terrain", "topo", "cm"))) fill <- get(paste0(fill[1], ".colors"))(ngroups)
else if(getRversion() >= "3.6.0" & length(fill)==1) {
- # choose an HCL palette 20190411
- # matching adapted from hcl.colors()
+ # Choose an HCL palette 20190411
+ # Matching adapted from hcl.colors()
fx <- function(x) tolower(gsub("[-, _, \\,, (, ), \\ , \\.]", "", x))
p <- charmatch(fx(fill), fx(hcl.pals()))
if(!is.na(p)) {
@@ -678,13 +677,13 @@
}
}
fill <- rep(fill, length.out=ngroups)
- # the x and y values
+ # The x and y values
xs <- eout$vals[[1]]
ys <- eout$vals[[2]]
- # the limits of the calculation; they aren't necessarily increasing, so don't use range()
+ # The limits of the calculation; they aren't necessarily increasing, so don't use range()
xlim.calc <- c(xs[1], tail(xs, 1))
ylim.calc <- c(ys[1], tail(ys, 1))
- # add if(is.null) to allow user-specified limits 20190223
+ # Add if(is.null) to allow user-specified limits 20190223
if(is.null(xlim)) {
if(add) xlim <- par("usr")[1:2]
else xlim <- xlim.calc
@@ -693,7 +692,7 @@
if(add) ylim <- par("usr")[3:4]
else ylim <- ylim.calc
}
- # initialize the plot
+ # Initialize the plot
if(!add) {
if(is.null(xlab)) xlab <- axis.label(eout$vars[1], basis=eout$basis, molality=molality)
if(is.null(ylab)) ylab <- axis.label(eout$vars[2], basis=eout$basis, molality=molality)
@@ -700,17 +699,18 @@
if(tplot) thermo.plot.new(xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab,
cex=cex, cex.axis=cex.axis, mar=mar, yline=yline, side=side, ...)
else plot(0, 0, type="n", xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, ...)
- # add a title
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/chnosz -r 746
More information about the CHNOSZ-commits
mailing list