Introduction

This document outlines the basic usage of the textureTriangleSummary() function from the aqp package (version >= 1.10-4). This function displays approximate "ranges" (based on user-defined percentiles) in sand, silt, and clay percentages on the soil texture triangle in the form of a shaded polygon. The outline of the polygon is created by connecting percentiles (typically 10th and 90th) of sand, silt, and clay fractions. The "Sample RV" is based on the univariate 50th percentiles (e.g. median sand, silt, and clay fractions). A reference, "normal composition" range and RV can be generated, based on the mean and (compositional) covariance of the source data. These calculations require the compositions package.

Prepare Example Data

library(aqp)
library(soilDB)
library(compositions)
library(reshape2)
library(latticeExtra)
library(plyr)

# get sand, silt, clay data from the named series
# return weighted mean of top 10cm
get.data <- function(series) {
  x <- fetchKSSL(series)
  s <- slab(x, pedon_key ~ clay + sand + silt, slab.structure=c(0, 10), strict = FALSE, slab.fun = mean, na.rm=TRUE)
  s <- dcast(s, pedon_key ~ variable, value.var = 'value')
  ssc <- na.omit(s[, c('sand', 'silt', 'clay')])
  return(ssc)
}

# combine original texture + simulated textures from normal composition
sim.data <- function(ssc, n.sim=1000) {
  # convert to compositional class, note range is now [0,1]
  ssc.acomp <- acomp(ssc)
  # simulate normally-distributed composition based on data
  # note that var() is dispatched to var.acomp()
  ssc.sim <- rnorm.acomp(n=n.sim, mean=meanCol(ssc.acomp), var=var(ssc.acomp))
  # convert to data.frame and original range of [0,100]
  ssc.sim <- as.data.frame(unclass(ssc.sim) * 100)
  # return stacked original + simulated data
  return(make.groups(original=ssc, simulated=ssc.sim))
}


# get data for named series
s.list <- c('clarksville', 'vista', 'auburn', 'cecil', 'drummer', 'capay')
d <- lapply(s.list, get.data)
names(d) <- s.list
d <- ldply(d)

# simulate from a normal composition, using mean and var-cov matrix from original data
d.sim <- ddply(d, '.id', function(i) sim.data(i[, -1]))
m <- melt(d.sim, id.vars=c('.id', 'which'))

Compare Univariate Distributions

# plotting style
cols <- brewer.pal(n=3, 'Set1')
tps <- list(superpose.line=list(col=cols), superpose.symbol=list(col=cols, cex=0.5, pch='|'))

# plot: panel according to original | simulated 
p.1 <- densityplot(~ value | which + .id, groups=variable, data=m, xlab='Sand, Silt, Clay (%)', auto.key=list(columns=3), par.settings=tps, scales=list(alternating=3), panel=function(...) {
  panel.grid(-1, -1)
  panel.densityplot(...)
})
p.1 <- useOuterStrips(p.1)

# plot: panel according to texture component
p.2 <- densityplot(~ value | variable + .id, groups=which, data=m, xlab='Sand, Silt, Clay (%)', auto.key=list(columns=2), par.settings=tps, scales=list(alternating=3), panel=function(...) {
  panel.grid(-1, -1)
  panel.densityplot(...)
})
p.2 <- useOuterStrips(p.2)

# display plots
print(p.1)

print(p.2)

Display on Soil Texture Triangle

See the manual page for details on the various arguments to textureTriangleSummary().

textureTriangleSummary(d[d$.id == 'clarksville', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Clarksville\n0-10cm Weighted Mean')

textureTriangleSummary(d[d$.id == 'drummer', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Drummer\n0-10cm Weighted Mean')

textureTriangleSummary(d[d$.id == 'auburn', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Auburn\n0-10cm Weighted Mean')

textureTriangleSummary(d[d$.id == 'cecil', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Cecil\n0-10cm Weighted Mean')

textureTriangleSummary(d[d$.id == 'vista', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Vista\n0-10cm Weighted Mean')

textureTriangleSummary(d[d$.id == 'capay', -1], sim=TRUE, cex=0.5, p = c(0.10, 0.5, 0.9)); title('Capay\n0-10cm Weighted Mean')

This document is based on aqp version 1.10-3.