R/lith.R

In abarbour/poel: Tools to facilitate working with simulations from poel

#' @title Poel's lithology structure
#' @description Manipulate or generate poel-type lithology specifications
#' @export
#' @aliases lith
#' @examples
#'
#' # just do something random to demonstrate
#' set.seed(1234)
#'
#' n <- 5 # five layers
#' deps <- seq_len(n)-1 # at these depths
#' rand <- runif(n)
#'
#' (l <- lithology(deps, rand, rand, rand))
#' is.lith(l)
#'
#' # plot values
#' plot(l)
#'
#' # Functions for plotting lith params at any depth:
#' funs <- depth_param_functions(l)
#'
#' HD <- funs[['Diffusivity']]
#' # get many more values than there were before
#' new.depths <- sort(jitter(rep(deps, each=n)))
#' new.hd <- HD(new.depths)
#' plot(new.depths, new.hd, type='l')
#'
lithology <- function(Depth,
                      ShearModulus=1.e9, Skempton=0.85, Diffusivity=0.1,
                      Nu=0.25, Nuu=0.35, Number=NULL){
  if (is.null(Number)) Number <- seq_along(Depth)
  Depth <- as.numeric(Depth)
  stopifnot(all(Depth>=0))
  ord <- order(Depth)
  L <- cbind(Number, Depth, ShearModulus, Nu, Nuu, Skempton, Diffusivity)[ord, , drop=FALSE]
  Lp <- ._lith_()
  colnames(L) <- Lp$Params
  class(L) <- 'lith'
  attr(L, "Units") <- Lp$Units
  attr(L, 'Nlith') <- length(Depth)
  return(L)
}

._lith_ <- function(){
    data.frame(Params=c('Number','Depth','ShearModulus','Nu','Nuu','Skempton','Diffusivity'),
               Units=c("","m","Pa","","","","m^2/s"))
}

#' @rdname lithology
#' @export
is.lith <- function(x, ...){
    mat <- is.matrix(x)
    lithn <- colnames(x)
    np <- ncol(x)
    nl <- nrow(x)
    nms <- np == 7
    unts <- length(attr(x, "Units")) == 7
    lays <- attr(x, 'Nlith') == nl
    tests <- c(Matrix=mat, N.Params=nms, N.Units=unts, N.Layers=lays)
    status <- all(tests)
    attr(status, 'tested') <- tests
    return(status)
}

#' @rdname lithology
#' @export
plot.lith <- function(L, no.layout=TRUE, add=FALSE, plot.inds=c(3,4,5,6,7), ...){
  Units <- attr(L,"Units")
  Deps <- L[,"Depth"]
  ldim <- dim(L)
  nms <- colnames(L)
  #
  nlay <- ldim[1]
  if (nlay==1){
    stop("Cannot plot this lithology -- only a single layer")
  }
  #
  layers <- seq_len(ldim[1])
  Deps.t <- Deps[layers[-1]]
  drng <- range(c(0,pretty(1.05*Deps[Deps>0])))
  #
  PLT <- function(n, add.=FALSE, ...){
    D <- L[,n]
    Du <- Units[n]
    Dnm <- nms[n]
    sf <- stepfun(Deps.t, D)
    if (!add.) plot(sf, lty=2, pch=NA, xaxs="i", xlab="Depth", main="", ylab=parse(text=Du), ...)
    lines(sf, vertical=FALSE, lwd=3, pch=NA, xaxs="i", ...)
    lims <- par("usr")
    text(lims[2]*0.98, lims[4]*0.95, Dnm, pos=2, font=2)
    return(sf)
  }
  toPlot <- plot.inds
  if (!no.layout){
    layout(matrix(seq_along(toPlot),ncol=1))
  }
  op <- par(mar=c(3.8,4.1,1,1), mgp=c(2.5,1,0), las=1, lend="butt")
  on.exit(par(op))
  #
  invisible(lapply(toPlot, PLT, xlim=drng, add.=add, ...))
}

#' @rdname lithology
#' @export
depth_param_functions <- function(x, ...) UseMethod('depth_param_functions')

#' @rdname lithology
#' @export
depth_param_functions.lith <- function(x, ...){
    stopifnot(is.lith(x))
    all.params <- colnames(x)
    params <- all.params[!(all.params %in% c('Number','Depth'))]
    .depths. <- x[,'Depth']
    .GetDepFun <- function(param.name){
        X <- .depths.
        Y <- x[, param.name]
        stats::approxfun(X, Y, method='constant', rule=2, f=0, ties = function(z){z[length(z)]})
    }
    names(params) <- params
    plyr::llply(params, .GetDepFun)
}