Nothing
R/texture.R
In aqp: Algorithms for Quantitative Pedology
Defines functions
texcl_to_ssc
Documented in
texcl_to_ssc
#' Textural conversions
#'
#' These functions consist of several conversions between sand, silt and clay
#' to texture class and visa versa, textural modifiers to rock fragments, and
#' grain size composition to the family particle size class.
#'
#' These functions are intended to estimate missing values or allocate particle
#' size fractions to classes. The \code{ssc_to_texcl()} function uses the same
#' logic as the particle size estimator calculation in NASIS to classify sand
#' and clay into texture class. The results are stored in \code{soiltexture}
#' and used by \code{texcl_to_ssc()} as a lookup table to convert texture class
#' to sand, silt and clay. The function \code{texcl_to_ssc()} replicates the
#' functionality described by Levi (2017). The \code{texmod_to_fragvol()}
#' function similarly uses the logical from the
#' Exhibit618-11_texture_modifier.xls spreadsheet to determine the textural
#' modifier from the various combinations of rock and pararock fragments
#' (e.g. GR and PGR).
#'
#' When \code{sample = TRUE}, the results can be used to estimate within-class,
#' marginal distributions of sand, silt, and clay fractions. It is recommended
#' that at least 10 samples be drawn for reasonable estimates.
#'
#' The function \code{texmod_to_fragvoltot} returns a data.frame with multiple
#' fragvoltot columns differentiated by tailing abbreviations (e.g. _r) which
#' refer to the following:
#' 1. l = low
#' 2. r = representative
#' 3. h = high
#' 4. nopf = no pararock fragments (i.e. total fragments - pararock fragments)
#'
#' The function \code{texture_to_texmod()} parses texture (e.g. GR-CL) to extract the texmod values from it in the scenario where it is missing from
#' texmod column. If multiple texmod values are present (for example in the
#' case of stratified textures) and \code{duplicates = "combine"} they will be combined in the output (e.g. GR & CBV). Otherwise if \code{duplicates = "max"}
#' the texmod with the highest rock fragment (e.g. CBV) will be returned.
#'
#' Unlike the other functions, \code{texture_to_taxpartsize()} is intended to
#' be computed on weighted averages within the family particle size control
#' section. Also recall from the criteria that carbonate clay should be
#' subtracted from clay content and added to silt content. Similarly, if the
#' percent of very fine sand is known it should be subtracted from the sand,
#' and added to the silt content. Unlike the other functions,
#' \code{texture_to_taxpartsize()} is intended to be computed on weighted
#' averages within the family particle size control section. Also recall from
#' the criteria that carbonate clay should be subtracted from clay content and
#' added to silt content. Similarly, if the percent of very fine sand is known
#' it should be subtracted from the sand, and added to the silt content.
#' @param texcl vector of texture classes than conform to the USDA code
#' conventions (e.g. c|C, sil|SIL, sl|SL, cos|COS)
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#' @param lieutex vector of in lieu of texture terms that conform to the USDA
#' code convenctions (e.g. gr|GR, pg|PG), only used when fragments or artifacts
#' are > 90 percent by volume (default: NULL))
#' @param texture vector of combinations of texcl, texmod, and lieutex (e.g. CL, GR-CL, CBV-S, GR)
#'
#' @param clay vector of clay percentages
#' @param sand vector of sand percentages
#' @param fragvoltot vector of rock fragment percentages
#'
#' @param as.is logical: should character vectors be converted to factors?
#' (default: TRUE)
#'
#' @param droplevels logical: indicating whether to drop unused levels in
#' factors. This is useful when the results have a large number of unused
#' classes, which can waste space in tables and figures.
#' @param sample logical: should ssc be random sampled from the lookup table?
#' (default: FALSE)
#'
#' @return - `texcl_to_ssc`: A `data.frame` containing columns `"sand"`,`"silt"`, `"clay"`
#'
#' @seealso \code{\link{SoilTextureLevels}}
#'
#' @author Stephen Roecker
#'
#' @references Matthew R. Levi, Modified Centroid for Estimating Sand, Silt, and Clay from Soil Texture Class, Soil Science Society of America Journal, 2017, 81(3):578-588, ISSN 1435-0661, \doi{10.2136/sssaj2016.09.0301}.
#'
#' @rdname texture
#' @keywords manip
#' @export
#' @examples
#' \donttest{
#' # example of ssc_to_texcl()
#' tex <- expand.grid(sand = 0:100, clay = 0:100)
#' tex <- subset(tex, (sand + clay) < 101)
#' tex$texcl <- ssc_to_texcl(sand = tex$sand, clay = tex$clay)
#' head(tex)
#'
#' # example of texcl_to_ssc(texcl)
#' texcl <- c("cos", "s", "fs", "vfs", "lcos", "ls",
#' "lfs", "lvfs", "cosl", "sl", "fsl", "vfsl", "l",
#' "sil", "si", "scl", "cl", "sicl", "sc", "sic", "c"
#' )
#' test <- texcl_to_ssc(texcl)
#' head(test <- cbind(texcl, test), 10)
#'
#'
#' # example of texcl_to_ssc(texcl, clay)
#' data(soiltexture)
#' st <- soiltexture$values
#' idx <- sample(1:length(st$texcl), 10)
#' st <- st[idx, ]
#' ssc <- texcl_to_ssc(texcl = st$texcl)
#' head(cbind(texcl = st$texcl, clay = ssc$clay))
#'
#'
#' # example of texmod_to_fragvoltol
#' frags <- c("gr", "grv", "grx", "pgr", "pgrv", "pgrx")
#' head(texmod_to_fragvoltot(frags))
#'
#'
#' # example of texture_to_taxpartsize()
#' tex <- data.frame(texcl = c("c", "cl", "l", "ls", "s"),
#' clay = c(55, 33, 18, 6, 3),
#' sand = c(20, 33, 42, 82, 93),
#' fragvoltot = c(35, 15, 34, 60, 91))
#' tex$fpsc <- texture_to_taxpartsize(texcl = tex$texcl,
#' clay = tex$clay,
#' sand = tex$sand,
#' fragvoltot = tex$fragvoltot)
#' head(tex)
#'
#'
#' # example of texture_to_taxpartsize() with carbonate clay and very fine sand
#' carbclay <- rnorm(5, 2, 3)
#' vfs <- rnorm(5, 10, 3)
#' st$fpsc <- texture_to_taxpartsize(texcl = tex$texcl,
#' clay = tex$clay - carbclay,
#' sand = tex$sand - vfs,
#' fragvoltot = tex$fragvoltot)
#' head(tex)
#'
#'
#' # example of sample = TRUE
#' texcl <- rep(c("cl", "sil", "sl"), 10)
#' ssc1 <- cbind(texcl, texcl_to_ssc(texcl = texcl, sample = FALSE))
#' ssc2 <- cbind(texcl, texcl_to_ssc(texcl = texcl, sample = TRUE))
#' ssc1$sample <- FALSE
#' ssc2$sample <- TRUE
#' ssc <- rbind(ssc1, ssc2)
#' aggregate(clay ~ sample + texcl, data = ssc, summary)
#' }
texcl_to_ssc <- function(texcl = NULL, clay = NULL, sample = FALSE) {
# fix for R CMD check
# texcl_to_ssc: no visible binding for global variable ‘soiltexture’
soiltexture <- NULL
# clay is not NULL
clay_not_null <- all(!is.null(clay))
clay_is_null <- !clay_not_null
# standardize the inputs
df <- data.frame(texcl = tolower(as.character(texcl)),
stringsAsFactors = FALSE
)
if (clay_not_null) {
df$clay <- as.integer(round(clay))
}
df$rn <- row.names(df)
load(system.file("data/soiltexture.rda", package="aqp")[1])
# convert fine sand classes to their generic counterparts
df <- within(df, {
texcl = ifelse(texcl %in% c("cos", "fs", "vfs"), "s", texcl)
texcl = ifelse(texcl %in% c("lcos", "lfs", "lvfs"), "ls", texcl)
texcl = ifelse(texcl %in% c("cosl", "fsl", "vfsl"), "sl", texcl)
})
# check for texcl that don't match
idx <- ! df$texcl %in% unique(soiltexture$averages$texcl)
if (any(idx)) {
warning("not all the user supplied texcl values match the lookup table")
}
# check clay values within texcl
if (clay_not_null) {
clay_not_na <- !is.na(df$clay)
idx <- paste(df$texcl[clay_not_na], df$clay[clay_not_na]) %in% paste(soiltexture$values$texcl, soiltexture$values$clay)
if (any(!idx)) {
warning("not all the user supplied clay values fall within the texcl, so they will be set to NA")
df$clay[which(!idx)] <- NA
clay_not_null <- all(!is.na(df$clay))
clay_is_null <- !clay_not_null
}
}
# check clay ranges 0-100
idx <- clay_not_null & any(clay < 0, na.rm = TRUE) & any(clay > 100, na.rm = TRUE)
if (idx) {
warning("some clay records < 0 or > 100%")
}
# if clay is present
if (clay_not_null & sample == FALSE) {
st <- aggregate(sand ~ texcl + clay, data = soiltexture$values, function(x) as.integer(round(mean(x))))
st$silt <- 100 - st$clay - st$sand
# some clay present (compute clay weighted averages)
idx <- is.na(df$clay)
if (any(idx)) {
df_na <- merge(df[idx, c("texcl", "rn")], soiltexture$averages, by = "texcl", all.x = TRUE, sort = FALSE)[c("texcl", "clay", "rn")]
df[idx, ] <- df_na
}
df <- merge(df[c("texcl", "clay", "rn")], st, by = c("texcl", "clay"), all.x = TRUE, sort = FALSE)
} else {
# clay missing (use average)
df <- merge(df[c("texcl", "rn")], soiltexture$averages, by = "texcl", all.x = TRUE, sort = FALSE)
}
# randomly sample ssc from texcl lookup table
if (sample == TRUE) {
if (clay_is_null) df$clay <- NA
split(df, df$texcl, drop = TRUE) ->.;
lapply(., function(x) {
st <- soiltexture$values
# clay present
x$idx <- is.na(x$clay)
x1 <- x[x$idx == FALSE, ]
x2 <- x[x$idx == TRUE, ]
if (clay_not_null) {
temp1 <- st[st$texcl == x1$texcl[1] &
st$clay %in% unique(x1$clay),
]
temp1 <- temp1[sample(1:nrow(temp1), size = nrow(x1), replace = TRUE), ]
temp1 <- cbind(x1[x1$idx == FALSE, c("rn", "texcl")], temp1[c("clay", "sand")])
} else temp1 <- NULL
# clay missing
temp2 <- st[st$texcl == x2$texcl[1], ]
temp2 <- temp2[sample(1:nrow(temp2), size = nrow(x2), replace = TRUE), ]
temp2 <- cbind(x2[x2$idx == TRUE, c("rn", "texcl")], temp2[c("clay", "sand")])
return(rbind(temp1, temp2))
}) ->.;
do.call("rbind", .) -> df
df$silt <- 100 - df$clay - df$sand
}
# standardize outputs
vars <- c("sand", "silt", "clay")
df <- df[(order(as.integer(df$rn))), vars]
df$rn <- NULL
df$texcl <- NULL
rownames(df) <- NULL
return(df)
}
#' Convert sand, silt and clay to texture class
#' @param simplify Passed to `SoilTextureLevels()` to set the number of possible texture classes. If `TRUE`, the ordered factor has a maximum of 12 levels, if `FALSE` (default) the ordered factor has a maximum of 21 levels (including e.g. very fine/fine/coarse variants)
#' @rdname texture
#' @return - `ssc_to_texcl`: A `character` vector containing texture class
#' @export
#'
ssc_to_texcl <- function(sand = NULL, clay = NULL, simplify = FALSE, as.is = FALSE, droplevels = TRUE) {
# fix for R CMD check:
# ssc_to_texcl: no visible binding for global variable ‘silt’
silt <- NULL
# check lengths
idx <- length(clay) != length(sand)
if (idx) {
stop("length of inputs do not match")
}
# standardize inputs
df <- data.frame(sand = as.integer(round(sand)),
clay = as.integer(round(clay)),
stringsAsFactors = FALSE
)
df$silt <- 100 - df$clay - df$sand
## TODO: this needs some more work: sum will always be 100, but silt-by-difference may be illogical
# # check sand, silt and clay sum to 100
# idx <- (df$sand + df$silt + df$clay) > 100 | (df$sand + df$silt + df$clay) < 100
# if (any(idx) & any(complete.cases(df[c("sand", "clay")]))) {
# warning("some records sand, silt, and clay do not sum to 100 %")
# }
# logic from the particle size estimator calculation from NASIS
df <- within(df, {
texcl = NA
texcl[silt >= 79.99 & clay < 11.99] = "si"
texcl[silt >= 49.99 & clay < 26.99 & (silt < 79.99 | clay >= 11.99)] = "sil"
texcl[clay >= 26.99 & clay < 39.99 & sand <= 20.01] = "sicl"
texcl[clay >= 39.99 & silt >= 39.99] = "sic"
texcl[clay >= 39.99 & sand <= 45.01 & silt < 39.99] = "c"
texcl[clay >= 26.99 & clay < 39.99 & sand > 20.01 & sand <= 45.01] = "cl"
texcl[clay >= 6.99 & clay < 26.99 & silt >= 27.99 & silt < 49.99 & sand <= 52.01] = "l"
texcl[clay >= 19.99 & clay < 34.99 & silt < 27.99 & sand > 45.01] = "scl"
texcl[clay >= 34.99 & sand > 45.01] = "sc"
texcl[(silt + 1.5 * clay) < 15] = "s"
texcl[(silt + 1.5 * clay) >= 15 & (silt + 2 * clay) < 29.99] = "ls"
texcl[!is.na(sand) & !is.na(clay) & is.na(texcl)] = "sl"
})
# encoding according to approximate AWC, from Red Book version 3.0
if (!as.is) {
df$texcl <- factor(df$texcl, levels = SoilTextureLevels(which = 'codes', simplify = simplify), ordered = TRUE)
if (droplevels) {
df$texcl <- droplevels(df$texcl)
}
}
return(df$texcl)
}
#' Convert from fragment modifier to estimated total volume
#'
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#'
#' @param lieutex vector of in lieu of texture terms that conform to the USDA
#' code conventions (e.g. gr|GR, pg|PG), only used when fragments or artifacts
#' are > 90 percent by volume (default: NULL))
#' @rdname texture
#'
#' @return - `texmod_to_fragvoltot`: A `data.frame` containing columns `"fragvoltot_l"`,
#' `"fragvoltot_r"`, `"fragvoltot_h"`, `"fragvoltot_l_nopf"`, `"fragvoltot_r_nopf"`, `"fragvoltot_h_nopf"`
#'
#' @export
texmod_to_fragvoltot <- function(texmod = NULL, lieutex = NULL) {
# standardize inputs ----
vars_l <- list(texmod = texmod, lieutex = lieutex)
# vars_l <- list(texmod = c("GR", "CBV", NA), lieutex = c(NA, NA, "BY"))
df <- NULL
df <- .format_inputs(vars_l, names(vars_l), NA_character_, "character")
df$texmod <- tolower(df$texmod)
df$lieutex <- toupper(df$lieutex)
df$rn <- 1:nrow(df)
# load lookup table
# fix for R CMD check
# texcl_to_ssc: no visible binding for global variable ‘soiltexture’
soiltexture <- NULL
load(system.file("data/soiltexture.rda", package="aqp")[1])
# check inputs ----
## both present?
idx <- any(complete.cases(df))
if (idx) {
warning("texmod and lieutex should not both be present, they are mutually exclusive, only the results for texmod will be returned")
}
# texmod and lieutex don't match
idx <- ! df$texmod %in% soiltexture$texmod$texmod
if (any(idx & !is.na(df$texmod))) {
message("not all the texmod supplied match the lookup table")
}
idx <- ! df$lieutex %in% c("GR", "CB", "ST", "BY", "CN", "FL", "PG", "PCB", "PST", "PBY", "PCN", "PFL", "BR", "HMM", "MPM", "SPM", "MUCK", "PEAT", "ART", "CGM", "FGM", "ICE", "MAT", "W")
if (any(idx & !is.na(df$lieutex))) {
message("not all the lieutex supplied match the lookup table")
}
# merge
df <- merge(df, soiltexture$texmod, by = "texmod", all.x = TRUE, sort = FALSE)
df <- df[order(df$rn), ]
df$rn <- NULL
# lieutex
if (any(!is.na(df$lieutex))) {
idx1 <- !is.na(df$lieutex) & df$lieutex %in% c("GR", "CB", "ST", "BY", "CN", "FL")
idx2 <- !is.na(df$lieutex) & df$lieutex %in% c("PG", "PCB", "PST", "PBY", "PCN", "PFL")
# df$lieutex <- toupper(lieutex)
df <- within(df, {
fragvoltot_l = ifelse(idx1, 90, fragvoltot_l)
fragvoltot_r = ifelse(idx1, 95, fragvoltot_l)
fragvoltot_h = ifelse(idx1, 100, fragvoltot_l)
fragvoltot_l_nopf = ifelse(idx2, 0, fragvoltot_l)
fragvoltot_r_nopf = ifelse(idx2, 0, fragvoltot_l)
fragvoltot_h_nopf = ifelse(idx2, 0, fragvoltot_l)
})
# df$lieutex <- lieutex
}
df[c("texmod", "lieutex", "texmod_label")] <- NULL
return(df)
}
#' Convert sand, silt and clay to the family particle size class
#'
#' @param texcl vector of texture classes than conform to the USDA code
#' conventions (e.g. c|C, sil|SIL, sl|SL, cos|COS)
#'
#' @param clay vector of clay percentages
#' @param sand vector of sand percentages
#' @param sandvf vector of very fine sand percentages
#'
#' @param fragvoltot vector of total rock fragment percentages
#'
#' @return - `texture_to_taxpartsize`: a character vector containing `"taxpartsize"` classes
#'
#' @seealso [hz_to_taxpartsize()], [lookup_taxpartsize()]
#'
#' @rdname texture
#'
#' @export
#'
texture_to_taxpartsize <- function(texcl = NULL, clay = NULL, sand = NULL, sandvf = NULL, fragvoltot = NULL) {
# check lengths
idx <- length(texcl) == length(clay) & length(clay) == length(sand) & length(sand) == length(fragvoltot)
if (!idx) {
stop("length of inputs do not match")
}
# standarize inputs
if (is.null(sandvf)) sandvf <- NA
df <- data.frame(texcl = tolower(texcl),
clay = as.integer(round(clay)),
sand = as.integer(round(sand)),
sandvf = as.integer(round(sandvf)),
fragvoltot = as.integer(round(fragvoltot)),
fpsc = as.character(NA),
stringsAsFactors = FALSE
)
df$silt <- 100 - df$sand - df$clay
sandytextures <- c("cos", "s", "fs", "lcos", "ls", "lfs")
# check texcl lookup
idx <- any(! df$texcl[!is.na(df$texcl)] %in% SoilTextureLevels(which = 'codes'))
if (idx) {
warning("not all the texcl supplied match the lookup table")
}
# check percentages
idx <- df$silt > 100 | df$silt < 0 | df$clay > 100 | df$clay < 0 | df$sand > 100 | df$sand < 0 | df$fragvoltot > 100 | df$fragvoltot < 0
if (any(idx, na.rm = TRUE)) {
warning("some records are > 100% or < 0%, or the calcuated silt fraction is > 100% or < 0%")
}
if (any(sandvf > sand & all(!is.na(sandvf)))) {
warning("the sandvf values should not be greater than the sand values")
}
# check ssc_to_texcl() vs texcl
df$texcl_calc <- suppressMessages(ssc_to_texcl(sand = df$sand, clay = df$clay, as.is = TRUE))
df <- within(df, {
texcl_calc = ifelse(texcl_calc == "s" & grepl("^cos$|^fs$|^vfs$", texcl), texcl, texcl_calc)
texcl_calc = ifelse(texcl_calc == "ls" & grepl("^lcos$|^lfs$|^lvfs$", texcl), texcl, texcl_calc)
texcl_calc = ifelse(texcl_calc == "sl" & grepl("^cosl$|^fsl$|^vfsl$", texcl), texcl, texcl_calc)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("vfs", "lvfs"), 50, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("vfsl"), 40, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("fsl"), 15, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("sl", "l", "scl"), 10, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("fsl"), 7, sandvf)
sand = ifelse(!is.na(sandvf), sand - sandvf, sand)
silt = ifelse(!is.na(sandvf), silt + sandvf, silt)
})
idx <- any(df$texcl != df$texcl_calc, na.rm = TRUE)
if (idx) {
warning("some of the texcl records don't match the calculated texcl via ssc_to_texcl()")
}
# calculate family particle size control section
idx <- df$fragvoltot >= 35
if (any(idx)) {
df[idx,] <- within(df[idx,], {
fpsc[clay >= 35] = "clayey-skeletal"
fpsc[clay < 35] = "loamy-skeletal"
fpsc[texcl %in% sandytextures] = "sandy-skeletal"
})
}
idx <- df$fragvoltot < 35 & df$texcl %in% sandytextures
if (any(idx)) {
df[idx, ]$fpsc <- "sandy"
}
idx <- df$fragvoltot < 35 & ! df$texcl %in% sandytextures
if (any(idx)) {
df[idx, ] <- within(df[idx,], {
fpsc[clay < 18 & sand >= 15] = "coarse-loamy"
fpsc[clay < 18 & sand < 15] = "coarse-silty"
fpsc[clay >= 18 & clay < 35] = "fine-loamy"
fpsc[clay >= 18 & clay < 35 & sand < 15] = "fine-silty"
fpsc[clay >= 35 & clay < 60] = "fine"
fpsc[clay > 60] = "very-fine"
})
}
df$fpsc <- ifelse(df$fragvoltot > 90, "fragmental", df$fpsc)
return(df$fpsc)
}
#' Parse texmod from texture
#'
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#' @param texture vector of combinations of texcl, texmod, and lieutex (e.g. CL, GR-CL, CBV-S, GR)
#'
#' @param duplicates character: specifying how multiple values should be handled, options are `"combined"` (e.g. 'GR & GRV) or `"max"`(e.g. 'GRV')
#'
#' @return - `texture_to_texmod`: a character vector containing `"texmod"` classes
#'
#' @rdname texture
#'
#' @export
#'
#' @examples
#' \donttest{
#' # example of texture_to_texmod()
#' tex <- c("SL", "GR-SL", "CBV-L", "SR- GR-FS GRX-COS")
#' texture_to_texmod(tex)
#' texture_to_texmod(tex, duplicates = "max")
#' }
texture_to_texmod <- function(texture, duplicates = "combine") {
# test
is_char <- is.character(texture)
is_fact <- is.factor(texture)
if (!is_char & !is_fact) stop("texture must be a character or a factor")
if (is_fact) texture <- as.character(texture)
# load lookup tables
soiltexture <- NULL
load(system.file("data/soiltexture.rda", package = "aqp")[1])
texmod_df <- soiltexture$texmod
texmod <- texmod_df$texmod[grepl("gravel|cobbl|ston|boulder|channer|flag", texmod_df$texmod_label)]
# lieutex <- metadata[metadata$ColumnPhysicalName == "lieutex", ]$ChoiceName
tex <- tolower(texture)
tex_l <- strsplit(tex, "-| |_|\\/|&|,|;|\\.|\\*|\\+")
# find which records have a texmod
idx_mod <- which(sapply(tex_l, function(x) any(x[-length(x)] %in% texmod)))
# idx_lieu <- which(sapply(tex_l, function(x) any(x %in% lieutex)))
# create a logical matrix of texmod presence/absence
if (length(idx_mod > 0)) {
tex_m <- lapply(tex_l[idx_mod], function(x)
texmod %in% x)
tex_df <- as.data.frame(do.call("rbind", tex_m))
names(tex_df) <- texmod
# slightly slower simpler(/better?) alternative
# not better, this mistakenly matches words entered into this field
# I would rather have false negatives than false positives
# tex_m <- sapply(texmod, function(x) grepl(x, tex))
# duplicates ----
## combine texmod (if multiple exist)
if (duplicates == "combine") {
texmod_parse <- apply(tex_df, 1, function(x) {
paste0(names(x)[which(x)], collapse = " & ")
})
}
## find texmod with the highest rock fragment content (if multiple exist)
if (duplicates == "max") {
tex_m <- sapply(texmod, function(x) {
vals <- suppressMessages(texmod_to_fragvoltot(x))
rf <- ifelse(tex_df[x] == TRUE, vals$fragvoltot_r, 0)
})
idx_col <- max.col(tex_m)
texmod_parse <- names(tex_df)[idx_col]
}
}
# results ----
n <- length(texture)
df <- data.frame(texmod = rep(NA_character_, n)) #, lieutex = rep(NA, n))
if (length(idx_mod) > 0) df[idx_mod, "texmod"] <- texmod_parse
return(df$texmod)
}
# add feature to uncode
# texture_to_texcl <- function(texclLabel, droplevels = FALSE
# # , stringsAsFactors = FALSE
# ) {
#
# # standardize inputs
# texclLabel <- tolower(texclLabel)
#
# # load lookup table
# load(system.file("data/soiltexture.rda", package = "aqp")[1])
# sub <- soiltexture$averages
#
# # check
# idx <- ! texclLabel %in% sub$texcl_label
# if (any(idx)) {
# warning("not all the user supplied texture values match the lookup table")
# }
#
# # convert
# texcl <- factor(texclLabel, levels = tolower(sub$texcl_label), label = sub$texcl, ordered = TRUE)
#
# if (droplevels == TRUE) {
# texcl <- droplevels(texcl)
# }
#
# # if (stringsAsFactors == TRUE) {
# # texcl <- as.character(texcl)
# # }
#
# return(texcl)
# }
#' Convert frags to texmod
#'
#' @param object data.frame: containing the following column names: gravel, cobbles,
#' stones, boulders, channers, flagstones, paragravel, paracobbles, parastones,
#' paraboulders, parachanners, paraflagstones
#' @param gravel numeric: gravel volume %
#' @param cobbles numeric: cobble volume %
#' @param stones numeric: stone volume %
#' @param boulders numeric: boulder volume %
#' @param channers numeric: channer volume %
#' @param flagstones numeric: flagstone volume %
#' @param paragravel numeric: para gravel volume %
#' @param paracobbles numeric: para cobble volume %
#' @param parastones numeric: para stone volume %
#' @param paraboulders numeric: para boulder volume %
#' @param parachanners numeric: para channer volume %
#' @param paraflagstones numeric: para flagstone volume %
#'
#' @return - `texmod_to_fragvol`: a data.frame containing `"texmod"` and `"lieutex"` classes
#'
#' @rdname texture
#'
#' @export
#'
#' @examples
#' \donttest{
#' # example of fragvol_to_texmod()
#' df <- expand.grid(
#' gravel = seq(0, 100, 5),
#' cobbles = seq(0, 100, 5),
#' stones = seq(0, 100, 5),
#' boulders = seq(0, 100, 5)
#' )
#' df <- df[rowSums(df) < 100, ]
#'
#' # data.frame input
#' test <- fragvol_to_texmod(df)
#' table(test$texmod)
#' table(test$lieutex)
#'
#' # vector inputs
#' fragvol_to_texmod(gravel = 10, cobbles = 10)
#'
#' }
fragvol_to_texmod <- function(
object = NULL,
gravel = NULL,
cobbles = NULL,
stones = NULL,
boulders = NULL,
channers = NULL,
flagstones = NULL,
paragravel = NULL,
paracobbles = NULL,
parastones = NULL,
paraboulders = NULL,
parachanners = NULL,
paraflagstones = NULL,
as.is = TRUE,
droplevels = TRUE
) {
# standardize inputs ----
var_cols <- c("gravel", "cobbles", "stones", "boulders", "channers", "flagstones", "paragravel", "paracobbles", "parastones", "paraboulders", "parachanners", "paraflagstones")
var_mods <- c("gr", "cb", "st", "by", "cn", "fl", "pgr", "pcb", "pst", "pby", "pcn", "pfl")
# object = NULL; gravel = 10; cobbles = 10; stones = NULL; boulders = NULL; channers = NULL; flagstones = NULL; paragravel = NULL; paracobbles = NULL; parastones = NULL; paraboulders = NULL; parachanners = NULL; paraflagstones = NULL
# object <- expand.grid(gravel = seq(0, 105, 5), cobbles = seq(0, 100, 5), stones = seq(0, 100, 5), boulders = seq(0, 100, 5))
# object <- object[rowSums(object) <= 105, ]
# object <- rbind(object, rep(NA, ncol(object)))
vars_l <- list(gravel, cobbles, stones, boulders, channers, flagstones, paragravel, paracobbles, parastones, paraboulders, parachanners, paraflagstones)
names(vars_l) <- var_cols
# check inputs----
vars_null <- sapply(vars_l, is.null)
df_null <- is.null(object)
vars_len <- sapply(vars_l, length)
## overlapping inputs
if (!df_null & any(!vars_null)) {
warning("if object and any other rock fragment arguments are both not null, only object will be used")
}
# ## object matching column
# if (!df_null & all(!var_cols %in% names(object))) {
# stop(paste("object is missing columns matching any of the following columns",
# paste0(var_cols, collapse = ", "))
# )
# }
# ## vars_l length of inputs
# if (all(vars_len == max(vars_len)) & df_null) {stop("length of inputs do not match")}
# standardize inputs again ----
## subset and select correct inputs
df <- NULL
if (any(!vars_null) & is.null(object)) {
# df <- as.data.frame(vars_l[!vars_null])
df <- .format_inputs(vars_l, var_cols, 0, "numeric")
} else df <- .format_inputs(as.list(object), var_cols, 0, "numeric")
# idx <- var_cols %in% names(df)
# var_cols_sub <- var_cols[idx]
# # var_mods_sub <- var_mods[idx]
# df <- df[var_cols_sub]
#
#
# ## append missing columns
# df_mis <- as.data.frame(matrix(data = 0, nrow = nrow(df), ncol = sum(!idx)))
# names(df_mis) <- var_cols[!idx]
# df <- cbind(df_mis, df)[var_cols]
names(df) <- var_mods
# df[is.na(df)] <- 0
#
#
# # check inputs again ----
# ## correct datatype
# if (!all(sapply(df, function(x) inherits(x, "integer") | inherits(x, "numeric")))) {
# stop("all fragments must be numeric or integer")
# }
## calculate sums ----
gr <- NULL; cb <- NULL; st <- NULL; by <- NULL; cn <- NULL; fl <- NULL
pgr <- NULL; pcb <- NULL; pst <- NULL; pby <- NULL; pcn <- NULL; pfl <- NULL
gr_cn <- NULL; cb_fl <- NULL; pgr_pcn <- NULL; pcb_pfl <- NULL;
sum_nopf <- NULL; sum_pf <- NULL;
fgr <- NULL; fcb <- NULL; fst <- NULL; fby <- NULL; fcn <- NULL; ffl <- NULL
fgr_fcn <- NULL; fcb_ffl <- NULL;
sum_f <- NULL
df$sum_nopf <- rowSums(df[, 1:6], na.rm = TRUE)
df$sum_pf <- rowSums(df[, 7:12], na.rm = TRUE)
df$sum_f <- rowSums(df[, 1:12], na.rm = TRUE)
df$gr_cn <- rowSums(df[c("gr", "cn")], na.rm = TRUE)
df$cb_fl <- rowSums(df[c("cb", "fl")], na.rm = TRUE)
df$pgr_pcn <- rowSums(df[c("pgr", "pcn")], na.rm = TRUE)
df$pcb_pfl <- rowSums(df[c("pcb", "pfl")], na.rm = TRUE)
df$fgr_fcn <- rowSums(df[c("gr", "cn", "pgr", "pcn")], na.rm = TRUE)
df$fcb_ffl <- rowSums(df[c("cb", "fl", "pcb", "pfl")], na.rm = TRUE)
df$fgr <- rowSums(df[c("gr", "pgr")], na.rm = TRUE)
df$fcb <- rowSums(df[c("cb", "pcb")], na.rm = TRUE)
df$fst <- rowSums(df[c("st", "pst")], na.rm = TRUE)
df$fby <- rowSums(df[c("by", "pby")], na.rm = TRUE)
df$fcn <- rowSums(df[c("cn", "pcn")], na.rm = TRUE)
df$ffl <- rowSums(df[c("fl", "pfl")], na.rm = TRUE)
# # load lookup table
# soiltexture <- NULL
# load(system.file("data/soiltexture.rda", package="aqp")[1])
# texmod <- soiltexture$texmod
## filter rows with no rock fragments
idx_sum <- df$sum_nopf < 15 & df$sum_pf < 15 & df$sum_f < 15
df_sub <- df[!idx_sum, ]
# check inputs again ----
if (any(df$sum_nopf > 100)) {
warning("some rows total rock fragments > 100")
}
if (any(df$sum_pf > 100)) {
warning("some rows total pararock fragments > 100")
}
# calculate texmod & lieutex----
df_sub <- within(df_sub, {
# nopf
texmod = NA_character_
lieutex = NA_character_
x_gr_by = gr_cn >= ((1.5 * cb_fl) + (2 * st) + (2.5 * by))
x_cb_by = cb_fl >= ((1.5 * st) + (2 * by))
# 15-34%
x1534 = sum_nopf >= 15 & sum_nopf < 35
texmod = ifelse(x1534 & x_gr_by & gr >= cn & gr > 0 , "gr", texmod)
texmod = ifelse(x1534 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cn", texmod)
texmod = ifelse(x1534 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cb", texmod)
texmod = ifelse(x1534 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "fl", texmod)
texmod = ifelse(x1534 & st >= 1.5 * by & st > 0 & is.na(texmod), "st", texmod)
texmod = ifelse(x1534 & st < 1.5 * by & by > 0 & is.na(texmod), "by", texmod)
# 35-59%
x3559 = sum_nopf >= 35 & sum_nopf < 60
texmod = ifelse(x3559 & x_gr_by & gr >= cn & gr > 0, "grv", texmod)
texmod = ifelse(x3559 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cnv", texmod)
texmod = ifelse(x3559 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cbv", texmod)
texmod = ifelse(x3559 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "flv", texmod)
texmod = ifelse(x3559 & st >= 1.5 * by & st > 0 & is.na(texmod), "stv", texmod)
texmod = ifelse(x3559 & st < 1.5 * by & by > 0 & is.na(texmod), "byv", texmod)
# 60-89%
x6089 = sum_nopf >= 60 & sum_nopf < 90
texmod = ifelse(x6089 & x_gr_by & gr >= cn & gr > 0 , "grx", texmod)
texmod = ifelse(x6089 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cnx", texmod)
texmod = ifelse(x6089 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cbx", texmod)
texmod = ifelse(x6089 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "flx", texmod)
texmod = ifelse(x6089 & st >= 1.5 * by & st > 0 & is.na(texmod), "stx", texmod)
texmod = ifelse(x6089 & st < 1.5 * by & by > 0 & is.na(texmod), "byx", texmod)
# 90-100%
x90 = sum_nopf >= 90
lieutex = ifelse(x90 & x_gr_by & gr >= cn & gr > 0, "gr", lieutex)
lieutex = ifelse(x90 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cn", lieutex)
lieutex = ifelse(x90 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cb", lieutex)
lieutex = ifelse(x90 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "fl", lieutex)
lieutex = ifelse(x90 & st >= 1.5 * by & st > 0 & is.na(texmod), "st", lieutex)
lieutex = ifelse(x90 & st < 1.5 * by & by > 0 & is.na(texmod), "by", lieutex)
# pf
texmod_pf = NA_character_
lieutex_pf = NA_character_
x_pgr_pby = pgr_pcn >= ((1.5 * pcb_pfl) + (2 * pst) + (2.5 * pby))
x_pcb_pby = pcb_pfl >= ((1.5 * pst) + (2 * pby))
# 15-34%
x1534 = sum_pf >= 15 & sum_pf < 35 & sum_nopf <= 0
texmod_pf = ifelse(x1534 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgr", texmod_pf)
texmod_pf = ifelse(x1534 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcn", texmod_pf)
texmod_pf = ifelse(x1534 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcb", texmod_pf)
texmod_pf = ifelse(x1534 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pfl", texmod_pf)
texmod_pf = ifelse(x1534 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pst", texmod_pf)
texmod_pf = ifelse(x1534 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pby", texmod_pf)
# 35-59%
x3559 = sum_pf >= 35 & sum_pf < 60 & sum_nopf <= 0
texmod_pf = ifelse(x3559 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgrv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcnv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcbv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pflv", texmod_pf)
texmod_pf = ifelse(x3559 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pstv", texmod_pf)
texmod_pf = ifelse(x3559 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pbyv", texmod_pf)
# 60-89%
x6089 = sum_pf >= 60 & sum_pf < 90 & sum_nopf <= 0
texmod_pf = ifelse(x6089 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgrx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcnx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcbx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pflx", texmod_pf)
texmod_pf = ifelse(x6089 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pstx", texmod_pf)
texmod_pf = ifelse(x6089 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pbyx", texmod_pf)
# 90-100%
x90 = sum_pf >= 90 & sum_nopf <= 0
lieutex_pf = ifelse(x90 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgr", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcn", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcb", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pfl", lieutex_pf)
lieutex_pf = ifelse(x90 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pst", lieutex_pf)
lieutex_pf = ifelse(x90 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pby", lieutex_pf)
# f
texmod_f = NA_character_
lieutex_f = NA_character_
x_fgr_fby = fgr_fcn >= ((1.5 * fcb_ffl) + (2 * fst) + (2.5 * fby))
x_fcb_fby = fcb_ffl >= ((1.5 * fst) + (2 * fby))
x1534 = sum_f >= 15 & sum_f < 35 & sum_nopf < 15
texmod_f = ifelse(x1534 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgr", texmod_f)
texmod_f = ifelse(x1534 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcn", texmod_f)
texmod_f = ifelse(x1534 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcb", texmod_f)
texmod_f = ifelse(x1534 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pfl", texmod_f)
texmod_f = ifelse(x1534 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pst", texmod_f)
texmod_f = ifelse(x1534 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pby", texmod_f)
# 35-59%
x3559 = sum_f >= 35 & sum_f < 60 & sum_nopf < 15
texmod_f = ifelse(x3559 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgrv", texmod_f)
texmod_f = ifelse(x3559 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcnv", texmod_f)
texmod_f = ifelse(x3559 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcbv", texmod_f)
texmod_f = ifelse(x3559 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pflv", texmod_f)
texmod_f = ifelse(x3559 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pstv", texmod_f)
texmod_f = ifelse(x3559 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pbyv", texmod_f)
# 60-89%
x6089 = sum_f >= 60 & sum_f < 90 & sum_nopf < 15
texmod_f = ifelse(x6089 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgrx", texmod_f)
texmod_f = ifelse(x6089 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcnx", texmod_f)
texmod_f = ifelse(x6089 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcbx", texmod_f)
texmod_f = ifelse(x6089 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pflx", texmod_f)
texmod_f = ifelse(x6089 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pstx", texmod_f)
texmod_f = ifelse(x6089 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pbyx", texmod_f)
# 90-100%
x90 = sum_f >= 90 & sum_nopf < 15
lieutex_f = ifelse(x90 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgr", lieutex_f)
lieutex_f = ifelse(x90 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcn", lieutex_f)
lieutex_f = ifelse(x90 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcb", lieutex_f)
lieutex_f = ifelse(x90 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pfl", lieutex_f)
lieutex_f = ifelse(x90 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pst", lieutex_f)
lieutex_f = ifelse(x90 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pby", lieutex_f)
# combine texmods and lieutex ----
texmod = ifelse(is.na(texmod), texmod_f, texmod)
texmod = ifelse(is.na(texmod), texmod_pf, texmod)
lieutex = ifelse(is.na(lieutex), lieutex_f, lieutex)
lieutex = ifelse(is.na(lieutex), lieutex_pf, lieutex)
})
# idx <- 1:ncol(df_sub)
# df_mod <- df_sub
# df_mod[idx] <- lapply(idx, function(i) {
# texmod <- rep(NA_character_, nrow(df_sub))
# texmod[df_sub[, i] >= 15 & df_sub[, i] < 35] <- var_mods_sub[i]
# texmod[df_sub[, i] >= 35 & df_sub[, i] < 60] <- paste0(var_mods_sub[i], "v")
# texmod[df_sub[, i] >= 60 & df_sub[, i] < 90] <- paste0(var_mods_sub[i], "x")
# return(texmod)
# })
#
# df_lieu <- df_sub
# df_lieu[idx] <- lapply(idx, function(i) {
# lieutex <- rep(NA_character_, nrow(df_sub))
# lieutex[df_sub[i] >= 90] <- var_mods_sub[i]
# return(lieutex)
# })
# standardize output----
if (nrow(df_sub) > 0) {
df$texmod[!idx_sum] <- df_sub$texmod
df$lieutex[!idx_sum] <- df_sub$lieutex
} else df[c("texmod", "lieutex")] <- NA_character_
df <- df[c("texmod", "lieutex")]
if (as.is == FALSE) {
lv <- paste0(
c(rep(var_mods[(7:12)], times = 3, length.out = 18),
rep(var_mods[(1:6)], times = 3, length.out = 18)
),
rep(c("", "v", "x"), each = 6, length.out = 36)
)
tn <- names(sort(table(df$texmod), decreasing = TRUE))
lv <- c(lv, tn[! tn %in% lv])
df$texmod <- factor(df$texmod, levels = lv)
df$lieutex <- factor(df$lieutex, levels = var_mods[1:6])
}
if (as.is == FALSE & droplevels == TRUE) {
df <- droplevels.data.frame(df)
}
return(df)
}
# format function inputs ----
# format vector argument inputs into a data.frame
.format_inputs <- function(vars_l = NULL, var_cols, mis_value = NA, datatype = NULL) {
# check inputs----
vars_null <- sapply(vars_l, is.null)
vars_len <- sapply(vars_l[!vars_null], length)
nm <- var_cols
n <- length(var_cols)
df <- NULL
## all inputs are NULL
if (all(vars_null)) {
stop("all inputs are NULL")
}
## length of inputs
if (!all(vars_len == max(vars_len))) {
stop("length of inputs do not match")
}
# standardize inputs
# create data.frame
df <- as.data.frame(vars_l[!vars_null])
# append missing columns
idx <- nm %in% names(df)
df_mis <- as.data.frame(matrix(data = mis_value, nrow = nrow(df), ncol = sum(!idx)))
names(df_mis) <- nm[!idx]
df <- cbind(df_mis, df)[nm]
## correct datatype
df[1:n] <- lapply(df, function(x) if (is.factor(x)) as.character(x) else x)
df[1:n] <- lapply(df, function(x) if (is.integer(x)) as.numeric(x) else x)
dt <- unlist(lapply(df, class))
## check datatype
idx <- dt == datatype | dt == "logical"
if (!all(idx)) {
stop(paste("unexpected input data types, must be", datatype))
}
return(df)
}
#' @title Ranking Systems for USDA Taxonomic Particle-Size and Substitute Classes of Mineral Soils
#'
#' @description Generate a lookup table of USDA Particle-Size and Substitute Classes names, ranked according to approximate particle size
#'
#' @references \href{https://nrcspad.sc.egov.usda.gov/DistributionCenter/product.aspx?ProductID=991}{Field Book for Describing and Sampling Soils, version 3.0}
#'
#' @return A data.frame with a rank column, taxonomic family particle size class, and a flag for contrasting.
#'
#' @author Stephen Roecker
#'
#' @seealso [hz_to_taxpartsize()], [texture_to_taxpartsize()], [SoilTextureLevels()]
#'
#' @export
#' @examples
#'
#' # class codes
#' lu <- lookup_taxpartsize()
#'
#' idx <- lu$contrasting == FALSE
#'
#' lu$taxpartsize[idx]
#'
#' lu$rank[as.integer(lu$taxpartsize)[idx]]
#'
lookup_taxpartsize <- function() {
fe <- c("diatomaceous", "very-fine", "clayey", "fine", "hydrous", "fine-silty",
"fine-gypseous", "fine-loamy", "medial", "loamy", "coarse-loamy",
"coarse-silty", "coarse-gypseous", "ashy", "sandy", "hydrous-pumiceous",
"medial-pumiceous", "ashy-pumiceous", "clayey-skeletal", "hydrous-skeletal",
"medial-skeletal", "loamy-skeletal", "gypseous-skeletal", "ashy-skeletal",
"sandy-skeletal", "pumiceous", "cindery", "fragmental")
rank <- c(84, 74, 60.02, 46.04, 44.04, 26, 25.8, 25.6, 24, 17.24, 8.88,
8.5, 7.5, 6.5, 4.67, -55.96, -76, -93.5, -43.33, -55.96, -76,
-83.23, -83.35, -93.5, -95.33, -95.83, -96.33, -98.94)
names(rank) <- fe
# cf <- c("fragmental", "sandy-skeletal", "loamy-skeletal", "clay-skeletal")
test <- strsplit(.pscs_sc, " over | or ")
names(test) <- .pscs_sc
idx <- lapply(test, function(x) {
idx <- sapply(x, function(y) rank[which(fe == y)]) |> unlist()
# select the 3rd value when "or" results in 3 values
if (length(idx) > 2) idx <- c(idx[1], idx[3])
dif <- diff(idx)
idx <- idx[1] + sqrt(abs(dif)) * sign(dif)
# l <- dplyr::lag(idx)
# l <- idx[c(NA, 1:(length(idx) - 1))]
# idx <- ifelse(idx < l & !is.na(l), idx * -1, idx * 1)
# n <- length(idx)
# idx <- sum(idx * c(1, 0.1, 0.01, 0.001)[1:n])
return(idx)
})
idx <- round(unlist(idx), 1)
# fe <- data.frame(rn = 1:length(fe), fe = fe)
fe_df <- data.frame(rank = unname(rank), taxpartsize = fe)
sc_df <- data.frame(rank = unname(idx), taxpartsize = .pscs_sc)
lu <- rbind(fe_df, sc_df)
lu <- lu[order(-lu$rank), ]
# lu$rn <- 1:nrow(lu)
lu$contrasting <- grepl(" over ", lu$taxpartsize)
lu$taxpartsize <- factor(lu$taxpartsize, levels = lu$taxpartsize, ordered = TRUE)
row.names(lu) <- NULL
return(lu)
}
.pscs_sc <- c("Ashy over clayey", "Ashy over clayey-skeletal", "Ashy over loamy", "Ashy over loamy-skeletal", "Ashy over medial", "Ashy over medial-skeletal", "Ashy over pumiceous or cindery", "Ashy over sandy or sandy-skeletal", "Ashy-skeletal over clayey", "Ashy-skeletal over fragmental or cindery", "Ashy-skeletal over loamy-skeletal", "Ashy-skeletal over sandy or sandy-skeletal", "Cindery over loamy", "Cindery over medial", "Cindery over medial-skeletal", "Clayey over coarse-gypseous", "Clayey over fine-gypseous", "Clayey over fragmental", "Clayey over gypseous-skeletal", "Clayey over loamy", "Clayey over loamy-skeletal", "Clayey over sandy or sandy-skeletal", "Clayey-skeletal over sandy or sandy-skeletal", "Coarse-loamy over clayey", "Coarse-loamy over fragmental", "Coarse-loamy over sandy or sandy-skeletal", "Coarse-silty over clayey", "Coarse-silty over sandy or sandy-skeletal", "Fine-loamy over clayey", "Fine-loamy over fragmental", "Fine-loamy over sandy or sandy-skeletal", "Fine-silty over clayey", "Fine-silty over fragmental", "Fine-silty over sandy or sandy-skeletal", "Hydrous over clayey", "Hydrous over clayey-skeletal", "Hydrous over fragmental", "Hydrous over loamy", "Hydrous over loamy-skeletal", "Hydrous over sandy or sandy-skeletal", "Loamy over ashy or ashy-pumiceous", "Loamy over coarse-gypseous", "Loamy over fine-gypseous", "Loamy over pumiceous or cindery", "Loamy over sandy or sandy-skeletal", "Loamy-skeletal over cindery", "Loamy-skeletal over clayey", "Loamy-skeletal over fragmental", "Loamy-skeletal over gypseous-skeletal", "Loamy-skeletal over sandy or sandy-skeletal", "Medial over ashy", "Medial over ashy-pumiceous or ashy-skeletal", "Medial over clayey", "Medial over clayey-skeletal", "Medial over fragmental", "Medial over hydrous", "Medial over loamy", "Medial over loamy-skeletal", "Medial over pumiceous or cindery", "Medial over sandy or sandy-skeletal", "Medial-skeletal over fragmental or cindery", "Medial-skeletal over loamy-skeletal", "Medial-skeletal over sandy or sandy-skeletal", "Pumiceous or ashy-pumiceous over loamy", "Pumiceous or ashy-pumiceous over loamy-skeletal", "Pumiceous or ashy-pumiceous over medial", "Pumiceous or ashy-pumiceous over medial-skeletal", "Pumiceous or ashy-pumiceous over sandy or sandy-skeletal", "Sandy over clayey", "Sandy over loamy", "Sandy-skeletal over loamy") |>
tolower()
Try the aqp package in your browser
Any scripts or data that you put into this service are public.
aqp documentation built on Sept. 11, 2024, 7:11 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions texcl_to_ssc
Documented in texcl_to_ssc
#' Textural conversions
#'
#' These functions consist of several conversions between sand, silt and clay
#' to texture class and visa versa, textural modifiers to rock fragments, and
#' grain size composition to the family particle size class.
#'
#' These functions are intended to estimate missing values or allocate particle
#' size fractions to classes. The \code{ssc_to_texcl()} function uses the same
#' logic as the particle size estimator calculation in NASIS to classify sand
#' and clay into texture class. The results are stored in \code{soiltexture}
#' and used by \code{texcl_to_ssc()} as a lookup table to convert texture class
#' to sand, silt and clay. The function \code{texcl_to_ssc()} replicates the
#' functionality described by Levi (2017). The \code{texmod_to_fragvol()}
#' function similarly uses the logical from the
#' Exhibit618-11_texture_modifier.xls spreadsheet to determine the textural
#' modifier from the various combinations of rock and pararock fragments
#' (e.g. GR and PGR).
#'
#' When \code{sample = TRUE}, the results can be used to estimate within-class,
#' marginal distributions of sand, silt, and clay fractions. It is recommended
#' that at least 10 samples be drawn for reasonable estimates.
#'
#' The function \code{texmod_to_fragvoltot} returns a data.frame with multiple
#' fragvoltot columns differentiated by tailing abbreviations (e.g. _r) which
#' refer to the following:
#' 1. l = low
#' 2. r = representative
#' 3. h = high
#' 4. nopf = no pararock fragments (i.e. total fragments - pararock fragments)
#'
#' The function \code{texture_to_texmod()} parses texture (e.g. GR-CL) to extract the texmod values from it in the scenario where it is missing from
#' texmod column. If multiple texmod values are present (for example in the
#' case of stratified textures) and \code{duplicates = "combine"} they will be combined in the output (e.g. GR & CBV). Otherwise if \code{duplicates = "max"}
#' the texmod with the highest rock fragment (e.g. CBV) will be returned.
#'
#' Unlike the other functions, \code{texture_to_taxpartsize()} is intended to
#' be computed on weighted averages within the family particle size control
#' section. Also recall from the criteria that carbonate clay should be
#' subtracted from clay content and added to silt content. Similarly, if the
#' percent of very fine sand is known it should be subtracted from the sand,
#' and added to the silt content. Unlike the other functions,
#' \code{texture_to_taxpartsize()} is intended to be computed on weighted
#' averages within the family particle size control section. Also recall from
#' the criteria that carbonate clay should be subtracted from clay content and
#' added to silt content. Similarly, if the percent of very fine sand is known
#' it should be subtracted from the sand, and added to the silt content.
#' @param texcl vector of texture classes than conform to the USDA code
#' conventions (e.g. c|C, sil|SIL, sl|SL, cos|COS)
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#' @param lieutex vector of in lieu of texture terms that conform to the USDA
#' code convenctions (e.g. gr|GR, pg|PG), only used when fragments or artifacts
#' are > 90 percent by volume (default: NULL))
#' @param texture vector of combinations of texcl, texmod, and lieutex (e.g. CL, GR-CL, CBV-S, GR)
#'
#' @param clay vector of clay percentages
#' @param sand vector of sand percentages
#' @param fragvoltot vector of rock fragment percentages
#'
#' @param as.is logical: should character vectors be converted to factors?
#' (default: TRUE)
#'
#' @param droplevels logical: indicating whether to drop unused levels in
#' factors. This is useful when the results have a large number of unused
#' classes, which can waste space in tables and figures.
#' @param sample logical: should ssc be random sampled from the lookup table?
#' (default: FALSE)
#'
#' @return - `texcl_to_ssc`: A `data.frame` containing columns `"sand"`,`"silt"`, `"clay"`
#'
#' @seealso \code{\link{SoilTextureLevels}}
#'
#' @author Stephen Roecker
#'
#' @references Matthew R. Levi, Modified Centroid for Estimating Sand, Silt, and Clay from Soil Texture Class, Soil Science Society of America Journal, 2017, 81(3):578-588, ISSN 1435-0661, \doi{10.2136/sssaj2016.09.0301}.
#'
#' @rdname texture
#' @keywords manip
#' @export
#' @examples
#' \donttest{
#' # example of ssc_to_texcl()
#' tex <- expand.grid(sand = 0:100, clay = 0:100)
#' tex <- subset(tex, (sand + clay) < 101)
#' tex$texcl <- ssc_to_texcl(sand = tex$sand, clay = tex$clay)
#' head(tex)
#'
#' # example of texcl_to_ssc(texcl)
#' texcl <- c("cos", "s", "fs", "vfs", "lcos", "ls",
#' "lfs", "lvfs", "cosl", "sl", "fsl", "vfsl", "l",
#' "sil", "si", "scl", "cl", "sicl", "sc", "sic", "c"
#' )
#' test <- texcl_to_ssc(texcl)
#' head(test <- cbind(texcl, test), 10)
#'
#'
#' # example of texcl_to_ssc(texcl, clay)
#' data(soiltexture)
#' st <- soiltexture$values
#' idx <- sample(1:length(st$texcl), 10)
#' st <- st[idx, ]
#' ssc <- texcl_to_ssc(texcl = st$texcl)
#' head(cbind(texcl = st$texcl, clay = ssc$clay))
#'
#'
#' # example of texmod_to_fragvoltol
#' frags <- c("gr", "grv", "grx", "pgr", "pgrv", "pgrx")
#' head(texmod_to_fragvoltot(frags))
#'
#'
#' # example of texture_to_taxpartsize()
#' tex <- data.frame(texcl = c("c", "cl", "l", "ls", "s"),
#' clay = c(55, 33, 18, 6, 3),
#' sand = c(20, 33, 42, 82, 93),
#' fragvoltot = c(35, 15, 34, 60, 91))
#' tex$fpsc <- texture_to_taxpartsize(texcl = tex$texcl,
#' clay = tex$clay,
#' sand = tex$sand,
#' fragvoltot = tex$fragvoltot)
#' head(tex)
#'
#'
#' # example of texture_to_taxpartsize() with carbonate clay and very fine sand
#' carbclay <- rnorm(5, 2, 3)
#' vfs <- rnorm(5, 10, 3)
#' st$fpsc <- texture_to_taxpartsize(texcl = tex$texcl,
#' clay = tex$clay - carbclay,
#' sand = tex$sand - vfs,
#' fragvoltot = tex$fragvoltot)
#' head(tex)
#'
#'
#' # example of sample = TRUE
#' texcl <- rep(c("cl", "sil", "sl"), 10)
#' ssc1 <- cbind(texcl, texcl_to_ssc(texcl = texcl, sample = FALSE))
#' ssc2 <- cbind(texcl, texcl_to_ssc(texcl = texcl, sample = TRUE))
#' ssc1$sample <- FALSE
#' ssc2$sample <- TRUE
#' ssc <- rbind(ssc1, ssc2)
#' aggregate(clay ~ sample + texcl, data = ssc, summary)
#' }
texcl_to_ssc <- function(texcl = NULL, clay = NULL, sample = FALSE) {
# fix for R CMD check
# texcl_to_ssc: no visible binding for global variable ‘soiltexture’
soiltexture <- NULL
# clay is not NULL
clay_not_null <- all(!is.null(clay))
clay_is_null <- !clay_not_null
# standardize the inputs
df <- data.frame(texcl = tolower(as.character(texcl)),
stringsAsFactors = FALSE
)
if (clay_not_null) {
df$clay <- as.integer(round(clay))
}
df$rn <- row.names(df)
load(system.file("data/soiltexture.rda", package="aqp")[1])
# convert fine sand classes to their generic counterparts
df <- within(df, {
texcl = ifelse(texcl %in% c("cos", "fs", "vfs"), "s", texcl)
texcl = ifelse(texcl %in% c("lcos", "lfs", "lvfs"), "ls", texcl)
texcl = ifelse(texcl %in% c("cosl", "fsl", "vfsl"), "sl", texcl)
})
# check for texcl that don't match
idx <- ! df$texcl %in% unique(soiltexture$averages$texcl)
if (any(idx)) {
warning("not all the user supplied texcl values match the lookup table")
}
# check clay values within texcl
if (clay_not_null) {
clay_not_na <- !is.na(df$clay)
idx <- paste(df$texcl[clay_not_na], df$clay[clay_not_na]) %in% paste(soiltexture$values$texcl, soiltexture$values$clay)
if (any(!idx)) {
warning("not all the user supplied clay values fall within the texcl, so they will be set to NA")
df$clay[which(!idx)] <- NA
clay_not_null <- all(!is.na(df$clay))
clay_is_null <- !clay_not_null
}
}
# check clay ranges 0-100
idx <- clay_not_null & any(clay < 0, na.rm = TRUE) & any(clay > 100, na.rm = TRUE)
if (idx) {
warning("some clay records < 0 or > 100%")
}
# if clay is present
if (clay_not_null & sample == FALSE) {
st <- aggregate(sand ~ texcl + clay, data = soiltexture$values, function(x) as.integer(round(mean(x))))
st$silt <- 100 - st$clay - st$sand
# some clay present (compute clay weighted averages)
idx <- is.na(df$clay)
if (any(idx)) {
df_na <- merge(df[idx, c("texcl", "rn")], soiltexture$averages, by = "texcl", all.x = TRUE, sort = FALSE)[c("texcl", "clay", "rn")]
df[idx, ] <- df_na
}
df <- merge(df[c("texcl", "clay", "rn")], st, by = c("texcl", "clay"), all.x = TRUE, sort = FALSE)
} else {
# clay missing (use average)
df <- merge(df[c("texcl", "rn")], soiltexture$averages, by = "texcl", all.x = TRUE, sort = FALSE)
}
# randomly sample ssc from texcl lookup table
if (sample == TRUE) {
if (clay_is_null) df$clay <- NA
split(df, df$texcl, drop = TRUE) ->.;
lapply(., function(x) {
st <- soiltexture$values
# clay present
x$idx <- is.na(x$clay)
x1 <- x[x$idx == FALSE, ]
x2 <- x[x$idx == TRUE, ]
if (clay_not_null) {
temp1 <- st[st$texcl == x1$texcl[1] &
st$clay %in% unique(x1$clay),
]
temp1 <- temp1[sample(1:nrow(temp1), size = nrow(x1), replace = TRUE), ]
temp1 <- cbind(x1[x1$idx == FALSE, c("rn", "texcl")], temp1[c("clay", "sand")])
} else temp1 <- NULL
# clay missing
temp2 <- st[st$texcl == x2$texcl[1], ]
temp2 <- temp2[sample(1:nrow(temp2), size = nrow(x2), replace = TRUE), ]
temp2 <- cbind(x2[x2$idx == TRUE, c("rn", "texcl")], temp2[c("clay", "sand")])
return(rbind(temp1, temp2))
}) ->.;
do.call("rbind", .) -> df
df$silt <- 100 - df$clay - df$sand
}
# standardize outputs
vars <- c("sand", "silt", "clay")
df <- df[(order(as.integer(df$rn))), vars]
df$rn <- NULL
df$texcl <- NULL
rownames(df) <- NULL
return(df)
}
#' Convert sand, silt and clay to texture class
#' @param simplify Passed to `SoilTextureLevels()` to set the number of possible texture classes. If `TRUE`, the ordered factor has a maximum of 12 levels, if `FALSE` (default) the ordered factor has a maximum of 21 levels (including e.g. very fine/fine/coarse variants)
#' @rdname texture
#' @return - `ssc_to_texcl`: A `character` vector containing texture class
#' @export
#'
ssc_to_texcl <- function(sand = NULL, clay = NULL, simplify = FALSE, as.is = FALSE, droplevels = TRUE) {
# fix for R CMD check:
# ssc_to_texcl: no visible binding for global variable ‘silt’
silt <- NULL
# check lengths
idx <- length(clay) != length(sand)
if (idx) {
stop("length of inputs do not match")
}
# standardize inputs
df <- data.frame(sand = as.integer(round(sand)),
clay = as.integer(round(clay)),
stringsAsFactors = FALSE
)
df$silt <- 100 - df$clay - df$sand
## TODO: this needs some more work: sum will always be 100, but silt-by-difference may be illogical
# # check sand, silt and clay sum to 100
# idx <- (df$sand + df$silt + df$clay) > 100 | (df$sand + df$silt + df$clay) < 100
# if (any(idx) & any(complete.cases(df[c("sand", "clay")]))) {
# warning("some records sand, silt, and clay do not sum to 100 %")
# }
# logic from the particle size estimator calculation from NASIS
df <- within(df, {
texcl = NA
texcl[silt >= 79.99 & clay < 11.99] = "si"
texcl[silt >= 49.99 & clay < 26.99 & (silt < 79.99 | clay >= 11.99)] = "sil"
texcl[clay >= 26.99 & clay < 39.99 & sand <= 20.01] = "sicl"
texcl[clay >= 39.99 & silt >= 39.99] = "sic"
texcl[clay >= 39.99 & sand <= 45.01 & silt < 39.99] = "c"
texcl[clay >= 26.99 & clay < 39.99 & sand > 20.01 & sand <= 45.01] = "cl"
texcl[clay >= 6.99 & clay < 26.99 & silt >= 27.99 & silt < 49.99 & sand <= 52.01] = "l"
texcl[clay >= 19.99 & clay < 34.99 & silt < 27.99 & sand > 45.01] = "scl"
texcl[clay >= 34.99 & sand > 45.01] = "sc"
texcl[(silt + 1.5 * clay) < 15] = "s"
texcl[(silt + 1.5 * clay) >= 15 & (silt + 2 * clay) < 29.99] = "ls"
texcl[!is.na(sand) & !is.na(clay) & is.na(texcl)] = "sl"
})
# encoding according to approximate AWC, from Red Book version 3.0
if (!as.is) {
df$texcl <- factor(df$texcl, levels = SoilTextureLevels(which = 'codes', simplify = simplify), ordered = TRUE)
if (droplevels) {
df$texcl <- droplevels(df$texcl)
}
}
return(df$texcl)
}
#' Convert from fragment modifier to estimated total volume
#'
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#'
#' @param lieutex vector of in lieu of texture terms that conform to the USDA
#' code conventions (e.g. gr|GR, pg|PG), only used when fragments or artifacts
#' are > 90 percent by volume (default: NULL))
#' @rdname texture
#'
#' @return - `texmod_to_fragvoltot`: A `data.frame` containing columns `"fragvoltot_l"`,
#' `"fragvoltot_r"`, `"fragvoltot_h"`, `"fragvoltot_l_nopf"`, `"fragvoltot_r_nopf"`, `"fragvoltot_h_nopf"`
#'
#' @export
texmod_to_fragvoltot <- function(texmod = NULL, lieutex = NULL) {
# standardize inputs ----
vars_l <- list(texmod = texmod, lieutex = lieutex)
# vars_l <- list(texmod = c("GR", "CBV", NA), lieutex = c(NA, NA, "BY"))
df <- NULL
df <- .format_inputs(vars_l, names(vars_l), NA_character_, "character")
df$texmod <- tolower(df$texmod)
df$lieutex <- toupper(df$lieutex)
df$rn <- 1:nrow(df)
# load lookup table
# fix for R CMD check
# texcl_to_ssc: no visible binding for global variable ‘soiltexture’
soiltexture <- NULL
load(system.file("data/soiltexture.rda", package="aqp")[1])
# check inputs ----
## both present?
idx <- any(complete.cases(df))
if (idx) {
warning("texmod and lieutex should not both be present, they are mutually exclusive, only the results for texmod will be returned")
}
# texmod and lieutex don't match
idx <- ! df$texmod %in% soiltexture$texmod$texmod
if (any(idx & !is.na(df$texmod))) {
message("not all the texmod supplied match the lookup table")
}
idx <- ! df$lieutex %in% c("GR", "CB", "ST", "BY", "CN", "FL", "PG", "PCB", "PST", "PBY", "PCN", "PFL", "BR", "HMM", "MPM", "SPM", "MUCK", "PEAT", "ART", "CGM", "FGM", "ICE", "MAT", "W")
if (any(idx & !is.na(df$lieutex))) {
message("not all the lieutex supplied match the lookup table")
}
# merge
df <- merge(df, soiltexture$texmod, by = "texmod", all.x = TRUE, sort = FALSE)
df <- df[order(df$rn), ]
df$rn <- NULL
# lieutex
if (any(!is.na(df$lieutex))) {
idx1 <- !is.na(df$lieutex) & df$lieutex %in% c("GR", "CB", "ST", "BY", "CN", "FL")
idx2 <- !is.na(df$lieutex) & df$lieutex %in% c("PG", "PCB", "PST", "PBY", "PCN", "PFL")
# df$lieutex <- toupper(lieutex)
df <- within(df, {
fragvoltot_l = ifelse(idx1, 90, fragvoltot_l)
fragvoltot_r = ifelse(idx1, 95, fragvoltot_l)
fragvoltot_h = ifelse(idx1, 100, fragvoltot_l)
fragvoltot_l_nopf = ifelse(idx2, 0, fragvoltot_l)
fragvoltot_r_nopf = ifelse(idx2, 0, fragvoltot_l)
fragvoltot_h_nopf = ifelse(idx2, 0, fragvoltot_l)
})
# df$lieutex <- lieutex
}
df[c("texmod", "lieutex", "texmod_label")] <- NULL
return(df)
}
#' Convert sand, silt and clay to the family particle size class
#'
#' @param texcl vector of texture classes than conform to the USDA code
#' conventions (e.g. c|C, sil|SIL, sl|SL, cos|COS)
#'
#' @param clay vector of clay percentages
#' @param sand vector of sand percentages
#' @param sandvf vector of very fine sand percentages
#'
#' @param fragvoltot vector of total rock fragment percentages
#'
#' @return - `texture_to_taxpartsize`: a character vector containing `"taxpartsize"` classes
#'
#' @seealso [hz_to_taxpartsize()], [lookup_taxpartsize()]
#'
#' @rdname texture
#'
#' @export
#'
texture_to_taxpartsize <- function(texcl = NULL, clay = NULL, sand = NULL, sandvf = NULL, fragvoltot = NULL) {
# check lengths
idx <- length(texcl) == length(clay) & length(clay) == length(sand) & length(sand) == length(fragvoltot)
if (!idx) {
stop("length of inputs do not match")
}
# standarize inputs
if (is.null(sandvf)) sandvf <- NA
df <- data.frame(texcl = tolower(texcl),
clay = as.integer(round(clay)),
sand = as.integer(round(sand)),
sandvf = as.integer(round(sandvf)),
fragvoltot = as.integer(round(fragvoltot)),
fpsc = as.character(NA),
stringsAsFactors = FALSE
)
df$silt <- 100 - df$sand - df$clay
sandytextures <- c("cos", "s", "fs", "lcos", "ls", "lfs")
# check texcl lookup
idx <- any(! df$texcl[!is.na(df$texcl)] %in% SoilTextureLevels(which = 'codes'))
if (idx) {
warning("not all the texcl supplied match the lookup table")
}
# check percentages
idx <- df$silt > 100 | df$silt < 0 | df$clay > 100 | df$clay < 0 | df$sand > 100 | df$sand < 0 | df$fragvoltot > 100 | df$fragvoltot < 0
if (any(idx, na.rm = TRUE)) {
warning("some records are > 100% or < 0%, or the calcuated silt fraction is > 100% or < 0%")
}
if (any(sandvf > sand & all(!is.na(sandvf)))) {
warning("the sandvf values should not be greater than the sand values")
}
# check ssc_to_texcl() vs texcl
df$texcl_calc <- suppressMessages(ssc_to_texcl(sand = df$sand, clay = df$clay, as.is = TRUE))
df <- within(df, {
texcl_calc = ifelse(texcl_calc == "s" & grepl("^cos$|^fs$|^vfs$", texcl), texcl, texcl_calc)
texcl_calc = ifelse(texcl_calc == "ls" & grepl("^lcos$|^lfs$|^lvfs$", texcl), texcl, texcl_calc)
texcl_calc = ifelse(texcl_calc == "sl" & grepl("^cosl$|^fsl$|^vfsl$", texcl), texcl, texcl_calc)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("vfs", "lvfs"), 50, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("vfsl"), 40, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("fsl"), 15, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("sl", "l", "scl"), 10, sandvf)
sandvf = ifelse(is.na(sandvf) & texcl %in% c("fsl"), 7, sandvf)
sand = ifelse(!is.na(sandvf), sand - sandvf, sand)
silt = ifelse(!is.na(sandvf), silt + sandvf, silt)
})
idx <- any(df$texcl != df$texcl_calc, na.rm = TRUE)
if (idx) {
warning("some of the texcl records don't match the calculated texcl via ssc_to_texcl()")
}
# calculate family particle size control section
idx <- df$fragvoltot >= 35
if (any(idx)) {
df[idx,] <- within(df[idx,], {
fpsc[clay >= 35] = "clayey-skeletal"
fpsc[clay < 35] = "loamy-skeletal"
fpsc[texcl %in% sandytextures] = "sandy-skeletal"
})
}
idx <- df$fragvoltot < 35 & df$texcl %in% sandytextures
if (any(idx)) {
df[idx, ]$fpsc <- "sandy"
}
idx <- df$fragvoltot < 35 & ! df$texcl %in% sandytextures
if (any(idx)) {
df[idx, ] <- within(df[idx,], {
fpsc[clay < 18 & sand >= 15] = "coarse-loamy"
fpsc[clay < 18 & sand < 15] = "coarse-silty"
fpsc[clay >= 18 & clay < 35] = "fine-loamy"
fpsc[clay >= 18 & clay < 35 & sand < 15] = "fine-silty"
fpsc[clay >= 35 & clay < 60] = "fine"
fpsc[clay > 60] = "very-fine"
})
}
df$fpsc <- ifelse(df$fragvoltot > 90, "fragmental", df$fpsc)
return(df$fpsc)
}
#' Parse texmod from texture
#'
#' @param texmod vector of textural modifiers that conform to the USDA code
#' conventions (e.g. gr|GR, grv|GRV)
#' @param texture vector of combinations of texcl, texmod, and lieutex (e.g. CL, GR-CL, CBV-S, GR)
#'
#' @param duplicates character: specifying how multiple values should be handled, options are `"combined"` (e.g. 'GR & GRV) or `"max"`(e.g. 'GRV')
#'
#' @return - `texture_to_texmod`: a character vector containing `"texmod"` classes
#'
#' @rdname texture
#'
#' @export
#'
#' @examples
#' \donttest{
#' # example of texture_to_texmod()
#' tex <- c("SL", "GR-SL", "CBV-L", "SR- GR-FS GRX-COS")
#' texture_to_texmod(tex)
#' texture_to_texmod(tex, duplicates = "max")
#' }
texture_to_texmod <- function(texture, duplicates = "combine") {
# test
is_char <- is.character(texture)
is_fact <- is.factor(texture)
if (!is_char & !is_fact) stop("texture must be a character or a factor")
if (is_fact) texture <- as.character(texture)
# load lookup tables
soiltexture <- NULL
load(system.file("data/soiltexture.rda", package = "aqp")[1])
texmod_df <- soiltexture$texmod
texmod <- texmod_df$texmod[grepl("gravel|cobbl|ston|boulder|channer|flag", texmod_df$texmod_label)]
# lieutex <- metadata[metadata$ColumnPhysicalName == "lieutex", ]$ChoiceName
tex <- tolower(texture)
tex_l <- strsplit(tex, "-| |_|\\/|&|,|;|\\.|\\*|\\+")
# find which records have a texmod
idx_mod <- which(sapply(tex_l, function(x) any(x[-length(x)] %in% texmod)))
# idx_lieu <- which(sapply(tex_l, function(x) any(x %in% lieutex)))
# create a logical matrix of texmod presence/absence
if (length(idx_mod > 0)) {
tex_m <- lapply(tex_l[idx_mod], function(x)
texmod %in% x)
tex_df <- as.data.frame(do.call("rbind", tex_m))
names(tex_df) <- texmod
# slightly slower simpler(/better?) alternative
# not better, this mistakenly matches words entered into this field
# I would rather have false negatives than false positives
# tex_m <- sapply(texmod, function(x) grepl(x, tex))
# duplicates ----
## combine texmod (if multiple exist)
if (duplicates == "combine") {
texmod_parse <- apply(tex_df, 1, function(x) {
paste0(names(x)[which(x)], collapse = " & ")
})
}
## find texmod with the highest rock fragment content (if multiple exist)
if (duplicates == "max") {
tex_m <- sapply(texmod, function(x) {
vals <- suppressMessages(texmod_to_fragvoltot(x))
rf <- ifelse(tex_df[x] == TRUE, vals$fragvoltot_r, 0)
})
idx_col <- max.col(tex_m)
texmod_parse <- names(tex_df)[idx_col]
}
}
# results ----
n <- length(texture)
df <- data.frame(texmod = rep(NA_character_, n)) #, lieutex = rep(NA, n))
if (length(idx_mod) > 0) df[idx_mod, "texmod"] <- texmod_parse
return(df$texmod)
}
# add feature to uncode
# texture_to_texcl <- function(texclLabel, droplevels = FALSE
# # , stringsAsFactors = FALSE
# ) {
#
# # standardize inputs
# texclLabel <- tolower(texclLabel)
#
# # load lookup table
# load(system.file("data/soiltexture.rda", package = "aqp")[1])
# sub <- soiltexture$averages
#
# # check
# idx <- ! texclLabel %in% sub$texcl_label
# if (any(idx)) {
# warning("not all the user supplied texture values match the lookup table")
# }
#
# # convert
# texcl <- factor(texclLabel, levels = tolower(sub$texcl_label), label = sub$texcl, ordered = TRUE)
#
# if (droplevels == TRUE) {
# texcl <- droplevels(texcl)
# }
#
# # if (stringsAsFactors == TRUE) {
# # texcl <- as.character(texcl)
# # }
#
# return(texcl)
# }
#' Convert frags to texmod
#'
#' @param object data.frame: containing the following column names: gravel, cobbles,
#' stones, boulders, channers, flagstones, paragravel, paracobbles, parastones,
#' paraboulders, parachanners, paraflagstones
#' @param gravel numeric: gravel volume %
#' @param cobbles numeric: cobble volume %
#' @param stones numeric: stone volume %
#' @param boulders numeric: boulder volume %
#' @param channers numeric: channer volume %
#' @param flagstones numeric: flagstone volume %
#' @param paragravel numeric: para gravel volume %
#' @param paracobbles numeric: para cobble volume %
#' @param parastones numeric: para stone volume %
#' @param paraboulders numeric: para boulder volume %
#' @param parachanners numeric: para channer volume %
#' @param paraflagstones numeric: para flagstone volume %
#'
#' @return - `texmod_to_fragvol`: a data.frame containing `"texmod"` and `"lieutex"` classes
#'
#' @rdname texture
#'
#' @export
#'
#' @examples
#' \donttest{
#' # example of fragvol_to_texmod()
#' df <- expand.grid(
#' gravel = seq(0, 100, 5),
#' cobbles = seq(0, 100, 5),
#' stones = seq(0, 100, 5),
#' boulders = seq(0, 100, 5)
#' )
#' df <- df[rowSums(df) < 100, ]
#'
#' # data.frame input
#' test <- fragvol_to_texmod(df)
#' table(test$texmod)
#' table(test$lieutex)
#'
#' # vector inputs
#' fragvol_to_texmod(gravel = 10, cobbles = 10)
#'
#' }
fragvol_to_texmod <- function(
object = NULL,
gravel = NULL,
cobbles = NULL,
stones = NULL,
boulders = NULL,
channers = NULL,
flagstones = NULL,
paragravel = NULL,
paracobbles = NULL,
parastones = NULL,
paraboulders = NULL,
parachanners = NULL,
paraflagstones = NULL,
as.is = TRUE,
droplevels = TRUE
) {
# standardize inputs ----
var_cols <- c("gravel", "cobbles", "stones", "boulders", "channers", "flagstones", "paragravel", "paracobbles", "parastones", "paraboulders", "parachanners", "paraflagstones")
var_mods <- c("gr", "cb", "st", "by", "cn", "fl", "pgr", "pcb", "pst", "pby", "pcn", "pfl")
# object = NULL; gravel = 10; cobbles = 10; stones = NULL; boulders = NULL; channers = NULL; flagstones = NULL; paragravel = NULL; paracobbles = NULL; parastones = NULL; paraboulders = NULL; parachanners = NULL; paraflagstones = NULL
# object <- expand.grid(gravel = seq(0, 105, 5), cobbles = seq(0, 100, 5), stones = seq(0, 100, 5), boulders = seq(0, 100, 5))
# object <- object[rowSums(object) <= 105, ]
# object <- rbind(object, rep(NA, ncol(object)))
vars_l <- list(gravel, cobbles, stones, boulders, channers, flagstones, paragravel, paracobbles, parastones, paraboulders, parachanners, paraflagstones)
names(vars_l) <- var_cols
# check inputs----
vars_null <- sapply(vars_l, is.null)
df_null <- is.null(object)
vars_len <- sapply(vars_l, length)
## overlapping inputs
if (!df_null & any(!vars_null)) {
warning("if object and any other rock fragment arguments are both not null, only object will be used")
}
# ## object matching column
# if (!df_null & all(!var_cols %in% names(object))) {
# stop(paste("object is missing columns matching any of the following columns",
# paste0(var_cols, collapse = ", "))
# )
# }
# ## vars_l length of inputs
# if (all(vars_len == max(vars_len)) & df_null) {stop("length of inputs do not match")}
# standardize inputs again ----
## subset and select correct inputs
df <- NULL
if (any(!vars_null) & is.null(object)) {
# df <- as.data.frame(vars_l[!vars_null])
df <- .format_inputs(vars_l, var_cols, 0, "numeric")
} else df <- .format_inputs(as.list(object), var_cols, 0, "numeric")
# idx <- var_cols %in% names(df)
# var_cols_sub <- var_cols[idx]
# # var_mods_sub <- var_mods[idx]
# df <- df[var_cols_sub]
#
#
# ## append missing columns
# df_mis <- as.data.frame(matrix(data = 0, nrow = nrow(df), ncol = sum(!idx)))
# names(df_mis) <- var_cols[!idx]
# df <- cbind(df_mis, df)[var_cols]
names(df) <- var_mods
# df[is.na(df)] <- 0
#
#
# # check inputs again ----
# ## correct datatype
# if (!all(sapply(df, function(x) inherits(x, "integer") | inherits(x, "numeric")))) {
# stop("all fragments must be numeric or integer")
# }
## calculate sums ----
gr <- NULL; cb <- NULL; st <- NULL; by <- NULL; cn <- NULL; fl <- NULL
pgr <- NULL; pcb <- NULL; pst <- NULL; pby <- NULL; pcn <- NULL; pfl <- NULL
gr_cn <- NULL; cb_fl <- NULL; pgr_pcn <- NULL; pcb_pfl <- NULL;
sum_nopf <- NULL; sum_pf <- NULL;
fgr <- NULL; fcb <- NULL; fst <- NULL; fby <- NULL; fcn <- NULL; ffl <- NULL
fgr_fcn <- NULL; fcb_ffl <- NULL;
sum_f <- NULL
df$sum_nopf <- rowSums(df[, 1:6], na.rm = TRUE)
df$sum_pf <- rowSums(df[, 7:12], na.rm = TRUE)
df$sum_f <- rowSums(df[, 1:12], na.rm = TRUE)
df$gr_cn <- rowSums(df[c("gr", "cn")], na.rm = TRUE)
df$cb_fl <- rowSums(df[c("cb", "fl")], na.rm = TRUE)
df$pgr_pcn <- rowSums(df[c("pgr", "pcn")], na.rm = TRUE)
df$pcb_pfl <- rowSums(df[c("pcb", "pfl")], na.rm = TRUE)
df$fgr_fcn <- rowSums(df[c("gr", "cn", "pgr", "pcn")], na.rm = TRUE)
df$fcb_ffl <- rowSums(df[c("cb", "fl", "pcb", "pfl")], na.rm = TRUE)
df$fgr <- rowSums(df[c("gr", "pgr")], na.rm = TRUE)
df$fcb <- rowSums(df[c("cb", "pcb")], na.rm = TRUE)
df$fst <- rowSums(df[c("st", "pst")], na.rm = TRUE)
df$fby <- rowSums(df[c("by", "pby")], na.rm = TRUE)
df$fcn <- rowSums(df[c("cn", "pcn")], na.rm = TRUE)
df$ffl <- rowSums(df[c("fl", "pfl")], na.rm = TRUE)
# # load lookup table
# soiltexture <- NULL
# load(system.file("data/soiltexture.rda", package="aqp")[1])
# texmod <- soiltexture$texmod
## filter rows with no rock fragments
idx_sum <- df$sum_nopf < 15 & df$sum_pf < 15 & df$sum_f < 15
df_sub <- df[!idx_sum, ]
# check inputs again ----
if (any(df$sum_nopf > 100)) {
warning("some rows total rock fragments > 100")
}
if (any(df$sum_pf > 100)) {
warning("some rows total pararock fragments > 100")
}
# calculate texmod & lieutex----
df_sub <- within(df_sub, {
# nopf
texmod = NA_character_
lieutex = NA_character_
x_gr_by = gr_cn >= ((1.5 * cb_fl) + (2 * st) + (2.5 * by))
x_cb_by = cb_fl >= ((1.5 * st) + (2 * by))
# 15-34%
x1534 = sum_nopf >= 15 & sum_nopf < 35
texmod = ifelse(x1534 & x_gr_by & gr >= cn & gr > 0 , "gr", texmod)
texmod = ifelse(x1534 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cn", texmod)
texmod = ifelse(x1534 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cb", texmod)
texmod = ifelse(x1534 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "fl", texmod)
texmod = ifelse(x1534 & st >= 1.5 * by & st > 0 & is.na(texmod), "st", texmod)
texmod = ifelse(x1534 & st < 1.5 * by & by > 0 & is.na(texmod), "by", texmod)
# 35-59%
x3559 = sum_nopf >= 35 & sum_nopf < 60
texmod = ifelse(x3559 & x_gr_by & gr >= cn & gr > 0, "grv", texmod)
texmod = ifelse(x3559 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cnv", texmod)
texmod = ifelse(x3559 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cbv", texmod)
texmod = ifelse(x3559 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "flv", texmod)
texmod = ifelse(x3559 & st >= 1.5 * by & st > 0 & is.na(texmod), "stv", texmod)
texmod = ifelse(x3559 & st < 1.5 * by & by > 0 & is.na(texmod), "byv", texmod)
# 60-89%
x6089 = sum_nopf >= 60 & sum_nopf < 90
texmod = ifelse(x6089 & x_gr_by & gr >= cn & gr > 0 , "grx", texmod)
texmod = ifelse(x6089 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cnx", texmod)
texmod = ifelse(x6089 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cbx", texmod)
texmod = ifelse(x6089 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "flx", texmod)
texmod = ifelse(x6089 & st >= 1.5 * by & st > 0 & is.na(texmod), "stx", texmod)
texmod = ifelse(x6089 & st < 1.5 * by & by > 0 & is.na(texmod), "byx", texmod)
# 90-100%
x90 = sum_nopf >= 90
lieutex = ifelse(x90 & x_gr_by & gr >= cn & gr > 0, "gr", lieutex)
lieutex = ifelse(x90 & x_gr_by & gr < cn & cn > 0 & is.na(texmod), "cn", lieutex)
lieutex = ifelse(x90 & x_cb_by & cb >= fl & cb > 0 & is.na(texmod), "cb", lieutex)
lieutex = ifelse(x90 & x_cb_by & cb < fl & fl > 0 & is.na(texmod), "fl", lieutex)
lieutex = ifelse(x90 & st >= 1.5 * by & st > 0 & is.na(texmod), "st", lieutex)
lieutex = ifelse(x90 & st < 1.5 * by & by > 0 & is.na(texmod), "by", lieutex)
# pf
texmod_pf = NA_character_
lieutex_pf = NA_character_
x_pgr_pby = pgr_pcn >= ((1.5 * pcb_pfl) + (2 * pst) + (2.5 * pby))
x_pcb_pby = pcb_pfl >= ((1.5 * pst) + (2 * pby))
# 15-34%
x1534 = sum_pf >= 15 & sum_pf < 35 & sum_nopf <= 0
texmod_pf = ifelse(x1534 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgr", texmod_pf)
texmod_pf = ifelse(x1534 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcn", texmod_pf)
texmod_pf = ifelse(x1534 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcb", texmod_pf)
texmod_pf = ifelse(x1534 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pfl", texmod_pf)
texmod_pf = ifelse(x1534 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pst", texmod_pf)
texmod_pf = ifelse(x1534 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pby", texmod_pf)
# 35-59%
x3559 = sum_pf >= 35 & sum_pf < 60 & sum_nopf <= 0
texmod_pf = ifelse(x3559 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgrv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcnv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcbv", texmod_pf)
texmod_pf = ifelse(x3559 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pflv", texmod_pf)
texmod_pf = ifelse(x3559 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pstv", texmod_pf)
texmod_pf = ifelse(x3559 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pbyv", texmod_pf)
# 60-89%
x6089 = sum_pf >= 60 & sum_pf < 90 & sum_nopf <= 0
texmod_pf = ifelse(x6089 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgrx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcnx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcbx", texmod_pf)
texmod_pf = ifelse(x6089 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pflx", texmod_pf)
texmod_pf = ifelse(x6089 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pstx", texmod_pf)
texmod_pf = ifelse(x6089 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pbyx", texmod_pf)
# 90-100%
x90 = sum_pf >= 90 & sum_nopf <= 0
lieutex_pf = ifelse(x90 & x_pgr_pby & pgr >= pcn & pgr > 0, "pgr", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pgr_pby & pgr < pcn & pcn > 0 & is.na(texmod_pf), "pcn", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pcb_pby & pcb >= pfl & pcb > 0 & is.na(texmod_pf), "pcb", lieutex_pf)
lieutex_pf = ifelse(x90 & x_pcb_pby & pcb < pfl & pfl > 0 & is.na(texmod_pf), "pfl", lieutex_pf)
lieutex_pf = ifelse(x90 & pst >= 1.5 * pby & pst > 0 & is.na(texmod_pf), "pst", lieutex_pf)
lieutex_pf = ifelse(x90 & pst < 1.5 * pby & pby > 0 & is.na(texmod_pf), "pby", lieutex_pf)
# f
texmod_f = NA_character_
lieutex_f = NA_character_
x_fgr_fby = fgr_fcn >= ((1.5 * fcb_ffl) + (2 * fst) + (2.5 * fby))
x_fcb_fby = fcb_ffl >= ((1.5 * fst) + (2 * fby))
x1534 = sum_f >= 15 & sum_f < 35 & sum_nopf < 15
texmod_f = ifelse(x1534 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgr", texmod_f)
texmod_f = ifelse(x1534 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcn", texmod_f)
texmod_f = ifelse(x1534 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcb", texmod_f)
texmod_f = ifelse(x1534 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pfl", texmod_f)
texmod_f = ifelse(x1534 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pst", texmod_f)
texmod_f = ifelse(x1534 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pby", texmod_f)
# 35-59%
x3559 = sum_f >= 35 & sum_f < 60 & sum_nopf < 15
texmod_f = ifelse(x3559 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgrv", texmod_f)
texmod_f = ifelse(x3559 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcnv", texmod_f)
texmod_f = ifelse(x3559 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcbv", texmod_f)
texmod_f = ifelse(x3559 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pflv", texmod_f)
texmod_f = ifelse(x3559 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pstv", texmod_f)
texmod_f = ifelse(x3559 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pbyv", texmod_f)
# 60-89%
x6089 = sum_f >= 60 & sum_f < 90 & sum_nopf < 15
texmod_f = ifelse(x6089 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgrx", texmod_f)
texmod_f = ifelse(x6089 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcnx", texmod_f)
texmod_f = ifelse(x6089 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcbx", texmod_f)
texmod_f = ifelse(x6089 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pflx", texmod_f)
texmod_f = ifelse(x6089 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pstx", texmod_f)
texmod_f = ifelse(x6089 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pbyx", texmod_f)
# 90-100%
x90 = sum_f >= 90 & sum_nopf < 15
lieutex_f = ifelse(x90 & x_fgr_fby & fgr >= fcn & fgr > 0, "pgr", lieutex_f)
lieutex_f = ifelse(x90 & x_fgr_fby & fgr < fcn & fcn > 0 & is.na(texmod_f), "pcn", lieutex_f)
lieutex_f = ifelse(x90 & x_fcb_fby & fcb >= ffl & fcb > 0 & is.na(texmod_f), "pcb", lieutex_f)
lieutex_f = ifelse(x90 & x_fcb_fby & fcb < ffl & ffl > 0 & is.na(texmod_f), "pfl", lieutex_f)
lieutex_f = ifelse(x90 & fst >= 1.5 * fby & fst > 0 & is.na(texmod_f), "pst", lieutex_f)
lieutex_f = ifelse(x90 & fst < 1.5 * fby & fby > 0 & is.na(texmod_f), "pby", lieutex_f)
# combine texmods and lieutex ----
texmod = ifelse(is.na(texmod), texmod_f, texmod)
texmod = ifelse(is.na(texmod), texmod_pf, texmod)
lieutex = ifelse(is.na(lieutex), lieutex_f, lieutex)
lieutex = ifelse(is.na(lieutex), lieutex_pf, lieutex)
})
# idx <- 1:ncol(df_sub)
# df_mod <- df_sub
# df_mod[idx] <- lapply(idx, function(i) {
# texmod <- rep(NA_character_, nrow(df_sub))
# texmod[df_sub[, i] >= 15 & df_sub[, i] < 35] <- var_mods_sub[i]
# texmod[df_sub[, i] >= 35 & df_sub[, i] < 60] <- paste0(var_mods_sub[i], "v")
# texmod[df_sub[, i] >= 60 & df_sub[, i] < 90] <- paste0(var_mods_sub[i], "x")
# return(texmod)
# })
#
# df_lieu <- df_sub
# df_lieu[idx] <- lapply(idx, function(i) {
# lieutex <- rep(NA_character_, nrow(df_sub))
# lieutex[df_sub[i] >= 90] <- var_mods_sub[i]
# return(lieutex)
# })
# standardize output----
if (nrow(df_sub) > 0) {
df$texmod[!idx_sum] <- df_sub$texmod
df$lieutex[!idx_sum] <- df_sub$lieutex
} else df[c("texmod", "lieutex")] <- NA_character_
df <- df[c("texmod", "lieutex")]
if (as.is == FALSE) {
lv <- paste0(
c(rep(var_mods[(7:12)], times = 3, length.out = 18),
rep(var_mods[(1:6)], times = 3, length.out = 18)
),
rep(c("", "v", "x"), each = 6, length.out = 36)
)
tn <- names(sort(table(df$texmod), decreasing = TRUE))
lv <- c(lv, tn[! tn %in% lv])
df$texmod <- factor(df$texmod, levels = lv)
df$lieutex <- factor(df$lieutex, levels = var_mods[1:6])
}
if (as.is == FALSE & droplevels == TRUE) {
df <- droplevels.data.frame(df)
}
return(df)
}
# format function inputs ----
# format vector argument inputs into a data.frame
.format_inputs <- function(vars_l = NULL, var_cols, mis_value = NA, datatype = NULL) {
# check inputs----
vars_null <- sapply(vars_l, is.null)
vars_len <- sapply(vars_l[!vars_null], length)
nm <- var_cols
n <- length(var_cols)
df <- NULL
## all inputs are NULL
if (all(vars_null)) {
stop("all inputs are NULL")
}
## length of inputs
if (!all(vars_len == max(vars_len))) {
stop("length of inputs do not match")
}
# standardize inputs
# create data.frame
df <- as.data.frame(vars_l[!vars_null])
# append missing columns
idx <- nm %in% names(df)
df_mis <- as.data.frame(matrix(data = mis_value, nrow = nrow(df), ncol = sum(!idx)))
names(df_mis) <- nm[!idx]
df <- cbind(df_mis, df)[nm]
## correct datatype
df[1:n] <- lapply(df, function(x) if (is.factor(x)) as.character(x) else x)
df[1:n] <- lapply(df, function(x) if (is.integer(x)) as.numeric(x) else x)
dt <- unlist(lapply(df, class))
## check datatype
idx <- dt == datatype | dt == "logical"
if (!all(idx)) {
stop(paste("unexpected input data types, must be", datatype))
}
return(df)
}
#' @title Ranking Systems for USDA Taxonomic Particle-Size and Substitute Classes of Mineral Soils
#'
#' @description Generate a lookup table of USDA Particle-Size and Substitute Classes names, ranked according to approximate particle size
#'
#' @references \href{https://nrcspad.sc.egov.usda.gov/DistributionCenter/product.aspx?ProductID=991}{Field Book for Describing and Sampling Soils, version 3.0}
#'
#' @return A data.frame with a rank column, taxonomic family particle size class, and a flag for contrasting.
#'
#' @author Stephen Roecker
#'
#' @seealso [hz_to_taxpartsize()], [texture_to_taxpartsize()], [SoilTextureLevels()]
#'
#' @export
#' @examples
#'
#' # class codes
#' lu <- lookup_taxpartsize()
#'
#' idx <- lu$contrasting == FALSE
#'
#' lu$taxpartsize[idx]
#'
#' lu$rank[as.integer(lu$taxpartsize)[idx]]
#'
lookup_taxpartsize <- function() {
fe <- c("diatomaceous", "very-fine", "clayey", "fine", "hydrous", "fine-silty",
"fine-gypseous", "fine-loamy", "medial", "loamy", "coarse-loamy",
"coarse-silty", "coarse-gypseous", "ashy", "sandy", "hydrous-pumiceous",
"medial-pumiceous", "ashy-pumiceous", "clayey-skeletal", "hydrous-skeletal",
"medial-skeletal", "loamy-skeletal", "gypseous-skeletal", "ashy-skeletal",
"sandy-skeletal", "pumiceous", "cindery", "fragmental")
rank <- c(84, 74, 60.02, 46.04, 44.04, 26, 25.8, 25.6, 24, 17.24, 8.88,
8.5, 7.5, 6.5, 4.67, -55.96, -76, -93.5, -43.33, -55.96, -76,
-83.23, -83.35, -93.5, -95.33, -95.83, -96.33, -98.94)
names(rank) <- fe
# cf <- c("fragmental", "sandy-skeletal", "loamy-skeletal", "clay-skeletal")
test <- strsplit(.pscs_sc, " over | or ")
names(test) <- .pscs_sc
idx <- lapply(test, function(x) {
idx <- sapply(x, function(y) rank[which(fe == y)]) |> unlist()
# select the 3rd value when "or" results in 3 values
if (length(idx) > 2) idx <- c(idx[1], idx[3])
dif <- diff(idx)
idx <- idx[1] + sqrt(abs(dif)) * sign(dif)
# l <- dplyr::lag(idx)
# l <- idx[c(NA, 1:(length(idx) - 1))]
# idx <- ifelse(idx < l & !is.na(l), idx * -1, idx * 1)
# n <- length(idx)
# idx <- sum(idx * c(1, 0.1, 0.01, 0.001)[1:n])
return(idx)
})
idx <- round(unlist(idx), 1)
# fe <- data.frame(rn = 1:length(fe), fe = fe)
fe_df <- data.frame(rank = unname(rank), taxpartsize = fe)
sc_df <- data.frame(rank = unname(idx), taxpartsize = .pscs_sc)
lu <- rbind(fe_df, sc_df)
lu <- lu[order(-lu$rank), ]
# lu$rn <- 1:nrow(lu)
lu$contrasting <- grepl(" over ", lu$taxpartsize)
lu$taxpartsize <- factor(lu$taxpartsize, levels = lu$taxpartsize, ordered = TRUE)
row.names(lu) <- NULL
return(lu)
}
.pscs_sc <- c("Ashy over clayey", "Ashy over clayey-skeletal", "Ashy over loamy", "Ashy over loamy-skeletal", "Ashy over medial", "Ashy over medial-skeletal", "Ashy over pumiceous or cindery", "Ashy over sandy or sandy-skeletal", "Ashy-skeletal over clayey", "Ashy-skeletal over fragmental or cindery", "Ashy-skeletal over loamy-skeletal", "Ashy-skeletal over sandy or sandy-skeletal", "Cindery over loamy", "Cindery over medial", "Cindery over medial-skeletal", "Clayey over coarse-gypseous", "Clayey over fine-gypseous", "Clayey over fragmental", "Clayey over gypseous-skeletal", "Clayey over loamy", "Clayey over loamy-skeletal", "Clayey over sandy or sandy-skeletal", "Clayey-skeletal over sandy or sandy-skeletal", "Coarse-loamy over clayey", "Coarse-loamy over fragmental", "Coarse-loamy over sandy or sandy-skeletal", "Coarse-silty over clayey", "Coarse-silty over sandy or sandy-skeletal", "Fine-loamy over clayey", "Fine-loamy over fragmental", "Fine-loamy over sandy or sandy-skeletal", "Fine-silty over clayey", "Fine-silty over fragmental", "Fine-silty over sandy or sandy-skeletal", "Hydrous over clayey", "Hydrous over clayey-skeletal", "Hydrous over fragmental", "Hydrous over loamy", "Hydrous over loamy-skeletal", "Hydrous over sandy or sandy-skeletal", "Loamy over ashy or ashy-pumiceous", "Loamy over coarse-gypseous", "Loamy over fine-gypseous", "Loamy over pumiceous or cindery", "Loamy over sandy or sandy-skeletal", "Loamy-skeletal over cindery", "Loamy-skeletal over clayey", "Loamy-skeletal over fragmental", "Loamy-skeletal over gypseous-skeletal", "Loamy-skeletal over sandy or sandy-skeletal", "Medial over ashy", "Medial over ashy-pumiceous or ashy-skeletal", "Medial over clayey", "Medial over clayey-skeletal", "Medial over fragmental", "Medial over hydrous", "Medial over loamy", "Medial over loamy-skeletal", "Medial over pumiceous or cindery", "Medial over sandy or sandy-skeletal", "Medial-skeletal over fragmental or cindery", "Medial-skeletal over loamy-skeletal", "Medial-skeletal over sandy or sandy-skeletal", "Pumiceous or ashy-pumiceous over loamy", "Pumiceous or ashy-pumiceous over loamy-skeletal", "Pumiceous or ashy-pumiceous over medial", "Pumiceous or ashy-pumiceous over medial-skeletal", "Pumiceous or ashy-pumiceous over sandy or sandy-skeletal", "Sandy over clayey", "Sandy over loamy", "Sandy-skeletal over loamy") |>
tolower()
Try the aqp package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.