R/get.increase.matrix.R In aqp: Algorithms for Quantitative Pedology

Documented in get.increase.depthsget.increase.matrix

```## TODO: clarify, new title / description

#' Compute Pair-wise Distances of Soil Properties over Depth
#'
#' Computes pair-wise distance matrix to determine where an attribute increases within
#' a specified vertical distance threshold.
#'
#' Uses matrix outer product to determine all pair-wise differences in `attr`
#' for the horizons of `p`. Supplies `attr` to `threshold.fun` to determine the
#' minimum value criterion to return TRUE in output matrix for an "increase".
#' Also, computes all pair-wise distances in depth dimension to determine
#' whether the vertical distance criteria have been met simultaneously with
#' `attr` increase.
#'
#' This function assumes that the `threshold.fun` supplied by the user returns
#' either a constant or a vector of equal length to its input.
#'
#' Note that the `threshold.fun` result is allowed to contain NA, but that will
#' result in no output for affected cells.
#'
#' `get.increase.depths()` performs the conversion of the square matrix
#' output of `get.increase.matrix` back to horizon top depth for where
#' criteria were met.
#'
#'
#' @param p a `SoilProfileCollection`, containing a single profile
#' @param attr horizon attribute name to get the "increase" of
#' @param threshold.fun a function that returns the threshold (as a function of
#' attr); may return a constant single value
#' @param vertical.distance the vertical distance (determined from difference
#' SPC top depth variable) within which increase must be met
#' @return Returns a square logical matrix reflecting where the increase
#' criteria were met.
#'
#' `get.increase.depths` converts to horizon top depth by using above
#' matrix output to determine depths where increase is met.
#'
#' @author Andrew Gene Brown
#' @seealso [getArgillicBounds()], [crit.clay.argillic()]
#' @keywords manip
#' @export
#' @examples
#'
#' data(sp1, package = 'aqp')
#' depths(sp1) <- id ~ top + bottom
#' site(sp1) <- ~ group
#'
#' p <- sp1[1]
#' attr <- 'prop' # clay contents
#' foo <- get.increase.matrix(p, threshold.fun = crit.clay.argillic,
#'                            attr = attr, vertical.distance = 30)
#' foo
#'
get.increase.matrix <- function(p, attr, threshold.fun, vertical.distance) {
# p - a SoilProfileCollection, containing a single profile
# attr - attribute name to get the "increase" of
# threshold.fun - a function that returns the threshold (as a function of attr); may return a constant single value
# vertical distance - the vertical distance (determined from difference SPC top depth variable) within which increase must be met

# NOTE: should topdepths or midpoints be used? typically the same horizon index is chosen
#
#       with thick horizons/diffuse transitions -- the midpoint is probably more appropriate.
#
#       I originally chose top depth because it is more straightforward to validate against known data
#       whereas the midpoints required some more thought. now that I am pleased with stability of algo,
#       I think it is worth switching to midpoint given it probably gives a better estimate of
#       the thickness of transitional zone between horizons when it matters
#
h <- horizons(p)
depthz <- horizonDepths(p)
middepth <- h[[depthz[1]]] + (h[[depthz[2]]] - h[[depthz[1]]]) / 2 # TODO: evaluate

increase.var <- horizons(p)[[attr]]

threshold.vector <- threshold.fun(increase.var)

if(length(threshold.vector) == 1) {
# if threshold.fun() returns a constant value, expand it to match the length of the attribute
threshold.vector <- rep(threshold.vector, length(increase.var))

} else if(length(threshold.vector) != length(increase.var)) {
# this function assumes that the threshold.fun() supplied by the user returns either a constant,
# or a vector of equal length to `increase.var` when called above. otherwise, we cannot calculate the result.
# note that the threshold.fun() result _is_ allowed to contain NA, but that will result in no output for affected cells
stop(paste0('profile ID:', profile_id(site(p)),
" - threshold.fun() result should be length 1 or equal to length of attribute \'",attr,"\' (n=",length(increase.var),")."))
}

# repeat attr content by horizon in columns; number of columns = number of horizons
attr.mat <- outer(increase.var, rep(1, length(increase.var)))

# create a matrix of corresponding thrsholds
thresh.mat <- outer(threshold.vector, rep(1, length(increase.var)))

# calculate a matrix of attr differences (between all horizons, not just adjacent)
attr.inc.mat <- outer(increase.var, increase.var, `-`)

# calculate a vertical distance matrix (between all horizons)
vdist.mat <- outer(middepth, c(0, middepth[-length(middepth)]), `-`)

# crit1 "an increase of at least [thresh.mat]"
increase.met <- (attr.mat - thresh.mat) > (attr.inc.mat * upper.tri(attr.inc.mat))
increase.met[is.na(increase.met)] <- FALSE

# crit2 "within a vertical distance of [vertical.distance]"
vdist.met <- abs(vdist.mat) <= vertical.distance

# are crit 1 and crit 2 met?
criteria.met <- increase.met & vdist.met
return(criteria.met)
}

#' Return the horizon top depths from a call to `get.increase.matrix()`
#'
#' \code{get.increase.depths} performs the conversion of the square matrix
#' output of \code{get.increase.matrix} back to horizon top depth for where
#' criteria were met.
#'
#' Note that the `threshold.fun` result is allowed to contain `NA`, but that will
#' result in no output for affected cells.
#'
#'
#' @param p a SoilProfileCollection, containing a single profile
#' @param attr horizon attribute name to get the "increase" of
#' @param threshold.fun a function that returns the threshold (as a function of
#' attr); may return a constant single value
#' @param vertical.distance the vertical distance (determined from difference
#' SPC top depth variable) within which increase must be met
#' @return Returns a numeric vector of depths where the increase requirement is
#' met. For the argillic, the first is the one of interest.
#'
#' `get.increase.depths()` converts to horizon top depth by using above
#' matrix output to determine depths where increase is met.
#' @seealso [getArgillicBounds()] [crit.clay.argillic()]
#' @rdname get.increase.matrix
#' @keywords manip
#' @export
#' @examples
#'
#' data(sp1, package = 'aqp')
#' depths(sp1) <- id ~ top + bottom
#' site(sp1) <- ~ group
#'
#' p <- sp1[1]
#' attr <- 'prop' # clay contents
#' foo <- get.increase.depths(p, threshold.fun = crit.clay.argillic,
#'                            attr = attr, vertical.distance = 30)
#' foo
#'
get.increase.depths <- function(p, attr, threshold.fun, vertical.distance) {
topdepth <- horizons(p)[[horizonDepths(p)[1]]]
#TODO: is there interesting pedogenic information that can be derived
# from the lower triangle of the criteria.met matrix we omit from the
# "increase" matrix?) -- i.e. pale/haplo clay decrease at depth?
criteria.met <- get.increase.matrix(p, attr, threshold.fun, vertical.distance)
# get the index of  _first_ column from matrix criteria.met that has nonzero sum,
# this column index is the index of horizon where the attr increase is met
# then return the top depth.
return(topdepth[which(colSums(criteria.met) > 0)])
}
```

Try the aqp package in your browser

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

aqp documentation built on Sept. 11, 2024, 7:11 p.m.