Introduction

Plotting several SoilProfileCollection objects in the same context can be challenging. This document illustrates two methods for organizing profiles into groups within a single sketch.

Get Data

We are going to demonstrate some concepts using both KSSL (lab) and OSD (typical morphology) data for the Sobrante and Auburn soil series.

library(aqp)
library(soilDB)
library(scales)
library(RColorBrewer)

# get OSD data for select series
osd <- fetchOSD(c('sobrante', 'auburn'))

# get KSSL data for these series
kssl.1 <- fetchKSSL('sobrante')
kssl.2 <- fetchKSSL('auburn')

Note that the KSSL-derived SPC objects share the same structure, but the OSD-derived object does not. This will have implications on how the data can be displayed within the same sketch.

str(kss1.1)
str(kss1.2)
str(osd)

Method 1: Split Sketch into Groups

We can easily combine the two KSSL-derived objects using rbind(). The colors and legend can be automatically generated from the entire set of (combined) clay values.

# combine similar SPCs into a new SPC
g <- rbind(kssl.1, kssl.2)

# there is some variabilty in capitalization: need to normalize these names
table(g$taxonname)
## 
##   Auburn   AUBURN Sobrante SOBRANTE 
##       16        2        2        4
# normalize taxonname (soil series) via capitalization
g$taxonname <- toupper(g$taxonname)

# plot groups of profiles organized by normalized taxonname
par(mar=c(1,1,4,3))
groupedProfilePlot(g, groups = 'taxonname', name='hzn_desgn', color='clay', label='pedon_id', id.style='side')

Method 2: Plot SPC Objects of Variable Provenance

SoilProfileCollection objects from any source can be included into the same sketch using plotMultipleSPC(). However, the symbology of soil property values via color must be manually defined.

Create "Merged" Color Scheme

# color ramp function
cr <- colorRamp(rev(brewer.pal(10, 'Spectral')))

# get the full range of clay
combined.data <- c(kssl.1$clay, kssl.2$clay)
combined.range <- range(combined.data, na.rm = TRUE)

# NA-padded value -> color mapping for full range of some horizon attribute
mapColor <- function(x, r, col.ramp) {
  c.rgb <- cr(scales::rescale(x, from=r))
  cc <- which(complete.cases(c.rgb))
  cols <- rep(NA, times=nrow(c.rgb))
  cols[cc] <- rgb(c.rgb[cc, ], maxColorValue=255)
  return(cols)
}

# convert non-NA values into colors
kssl.1$.color <- mapColor(kssl.1$clay, combined.range, cr)
kssl.2$.color <- mapColor(kssl.2$clay, combined.range, cr)

# generate combined range / colors for legend
pretty.vals <- pretty(combined.data, n = 8)
legend.data <- list(legend=pretty.vals, col=rgb(cr(scales::rescale(pretty.vals)), maxColorValue=255))

# assemble into list
spc.list <- list(osd, kssl.1, kssl.2)

Plot

# assemble list of arguments to plotSPC for each object
# note the use of .color for pre-computed value->color 
args.list <- list(
  list(name='hzname', color='soil_color', id.style='side'), 
  list(name='hzn_desgn', color='.color', label='pedon_id', id.style='side'), 
  list(name='hzn_desgn', color='.color', label='pedon_id', id.style='side')
  )

par(mar=c(1,1,3,3))
plotMultipleSPC(spc.list, group.labels=c('OSD', 'Sobrante', 'Auburn'), args=args.list)

# legend is the combined legend
mtext(side=3, text='Clay Content (%)', font=2, line=1.6)
legend('bottom', legend=legend.data$legend, col=legend.data$col, bty='n', pch=15, horiz=TRUE, xpd=TRUE, inset=c(0, 0.99))


This document is based on aqp version 1.15.3 and soilDB version 2.0-1.