R/panel.depth_function.R

In aqp: Algorithms for Quantitative Pedology

## TODO: large gaps in data result in some strange-looking polygons

## TODO: 'id' isn't a very safe argument for defining profile IDs... (e.g. sp1)

## TODO: iterate over profile IDs instead of groups
## TODO: confidence bands not defined when NA is present


#' @title Lattice Panel Function for Soil Profiles
#'
#' @description Panel function for plotting grouped soil property data, along with upper and
#' lower estimates of uncertainty.
#'
#' This function can be used to replace \code{panel.superpose} when plotting
#' depth function data. When requested, contributing fraction data are printed
#' using colors the same color as corresponding depth function lines unless a
#' single color value is given via \code{cf.col}.
#' 
#' This function is not able to apply transformations (typically `log = 10`) applied in the `scales` argument to `xyplot` to  upper/lower bounds. These will have to be manually applied. See examples.
#'
#' @aliases panel.depth_function prepanel.depth_function
#' 
#' @param x x values (generated by calling lattice function)
#' 
#' @param y y values (generated by calling lattice function)
#' 
#' @param id vector of id labels, same length as x and y--only required when
#' plotting segments (see Details section)
#' 
#' @param upper vector of upper confidence envelope values
#' 
#' @param lower vector of lower confidence envelope values
#' 
#' @param subscripts paneling indices (generated by calling lattice function)
#' 
#' @param groups grouping data (generated by calling lattice function)
#' 
#' @param sync.colors optionally sync the fill color within the region bounded
#' by (lower--upper) with the line colors
#' 
#' @param cf optionally annotate contributing fraction data at regular depth
#' intervals see \code{\link{slab}}
#' 
#' @param cf.col optional color for contributing fraction values, typically
#' used to override the line color
#' 
#' @param cf.interval number of depth units to space printed contributing
#' fraction values
#' 
#' @param \dots further arguments to lower-level lattice plotting functions,
#' see below
#' 
#' @author D.E. Beaudette
#' 
#' @seealso \code{\link{sp1}}, \code{\link{slice}}, \code{\link{slab}}
#' @export
#' @keywords hplot
#' @examples
#'
#' library(lattice)
#' data(sp1)
#' 
#' # 1. plotting mechanism for step-functions derived from soil profile data
#' xyplot(
#'   cbind(top, bottom) ~ prop,
#'   data = sp1,
#'   id = sp1$id,
#'   panel = panel.depth_function,
#'   ylim = c(250, -10),
#'   scales = list(y = list(tick.number = 10)),
#'   xlab = 'Property',
#'   ylab = 'Depth (cm)',
#'   main = 'panel.depth_function() demo'
#' )
#' 
#' # 1.1 include groups argument to leverage lattice styling framework
#' sp1$group <- factor(sp1$group, labels = c('Group 1', 'Group2'))
#' 
#' xyplot(
#'   cbind(top, bottom) ~ prop,
#'   groups = group,
#'   data = sp1,
#'   id = sp1$id,
#'   panel = panel.depth_function,
#'   ylim = c(250, -10),
#'   scales = list(y = list(tick.number = 10)),
#'   xlab = 'Property',
#'   ylab = 'Depth (cm)',
#'   main = 'panel.depth_function() demo',
#'   auto.key = list(
#'     columns = 2,
#'     points = FALSE,
#'     lines = TRUE
#'   ),
#'   par.settings = list(superpose.line = list(col = c(
#'     'Orange', 'RoyalBlue'
#'   )))
#' )
#' 
#' 
#' # more complex examples, using step functions with grouped data
#' # better looking figures with less customization via tactile package
#' if(requireNamespace('tactile')) {
#'   
#'   library(data.table)
#'   library(lattice)
#'   library(tactile)
#'   
#'   # example data
#'   data(sp6)
#'   
#'   # a single profile
#'   x <- sp6[1:5, ]
#'   
#'   # wide -> long format
#'   x.long <- data.table::melt(
#'     data.table::data.table(x), 
#'     id.vars = c('id', 'top', 'bottom'), 
#'     measure.vars = c('sand', 'silt', 'clay')
#'   )
#'   
#'   # (optional) convert back to data.frame
#'   x.long <- as.data.frame(x.long)
#'   
#'   # three variables sharing a common axis
#'   # factor levels set by melt()
#'   xyplot(
#'     cbind(top, bottom) ~ value | id,
#'     groups = variable,
#'     data = x.long,
#'     id = x.long$id,
#'     ylim = c(200, -5), xlim = c(10, 60),
#'     scales = list(alternating = 1, y = list(tick.number = 10)),
#'     par.settings = tactile.theme(superpose.line = list(lwd = 2)),
#'     xlab = 'Sand, Silt, Clay (%)',
#'     ylab = 'Depth (cm)',
#'     panel = panel.depth_function,
#'     auto.key = list(columns = 3, lines = TRUE, points = FALSE),
#'     asp = 1.5
#'   )
#'   
#'   
#'   # all profiles
#'   x <- sp6
#'   
#'   # wide -> long format
#'   x.long <- data.table::melt(
#'     data.table::data.table(x), 
#'     id.vars = c('id', 'top', 'bottom'), 
#'     measure.vars = c('sand', 'silt', 'clay')
#'   )
#'   
#'   # (optional) convert back to data.frame
#'   x.long <- as.data.frame(x.long)
#'   
#'   # three variables sharing a common axis
#'   # factor levels set by melt()
#'   xyplot(
#'     cbind(top, bottom) ~ value | id,
#'     groups = variable,
#'     data = x.long,
#'     id = x.long$id,
#'     ylim = c(200, -5), xlim = c(0, 70),
#'     scales = list(alternating = 1, y = list(tick.number = 10)),
#'     par.settings = tactile.theme(superpose.line = list(lwd = 2)),
#'     xlab = 'Sand, Silt, Clay (%)',
#'     ylab = 'Depth (cm)',
#'     panel = panel.depth_function,
#'     auto.key = list(columns = 3, lines = TRUE, points = FALSE),
#'     as.table = TRUE
#'   )
#'   
#'   xyplot(
#'     cbind(top, bottom) ~ value,
#'     groups = variable,
#'     data = x.long,
#'     id = x.long$id,
#'     ylim = c(200, -5), xlim = c(0, 70),
#'     scales = list(alternating = 1, y = list(tick.number = 10)),
#'     par.settings = tactile.theme(superpose.line = list(lwd = 2)),
#'     xlab = 'Sand, Silt, Clay (%)',
#'     ylab = 'Depth (cm)',
#'     panel = panel.depth_function,
#'     auto.key = list(columns = 3, lines = TRUE, points = FALSE),
#'     as.table = TRUE
#'   )
#'   
#'   xyplot(
#'     cbind(top, bottom) ~ value | variable,
#'     groups = variable,
#'     data = x.long,
#'     id = x.long$id,
#'     ylim = c(200, -5), xlim = c(0, 70),
#'     scales = list(alternating = 1, y = list(tick.number = 10)),
#'     par.settings = tactile.theme(superpose.line = list(lwd = 2)),
#'     xlab = 'Sand, Silt, Clay (%)',
#'     ylab = 'Depth (cm)',
#'     panel = panel.depth_function,
#'     auto.key = list(columns = 3, lines = TRUE, points = FALSE),
#'     as.table = TRUE
#'   )
#'   
#'   
#'   xyplot(
#'     cbind(top, bottom) ~ value | variable,
#'     data = x.long,
#'     id = x.long$id,
#'     ylim = c(200, -5), xlim = c(0, 70),
#'     scales = list(alternating = 1, y = list(tick.number = 10)),
#'     par.settings = tactile.theme(superpose.line = list(lwd = 2)),
#'     xlab = 'Sand, Silt, Clay (%)',
#'     ylab = 'Depth (cm)',
#'     panel = panel.depth_function,
#'     auto.key = list(columns = 3, lines = TRUE, points = FALSE),
#'     as.table = TRUE
#'   )
#'   
#'   
#' }
panel.depth_function <- function(x, y, id, upper=NA, lower=NA, subscripts=NULL, groups=NULL, sync.colors=FALSE, cf=NA, cf.col=NA, cf.interval=20, ...) {
  
  # add grid
  panel.grid(h=-1, v=-1, lty=3, col=1)
  
  # load current line style
  superpose.line <- trellis.par.get("superpose.line")
  
  # when the length of 'y' is > 'x', we are plotting a step function
  if(length(y) > length(x)) {
    
    # sanity check
    if(missing(id)) {
      stop('must provide a profile id')
    }
      
    # re-make a nice data.frame, assuming that we have 'groups' defined
    if(!missing(groups)) {
      d <- data.frame(
        prop = x, 
        bnd = y, 
        upper = upper[subscripts], 
        lower = lower[subscripts], 
        groups = groups[subscripts], 
        id = id[subscripts]
      )
    } else {
      # 'groups' is missing, add a fake 'groups' column
      d <- data.frame(
        prop = x, 
        bnd = y, 
        upper = upper[subscripts], 
        lower = lower[subscripts], 
        groups = factor(1), 
        id = id[subscripts]
      )
    }
    
    
    # iterate over IDs creating step functions, by group
    by(d, d$id, .make.segments, ...)
  }
  
  # normal plot -- not a step function
  else {
    # if we have an upper and lower bound defined, plot them
    if(!missing(upper) & !missing(lower)) {
      # working with grouped data and paneled data
      if(!missing(groups) & !missing(subscripts)) {
        d <- data.frame(
          yhat = x, 
          top = y, 
          upper = upper[subscripts], 
          lower = lower[subscripts], 
          groups = groups[subscripts]
        )
        
        # levels in the groups, for color matching
        ll <- levels(d$groups)
        n_groups <- length(ll)
      }
      
      # no grouping, add a fake group for compatibility
      if(missing(groups)) {
        fake.groups <- factor(1)
        d <- data.frame(
          yhat = x, 
          top = y, 
          upper = upper[subscripts], 
          lower = lower[subscripts], 
          groups = fake.groups
        )
        
        # levels in the groups, for color matching
        ll <- levels(d$groups)
        n_groups <- length(ll)
      }
      
      # optionally sync region + main line colors
      if(sync.colors) {
        region.col <- rep(superpose.line$col, length.out = n_groups)
      } else {
        region.col <- rep(grey(0.7), length.out = n_groups)
      }
        
      
      # add conf. intervals / aggregation uncertainty
      by(d, d$groups, function(d_i) {
        
        # lookup color
        m <- match(as.character(d_i$group[1]), ll)
        
        # cannot have NA in values that define polygon boundaries
        d_i <- d_i[which(!is.na(d_i$upper) & !is.na(d_i$lower)), ]
        
        # make conf.int polygon
        x.coords <- c(d_i$lower, rev(d_i$upper))
        y.coords <- c(d_i$top, rev(d_i$top))
        poly.col <- region.col[m]
        
        # note: using tactile interfers with alpha
        panel.polygon(x = x.coords, y = y.coords, col = poly.col, border = NA, ...)
      })
    }
    
    # no upper/lower polygon boundaries defined
    else {
      if(missing(groups)) {
        fake.groups <- factor(1)
        d <- data.frame(yhat=x, top=y, groups=fake.groups)
      } else {
        d <- data.frame(yhat=x, top=y, groups=groups[subscripts])
      }
      
      # levels in the groups, for color matching
      ll <- levels(d$groups)
      n_groups <- length(ll)
    }
    
    # setup style parameters for main lines
    # repeat enough times for the current number of groups
    line.col <- rep(superpose.line$col, length.out=n_groups)
    line.lty <- rep(superpose.line$lty, length.out=n_groups)
    line.lwd <- rep(superpose.line$lwd, length.out=n_groups)
    
    # add main lines
    by(d, d$groups, function(d_i){
      # lookup color
      m <- match(unique(d_i$group), ll)
      
      # add line
      panel.lines(d_i$yhat, d_i$top, lwd=line.lwd[m], col=line.col[m], lty=line.lty[m])
    })
  }
  
  
  # annotate with contributing fraction
  if(!missing(cf)) {
    
    # add CF data to panel's worth of data
    d$cf <- cf[subscripts]
    
    # annotate with contributing fraction by group
    by(d, d$groups, function(d_i){
      
      # lookup linestyle by group
      m <- match(unique(d_i$group), ll)
      
      # if there is no user-specified CF color, then use the same as line style
      if(is.na(cf.col)) {
        cf.col <- line.col[m]
      }
      
      # make a function for linear interpolation of CF values based on depth
      cf.approx.fun <- approxfun(d_i$top, d_i$cf, method='linear')
      
      # generate annotated depths: 5 cm to 95th percentile of max depth
      y.q95 <- quantile(d_i$top, probs=c(0.95), na.rm=TRUE)
      a.seq <- seq(from=2, to=y.q95, by=cf.interval)
      
      # offset CF sequence according to group index
      a.seq <- a.seq + ((m-1) * cf.interval/4)
      
      # interpolate CF at annotated depths
      a.CF <- cf.approx.fun(a.seq)
      a.text <- paste(round(a.CF * 100), '%')
      
      # remove any NAs in CF sequence and text
      not.na.idx <- which(! is.na(a.CF))
      a.seq <- a.seq[not.na.idx]
      a.text <- a.text[not.na.idx]
      
      # add to right-hand side of the panel
      unit <- gpar <- NULL
      grid::grid.text(
        a.text,
        x = grid::unit(0.99, 'npc'),
        y = grid::unit(a.seq, 'native'),
        just = 'right',
        gp = grid::gpar(font = 3, cex = 0.8, col = cf.col)
      )
    })
    
  }
  
}

# internally used function, for panel.depth_function() step-function mode
.make.segments <- function(df) {
  
  # group data
  ll <- levels(df$groups)	
  n_groups = length(ll)
  
  # get the current group index
  m <- match(unique(df$group), ll)
  
  # load current line style
  superpose.line <- trellis.par.get("superpose.line")
  
  # repeat enough times for the current number of groups
  line.col <- rep(superpose.line$col, length.out = n_groups)
  line.lty <- rep(superpose.line$lty, length.out = n_groups)
  line.lwd <- rep(superpose.line$lwd, length.out = n_groups)
  
  for(i in m) {
    
    this.group <- df[df$groups == ll[i], ]
    
    # get the number of depth slices
    n_hz <- length(this.group$prop) / 2
    
    # need at least 2 horizons
    if(n_hz > 1) {
      
      .top <- this.group$bnd[1:n_hz]
      .bottom <- this.group$bnd[(n_hz+1):length(this.group$prop)]
      .prop <- this.group$prop[1:n_hz]
      
      # vertical segments
      panel.segments(x0 = .prop, y0 = .top, x1 = .prop, y1 = .bottom, 
                     lwd=line.lwd[i], col=line.col[i], lty=line.lty[i]) 
      
      # horizontal segments
      panel.segments(x0 = .prop[-n_hz], y0 = .bottom[-n_hz], x1 = .prop[-1], y1 = .top[-1], 
                     lwd=line.lwd[i], col=line.col[i], lty=line.lty[i])
    }
    
    else {
      message(paste('only 1 horizon, skipping!', df$groups[1]))
    }
    
  }
  
}