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')