It is possible to have multiple "densic" in a profile, with only the bottom-most horizon meeting the criteria for "contact". One solution is to integrate diagnostic features into the calculation.
See the attached data from TX for reference.
pedons.zip
library(aqp)
library(soilDB)
library(plyr)
library(lattice)
library(reshape2)
# load Rich's data
load('pedons.rda')
# get selected set of pedons, no removal of pedons with inconsistent horizonation
# x <- fetchNASIS(rmHzErrors = FALSE)
# compute soil depth and depth class for all pedons, using pattern matching of horizon names
sdc <- getSoilDepthClass(x, p = 'Cr|R|Cd')
head(sdc)
# join depth and depth class data to pedons
site(x) <- sdc
# normalize taxonname
soils <- c('Houston Black', 'Heiden', 'Ferris', 'Vertel')
for(i in soils)
x$taxonname[grep(i, x$taxonname, ignore.case = TRUE)] <- i
table(x$taxonname)
# access diagnostict HZ data, 1:many per pedon
d <- diagnostic_hz(x)
# keep just records of "densic" anything
idx <- grep('densic', d$diag_kind, ignore.case = TRUE)
d <- d[idx, ]
table(d$diag_kind)
# convert long -> wide format: top depths
d.wide <- dcast(d, peiid ~ diag_kind, value.var='featdept')
# fix names
names(d.wide) <- c('peiid', 'densic.contact.top', 'densic.materials.top')
# join to SPC
site(x) <- d.wide
# convert long -> wide format: bottom depths
d.wide <- dcast(d, peiid ~ diag_kind, value.var='featdepb')
# fix names
names(d.wide) <- c('peiid', 'densic.contact.bottom', 'densic.materials.bottom')
# join to SPC
site(x) <- d.wide
# check depth via hz name pattern matching vs. diagnostic features
head(site(x)[, c('pedon_id', 'depth', 'densic.contact.top', 'densic.materials.top', 'densic.materials.bottom')], 10)
# rules for refinement of "soil depth"
x$depth_revised <- rep(NA, times=length(x))
# 1. "densic" diagnostic features not present or populated: use depth via hzname
idx <- which(is.na(x$densic.materials.top) & is.na(x$densic.contact.top))
x$depth_revised[idx] <- x$depth[idx]
# 2. densic contact present: use top depth of densic contact
idx <- which(!is.na(x$densic.contact.top))
x$depth_revised[idx] <- x$densic.contact.top[idx]
# 3. densic contact not present, but densic materials present: use bottom of densic materials
idx <- which(is.na(x$densic.contact.top) & !is.na(x$densic.materials.bottom))
x$depth_revised[idx] <- x$densic.materials.bottom[idx]
# what is left?
# errors in data population: missing "densic materials" bottom depth
# TODO: fix these in NASIS
site(x)[is.na(x$depth_revised), c('pedon_id', 'depth', 'depth_revised', 'densic.contact.top', 'densic.materials.top', 'densic.materials.bottom')]
## Note: use depth_revised for subsequent calculations
## Note: previously estimated depth classes are no longer correct
# what fraction of pedons had the "wrong" soil depth?
# about 33% ! wow !
prop.table(table(x$depth_revised > x$depth))
It is possible to have multiple "densic" in a profile, with only the bottom-most horizon meeting the criteria for "contact". One solution is to integrate diagnostic features into the calculation.
See the attached data from TX for reference.
pedons.zip