R/piper.R

.packageName <- "hydrogeo"
# Copyright Myles English 2008-2016
# Contact: myles@rockhead.biz
# License: BSD_2_clause

#' @import methods
NULL

#' @importFrom graphics plot points Axis
NULL

#' Major ions as a percentage of total major ions
#'
#' Expects certain column names
#'
#' @param d list or data.frame with the following columns:
#' Ca, Mg, Na, K and Cl, SO4, CO3, HCO3
#' @examples
#' library(hydrogeo)
#' l <- list( Ca = c(43,10,73,26,32),
#'           Mg = c(30,50,83,14,62),
#'           Na = c(54,76,3,14,12),
#'           K = c(31,22,32,22,11),
#'           Cl = c(24,10,12,30,43),
#'           SO4 = c(24,10,12,30,43),
#' 	  CO3 = c(24,10,12,30,43),
#' 	  HCO3 = c(42,110,12,3,4),
#'           IDs = c("A","B","C","D","E") )
#' d <- toPercent(l)
#' # check, should add up to 100%
#' z <- as.data.frame(d)
#' for(i in 1:length(z[[1]])) { print(sum(z[i,5:8])) }
#' for(i in 1:length(z[[1]])) { print(sum(z[i,1:4])) }
#' @export
toPercent <- function(d){
    totalCations <- d$Ca + d$Mg + d$Na + d$K
    d$Ca <- 100 * (d$Ca / totalCations)
    d$Mg <- 100 * (d$Mg / totalCations)
    d$Na <- 100 * (d$Na / totalCations)
    d$K <- 100 * (d$K / totalCations)

    totalAnions<- d$Cl + d$SO4 + d$CO3 + d$HCO3
    d$Cl <- 100 * (d$Cl / totalAnions)
    d$SO4 <- 100 * (d$SO4 / totalAnions)
    d$CO3 <- 100 * (d$CO3 / totalAnions)
    d$HCO3 <- 100 * (d$HCO3 / totalAnions)
    return(d)
}

#' Major ions as a percentage of total major ions - Test Data
#'
#' @param n Number of test samples to be generated.
#' @importFrom stats runif
#' @examples
#' library(hydrogeo)
#' lp <- piper( testData(26) )
#' @export
testData <- function(n) {
      Ca <- runif(n, 0, 100)
      Mg <- runif(n, 0, 100-Ca)
      Cl <- runif(n, 0, 100)
      SO4 <- runif(n, 0, 100-Cl)
      list(Ca=Ca, Mg=Mg, Cl=Cl, SO4=SO4)
}

#' Class \code{piperplot}
#'
#' Objects of this class are plottable as empty (i.e. no points) Piper-Hill
#' diagrams
#'
#' @slot size Object of class \code{numeric} --- Length of the (square) plot
#' area, defaults to 300
#' @slot call R call that created it
setClass("piperplot",
         slots = list(size="numeric",
            ## plotted="logical",
             call="call"),
         prototype = list(size=300)##,plotted=FALSE)
         )

#' Class \code{piper}
#'
#' Objects of this class are plotable as Piper-Hill diagrams.  A dataframe of
#' major ions as percentages can be used to initialise a piper object.
#'
#' @slot Ca Object of class \code{vector} --- Calcium
#' @slot Mg Object of class \code{vector} --- Magnesium
#' @slot Cl Object of class \code{vector} --- Chloride
#' @slot SO4 Object of class \code{vector} --- Sulphate
#' @slot anion.x x coordinate of the point on the anion triangle (internal)
#' @slot anion.y y coordinate of the point on the anion triangle (internal)
#' @slot cation.x x coordinate of the point on the cation triangle (internal)
#' @slot cation.y y coordinate of the point on the cation triangle (internal)
#' @slot diamond.x x coordinate of the point on the diamond (internal)
#' @slot diamond.y y coordinate of the point on the anion diamond (internal)
#' @slot IDs Object of class \code{vector} of sample identifiers
#' @slot pt.col Object of class \code{vector} of colours for points
#' @slot pt.pch Object of class \code{vector} of symbols for  points
#' @slot call Object of class \code{character} --- call that created it
#' @references A. Zaporozec, ``Graphical interpretation of water quality
#' data,'' Ground Water 10, no. 2 (1972): 32--43.
#' @author Myles English \email{myles@@rockhead.biz}
#' @keywords classes
#' @examples
#' showClass("piper")
#'
#' l <- list( Ca = c(43,10,73,26,32),
#'            Mg = c(30,50,3,14,12),
#'            Cl = c(24,10,12,30,43),
#'            SO4 = c(24,10,12,30,43))
#'
#' lp <- piper(l)
#' plot( lp, main="Piper-Hill Diagram of Water Quality" )
# \dontrun{
# x = identify( c(p@@anion.x,p@@cation.x,p@@diamond.x),
#               c(p@@anion.y,p@@cation.y,p@@diamond.y),
#               p@@IDs)
#
# Now, click on the plot near some points with the LH mouse button, press the
# RH mouse button when finished. Vector x now contains the indicies of the
# points, however there are three points for each water sample, so, to look
# up the IDs in the piper object, the indicies are modulo the number of
# samples, i.e.:
#
# p@@IDs[ x %% length(p@@IDs) ]
# }
# @export piper
setClass("piper",
         slots = list(  Ca="vector",
                        Mg="vector",
                        Cl="vector",
                        SO4="vector",
                        IDs="vector",       # NOTE: row.names changed to IDs
                        pt.col="vector",
                        pt.pch="vector",
             anion.x="vector",
             anion.y="vector",
             cation.x="vector",
             cation.y="vector",
             diamond.x="vector",
             diamond.y="vector"
                      ),
         contains = "piperplot"
         )

cationy = function (Mg, size) { c(Mg * 5 * size/1100) }
cationx = function (Ca, Mg, size) { c((5 * size/11) * (1 - (Ca/100) - (Mg/200))) }

#' @describeIn piper Initialiser
# Sample ID is shown by colour and pch via legend with IDs
#' @param .Object object of class piper
#' @param l list of data, see 'Examples' below
# as per \code{\linkS4class{piper}}
#' @param ... additional arguments, as for \code{\linkS4class{piper}}
#' @param call the call that asked for the new piper object
#' @param pt.col Object of class \code{vector} of colours for points
setMethod(f="initialize",
          signature(.Object = "piper"),
          function(.Object, l, ..., call=NULL, pt.col=NULL)
          {
            #TODO: test these checks work
            if ( any((l$Ca + l$Mg)>100) ) {
                stop("Ca + Mg must be <= 100")
            }
            if ( any((l$Cl + l$SO4)>100) ) {
                stop("Cl + SO4 must be <= 100")
            }
            .Object@Ca <- l$Ca
            .Object@Mg <- l$Mg
            .Object@Cl <- l$Cl
            .Object@SO4 <- l$SO4

            if ( !is.null(l$pt.pch) ) {
                .Object@pt.pch <- l$pt.pch
            } else {
                .Object@pt.pch <- rep_len( 0:25, length(l$Ca) )
            }
            if ( !is.null(l$pt.col) ) {
                .Object@pt.col <- l$pt.col
            } else {
                .Object@pt.col <- rep_len( 0:25, length(l$Ca) )
            }

            if ( !is.null(l$IDs) ) {
                .Object@IDs <- l$IDs
            } else {
                if ( !is.null(row.names(l)) ) {
                    .Object@IDs <- row.names(l)
                } else {
                    .Object@IDs <- seq(1,length(l$Ca))
                    }
            }

            .Object@call <- call

            callNextMethod(.Object, ...)   # to fill 'size'

            .Object@anion.x = c((6 * .Object@size/11)
                           + ((5 * .Object@size/11) * (.Object@Cl/100))
                           + (1/2) * (5 * .Object@size/11) * (.Object@SO4/100))
            .Object@anion.y = c(.Object@SO4 * 5 * .Object@size/1100)
            .Object@cation.x = cationx(.Object@Ca,.Object@Mg,.Object@size)
            .Object@cation.y = cationy(.Object@Mg,.Object@size)
            .Object@diamond.x = ( (1/2) * (.Object@cation.x + .Object@anion.x)
                            + ( (1/4) * (.Object@anion.y - .Object@cation.y) ) )
            .Object@diamond.y = ( .Object@anion.x - .Object@cation.x
                            + ( (1/2) * (.Object@anion.y + .Object@cation.y) ) )
          return(.Object)
          }
          )

#' Create a new piper object
#'
#' @param d list passed to class piper, \code{\linkS4class{piper}}
#' @param ... additional arguments, as for \code{\linkS4class{piper}}
#' @seealso \code{\link{piper-class}} and \code{\link{toPercent}}
#' @export
piper <- function(d, ...){
    new("piper", d, call=sys.call(), ...)
}

#' Create a new piperplot object
#'
#' @param size integer related to the size of the plot area
#' @param ... additional arguments, as for \code{\linkS4class{piperplot}}
#' @examples
#' library(hydrogeo)
#' p = piperPaper(size=1)
#' plot(p)
#' @export
piperPaper <- function(size=NULL, ...){
    new("piperplot", call=sys.call(), ...)
}

#' Plot the diagram area with two triangles and a diamond
#'
#' @param x object of class piperplot
#' @param axes logical saying whether to draw the axes or not, defaults to TRUE
#' @param ... further arguments to plot.default
#' @export plot
#' @importFrom graphics plot.default segments par
setMethod(
    f="plot",
    signature(x = "piperplot"),
        function(x, axes = TRUE, ...)
        {
            p <- (x@size/11)
            r <- (x@size/22)

            plot.default(0, 0, type="n", axes=FALSE, lty=1, lwd=1,
                         xlim=c(0, x@size + p), ylim=c(-p, x@size),
                         frame.plot=FALSE, ann=TRUE, ylab="", xlab="", ...)
            

    # draws grid lines
    thickxf <- c(0, (10 * r), (5 * r), (12 * r), x@size, (17 * r),
                 (x@size/2), (16 * r), (x@size/2), (6 * r))
    thickxt <- c((10 * r), (5 * r), 0, x@size, (17 * r), (12 * r),
                 (16 * r), (x@size/2), (6 * r), (x@size/2))
    thickyf <- c(0, 0, (10 * r), 0, 0, (10 * r), (2 * r), (12 * r), 
                 x@size, (12 * r))
    thickyt <- c(0, (10 * r), 0, 0, (10 * r), 0, (12 * r), x@size, 
                 (12 * r), (2 * r) )
    
    xf <- c(thickxf, (2 * r), (4 * r), (6 * r), (8 * r), (14 * r), 
            (16 * r), (18 * r), (20 * r), (21 * r), (20 * r), (19 * r),
            (18 * r), (21 * r), (20 * r), (19 * r), (18 * r), (9 * r), (8 * r),
            (7 * r), (6 * r), (9 * r), (8 * r), (7 * r), (6 * r), (7 * r),
            (8 * r), (9 * r), (10 * r), (7 * r), (8 * r), (9 * r), (10 * r))
    
    xt <- c(thickxt, r, (2 * r), (3 * r), (4 * r), (13 * r), (14 * r), (15 * r),
            (16 * r), (13 * r), (14 * r), (15 * r), (16 * r), (20 * r),
            (18 * r), (16 * r), (14 * r), (8 * r), (6 * r), (4 * r), (2 * r),
            r, (2 * r), (3 * r), (4 * r), (12 * r), (13 * r), (14 * r),
            (15 * r), (12 * r), (13 * r), (14 * r), (15 * r))
    
    yf <- c(thickyf, 0, 0, 0, 0, 0, 0, 0, 0, (2 * r), (4 * r), (6 * r), (8 * r),
            (2 * r), (4 * r), (6 * r), (8 * r), (2 * r), (4 * r), (6 * r),
            (8 * r), (2 * r), (4 * r), (6 * r), (8 * r), (14 * r), (16 * r),
            (18 * r), (20 * r), (10 * r), (8 * r), (6 * r), (4 * r))
    
    yt <- c(thickyt, (2 * r), (4 * r), (6 * r), (8 * r), (2 * r), (4 * r),
            (6 * r), (8 * r), (2 * r), (4 * r), (6 * r), (8 * r), 0, 0, 0, 0, 0,
            0, 0, 0, (2 * r), (4 * r), (6 * r), (8 * r), (4 * r), (6 * r),
            (8 * r), (10 * r), (20 * r), (18 * r), (16 * r), (14 * r))

    segments(xf, yf, xt, yt, col=par("fg"), lty=1, lwd=par("lwd"))

            if ( axes ) {
                Axis(x) }

##    x@plotted <- TRUE
}
    )


#' @describeIn piperplot Add axes to a piperplot
#' @param x an object of class piperplot
#' @importFrom graphics text
setMethod(
    f="Axis",
    signature(x="piperplot"),
    function(x = NULL) #, at = NULL, ..., side, labels = NULL)
        {
    p <- (x@size/11)      # for scaling
    r <- (x@size/22)      # for scaling

    cex.axis = 0.7 # FIXME
    # add axis titles
    vfont <-  c("serif", "italic")
    
    # label bottom axes (the two horizontal axes)
    xstr <- c(5 * r, 17 * r)
    ystr <- c(-r, -r)

    text(xstr, ystr, labels=c( expression(Ca^{2+''}), expression(paste(Cl,
                               ''^{symbol('-')})) ),
         vfont=vfont, cex=cex.axis)

    # label axes parallel with a line angled at 60% from horizontal
    # in a clockwise direction
    xgh <- c(14 * r, 8 * r, 20 * r)
    ygh <- c(18 * r, 6 * r, 6 * r)

    text(xgh, ygh,
         labels=c( expression( paste(Ca^{2+''}, " ", symbol('&'), " ",
                                       Mg^{2+''})),
                     expression( paste(Na^{symbol('+')}, " ", symbol('&'), " ",
                                       K^{symbol('+')})),
                     expression( paste(SO[4],''^{2-''}))),
         srt=300, vfont=vfont, cex=cex.axis)

    # label axes parallel with a line angled at 60% (approx.) from horizontal
    # in an anti clockwise direction    
    xla <- c(14 * r, 8 * r, 2 * r)
    yla <- c(6 * r, 18 * r, 6 * r)
    
    text(xla, yla,
         labels=c( expression( paste(CO[3],''^{2-''}, " ", symbol('&'), " ",
                                 HCO[3],''^{-''})),
                     expression( paste(SO[4],''^{2-''}, " ", symbol('&'), " ",
                                 Cl, ''^{symbol('-')})),
                     expression( paste(Mg^{2+''}))),
         srt=60, vfont=vfont, cex=cex.axis)

    labelAxes(x)
   })


setGeneric (
    "labelAxes",
    function(x, cex.axis=0.1, side=-1, ...)
    standardGeneric("labelAxes")
    )

#' @describeIn piper Label the axes
#' @param x an object of class piperplot
#' @param cex.axis magnification to be used for axis annotation relative
#' to the current setting of 'cex', see help("par")
#' @param side integer between 1 and 10 specifying which side to lable, the default is to label all
#' @importFrom graphics text
setMethod(
          f="labelAxes",
          signature(x="piperplot"),
          definition=function(x, cex.axis=0.35, side=-1, ...)
    {
    p <- (x@size/11)      # for scaling
    r <- (x@size/22)      # for scaling

    #TODO: Most of these don't work

    ## Tick mark labels:
    # add axis tick mark labels: 20, 40, 60, 80 percent labels
    adj <- 0.5     # to centre text on the point
    labels <- c(20 * (1:4))  # too clever by half: 20, 40, 60, 80 percent
    vfont <- c("serif", "italic")

    # for calculating offsets from line intersections
    ofs <- function(deg) {
      dd <- tan( (deg/360) * 2 * pi )
      offset <- cex.axis * dd * (x@size/50)
      return( offset )
    }

    if(side==-1 || side==1){
    # flat labels LHS - Mg2+
    xe <- c(r * 1:4)
    ye <- c(p * 1:4) 
    text(xe, ye, labels=labels, pos=2, offset=0.1, vfont=vfont, cex=cex.axis)
}
    if(side==-1 || side==2){
    # flat labels RHS - SO42-
    xe <- c(r * 21:18)
    ye <- c(p * 1:4) 
    text(xe, ye, labels=labels, pos=4, offset=0.1, vfont=vfont, cex=cex.axis)
}
    if(side==-1  || side==3){
    # labels rotated 60degrees
    # Ca2+ & Mg2+
    srt <- 60
    delta <- ofs(srt)
    xd <- c(r * 15:12) + delta
    yd <- c(p * 7:10) + delta
    text(xd, yd, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    if(side==-1 || side==4){
    # labels rotated 120degrees
    # Ca2+
    srt <- 120
    delta <- ofs(srt)
    xa <- c(p * 4:1) - delta
    ya <- c(0) + delta
    text(xa, ya, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}

    if(side==-1 || side==5){
    # CO32- & HCO3-, triangle
    srt <- 120
    delta <- ofs(srt)
    xa <- c(r * 16:13) + delta
    ya <- c(p * 4:1) - delta
    text(xa, ya, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    
    if(side==-1 || side==6){
    # CO32- & HCO3-, diamond
    srt <- 120
    delta <- ofs(srt)
    xa <- c(r * 15:12) - delta
    ya <- c(p * 5:2) + delta
    text(xa, ya, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    if(side==-1 || side==7){
    # labels rotated 240degrees
    # Cl-
    srt <- 240
    delta <- ofs(srt)
    xb <- c(p * 7:10) - delta
    yb <- c(0, 0, 0, 0) - delta
    text(xb, yb, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    if(side==-1 || side==8){
    # Na+ & K+, triangle
    srt <- 240
    delta <- ofs(srt)
    xb <- c(r * 6:9) + delta
    yb <- c(p * 4:1) + delta
    text(xb, yb, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    if(side==-1 || side==9){
    # Na+ & K+, diamond
    srt <- 240
    delta <- ofs(srt)
    xb <- c(r * 7:10) - delta
    yb <- c(p * 5:2) - delta
    text(xb, yb, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
    if(side==-1 || side==10){
    # labels rotated 300degrees
    # SO42- & Cl-
    srt <- 300
    delta <- ofs(srt)
    xc <- c(r * 7:10) + delta
    yc <- c(p * 7:10) - delta
    text(xc, yc, labels=labels, vfont=vfont, srt=srt, cex=cex.axis, adj=adj)
}
}
  )

#' @describeIn piper Plot an object of class \code{piper}
#' @param type what type of plot should be drawn, only "p" for *p*oints is useful
#' @param cex magnification to be used for symbols relative
#' to the current setting of 'cex', see help("par")
#' @importFrom graphics points.default
#' @examples
#' 
#' # change symbols and colours to differentiate water type groups
#' lp@pt.pch = c(2,2,4,4,4)
#' lp@pt.col = c(0,1,0,1,2)
#' 
#' # use larger symbols
#' plot( lp, main="Piper-Hill Diagram of Water Quality", cex=1.4 )
setMethod(
    f="plot",
    signature(x = "piper"),
    function(x, type = "p", cex=0.75, ...)
    {
        ##if( ! x@plotted ){ callNextMethod() }
        callNextMethod() # calls piperplot.plot()

        if ( type == "p" ) { # i.e. 'points'
            by(x, x@pt.col,
               function(j, ...) {
                   px <- c(j$anion.x, j$cation.x, j$diamond.x)
                   py <- c(j$anion.y, j$cation.y, j$diamond.y)
                   points.default(px, py, type="p", lty=1, lwd=1, pch=j$pt.pch,
                                  col=j$pt.col, bg=NA, cex=cex, ...)
               },
               x, ... )
        }
        invisible()
        return(TRUE) # for tests
    }
)

# Coercion functions
#' @export
as.data.frame.piper =
  function (x, ...)
  {
    as.data.frame.list( list( IDs=x@IDs, Ca=x@Ca, Mg=x@Mg, Cl=x@Cl,
                             SO4=x@SO4, pt.pch=x@pt.pch, pt.col=x@pt.col,
                             anion.x=x@anion.x, anion.y=x@anion.y,
                             cation.x=x@cation.x, cation.y=x@cation.y,
                             diamond.x=x@diamond.x, diamond.y=x@diamond.y ) )
  }


# TODO: work on this

# describeIn piper For checking the validity of an object of class \code{piper}
# param object an object of class piper
# name setValidity

#' @noRd
setValidity("piper",
            function(object)
            {
              n <- sort(names(object))
              f <- c("Ca", "Cl", "Mg", "SO4")
              if ( identical( intersect(n,f) , f ) ) TRUE
              else paste("ERROR: missing items: ", setdiff(f,n), sep="")

              a <- length(object[[1]])
              nms <- names(object)
              x <- 0
              for (i in object){
                x <- x + 1
                nme <- nms[x]
              
                if (length(i) != object@Ca) {
                  stop("ERROR: lengths of items differ") }
                if (nme != "IDs") {
                  for (j in i) {
                    if ( ! is.numeric(j) || j < 1)
                      stop("ERROR: there is a non-numeric or negative number")
                  }
                }
              }

              over100 <- object$row.name[ ((object$Ca + object$Mg) > 100 |
                                           (object$Cl + object$SO4) > 100) ]
              if (length(over100) != 0) {
                print(paste("ERROR: row.name == ",over100,
                       ", Either cations or anions add up to more than 100%.",
                       sep=""))
                return(invisible())      # represses echo
              }
            }
            )

#' @describeIn piper Show an object of class \code{piper}
#' @param object an object of class piper
#' @importFrom utils str
setMethod("show","piper",
          function(object){
              cat("Piper object: ")
              print(object@call)
              str(object)
          })

Try the hydrogeo package in your browser

Any scripts or data that you put into this service are public.

hydrogeo documentation built on May 2, 2019, 6:34 a.m.