Setup
# stable packages from CRAN
install.packages('aqp')
install.packages('soilDB')
install.packages('sharpshootR')
install.packages('FactoMineR')
install.packages('cowplot')
# latest versions from GitHub
remotes::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade_dependencies=FALSE)
remotes::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)
remotes::install_github("ncss-tech/sharpshootR", dependencies=FALSE, upgrade_dependencies=FALSE)
Hillslope Position / Soil Morphology
library(aqp)
library(soilDB)
library(sharpshootR)
# soil of interest
# https://casoilresource.lawr.ucdavis.edu/gmap/?loc=44.10219,-122.38589,z14
s <- 'peavine'
# s <- 'clarion'
# all siblings
sib <- siblings(s)
# all soils of interest
soils <- c(s, sib$sib$sibling)
# get morphology + extended summaries
osds <- fetchOSD(soils, extended = TRUE)
There may be missing series from either / both SPC and geomorphic
summaries. TBC.
# there may be records missing from SPC / hill slope position
nm <- intersect(profile_id(osd$SPC), osd$hillpos$series)
# keep only those series that exist in both
sub <- subset(osd$SPC, profile_id(osd$SPC) %in% nm)
## inverse problem: extra records in hill slope summaries
# subset hillpos
hillpos.sub <- subset(osd$hillpos, subset = series %in% profile_id(sub))
## TODO: use subset objects below
res <- vizHillslopePosition(osds$hillpos)
# re-order hillslope proportions according to clustering
hp <- osds$hillpos[res$order, ]
nm <- names(hp[, 2:6])
series
|
Toeslope
|
Footslope
|
Backslope
|
Shoulder
|
Summit
|
n
|
shannon_entropy
|
REMOTE
|
0.00
|
0.00
|
0.53
|
0.13
|
0.33
|
30
|
1.40
|
DIGGER
|
0.00
|
0.00
|
0.59
|
0.19
|
0.22
|
64
|
1.38
|
FORMADER
|
0.00
|
0.00
|
0.43
|
0.36
|
0.21
|
14
|
1.53
|
BOHANNON
|
0.00
|
0.00
|
0.46
|
0.37
|
0.18
|
136
|
1.49
|
SHIVIGNY
|
0.00
|
0.00
|
0.59
|
0.34
|
0.06
|
32
|
1.23
|
MEDA
|
0.00
|
0.25
|
0.25
|
0.25
|
0.25
|
4
|
2.00
|
PEAVINE
|
0.13
|
0.13
|
0.31
|
0.23
|
0.21
|
39
|
2.24
|
SLICKROCK
|
0.03
|
0.25
|
0.33
|
0.15
|
0.24
|
72
|
2.08
|
PREACHER
|
0.10
|
0.21
|
0.40
|
0.06
|
0.23
|
113
|
2.07
|
HONEYGROVE
|
0.19
|
0.29
|
0.26
|
0.00
|
0.26
|
80
|
1.98
|
HEMCROSS
|
0.18
|
0.38
|
0.21
|
0.00
|
0.23
|
39
|
1.93
|
PANTHER
|
0.18
|
0.82
|
0.00
|
0.00
|
0.00
|
146
|
0.68
|
MINNIECE
|
0.57
|
0.43
|
0.00
|
0.00
|
0.00
|
7
|
0.99
|
print(res$fig)
par(mar=c(0,0,1,1))
SoilTaxonomyDendrogram(osds$SPC, scaling.factor = 0.015, width = 0.28, name.style = 'center-center')
# colors
par(mar = c(1, 1, 1, 1))
hp.cols <- RColorBrewer::brewer.pal(n = 5, name = 'Set1')[c(2, 3, 4, 5, 1)]
soilPalette(hp.cols, lab = nm)
Composite Figures
Version 1: Lines
par(mar = c(0.5, 0, 0, 2))
layout(matrix(c(1,2)), widths = c(1,1), heights = c(2,1))
plotProfileDendrogram(osds$SPC, res$clust, dend.y.scale = 3, scaling.factor = 0.012, y.offset = 0.2, width = 0.32, name.style = 'center-center', cex.names = 0.7, shrink = TRUE, cex.id = 0.55)
matplot(y = hp[, 2:6], type = 'b', lty = 1, pch = 16, axes = FALSE, col = hp.cols, xlab = '', ylab = '', xlim = c(0.5, length(osds$SPC) + 1))
# grid(nx = 0, ny = NULL)
axis(side = 4, line = -1, las = 1, cex.axis = 0.7)
# axis(side = 2, line = -3, las = 1, cex.axis = 0.7)
legend('topleft', legend = rev(nm), col = rev(hp.cols), pch = 16, bty = 'n', cex = 0.8, pt.cex = 2, horiz = TRUE, inset = c(0.01, 0.01))
mtext('Probability', side = 2, line = -2, font = 2)
Version 2: Stacked Bars
par(mar = c(0.5, 0, 0, 2))
layout(matrix(c(1,2)), widths = c(1,1), heights = c(2,1))
plotProfileDendrogram(osds$SPC, res$clust, dend.y.scale = 3, scaling.factor = 0.012, y.offset = 0.2, width = 0.32, name.style = 'center-center', cex.names = 0.7, shrink = TRUE, cex.id = 0.55)
sp <- c(1.5, rep(1, times = length(osds$SPC) - 1))
barplot(height = t(as.matrix(hp[, 2:6])), beside = FALSE, width = 0.5, space = sp, col = hp.cols, axes = FALSE, xlab = '', ylab = '', xlim = c(0.5, length(osds$SPC) + 1), ylim = c(0, 1.2))
legend(x = 0.75, y = 1.2, legend = rev(nm), col = rev(hp.cols), pch = 15, bty = 'n', cex = 0.8, pt.cex = 1.25, horiz = TRUE)
mtext('Probability', side = 2, line = -2, font = 2)
3D Surface Shape
res <- vizSurfaceShape(osds$shape_across, title = 'Surface Shape (Across)')
# re-order proportions according to clustering
ss <- osds$shape_across[res$order, ]
nm <- names(ss[, 2:6])
series
|
Concave
|
Linear
|
Convex
|
Complex
|
Undulating
|
n
|
shannon_entropy
|
SHIVIGNY
|
0.00
|
0.00
|
1.00
|
0
|
0
|
32
|
0.00
|
FORMADER
|
0.12
|
0.12
|
0.76
|
0
|
0
|
17
|
1.02
|
DIGGER
|
|