This is a work in progress, please see the in-line comments for a description of what is going on here.
# load required libraries
library(soilDB)
library(aqp)
library(sharpshootR)
library(sf)
library(scales)
# sample data from soilDB packaeg
data("loafercreek")
# keep first 10 profiles
x <- loafercreek[1:10, ]
# default margins are too wide, make smaller
par(mar=c(1,1,1,1))
# profile sketches
plotSPC(x, label='pedon_id', id.style='side')# group profiles by hillslope position
groupedProfilePlot(x, 'hillslopeprof', label='pedon_id', id.style='side')# extract site data so that we can annotate wit PSCS depths
s <- site(x)
# reminder of ID used by SPC
idname(x)## [1] "peiid"# subset site data into ID and PSCS top / bottom
s <- s[, c('peiid', 'psctopdepth', 'pscbotdepth')]
# re-name for addBrackets() function
names(s)[2:3] <- c('top', 'bottom') # plot sketches and annotate with PSCS
par(mar = c(1, 1, 1, 1))
# omit pedon IDs for now
plotSPC(x, print.id = TRUE, max.depth = 125, name.style = 'center-center', depth.axis = list(line = -3))
# add PSCS brackets as vertical var
addBracket(s, hzDepths = c('top', 'bottom'), offset = -0.38, tick.length = 0, lwd = 10, col = rgb(red=0, green=0, blue=1, alpha=0.25))
# check structure of diagnostic horizon data
head(diagnostic_hz(x))## peiid featkind featdept featdepb
## 1 115595 ochric epipedon 0 10
## 2 115595 argillic horizon 10 74
## 3 115595 paralithic contact 74 81
## 4 115819 ochric epipedon 0 10
## 5 115819 argillic horizon 10 74
## 6 115819 paralithic contact 74 81# tabulate and sort
# sort(table(diagnostic_hz(x)$featkind))
# add brackets for some diagnostic features
addDiagnosticBracket(x, offset = -0.38, kind='argillic horizon', col='red')
addDiagnosticBracket(x, offset = -0.38, kind='paralithic contact', col='black')Align profiles to an elevation gradient, regular spacing of sketches.
# experiment with re-ordering sketches
# order according to elevation as described in the field
new.order <- order(x$elev_field)
new.order## [1] 8 1 2 9 10 7 3 6 4 5# plot profiles in the new order
par(mar=c(4.5,1,1,1))
plotSPC(x, plot.order=new.order, print.id=TRUE)
# the "brackets" will automatically follow the new ordering
addBracket(s, hzDepths = c('top', 'bottom'), tick.length = 0, lwd=10, col=rgb(red=0, green=0, blue=1, alpha=0.25))
addDiagnosticBracket(x, kind='argillic horizon', col='red')
addDiagnosticBracket(x, kind='paralithic contact', col='black')
# add an axis with elevation associated with each profile
axis(side=1, at = 1:length(x), labels = x$elev_field[new.order], line = 1)
mtext('Elevation (m)', side = 1, line = 3.25)Align profiles to an elevation gradient, relative spacing of
sketches. Note that figure elements need to be explicitly ordered with
[new.order].
pos <- rescale(x$elev_field, to = c(1, length(x)))
pos <- fixOverlap(pos, thresh = 0.65)
new.order <- order(x$elev_field)
# plot profiles in the new order
par(mar=c(4.5,1,1,1))
plotSPC(x, plot.order = new.order, print.id = TRUE, relative.pos = pos[new.order], , max.depth = 125, name.style = 'center-center', depth.axis = list(line = -3))
# the "brackets" will automatically follow the new ordering
addBracket(s, hzDepths = c('top', 'bottom'), offset = -0.33, tick.length = 0, lwd=8, col=rgb(red=0, green=0, blue=1, alpha=0.25))
addDiagnosticBracket(x, offset = -0.33, kind='argillic horizon', col='red')
# add an axis with elevation associated with each profile
axis(side=1, at = pos[new.order], labels = x$elev_field[new.order], line = 1)
mtext('Elevation (m) | Not to Scale', side = 1, line = 3.25)Test out some new features in aqp.
par(mar=c(4.5,1,3,1))
plotSPC(x, color='total_frags_pct', plot.order=new.order, print.id=FALSE, max.depth = 125, name.style = 'center-center', depth.axis = list(line = -3))
# the "brackets" will automatically follow the new ordering
addBracket(s, hzDepths = c('top', 'bottom'), offset = -0.38, tick.length = 0, lwd=10, col=rgb(red=0, green=0, blue=1, alpha=0.25))
addDiagnosticBracket(x, offset = -0.38, kind='argillic horizon', col='red')
addDiagnosticBracket(x, offset = -0.38, kind='paralithic contact', col='black')
addVolumeFraction(x, 'total_frags_pct')plotSPC(x, color='total_frags_pct', plot.order=new.order, print.id=FALSE, relative.pos=jitter(1:length(x)), max.depth = 125, name.style = 'center-center', depth.axis = list(line = -3))
# the "brackets" will automatically follow the new ordering
addBracket(s, hzDepths = c('top', 'bottom'), offset = -0.38, tick.length = 0, lwd=10, col=rgb(red=0, green=0, blue=1, alpha=0.25))
addDiagnosticBracket(x, offset = -0.38, kind='argillic horizon', col='red')
addDiagnosticBracket(x, offset = -0.38, kind='paralithic contact', col='black')
addVolumeFraction(x, 'total_frags_pct')plotSPC(x, color='total_frags_pct', print.id=FALSE, relative.pos=jitter(1:length(x)), max.depth = 125, name.style = 'center-center', depth.axis = list(line = -3))
# the "brackets" will automatically follow the new ordering
addBracket(s, hzDepths = c('top', 'bottom'), offset = -0.38, tick.length = 0, lwd=10, col=rgb(red=0, green=0, blue=1, alpha=0.25))
addDiagnosticBracket(x, offset = -0.38, kind='argillic horizon', col='red')
addDiagnosticBracket(x, offset = -0.38, kind='paralithic contact', col='black')
addVolumeFraction(x, 'total_frags_pct')Grouped soil profile sketches.
# example data
data("mineralKing")
# SPC sketch style for next couple of plots
options(.aqp.plotSPC.args = list(max.depth = 175, name.style = 'center-center', depth.axis = list(line = -3)))
par(mar = c(1,1,2,1))
# crowded
groupedProfilePlot(
mineralKing,
groups = 'taxonname'
)# vertical offsets for group labels
groupedProfilePlot(
mineralKing,
groups = 'taxonname',
group.name.offset = c(-10, -15),
id.style = 'side'
)# suppress profile ID labeling
groupedProfilePlot(
mineralKing,
groups = 'taxonname',
print.id = FALSE,
group.name.offset = c(-10, -15),
)# shift group lines to the right
groupedProfilePlot(
mineralKing,
groups = 'taxonname',
print.id = FALSE,
group.name.offset = c(-10, -15),
break.offset = 0.66
)# move group breaks
groupedProfilePlot(
mineralKing,
groups = 'taxonname',
print.id = FALSE,
group.name.offset = c(-10, -15),
break.offset = 0.6,
)# setup point locations
s <- site(mineralKing)
xy <- st_as_sf(s, coords = c('longstddecimaldegrees', 'latstddecimaldegrees'))
st_crs(xy) <- 4326
# convert to suitable projected cRS
# projected CRS, UTM z11 NAD83 (https://epsg.io/26911)
xy <- st_transform(xy, 26911)
par(mar=c(4.5,4,4,1))
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', scaling.factor = 1)plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', color='clay', scaling.factor = 1)Tinker with relative sketch positioning.
# standard transect plot, profile sketches arranged along integer sequence
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', name='hzname', scaling.factor = 1)
title('Regular Spacing')# attempt relative positioning based on scaled distances, no corrections for overlap
# profiles are clustered in space and therefore over-plot
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', name='hzname', spacing = 'relative', fix.relative.pos = FALSE, scaling.factor = 1)
title('Relative Spacing | No Overlap Fix')# default behavior, attempt adjustments to prevent over-plot and preserve relative spacing
# use set.seed() to fix outcome
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', name='hzname', spacing = 'relative', scaling.factor = 1)
title('Relative Spacing | Default Overlap Fix')# customize arguments to aqp::fixOverlap()
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', name='hzname', spacing = 'relative', fix.relative.pos = list(maxIter=6000, adj=0.2, thresh=0.7), scaling.factor = 1)
title('Relative Spacing | Custom Overlap Fix')# electrostatic simulation
plotTransect(mineralKing, xy, grad.var.name='elev_field', grad.axis.title='Elevation (m)', label='pedon_id', name='hzname', spacing = 'relative', fix.relative.pos = list(method = 'E', q = 1.85, thresh = 0.65), scaling.factor = 1)
title('Relative Spacing | Electrostatic Simulation')# load example dataset
data(gopheridge)
# randomly select 20 profiles
gopheridge <- sample(gopheridge, 20)
# define a new function
# get hz mid-point of hz with max clay content
f.max.clay.depth <- function(x) {
h <- horizons(x)
max.clay <- which.max(h$clay)
res <- with(h[max.clay, ], (hzdept+hzdepb) / 2)
return(res)
}
# apply function to each profile
# results are put into @site
gopheridge$max.clay.depth <- profileApply(gopheridge, f.max.clay.depth)
# reset style
options(.aqp.plotSPC.args = NULL)
# plot and annotate with max clay depths
par(mar=c(0,0,0,0))
plotSPC(gopheridge, name='hzname', print.id=FALSE, max.depth = 175, depth.axis = list(line = -3), width = 0.25)
x.pos <- 1:length(gopheridge)
points(x.pos, gopheridge$max.clay.depth, col='white', pch=21, cex=0.75, bg='black')This document is based on aqp version 2.3.1,
soilDB version 2.9.0, and sharpshootR version
2.4.1.