mollic thickness requirement, mollic epipedon, and more

[brownag]

See branch aqp/mollic

Here is a demonstration of the new functionality I am working on for evaluating taxonomic criteria with aqp.

mollic.thickness.requirement calculates the minimum thickness of the materials meeting other requirements for a mollic epipedon, per criterion 6 in U.S. Soil Taxonomy (12th edition). I have used it for QC of pedon data, and also to assess questions about how frequently e.g. the sliding scale portion of the thickness requirement is invoked.

I'll have more fun graphs to show -- using more data -- but for now here is a demo showing how thick the mollic would have to be, and the boundaries of PSCS, for a subset of loafercreek.

UPDATE (05/18/2020): to include better cambic finding via getCambicBounds and to portray "untruncated" sliding scale thickness requirements for QC.

library(aqp)
library(soilDB)

data("loafercreek")

# subset loafercreek -- profiles with all horizons with clay between 20 and 35 
#   mixture of argillic/cambics
foo <- mutate_profile(loafercreek, in_range = all(na.omit(clay > 20 & clay < 35)) )
s <- filter(foo, in_range)

# make a plot
par(mar=c(0,0,3,0))
plotSPC(s, print.id = F, axis.line.offset = -3.5)
mtext("aqp calculated\nMollic/Umbric Minimum Thickness Requirement\n& Particle Size Control Section Boundaries", 3)
legend("bottomleft", c("particle size control section",
                       "minimum mollic/umbric thickness",
                       "minimum mollic/umbric thickness (no truncate)",
                       "absolute depths: 18, 25, 54, 75 cm"), lwd=c(3,3,1,1), 
       col=c("black","green","blue","black"))

# calculate PSCS boundaries and minimum mollic thickness
s.diag <- profileApply(s, FUN = function(p) {

  # growing suite of aqp functions related to application of US soil taxonomic criteria
  ep <- estimatePSCS(p)
  mss <- getMineralSoilSurfaceDepth(p)
  mtt <- mollic.thickness.requirement(p)
  mt <- mollic.thickness.requirement(p, truncate = FALSE)

  # construct result for frameify/merge into @diagnostic slot
  df1 <- data.frame(feat="particle size control section", top=ep[1], bot=ep[2])
  df2 <- data.frame(feat="minimum mollic/umbric thickness", top=mss, bot=mss+mtt)
  df3 <- data.frame(feat="minimum mollic/umbric thickness (no truncate)", top=mss, bot=mss+mt)
  df <- data.frame(id=profile_id(p), rbind(df1, df2, df3))
  colnames(df) <- c(idname(p), "featkind", "featdept","featdepb")

  return(df) 

}, frameify=TRUE)

# co-opting the @diagnostics slot and related plotting functionality
#   these are diagnostic criteria/depths/control sections, rather than "features" per se
diagnostic_hz(s) <- rbind(diagnostic_hz(s), s.diag)

addDiagnosticBracket(s, kind = 'particle size control section', lwd=3)
addDiagnosticBracket(s, kind = 'minimum mollic/umbric thickness', lwd=3, 
                     offset=0, tick.length=0, col="green")
addDiagnosticBracket(s, kind = 'minimum mollic/umbric thickness (no truncate)', 
                     lwd=1, offset=0, tick.length=0, col="blue")
abline(h=c(18,25,54,75))

Obligatory profile plot

On a related note, along with this branch will come improvements to the hzdesgnname/hztexclname slot usage and better / consistent methods for guessing horizon attributes of interest to analyses like these (argillic, mollic, PSCS) that are used in the default arguments to various functions that take a SPC as an argument.

• guessHzDesgnName
• guessHzTexClName
• guessHzAttrName

There are many notes in the mollic.thickness.requirement code comments about things that need to be expanded on pertaining to other diagnostic feature identification. Currently, this relies on variants of a new set of aqp methods called depthOf.

EDIT (5/24/2020): there have been significant updates to getCambicBounds, the depthOf family and other critical functions that mollic.thickness.requirement relies on. The results are much more stable with the Region 2 pedon dataset. There is still room to expand on natric, oxic, spodic horizon identification beyond horizon designations -- so I'll be seeking additional test cases for those next up.

[brownag]

I was curious about answering the question “How often is Mollic criterion 6c used to determine the minimum thickness requirement?” so I wrote up a basic calculation that applies the thickness criteria.

I have ran mollic.thickness.requirement on over 31,000 NASIS pedons from Region 2 (thanks to @smroecker for soilDB::fetchNASIS('pedon_report')) Basically, it is using some assumptions to identify the relevant diagnostic features and properties for minimum thickness.

EDIT: significant updates to functions that estimate soil diagnostic feature boundaries have rendered the values previously reported here obsolete. I previously reported that 70% of pedons would see an increase in minimum thickness requirement (under proposed 25cm criterion), when in fact it is only approximately 50% (6d+6c v.s. 6a)

         criterion  Freq total        prop
1               6a 13790 31573 0.436765591
4               6d 12521 31573 0.396573021
3               6c  3734 31573 0.118265607
2               6b  1476 31573 0.046748804
5 cannot calculate    52 31573 0.001646977

[MollicMeyer]

I'm not sure if this is the correct place to post, but I'm very interested in translating different soil properties maps created from spatial prediction methods back into something interpretable, i.e., taxonomic classification. So that requires classifying depths and thicknesses of different diagnostic horizons.

The example property I'm working with is the depth of mollic colors, or a quasi proxy for mollic epipedon thickness. I've found it challenging to use horizon names for criteria, because naming is very biased, but color is a bit more quantitative (hopefully more reliable).

I was trying to create a function that could get the depth of contiguous horizons with mollic colors, as to avoid buried A's being included in the calculation. I see that getSurfaceHorizonDepth does this, but doesn't account for situations where the property depth coincides with the maximum depth of exploration. So I addressed it as such:

Here's the data: NASIS.zip

##creating new SPC object
nas <- read.csv('./tables/NASIS/Nasis_DSMlobe_hzs.csv')

depths(nas) <- peiid ~ hzdept + hzdepb

#append site level information
nas.site <- read.csv('./tables/Nasis/IA_DSMlobe_site.csv')

site(nas) <- nas.site
#I think that worked???
##now let's clean this up
n.chk <- checkHzDepthLogic(nas)
idx <- which(n.chk$valid)

# keep the good ones
nas2 <- nas[idx, ]

##get mollic colors
nas2@horizons$mollic <- ifelse(nas2$m_value <= 3 & nas2$m_chroma <= 3, "TRUE", "FALSE")

##determine mollic depth
nas2@site$moldepth2 <- profileApply(nas2, getSurfaceHorizonDepth, pattern = 'TRUE', hzdesgn = 'mollic')

#get bottom depth of exploration
nas.bottom <- function(i) {
  h <- horizons(i)
  bot.depth <- max(h$hzdepb)
  return(bot.depth)
}
##appended bottom depth at site level
nas2@site$total.depth <- profileApply(nas2, nas.bottom)
sort(nas2@site$total.depth)
nas2@site$moldepth <- profileApply(nas2, nas.mollic)

####now to make brackets for viz check
nmol.tops <- rep(0, times = length(nas2))
nas2@site$surface <- nmol.tops
nmols <- nas2@site[, c('peiid', 'surface', 'moldepth2')]
# re-name for addBrackets() function
names(nmols)[2:3] <- c('top', 'bottom') 

##plot first 10 with mollic brackets
plot(nas2[0:10], color='moist_soil_color', name='hzname', id.style='top')
addBracket(nmols, col = 'red')
#looks good

#identify pedons with same bottom depth as mollic depth
nas2@site$molsame2 <- ifelse(nas2$moldepth2 == nas2$total.depth, "TRUE", "FALSE")
prob3 <- subset(nas2, nas2$molsame2== "TRUE")

#Make brackets
nmolsp3 <- prob3@site[, c('peiid', 'surface', 'moldepth2')]

# re-name for addBrackets() function
names(nmolsp3)[2:3] <- c('top', 'bottom')

##Check out 10 with errors - looks like noncontiguous colors are still being included like pedon 58426
plot(prob3[30:40], color = 'moist_soil_color', name = 'hzname', id.style = 'top')
addBracket(nmolsp3, col = 'red')

Long story short, it would be helpful to alter the code to account for max depth of exploration. It would also be great to make a 'getDepthtoFeature' where you could find the top depth of a reduced matrix, for example.

Obviously, I'm crudely altering the code to fit my applications. I looked at the source code for 'getSurfaceHorizonDepth' but I had no clue what I was looking at 😆

[MollicMeyer]

I now see that my getDepthtoFeature is possible with EstimateSoilDepth

[dylanbeaudette]

@MollicMeyer thanks for asking, this is approximately the right place. A lot to unpack here, and as a side-note, makes me think we need better documentation.

[brownag]

Hey @MollicMeyer I have created a whole slew of new replacement functions in that category that generalize the behavior of estimateSoilDepth and related functionality. You might not have my latest branch. These functions are all well documented, but are not part of the master branch yet. Give me a moment and I will see if I can get you an example -- as I designed getSurfaceHorizonDepth and hasDarkColors for exactly the case you are describing.

[brownag]

@MollicMeyer

Thanks for reporting this. Very clever and good for you for trawling through the docs that we do have; it should have worked!

Turns out in that uncommented, impenetrable code of getSurfaceHorizonDepth there was a nasty bug that wasn't properly handling isolated (non-contiguous) instances of horizons that didn't meet the match criteria. It just skipped over them if there were contiguous matching horizons below.

This was never my intention. The test cases I used for this function are simply not robust enough -- as they are relatively simple cases all with traditional horizon designation field. So, thanks for some new tests :)

I fixed the code -- get it with devtools::install_github('ncss-tech/aqp@mollic').

As far as use of estimateSoilDepth related functions -- as you found out that will do what you want either based on g suffix or a calculated field like you did here.

You can also check out the new depthOf/maxDepthOf/minDepthOf in the mollic branch -- they have similar arguments as estimateSoilDepth but are a bit more flexible for criteria that may occur many times within a single profile.

[dylanbeaudette]

A couple of tips:

Don't access SPC slots directly: you can create/edit/remove site or horizon level attributes with the $ or [[ methods. For example:

# note that I removed the quotes around TRUE and FALSE so that they are interpreted as logical vectors
# this works but side-steps all of the consistency enforced by the SPC class
nas2@site$molsame2 <- ifelse(nas2$moldepth2 == nas2$total.depth, TRUE, FALSE)

# safer, consistency checks fire automatically
nas2$molsame2 <-  ifelse(nas2$moldepth2 == nas2$total.depth, TRUE, FALSE)

[brownag]

@MollicMeyer With latest version of mollic branch, the following code works as expected

devtools::install_github('ncss-tech/aqp@mollic')
library(aqp)
setwd("E:/workspace/MollicMeyer")

##creating new SPC object
nas <- read.csv('NASIS/Nasis_DSMlobe_hzs.csv')

depths(nas) <- peiid ~ hzdept + hzdepb

#append site level information
nas.site <- read.csv('NASIS/IA_DSMlobe_site.csv')

site(nas) <- nas.site
n.chk <- checkHzDepthLogic(nas)
idx <- which(n.chk$valid)

# keep the good ones
nas2 <- nas[idx, ]

# use hasDarkColors, ignoring dry value requirement (based on example)
nas2$mollic <- ifelse(hasDarkColors(nas2, d_value = NA),"dark","light")
nas2$mollic[is.na(nas2$mollic)] <- "light"
hzdesgnname(nas2) <- "mollic"
nas2$mollic_depth <- profileApply(nas2, function(p) {
  getSurfaceHorizonDepth(p, pattern="dark")
})
hzdesgnname(nas2) <- "hzname"
nas2$max_depth <- profileApply(nas2, estimateSoilDepth)

# check ones without mollic
subs <- filter(nas2, mollic_depth == 0)
subs$mollic <- factor(subs$mollic)
plot(subs[1:10], color="mollic", label="mollic_depth")

# check ones with mollic equal to max depth
subs <- filter(nas2, mollic_depth == max_depth)
subs$mollic <- factor(subs$mollic)
plot(subs[1:10], color="mollic", label="mollic_depth")

[MollicMeyer]

Thank you both so much! I've been running into a lot of issues getting this training data ready for spatial prediction. NRCS RFOs are due June 11, so i'm only slightly panicking 😅 . Is there a proper forum to post in should I run into more issues along the way? (for example, filtering and/or converting horizon designations from older taxonomy editions)

@brownag I will try this updated code out and get back to you. Once again, greatly appreciated!

and @dylanbeaudette thanks for the tips as well!

[brownag]

@MollicMeyer

The GitHub issue pages are as close as you will get to a "forum" for these aqp R package related topics for now

I would say that items related to the NASIS/SSURGO/LIMS data model specifically might be best discussed over on the https://github.com/ncss-tech/soilDB/issues/ page. If it pertains to these fundamental pattern-matching type functions in aqp, then this is the right place.

Your finding of unexpected/undesired behavior with getSurfaceHorizonDepth could easily have been a new issue here in https://github.com/ncss-tech/aqp/issues/, but your work was also very relevant to this issue (#132) so it was fine that you posted here.

You could post even if there wasn't an existing issue related to your topic. This issue is opened for discussion around the mollic branch, as I implement a variety of new stuff that will be fully refined soon for a CRAN release. It will be closed when the mollic branch gets integrated into master.

We are excited to have people applying these tools. As they get used outside their original scope, things will come up. It is fine to create "Issues" even if they don't turn out to be a bug -- there are enhancement, help wanted labels etc. None of us work on these packages as our primary duty, but generally will help if/when we can!

[MollicMeyer]

@brownag side note for the hasDarkColors function. It works great with EstimateSoilDepth, but will produce an unflagged error if there are any NAs in the color data.

[brownag]

Yes sorry I glossed over that in my example. I filled in any NAs with "light" -- I'll see what I can do to make behavior a bit more stable

[brownag]

@brownag side note for the hasDarkColors function. It works great with EstimateSoilDepth, but will produce an unflagged error if there are any NAs in the color data.

@MollicMeyer I wanted to verify that this function was working as I intended -- it seems like it does. You can make different color criteria optional in hasDarkColors by setting relevant argument/threshold to NA. I haven't been able to produce an actual error from sticking NA values anywhere -- do you have an example?

library(aqp)
library(soilDB)

data(loafercreek)

# default hasDarkColors will apply the "5-3-3" rule
# apply d_value, m_value, m_chroma rules -- returns NA if any of the 3 aren't populated
loafercreek$mollic533 <- hasDarkColors(loafercreek)
table(loafercreek$mollic533, useNA = "ifany")

# apply m_value and m_chroma, ignoring default d_value (some with NA dry value above now return T/F)
loafercreek$mollic33 <- hasDarkColors(loafercreek, d_value = NA)
table(loafercreek$mollic33, useNA = "ifany")

# apply no criteria at all (all horizons return TRUE)
loafercreek$mollicnone <- hasDarkColors(loafercreek, d_value = NA, m_value=NA, m_chroma=NA)
table(loafercreek$mollicnone, useNA = "ifany")

# convert T/F/NA into labels
loafercreek$mollic33 <- ifelse(as.logical(loafercreek$mollic33),"dark","light")

# inspect
plot(loafercreek[30:45,], color="mollic33")

# depth to first light color
estimateSoilDepth(loafercreek[31,], name = "mollic33", p="light")

# change around color assignments a bit
loaf31_modified <- mutate_profile(loafercreek[31, ], 
                                  mollic33_alt = c(NA,"dark","dark","dark","light","light","light",NA))

# now light colors are deeper
estimateSoilDepth(loaf31_modified, name = "mollic33_alt", p="light")

# now turn all attribute values to NA
loaf31_modified2 <- mutate_profile(loafercreek[31, ], 
                                  mollic33_alt = c(NA,NA,NA,NA,NA,NA,NA,NA))

# now bottom of profile is returned (can adjust no.contact.depth/no.contact.assigned as needed)
estimateSoilDepth(loaf31_modified2, name = "mollic33_alt", p="light")

# try new depthOf functions

# all top depths of light colors
depthOf(loafercreek[31,], hzdesgn="mollic33", pattern="light")

# all bottom depths of light colors
depthOf(loafercreek[31,], hzdesgn="mollic33", pattern="light", top=FALSE)

# max depth of dark color
maxDepthOf(loafercreek[31,], hzdesgn="mollic33", pattern="dark", top=FALSE)

[MollicMeyer]

Yep. Check this out if you append to the previous script using my data.

######Test Example for Andrew#####################
idx <- which(is.na(horizons(nas2)$m_hue) | is.na(horizons(nas2)$m_chroma))
ids <- nas2@horizons$peiid[idx]
bad.color <- nas2[site(nas2)$peiid %in% ids, ]

mols.bad <- site(bad.color)[, c('peiid', 'surface', 'mollic_depth')]
# re-name for addBrackets() function
names(mols.bad)[2:3] <- c('top', 'bottom')
plot(bad.color[11:20], color = 'moist_soil_color', name = 'mollic', id.style = 'top')
addBracket(mols.bad, col = 'red')
axis(side=1, at = 1:length(bad.color@site[11:20]), labels =bad.color$mollic_depth[11:20], line = 1)
###################################################

[brownag]

Ok, so here is what I get:

Pedon 1: bracket seems to be in right place Pedon 2: bracket stops at first NA value in mollic -- clearly there are dark colors below first horizon Pedons 3-10: bracket stops at first NA value in mollic -- dark colors are deeper than zero.

I see what the problem is though -- look at the data in the second profile (1092093)

It has a neutral (N) hue and NA for chroma -- so it is right to return NA, since chroma is NA. To work around that, you could try something like this:

EDIT: complete example, fix axis

library(aqp)
setwd("E:/workspace/MollicMeyer")

##creating new SPC object
nas <- read.csv('NASIS/Nasis_DSMlobe_hzs.csv', stringsAsFactors = FALSE)

depths(nas) <- peiid ~ hzdept + hzdepb

#append site level information
nas.site <- read.csv('NASIS/IA_DSMlobe_site.csv', stringsAsFactors = FALSE)

site(nas) <- nas.site
n.chk <- checkHzDepthLogic(nas)
idx <- which(n.chk$valid)

# keep the good ones
nas2 <- nas[idx, ]

# calculate surface
nas2$surface <- profileApply(nas2, getMineralSoilSurfaceDepth)

# get pedons with one or more horizon missing hue or chroma
bad.color.v2 <- filter(nas2, is.na(m_hue) | is.na(m_chroma))

# get colors with m_chroma and m_value <= 3
bad.color.v2$dark_color <- hasDarkColors(bad.color.v2, d_value = NA)

# get neutral hues with m_value <= 3
bad.color.v2$neutral_hue <- hasDarkColors(bad.color.v2, m_hue="N", d_value = NA, m_chroma=NA)

# create combined attribute reflecting dark and neutral colors
bad.color.v2$mollic <- bad.color.v2$dark_color | bad.color.v2$neutral_hue

# calculate mollic color thickness from surface
bad.color.v2$mollic_depth_new <- profileApply(bad.color.v2, getSurfaceHorizonDepth, pattern="TRUE", hzdesgn = "mollic")

# inspect
bad.sub <- bad.color.v2[11:20,]
plot(bad.sub, color = 'moist_soil_color', name = 'mollic', id.style = 'top')

mols.bad <- site(bad.sub)[, c('peiid', 'surface', 'mollic_depth_new')]
names(mols.bad) <- c("peiid","top","bottom")
addBracket(mols.bad, col = 'red')
axis(side=1, at = 1:length(bad.sub), labels =bad.sub$mollic_depth_new, line = 1)

Which yields this lovely specimen:

[MollicMeyer]

This is gooooood. I totally blanked on the Neutral hue being a potential troublemaker! I will work on this tomorrow.

Regex expressions are also troublesome lately. I wanted to filter the E horizons out of my column indicating a reduced matrix, but keep those that included Eg. That involved....

Idx <- grep(“e[^g]”, hzs(spc)$hzname)

Which worked... (kind of), but I later found out that plain “E” was not filtered with this expression. Only “E1, E2, EBt... etc.,). Any idea why that may be?

Sorry about format, I’m on mobile currently.

[brownag]

@MollicMeyer

Add a * at the end of the expression [zero or more quantifier] to make the "not g" character group / token optional. Regex is looking for any non-g character in that position for a match -- which is why it wont match plain "E"

idx <- grep("E[^g]*", horizons(spc)$hzname)

EDIT: formatting.

Also, this makes me think I should clean up those examples in getSurfaceHorizonDepth -- a couple "work" for the wrong reason.

[brownag]

Significant updates to functions that estimate soil diagnostic feature boundaries have rendered the values previously reported here obsolete. I previously reported that 70% of R02 pedons would see an increase in minimum thickness requirement (under proposed 25cm criterion), when in fact it is only approximately 50% (6d+6c v.s. 6a)

        criterion    Freq     total        prop
1               6a 13790 31573 0.436765591
4               6d 12521 31573 0.396573021
3               6c  3734 31573 0.118265607
2               6b  1476 31573 0.046748804
5 cannot calculate    52 31573 0.001646977

I had originally expected that the sliding scale would be invoked a lot. I was surprised originally at how few diagnostic features were being detected. The 6c (untruncated) usage is even higher still. 65% (two-thirds) of pedons use the sliding scale: ~20% (absolute) of which are truncated to 18cm and ~35% of which are truncated to 25cm.

criterion  Freq  total    prop
3     6c 20340 31573 0.64422133
4     6d  6276 31573 0.19877744
1     6a  3062 31573 0.09698160
2     6b  1895 31573 0.06001964

[brownag]

Quick demo graphic for Curtis Monger -- didnt have anywhere else to put this for now.

# example graphic
par(mar=c(0,0,3,0))
s$mollicthk <- profileApply(s, mollic.thickness.requirement)
s.sub <- union(list(filter(s, mollicthk == 25)[3], 
               filter(s, mollicthk == 18)[3],
               filter(s, mollicthk > 18 & mollicthk < 23)[1]))

s.sub$criteria <- c("25 cm Requirement", "18 cm Requirement", "Between 18 and 25 cm\nRequirement")

plotSPC(s.sub, label = "criteria", axis.line.offset = -3.5, cex.names = 0.9, id.style = "top")
addDiagnosticBracket(s.sub, kind = 'particle size control section', lwd=3)
addDiagnosticBracket(s.sub, kind = 'minimum mollic/umbric thickness', lwd=3, 
                     offset=0, tick.length=0, col="green")
abline(h=c(18+0.5,25+0.5,54+0.5,75+0.5), lty=2)

mtext("aqp calculated\nMollic/Umbric Minimum Thickness Requirement\n& Particle Size Control Section Boundaries", 3)
legend("bottomleft", c("particle size control section",
                       "minimum mollic/umbric thickness",
                       "absolute depths: 18, 25, 54, 75 cm"), lwd=c(3,3,1), lty=c(1,1,2),
       col=c("black","green","blue","black"))

[brownag]