1 Examples with the Loafercreek sample data

Load required packages, sample data, and subset of horizon attributes.

# SoilProfileCollection methods
# you will need the latest version of aqp
# https://github.com/ncss-tech/aqp/
library(aqp)
# data-getting / sample data
library(soilDB)
# soil texture triangle
library(soiltexture)
# convenient alpha() function
library(scales)

# load sample data from soilDB package
data("loafercreek")

# extract just horizon data
h <- horizons(loafercreek)

Subset sand, silt, clay percentages and filter missing records or inconsistent fractions (e.g. add to more or less than 100%). Column names are adjusted to accommodate subsequent usage.

# extract components of field texture, removing rows with missing data
ssc <- h[, c('genhz', 'sand', 'silt', 'clay')]
ssc <- na.omit(ssc)

# adjust names for plotting with TT.plot()
# names must be SAND, SILT, CLAY
names(ssc) <- toupper(names(ssc))

# test of bogus data
ssc$sum <- rowSums(ssc[, c('SAND', 'SILT', 'CLAY')])
# > 5% deviation from 100%
idx <- which(abs(ssc$sum - 100) > 5)

# check errors
ssc[idx, ]
##     GENHZ SAND SILT CLAY sum
## 623     A   45   62   17 124
# ommit errors
# maybe fix these later?
ssc <- ssc[-idx, ]

# check: OK
head(ssc)
##    GENHZ SAND SILT CLAY sum
## 23     A   40   48   12 100
## 24   Bt1   40   48   12 100
## 25   Bt1   45   41   18 104
## 26   Bt2   50   32   18 100
## 27   BCt   50   34   16 100
## 29     A   15   66   19 100
# split by generalized horizon name into a list of data.frames
ssc.list <- split(ssc, ssc$GENHZ)

record.count <- rev(sapply(ssc.list, nrow))

par(mar = c(4.5, 3, 3, 1))
dotchart(record.count, pch = 21, bg = 'royalblue', pt.cex = 1.5, xlab = 'Number of records', main = 'Loafercreek\nRecords per Generalized Horizon Label')