texcl_to_ssc() option to derive low/high values from input texture class

[brownag]

It would be useful if texcl_to_ssc() could optionally return a data.frame with class limits for sand, silt, clay given a texture class. This would help provide context for the estimated central / "RV" value, and the effect of optional specification of a custom clay content.

I have not investigated the easiest way to bolt this into the current texcl_to_ssc() code, but I note that logic from the NASIS validation used for ssc_to_texcl() has the required information, but not yet in a "reusable" format:

aqp/R/texture.R

Lines 298 to 312 in 37d6e81

	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"
	})

[brownag]

Some time before these functions were in AQP, I used this type of logic in some of my QC/validation scripts, which were a precursor to a bunch of stuff that ultimately was added to this package. Logic derived directly from NASIS calculation "Textural Class versus Particle Size Separates" written by Cathy Seybold (version updated 4/07/14)

For instance:
https://github.com/brownag/sANDREWbox/blob/27e25464921bfdcd2c45f3df4ebbb74056f1da98/PedonSummaryFunctions/fine_earth_fractions.R#L215-L289

[brownag]

Perhaps this would be better suited to a new function.
Here is a prototype that doesn't quite generalize all of the above referenced logic for re-use, but gives what I have in mind as desired output. I think that perhaps there are some opportunities to improve efficiency, but I have a few different use cases that I want covered, so starting with expected output then can optimize further.

#' Convert Texture Class to Class Limits
#'
#' @param x _character_ or _list_. A vector of texture class codes (e.g. `"l"` for "loam", `"sicl"` for silty clay loam) without texture class modifiers. If the input is a list, multiple texture classes within each list element are aggregated to create combined class limits.
#' @details Logic derived from NASIS validation NASIS calculation "Textural Class versus Particle Size Separates" written by Cathy Seybold (last updated 4/07/14)
#' 
#' @return A _data.frame_ with column names "texcl", "clay_l", "clay_m", "clay_h", "sand_l", "sand_m", "sand_h", "silt_l", "silt_m", "silt_h"
#' 
#' @export
#'
#' @examples
#' 
#' texcl_to_classlimit(c("l", "sicl", "cl"))
#' 
#' texcl_to_classlimit(list(c("l", "sicl", "cl")))
#' 
texcl_to_classlimit <- function(x) {
  x <- lapply(x, function(y) tolower(trimws(y)))
  xout <- sapply(x, paste0, collapse = ",")
  res <- as.data.frame(t(sapply(x, function(y) {
    xcl <- as.integer(factor(y, levels = SoilTextureLevels()))
    data.frame(
      clay_l = .getClayLow(xcl),
      clay_h = .getClayHigh(xcl),
      sand_l = .getSandLow(xcl),
      sand_h = .getSandHigh(xcl),
      silt_l = .getSiltLow(xcl),
      silt_h = .getSiltHigh(xcl)
    )
  })))
  res[] <- lapply(res, unlist)
  res$clay_m <- apply(res[,c("clay_l", "clay_h")], MARGIN = 1, mean)
  res$sand_m <- apply(res[,c("sand_l", "sand_h")], MARGIN = 1, mean)
  res$silt_m <- apply(res[,c("silt_l", "silt_h")], MARGIN = 1, mean)
  res <- res[c("clay_l", "clay_m", "clay_h",
               "sand_l", "sand_m", "sand_h",
               "silt_l", "silt_m", "silt_h")]
  res <- cbind(data.frame(texcl = xout), res)
  rownames(res) <- NULL
  res
}

####
# LOGIC DERIVED FROM TEXTURAL CLASS VERSUS PARTICLE SIZE SEPARATES
####
.getClayHigh <- function(texcl) {
  if (any(texcl == 21)) return(100) #clay
  else if (any(texcl == 20)) return(60)  #sic
  else if (any(texcl == 19)) return(55)  #sc
  else if (any(texcl == 18 | texcl == 17)) return(40)  #cl, sicl
  else if (any(texcl == 16)) return(35)  #scl
  else if (any(texcl == 14 | texcl == 13)) return(27)  #l, sil
  else if (any(texcl == 12 | texcl == 11 | texcl == 10 | texcl == 9)) return(20)  #sl, fsl, vfsl, cosl
  else if (any(texcl == 8 | texcl == 7 | texcl == 6 | texcl == 5)) return(15)  #ls, lvfs, lfs, lcos
  else if (any(texcl == 15)) return(12) #si
  else if (any(texcl == 4 | texcl == 3 | texcl == 2 | texcl == 1)) return(10)
  return(NA)
}

.getClayLow <- function(texcl) {
  if (any(texcl <= 12 & texcl >= 1) | any(texcl == 14 | texcl == 15)) return(0)
  else if (any(texcl == 13)) return(7)
  else if (any(texcl == 16)) return(20)
  else if (any(texcl == 17 | texcl == 18)) return(27)
  else if (any(texcl == 19)) return(35)
  else if (any(texcl == 20 | texcl == 21)) return(40)
  return(NA)
}

.getSiltHigh <- function(texcl) {
  if (any(texcl == 15)) return(100)
  else if (any(texcl == 14)) return(88)
  else if (any(texcl == 18)) return(73)
  else if (any(texcl == 20)) return(60)
  else if (any(texcl == 17)) return(53)
  else if (any(texcl <= 13 & texcl >= 9)) return(50)
  else if (any(texcl == 21)) return(40)
  else if (any(texcl <= 8 & texcl >= 5)) return(30)
  else if (any(texcl == 16)) return(28)
  else if (any(texcl == 19)) return(20)
  else if (any(texcl == 4 | texcl == 3 | texcl == 2 | texcl == 1)) return(15)
  return(NA)
}

.getSiltLow <- function(texcl) {
  if (any(texcl >= 1 & texcl <= 12) | any(texcl == 16 | texcl == 19 | texcl == 21)) return(0)
  else if (any(texcl == 17)) return(15)
  else if (any(texcl == 13)) return(28)
  else if (any(texcl == 20 | texcl == 18)) return(40)
  else if (any(texcl == 14)) return(50)
  else if (any(texcl == 15)) return(80)
  return(NA)
}

.getSandHigh <- function(texcl) {
  if (any(texcl == 4 | texcl == 3 | texcl == 2 | texcl == 1)) return(100)
  else if (any(texcl >= 5 & texcl <= 8)) return(90)
  else if (any(texcl >= 9 & texcl <= 12)) return(85)
  else if (any(texcl == 16)) return(80)
  else if (any(texcl == 19)) return(65)
  else if (any(texcl == 13)) return(52)
  else if (any(texcl == 14)) return(50)
  else if (any(texcl == 17 | texcl == 21)) return(45)
  else if (any(texcl == 15 | texcl == 18 | texcl == 20)) return(20)
  return(NA)
}

.getSandLow <- function(texcl) {
  if (any(texcl == 21 | texcl == 20 | texcl == 18 | texcl == 14 | texcl == 15)) return(0)
  else if (any(texcl == 17)) return(20)
  else if (any(texcl == 13)) return(23)
  else if (any(texcl >= 9 & texcl <= 12)) return(43)
  else if (any(texcl == 16 | texcl == 19)) return(45)
  else if (any(texcl >= 5 & texcl <= 8)) return(70)
  else if (any(texcl >= 1 & texcl <= 4)) return(85)
  return(NA)
}

library(aqp)
#> This is aqp 2.0.3

texcl_to_classlimit(c("l", "sicl", "cl"))
#>   texcl clay_l clay_m clay_h sand_l sand_m sand_h silt_l silt_m silt_h
#> 1     l      7   17.0     27     23   37.5     52     28   39.0     50
#> 2  sicl     27   33.5     40      0   10.0     20     40   56.5     73
#> 3    cl     27   33.5     40     20   32.5     45     15   34.0     53

texcl_to_classlimit(list(c("l", "sicl", "cl")))
#>       texcl clay_l clay_m clay_h sand_l sand_m sand_h silt_l silt_m silt_h
#> 1 l,sicl,cl      7   23.5     40      0     26     52     15     44     73

EDIT: fixed silt_m calculation

Note that the input could be a character vector or single texture classes, or a list of character vectors (~texture groups). In the latter case an aggregate range is given for the combined set of textures in each list element.

[brownag]

The above process can likely be made more efficient using a lookup table similar to the one generated below:

library(aqp)
#> This is aqp 2.0.3
texcl_to_classlimit(SoilTextureLevels())
#>    texcl clay_l clay_m clay_h sand_l sand_m sand_h silt_l silt_m silt_h
#> 1    cos      0    5.0     10     85   92.5    100      0    7.5     15
#> 2      s      0    5.0     10     85   92.5    100      0    7.5     15
#> 3     fs      0    5.0     10     85   92.5    100      0    7.5     15
#> 4    vfs      0    5.0     10     85   92.5    100      0    7.5     15
#> 5   lcos      0    7.5     15     70   80.0     90      0   15.0     30
#> 6     ls      0    7.5     15     70   80.0     90      0   15.0     30
#> 7    lfs      0    7.5     15     70   80.0     90      0   15.0     30
#> 8   lvfs      0    7.5     15     70   80.0     90      0   15.0     30
#> 9   cosl      0   10.0     20     43   64.0     85      0   25.0     50
#> 10    sl      0   10.0     20     43   64.0     85      0   25.0     50
#> 11   fsl      0   10.0     20     43   64.0     85      0   25.0     50
#> 12  vfsl      0   10.0     20     43   64.0     85      0   25.0     50
#> 13     l      7   17.0     27     23   37.5     52     28   39.0     50
#> 14   sil      0   13.5     27      0   25.0     50     50   69.0     88
#> 15    si      0    6.0     12      0   10.0     20     80   90.0    100
#> 16   scl     20   27.5     35     45   62.5     80      0   14.0     28
#> 17    cl     27   33.5     40     20   32.5     45     15   34.0     53
#> 18  sicl     27   33.5     40      0   10.0     20     40   56.5     73
#> 19    sc     35   45.0     55     45   55.0     65      0   10.0     20
#> 20   sic     40   50.0     60      0   10.0     20     40   50.0     60
#> 21     c     40   70.0    100      0   22.5     45      0   20.0     40

Another possible extension of the above function (and standard lookup table) would be to optionally generate limits for sand fractions, to differentiate the more detailed texture classes e.g. fs, vfs, cos

[dylanbeaudette]

This looks great, especially the possibility to pre-generate the results for more efficient operation. I'll save the simulation stuff for another time. There are still some issues that I'd like to better understand before bolting any more functionality onto what is here and working.

[brownag]

I am going to take a crack at pre-generating the results, and will consider incorporating into existing data structures rather than making something new. PR imminent