Soil Profile Summaries
D.E. Beaudette
2021-02-01
Introduction
# load required libaries
library(aqp)
library(soilDB)
library(sharpshootR)
library(cluster)
library(lattice)
library(reshape2)
library(colorspace)
# init vector of taxonnames to summarize
soils <- c('Whiterock', 'Copperopolis', 'Amador', 'Auburn', 'Loafercreek', 'Argonaut', 'Crimeahouse')Optionally re-make cached data. Requires NASIS and relevant soils in the selected set.
# get all pedons from the selected set
x <- fetchNASIS(from='pedons', rmHzErrors = FALSE, nullFragsAreZero = FALSE)
# convert vector of taxonnames into REGEX pattern for matching
pat <- paste0(soils, collapse = '|')
# subset pedons that match our REGEX pattern
idx <- grep(pat, x$taxonname, ignore.case = TRUE)
x <- x[idx, ]
# save for later
save(x, file='profile-summary-SPC.rda')Get the cached copy of the SoilProfileCollection object used in this tutorial.
# load cached copy
# load(url('http://ncss-tech.github.io/AQP/aqp/profile-summary-SPC.rda'))
load(file='profile-summary-SPC.rda')
# normalize taxonname via REGEX
for(i in soils) {
x$taxonname[grep(i, x$taxonname, ignore.case = TRUE)] <- i
}
# tabulate number of profiles
table(x$taxonname)##
## Amador Argonaut Auburn Copperopolis Crimeahouse Loafercreek Whiterock
## 15 31 63 27 29 109 34# convert soil colors to CIE LAB colorspace
# moist colors
cols.lab <- munsell2rgb(x$m_hue, x$m_value, x$m_chroma, returnLAB = TRUE)
x$m_L <- cols.lab$L
x$m_A <- cols.lab$A
x$m_B <- cols.lab$B
# dry colors
cols.lab <- munsell2rgb(x$d_hue, x$d_value, x$d_chroma, returnLAB = TRUE)
x$d_L <- cols.lab$L
x$d_A <- cols.lab$A
x$d_B <- cols.lab$B
# aggregate data by normalized taxonname, via slice-wise median
a.colors <- slab(x, taxonname ~ d_r + d_g + d_b + d_L + d_A + d_B + clay + phfield + total_frags_pct)
# throw out aggregate data that are deeper than 150cm
a.colors <- subset(a.colors, subset=bottom < 150)
# convert long -> wide format
x.colors <- dcast(a.colors, taxonname + top + bottom ~ variable, value.var = 'p.q50')
# check
head(a.colors)## variable taxonname p.q5 p.q25 p.q50 p.q75 p.q95 top bottom
## 1 d_r Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251 0 1
## 2 d_r Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251 1 2
## 3 d_r Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251 2 3
## 4 d_r Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251 3 4
## 5 d_r Amador 0.5087045 0.5695402 0.6617205 0.6617205 0.7136251 4 5
## 6 d_r Amador 0.4205719 0.5589675 0.6617205 0.6617205 0.7136251 5 6
## contributing_fraction
## 1 0.7333333
## 2 0.7333333
## 3 0.7333333
## 4 0.7333333
## 5 0.7333333
## 6 0.7333333head(x.colors)## taxonname top bottom d_r d_g d_b d_L d_A d_B clay phfield
## 1 Amador 0 1 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## 2 Amador 1 2 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## 3 Amador 2 3 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## 4 Amador 3 4 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## 5 Amador 4 5 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## 6 Amador 5 6 0.6617205 0.5656823 0.4098928 61.72703 4.870262 19.10015 17 5.6
## total_frags_pct
## 1 7
## 2 7
## 3 7
## 4 7
## 5 7
## 6 7# composite sRGB triplets into an R-compatible color
# note that missing colors must be padded with NA
x.colors$soil_color <- NA
not.na <- which(complete.cases(x.colors[, c('d_L', 'd_A', 'd_B')]))
cols.srgb <- data.frame(coords(as(LAB(x.colors$d_L, x.colors$d_A, x.colors$d_B), 'sRGB')))
x.colors$soil_color[not.na] <- with(cols.srgb[not.na, ], rgb(R, G, B, maxColorValue = 1))
# aggregate bedrock depth probability by taxonname
# at 90% level of confidence
dp <- aggregateSoilDepth(x, groups = 'taxonname', crit.prob=0.9)
# init a new SoilProfileCollection from aggregate data
depths(x.colors) <- taxonname ~ top + bottom
# join-in our depth to contact data
site(x.colors) <- dp
# generate index for new ordering
new.order <- match(c('Whiterock', 'Copperopolis', 'Amador', 'Auburn', 'Loafercreek', 'Argonaut', 'Crimeahouse'), profile_id(x.colors))
# make a table of data for adding depth brackets to profile sketches
# see ?addBracket for details
depth.brackets <- site(x.colors)
names(depth.brackets)[2:3] <- c('top', 'bottom')
# add a label
depth.brackets$label <- 'P(soil >= 90%)'par(mar=c(1, 0, 3, 1.5))
plotSPC(
x.colors,
divide.hz = FALSE,
name = NA,
plot.order = new.order,
col.label = 'Soil Color',
lwd = 1.25,
axis.line.offset = -3,
cex.depth.axis = 1,
cex.id = 1
)
addBracket(depth.brackets, col='black', label.cex=0.85)
title('Aggregate Soil Properties')
d <- profile_compare(x.colors, vars = c('d_L', 'd_A', 'd_B'), max_d = 150, k = 0)
plotProfileDendrogram(
x.colors,
clust = diana(d),
scaling.factor = 1,
name = NA,
divide.hz = FALSE,
y.offset = 10,
width = 0.2
)
addBracket(depth.brackets, col='black', label.cex=0.6)
title('Aggregate Soil Properties')
par(mar=c(1,0,3,0))
plotSPC(x.colors, divide.hz=FALSE, color='clay', name=NA, plot.order=new.order, col.label='Clay Content (%)', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addBracket(depth.brackets, col='black', label.cex=0.85)
plot(x.colors, divide.hz=FALSE, color='phfield', name=NA, plot.order=new.order, col.label='pH', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addBracket(depth.brackets, col='black', label.cex=0.85)
plot(x.colors, divide.hz=FALSE, color='total_frags_pct', name=NA, plot.order=new.order, col.label='Total Rock Fragment Volume (%)', lwd=1.25, axis.line.offset=-3, cex.depth.axis=1, cex.id=1)
addVolumeFraction(x.colors, 'total_frags_pct')
addBracket(depth.brackets, col='black', label.cex=0.85)
d <- profile_compare(x.colors, vars=c('clay', 'phfield', 'total_frags_pct'), max_d=150, k=0)
plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25)
plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25, color='clay')
addBracket(depth.brackets, col='black', label.cex=0.6)
addVolumeFraction(x.colors, 'total_frags_pct')
plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25, color='phfield')
plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25, color='total_frags_pct')
d <- profile_compare(x.colors, vars=c('d_L', 'clay', 'phfield', 'total_frags_pct'), max_d=150, k=0)
plotProfileDendrogram(x.colors, diana(d), scaling.factor = 1, name = NA, divide.hz=FALSE, y.offset=10, width = 0.25)
addBracket(depth.brackets, col='black', label.cex=0.6, offset = -0.35)
pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'colorBucket')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname
head(pig)## taxonname .white.pigment .red.pigment .green.pigment .yellow.pigment .blue.pigment
## Amador Amador 0.7534900 0.03788857 0.000000000 0.2086214 0
## Argonaut Argonaut 0.6021182 0.10820262 0.000000000 0.2896792 0
## Auburn Auburn 0.7076987 0.07564808 0.015961548 0.2006917 0
## Copperopolis Copperopolis 0.7264037 0.03890978 0.000000000 0.2346865 0
## Crimeahouse Crimeahouse 0.5267660 0.17433608 0.000000000 0.2988979 0
## Loafercreek Loafercreek 0.6496804 0.10102108 0.009483642 0.2398149 0d <- daisy(pig[, -1])
dd <- diana(d)
par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.002, y.offset = 0.04, width = 0.25, name = NA, divide.hz=FALSE)
addBracket(depth.brackets, col='black', label.cex=0.6, offset = -0.35)
pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'depthSlices')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname
d <- daisy(pig[, -1])
dd <- diana(d)
par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.5, y.offset = 7, width=0.25, name = NA, divide.hz=FALSE)
addBracket(depth.brackets, col='black', label.cex=0.6, offset = -0.35)
pig <- soilColorSignature(x.colors, r='d_r', g='d_g', b='d_b', RescaleLightnessBy = 1, method = 'pam')
# move row names over for distance matrix
row.names(pig) <- pig$taxonname
d <- daisy(pig[, -1])
dd <- diana(d)
par(mar=c(1,1,1,1))
plotProfileDendrogram(x.colors, dd, dend.y.scale = max(d) * 2, scaling.factor = 0.5, y.offset = 7, width=0.25, name = NA, divide.hz=FALSE)
addBracket(depth.brackets, col='black', label.cex=0.6, offset = -0.35)
This document is based on aqp version 1.27 and soilDB version 2.5.9.