From 5960183aa61cb4e6e20e47fd623b5fec89b755f3 Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Thu, 10 Oct 2024 15:24:59 -0700 Subject: [PATCH] adjustments to criteria, fully implement `'lieutex'` column values and metadata, add permafrost organic soil example --- NEWS.md | 3 +- R/estimatePSCS.R | 111 +++++++++++++++++++---------- man/estimatePSCS.Rd | 59 ++++++++++----- tests/testthat/test-estimatePSCS.R | 4 +- 4 files changed, 120 insertions(+), 57 deletions(-) diff --git a/NEWS.md b/NEWS.md index 1f8ab529..4f4fcef5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,7 +3,8 @@ * `munsell2rgb()` now safely selects the closest Munsell value and chroma to those available in the package LUT * new function `soilTextureColorPal()` for suggesting a color palette suitable for soil texture class * **Breaking Change**: `@sp` slot of the SoilProfileCollection object, and dependency on sp package, has been removed. - * Any SoilProfileCollection objects previously written to file (.rda, .rds) with aqp <2.1.x will need to be rebuilt using `rebuildSPC()` due to changes to S4 object structure + * any SoilProfileCollection objects previously written to file (.rda, .rds) with aqp <2.1.x will need to be rebuilt using `rebuildSPC()` due to changes to S4 object structure + * `estimatePSCS()` gains argument `"lieutex"` for in lieu textures which are used in the new routine for identification of the particle size control section of organic soils # aqp 2.0.4 (2024-07-30) * CRAN release diff --git a/R/estimatePSCS.R b/R/estimatePSCS.R index 2ebb3e92..6c2c2dab 100644 --- a/R/estimatePSCS.R +++ b/R/estimatePSCS.R @@ -1,44 +1,58 @@ #estimatePSCS() -#' Estimate boundaries of the particle size control section (U.S Soil Taxonomy; -#' 12th edition) -#' +#' Estimate boundaries of the U.S Soil Taxonomy Particle Size Control Section +#' #' Estimates the upper and lower boundary of the particle size control section -#' by applying a programmatic version of the particle size control section key -#' from the Keys to Soil Taxonomy (12th edition). +#' for Mineral or Organic soilsby applying a programmatic version of the +#' particle size control section key from the Keys to Soil Taxonomy (13th edition). +#' See details for assumptions and required data elements. +#' +#' @details #' #' Requires information to identify argillic horizons (clay contents, horizon #' designations) with \code{getArgillicBounds()} as well as the presence of #' plow layers and surface organic soil material. Any #' \code{getArgillicBounds()} arguments may be passed to \code{estimatePSCS}. #' -#' Requires information on taxonomic order (to handle Andisols). +#' Also, requires information on taxonomic order to handle Andisols. +#' +#' In aqp 2.1, particle size control sections of organic soils Histosols and +#' Histels are supported. This requires setting the `"lieutex"` horizon metadata +#' column using `hzmetaname<-()` Horizon designations containing `"f"` or `"W"` +#' are recognized as permafrost and water layers, respectively, for application +#' of the organic soils key to control sections. In lieu textures `"SPM"` and +#' `"PEAT"` are used to identify low density organic materials used for surface +#' tier thickness. To avoid using the 160 cm surface tier, set the `"lieutex"` +#' column to any column that does not contain `"SPM"` or `"PEAT"` values. #' #' WARNING: Soils in arenic or grossarenic subgroups, with fragipans, or with #' strongly contrasting PSCs may not be classified correctly. The author would #' welcome a dataset to develop this functionality for. #' -#' In aqp 2.0, control sections of Histosols and Histels are supported. -#' #' @param p A SoilProfileCollection +#' @param hzdesgn Name of the horizon attribute containing the horizon +#' designation. Default `hzdesgnname(p, required = TRUE)` #' @param clay.attr Name of the horizon attribute containing clay contents. -#' Default `"clay"` +#' Default `hzmetaname(p, "clay", required = TRUE)` #' @param texcl.attr Name of the horizon attribute containing textural class -#' (used for finding sandy textures). Default `"texcl"` -#' @param hzdesgn Name of the horizon attribute containing the horizon -#' designation. Default `"hzname"` +#' (used for finding sandy textures). Default `hztexclname(p, required = TRUE)` #' @param tax_order_field Name of the site attribute containing taxonomic #' order; for handling PSCS rules for Andisols in lieu of lab data. May be NA #' or column missing altogether, in which case Andisol PSC possibility is #' ignored. -#' @param lieutex Optional data element used in addition to the horizon designation to identify kinds of organic soil material for soils with organic surfaces. Default: `"lieutex"` +#' @param lieutex Optional data element used in addition to the horizon +#' designation to identify kinds of organic soil material for soils with organic +#' surfaces. Default: `hzmetaname(p, "lieutex")` #' @param bottom.pattern Regular expression pattern to match a root-restrictive #' contact. Default matches Cr, R, Cd or m. This argument is passed to both #' `minDepthOf()` and `getArgillicBounds()`. -#' @param simplify Return a length 2 vector with upper and lower boundary when p has length 1? Default TRUE. +#' @param simplify Return a length 2 vector with upper and lower boundary when +#' p has length 1? Default TRUE. #' @param ... additional arguments are passed to getArgillicBounds() -#' @return A numeric vector (when `simplify=TRUE`) containing the top and bottom depth of the particle +#' @return A numeric vector (when `simplify=TRUE`) containing the top and bottom +#' depth of the particle #' size control section. First value is top, second value is bottom. -#' If `p` contains more than one profile, the result is a data.frame with profile ID plus PSCS top and bottom depths. +#' If `p` contains more than one profile, the result is a data.frame with +#' profile ID plus PSCS top and bottom depths. #' @author Andrew Gene Brown #' @seealso [`getArgillicBounds()`], [`getSurfaceHorizonDepth()`] #' @references Soil Survey Staff. 2014. Keys to Soil Taxonomy, 12th ed. @@ -58,11 +72,27 @@ #' #' x <- estimatePSCS(sp1) #' x +#' +#' # change horizon texture and set inlieu texture column to turn +#' # first profile into an organic soil +#' sp1$name[1:6] <- c("Oi1", "Oi2", "Oi3", "Oaf", "Cf1", "Cf2") +#' sp1$texture <- as.character(sp1$texture) +#' sp1$texture[1:6] <- c("PEAT", "PEAT", "PEAT", "MUCK", "GRVLS", "GRVLS") +#' sp1$bottom[6] <- 200 +#' hzmetaname(sp1, 'lieutex') <- 'texture' +#' +#' y <- estimatePSCS(sp1[1, ], simplify = FALSE) +#' +#' # 74cm lower boundary is 25cm past the upper boundary of permafrost (49cm) +#' # but minimum depth is 100cm unless there is a root-limiting layer +#' y +#' estimatePSCS <- function( p, hzdesgn = hzdesgnname(p, required = TRUE), clay.attr = hzmetaname(p, "clay", required = TRUE), texcl.attr = hztexclname(p, required = TRUE), + lieutex = hzmetaname(p, "lieutex"), tax_order_field = "tax_order", bottom.pattern = 'Cr|R|Cd|m', simplify = TRUE, @@ -101,7 +131,8 @@ estimatePSCS <- function( default_t <- rep(25, length(soildepth)) default_b <- rep(100, length(soildepth)) - odepth <- getMineralSoilSurfaceDepth(p, hzdesgn, simplify = FALSE)[[hz.depths[2]]] + odepth <- getMineralSoilSurfaceDepth(p, hzdesgn = hzdesgn, simplify = FALSE)[[hz.depths[2]]] + othick <- thicknessOf(p, pattern = "O", hzdesgn = hzdesgn)$thickness ohzidx <- which(odepth > 0) # Key part A (soils with restriction in shallow depth) @@ -181,11 +212,12 @@ estimatePSCS <- function( } # handle thick o horizons that cause the mineral PSCS to be greater than soil depth - bad.idx <- which(default_t[ohzidx] >= soildepth[ohzidx] | odepth[ohzidx] > 15) - if (length(bad.idx) > 0) { - organic_cs <- histosol.control.section(p[ohzidx[bad.idx],], hzdesgn = hzdesgn, lieutex = lieutex) - default_t[ohzidx[bad.idx]] <- organic_cs[[2]] - default_b[ohzidx[bad.idx]] <- organic_cs[[3]] + organic.idx <- which(default_t >= soildepth | odepth >= 40 | + (othick >= (soildepth * 2 / 3) & (soildepth - othick) <= 10)) + if (length(organic.idx) > 0) { + organic_cs <- histosol.control.section(p[organic.idx,], hzdesgn = hzdesgn, lieutex = lieutex) + default_t[organic.idx] <- organic_cs[[2]] + default_b[organic.idx] <- organic_cs[[3]] } # Adjust PSCS bottom depth to restriction depth, if appropriate @@ -204,18 +236,19 @@ estimatePSCS <- function( } # calculate control section based on materials in surface tier, root limiting layers, permafrost and water -histosol.control.section <- function(p, hzdesgn = hzdesgnname(p), lieutex = "lieutex") { - +histosol.control.section <- function(p, + hzdesgn = hzdesgnname(p, required = TRUE), + lieutex = hzmetaname(p, "lieutex", required = TRUE)) { depthz <- horizonDepths(p) # control section stops at Cr|R|Cd - rll.depth <- minDepthOf(p, pattern = "Cr|R|Cd|m", no.contact.depth = 0, no.contact.assigned = 200, simplify = FALSE)[[depthz[1]]] + rll.depth <- minDepthOf(p, pattern = "Cr|R|Cd|m", hzdesgn = hzdesgn, no.contact.depth = 0, no.contact.assigned = 200, simplify = FALSE)[[depthz[1]]] - # control section stops at permafrost + 25cm - permafrost.depth <- minDepthOf(p, "f|ff", no.contact.assigned = Inf, simplify = FALSE)[[depthz[1]]] + # control section stops at permafrost + 25cm (must be at least 100cm) + permafrost.depth <- minDepthOf(p, "f|ff", hzdesgn = hzdesgn, no.contact.assigned = Inf, simplify = FALSE)[[depthz[1]]] # control section stops at water bottom depth greater than either 130 or 160 cm, as applicable - water.depth <- minDepthOf(p, "W", no.contact.assigned = Inf, simplify = FALSE)[[depthz[1]]] + water.depth <- minDepthOf(p, "W", hzdesgn = hzdesgn, no.contact.assigned = Inf, simplify = FALSE)[[depthz[1]]] # get surface tier thickness for p hst <- histosol.surface.tier(p, hzdesgn = hzdesgn, lieutex = lieutex, rll = rll.depth)$hst_b @@ -223,16 +256,15 @@ histosol.control.section <- function(p, hzdesgn = hzdesgnname(p), lieutex = "lie data.frame( peiid = profile_id(p), hcs_t = 0, - hcs_b = pmin(rll.depth, permafrost.depth + 25, pmin(water.depth, ifelse(hst == 30, 130, 160))) + hcs_b = pmin(rll.depth, pmax(permafrost.depth + 25, 100), pmin(water.depth, ifelse(hst == 30, 130, 160))) ) } -histosol.surface.tier <- - function(p, - hzdesgn = hzdesgnname(p), - lieutex = "lieutex", - rll = NULL) { +histosol.surface.tier <- function(p, + hzdesgn = hzdesgnname(p, required = TRUE), + lieutex = hzmetaname(p, "lieutex", required = TRUE), + rll = NULL) { .SD <- NULL; .LAST <- NULL; .BOTTOM <- NULL; .TOP <- NULL; .thk <- NULL depthz <- horizonDepths(p) @@ -245,10 +277,13 @@ histosol.surface.tier <- # column "lieutex" used in lieu of texture, and in lieu of horizon designations (Oa, Oe, Oi) if (!lieutex %in% names(p)) { - lieutex <- guessHzAttrName(p, "lieu", c("tex", "in"), required = FALSE) + # this will only trigger if one or pedons appears to meet the organic layer thickness requirements + stop("The in lieu texture ('lieutex') is required to apply rules for organic soils. + Set the 'lieutex' horizon designation metadata with `hzmetaname(p, 'lieutex') <- 'lieutex'` + Soils with low density fibric materials are indicated using lieutex column values `'PEAT'` or `'SPM'`.") } - mss <- getMineralSoilSurfaceDepth(p, simplify = FALSE) + mss <- getMineralSoilSurfaceDepth(p, hzdesgn = hzdesgn, simplify = FALSE) # if there is a mineral surface, no surface tier noo.idx <- which(mss[[depthz[2]]] == 0) @@ -259,7 +294,7 @@ histosol.surface.tier <- } if (is.null(rll)) { - rll <- minDepthOf(psub, pattern = "Cr|Cd|R|m", no.contact.depth = 0, no.contact.assigned = 200, simplify = FALSE)[[depthz[1]]] + rll <- minDepthOf(psub, hzdesgn = hzdesgn, pattern = "Cr|Cd|R|m", no.contact.depth = 0, no.contact.assigned = 200, simplify = FALSE)[[depthz[1]]] } p60 <- glom(psub, 0, pmin(rll, 60), truncate = TRUE, drop = FALSE) @@ -268,7 +303,7 @@ histosol.surface.tier <- # either lieutex of peat or "i" suffix throughout all horizons in upper 60 hz <- data.table::data.table(horizons(p60)) - idx <- which(grepl("i", hz[[hzdesgn]]) | hz[[lieutex]] %in% c("peat", "spm")) + idx <- which(grepl("i", hz[[hzdesgn]]) | tolower(hz[[lieutex]]) %in% c("peat", "spm")) hz$.thk <- 0 hz$.thk[idx] <- (hz[[depthz[2]]] - hz[[depthz[1]]])[idx] diff --git a/man/estimatePSCS.Rd b/man/estimatePSCS.Rd index ba3dfe3c..d6479878 100644 --- a/man/estimatePSCS.Rd +++ b/man/estimatePSCS.Rd @@ -2,16 +2,15 @@ % Please edit documentation in R/estimatePSCS.R \name{estimatePSCS} \alias{estimatePSCS} -\title{Estimate boundaries of the particle size control section (U.S Soil Taxonomy; -12th edition)} +\title{Estimate boundaries of the U.S Soil Taxonomy Particle Size Control Section} \usage{ estimatePSCS( p, hzdesgn = hzdesgnname(p, required = TRUE), clay.attr = hzmetaname(p, "clay", required = TRUE), texcl.attr = hztexclname(p, required = TRUE), + lieutex = hzmetaname(p, "lieutex"), tax_order_field = "tax_order", - lieutex = "lieutex", bottom.pattern = "Cr|R|Cd|m", simplify = TRUE, ... @@ -21,38 +20,44 @@ estimatePSCS( \item{p}{A SoilProfileCollection} \item{hzdesgn}{Name of the horizon attribute containing the horizon -designation. Default \code{"hzname"}} +designation. Default \code{hzdesgnname(p, required = TRUE)}} \item{clay.attr}{Name of the horizon attribute containing clay contents. -Default \code{"clay"}} +Default \code{hzmetaname(p, "clay", required = TRUE)}} \item{texcl.attr}{Name of the horizon attribute containing textural class -(used for finding sandy textures). Default \code{"texcl"}} +(used for finding sandy textures). Default \code{hztexclname(p, required = TRUE)}} + +\item{lieutex}{Optional data element used in addition to the horizon +designation to identify kinds of organic soil material for soils with organic +surfaces. Default: \code{hzmetaname(p, "lieutex")}} \item{tax_order_field}{Name of the site attribute containing taxonomic order; for handling PSCS rules for Andisols in lieu of lab data. May be NA or column missing altogether, in which case Andisol PSC possibility is ignored.} -\item{lieutex}{Optional data element used in addition to the horizon designation to identify kinds of organic soil material for soils with organic surfaces. Default: \code{"lieutex"}} - \item{bottom.pattern}{Regular expression pattern to match a root-restrictive contact. Default matches Cr, R, Cd or m. This argument is passed to both \code{minDepthOf()} and \code{getArgillicBounds()}.} -\item{simplify}{Return a length 2 vector with upper and lower boundary when p has length 1? Default TRUE.} +\item{simplify}{Return a length 2 vector with upper and lower boundary when +p has length 1? Default TRUE.} \item{...}{additional arguments are passed to getArgillicBounds()} } \value{ -A numeric vector (when \code{simplify=TRUE}) containing the top and bottom depth of the particle +A numeric vector (when \code{simplify=TRUE}) containing the top and bottom +depth of the particle size control section. First value is top, second value is bottom. -If \code{p} contains more than one profile, the result is a data.frame with profile ID plus PSCS top and bottom depths. +If \code{p} contains more than one profile, the result is a data.frame with +profile ID plus PSCS top and bottom depths. } \description{ Estimates the upper and lower boundary of the particle size control section -by applying a programmatic version of the particle size control section key -from the Keys to Soil Taxonomy (12th edition). +for Mineral or Organic soilsby applying a programmatic version of the +particle size control section key from the Keys to Soil Taxonomy (13th edition). +See details for assumptions and required data elements. } \details{ Requires information to identify argillic horizons (clay contents, horizon @@ -60,13 +65,20 @@ designations) with \code{getArgillicBounds()} as well as the presence of plow layers and surface organic soil material. Any \code{getArgillicBounds()} arguments may be passed to \code{estimatePSCS}. -Requires information on taxonomic order (to handle Andisols). +Also, requires information on taxonomic order to handle Andisols. + +In aqp 2.1, particle size control sections of organic soils Histosols and +Histels are supported. This requires setting the \code{"lieutex"} horizon metadata +column using \verb{hzmetaname<-()} Horizon designations containing \code{"f"} or \code{"W"} +are recognized as permafrost and water layers, respectively, for application +of the organic soils key to control sections. In lieu textures \code{"SPM"} and +\code{"PEAT"} are used to identify low density organic materials used for surface +tier thickness. To avoid using the 160 cm surface tier, set the \code{"lieutex"} +column to any column that does not contain \code{"SPM"} or \code{"PEAT"} values. WARNING: Soils in arenic or grossarenic subgroups, with fragipans, or with strongly contrasting PSCs may not be classified correctly. The author would welcome a dataset to develop this functionality for. - -In aqp 2.0, control sections of Histosols and Histels are supported. } \examples{ @@ -81,6 +93,21 @@ hzmetaname(sp1, 'clay') <- 'prop' x <- estimatePSCS(sp1) x + +# change horizon texture and set inlieu texture column to turn +# first profile into an organic soil +sp1$name[1:6] <- c("Oi1", "Oi2", "Oi3", "Oaf", "Cf1", "Cf2") +sp1$texture <- as.character(sp1$texture) +sp1$texture[1:6] <- c("PEAT", "PEAT", "PEAT", "MUCK", "GRVLS", "GRVLS") +sp1$bottom[6] <- 200 +hzmetaname(sp1, 'lieutex') <- 'texture' + +y <- estimatePSCS(sp1[1, ], simplify = FALSE) + +# 74cm lower boundary is 25cm past the upper boundary of permafrost (49cm) +# but minimum depth is 100cm unless there is a root-limiting layer +y + } \references{ Soil Survey Staff. 2014. Keys to Soil Taxonomy, 12th ed. diff --git a/tests/testthat/test-estimatePSCS.R b/tests/testthat/test-estimatePSCS.R index f48561c5..1715b7fe 100644 --- a/tests/testthat/test-estimatePSCS.R +++ b/tests/testthat/test-estimatePSCS.R @@ -93,9 +93,9 @@ test_that("estimatePSCS() multiple profiles",{ test_that("estimatePSCS() organic profiles",{ sp1@horizons$name[1:60] <- "Oi" sp1@horizons$name[48:51] <- "W" - sp1@horizons$name[57:60] <- "Oiff" + sp1@horizons$name[57:60] <- "Oif" horizons(sp1)$lieutex <- NA - oo <- estimatePSCS(sp1, clay.attr = 'prop', texcl.attr = "texture") + oo <- estimatePSCS(sp1, hzdesgn = "name", clay.attr = "prop", texcl.attr = "texture", lieutex = "texture") expect_equal(oo$pscs_top, rep(0, 9)) expect_equal(oo$pscs_bottom, c(89, 59, 67, 62, 68, 160, 160, 68, 133)) })