In this demo, we use the aqp functions glom() and profileApply() to analyze a depth-subset of horizons from each profile in a SoilProfileCollection.

glom() returns the horizon(s) that intersect the specified depth [z1] or depth interval [z1, z2].

Checking for cambic horizon criteria

First, the script queries components by component name, then gets records populated in the Component Diagnostic Features (codiagfeatures) table.

glom() is used to get all horizons that are within the cambic depth range. A custom function has_sandy_cambic() is defined to use glom() on individual profiles (p).

A regular expression pattern is used to check the glom()-ed horizon textures field (p$texture) for stratified and/or sandy texture classes that are not allowed in cambic horizons.

If either stratified or sandy textures are found, the function returns TRUE for that profile. Otherwise, the function returns FALSE.

This demonstration script loads some components from Yosemite National Park (CA790). In this simple example, we look at some components correlated to Glacierpoint.

The Glacierpoint series (Sandy-skeletal, isotic Xeric Dystrocryepts) is defined as having either an Umbric epipedon OR a cambic horizon. It appears that cambic horizons were populated even for some of the components whose RV textures do not meet cambic criteria.

library(aqp)
library(soilDB)

# get components from Yosemite National Park for a sandy-skeletal, cryic, inceptisol
my.spc <- fetchSDA(WHERE="compname = 'Glacierpoint'")

Now, access the SPC @diagnostics slot to figure out which components have cambics, and then make a subset SPC my.spc.sub with just those components.

# copy contents of diagnostics slot
diagz <- diagnostic_hz(my.spc)

# create a subset SPC with just profiles that have a cambic in their diagnostics
my.spc.sub <- my.spc[profile_id(my.spc) %in% 
                       diagz[grepl(diagz$featkind, pattern="[Cc]ambic"), idname(my.spc)],]

A total of 3 components out of 8 have cambic horizons.

Let’s inspect those 3 components visually.

Use plotSPC() to show sand fraction percentage by horizon – with cambic top and bottom depths shown with brackets using aqp::addDiagnosticBracket().

# set margins
par(mar=c(0,1,3,1))

# plot spc using total sand RV % as color
plotSPC(my.spc.sub, 
        color='sandtotal_r', 
        width=0.15, 
        cex.names = 1)

# add red brackets showing cambic horizon top and bottom depth
addDiagnosticBracket(my.spc.sub, 
                     kind='Cambic horizon', 
                     top='featdept_r',
                     bottom='featdepb_r',
                     col='red', 
                     offset=0,
                     lwd=10)

Now that we have the data and have done some preliminary inspection, let’s define a regular expression pattern and a function to check our cambic criteria.

Define a custom function has_sandy_cambic() to test horizons within the cambic depth range for each profile. The default argument sandy.texture.pattern matches stratified or sandy textures.

# has_sandy_cambic() has one required argument p, which is a single profile SPC 
has_sandy_cambic <- function(p, sandy.texture.pattern = "SR|-S$|^S$|COS$|L[^V]FS$|[^L]VFS$|LS$|LFS$") {
  d <- diagnostic_hz(p)
  
  # match the pattern cambic or Cambic
  d.cambic.idx <- grep(d$featkind, pattern="[Cc]ambic")
  cambic_top <- d[d.cambic.idx, 'featdept_r']
  cambic_bot <- d[d.cambic.idx, 'featdepb_r']
  
  # if top or bottom depth are NULL there is no cambic
  if(!is.null(cambic_top) | !is.null(cambic_bot)) {
    
    # glom horizons in p from cambic_top to cambic_bot 
    cambic_horizons <- glom(p, cambic_top, cambic_bot, as.data.frame = TRUE)
    
    print(paste0("COKEY: ",unique(p$cokey)," has cambic (",
                 cambic_top," - ",cambic_bot," cm) with ", paste0(cambic_horizons$hzname, collapse=","), 
                 " horizons of texture class: ", paste0(cambic_horizons$texture, collapse=", ")))
  
    # use a regular expression to match the sandy textures
    if(any(grepl(cambic_horizons$texture, 
                 pattern = sandy.texture.pattern, 
                 ignore.case = TRUE))) {
      return(TRUE) 
    }  
  }
  return(FALSE)
}

Apply the function to each profile.

# apply function has_sandy_cambic() to each profile
my.spc.sub$sandy_cambic <- profileApply(my.spc.sub, has_sandy_cambic)
## [1] "COKEY: 17639151 has cambic (14 - 45 cm) with Bw1,Bw2 horizons of texture class: cb-ls, cbv-lfs"
## [1] "COKEY: 17639343 has cambic (28 - 49 cm) with Bw horizons of texture class: byx-cosl"
## [1] "COKEY: 17639416 has cambic (41 - 69 cm) with Bw horizons of texture class: stv-ls"

Inspect results. Print out the offending unique profile IDs (cokey).

# do any components have a sandy cambic?
if(any(my.spc.sub$sandy_cambic)) {
  print(paste0("The following components have sandy cambics: ",
        paste0(profile_id(my.spc.sub)[my.spc.sub$sandy_cambic], collapse=",")))
} else {
  print("All components with cambics have appropriate textures.")
}
## [1] "The following components have sandy cambics: 17639151,17639416"

This document is based on aqp version 1.19, soilDB version 2.5, and sharpshootR version 1.5.