Skip to content

Instantly share code, notes, and snippets.

@brownag
Last active March 14, 2023 18:36
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save brownag/c14b66cc7b30819e8ab12b7803a1cce5 to your computer and use it in GitHub Desktop.
Save brownag/c14b66cc7b30819e8ab12b7803a1cce5 to your computer and use it in GitHub Desktop.
Example using `dplyr::case_match()` for assigning codes to groups of `taxpartsize` for a "similar soils" style key
library(aqp)
library(soilDB)
f <- fetchNASIS()
#> NOTE: some records are missing surface fragment cover
#> mixing dry colors ... [1 of 43 horizons]
#> mixing moist colors ... [1 of 43 horizons]
#> NOTE: some records are missing rock fragment volume
#> NOTE: all records are missing artifact volume

rate_taxpartsize <- function(f, 
                             taxpartsize = "taxpartsize", 
                             claytotest = "clay") {
  
  # encode PSCS to a factor using NASIS domain values (not necessary)
  # y <- soilDB::NASISChoiceList(f[[taxpartsize]], taxpartsize)
  
  y <- as.character(f[[taxpartsize]])
  z <- dplyr::case_match(y,
    "fragmental" ~ 1,
    c("sandy-skeletal", "ashy-skeletal") ~ 2,
    c("sandy", "ashy") ~ 3,
    c("coarse-loamy", "coarse-silty", "medial", "loamy") ~ 4,
    c("medial-skeletal", "loamy-skeletal") ~ 5,
    c("fine-loamy", "fine-silty") ~ 6,
    c("fine", "clayey") ~ 7,
    c("very fine") ~ 8
  )
  
  ## Custom handling of shallow soil PSCS
  p <- aqp::estimatePSCS(f)
  f$.claytotest <- f[[claytotest]]
  f$.pscs_clay <- aqp::mutate_profile(trunc(f, p$pscs_top, p$pscs_bottom),
                                 V1 = weighted.mean(.claytotest, hzdepb - hzdept))$V1
  
  # add 2 to loamy textures with >18% clay
  idx1 <- (f$.pscs_clay >= 18 & y == "loamy")
  z[idx1] <- z[idx1] + 2
  
  # add 1 to clayey textures with >60% clay
  idx2 <- (f$.pscs_clay >= 60 & y == "clayey")
  z[idx2] <- z[idx2] + 1
  z
}

f$taxpartsize_rating <- rate_taxpartsize(f)

data.frame(taxpartsize=f$taxpartsize, rating=f$taxpartsize_rating)
#>                                  taxpartsize rating
#> 1                             loamy-skeletal      5
#> 2                               coarse-loamy      4
#> 3  coarse-loamy over sandy or sandy-skeletal     NA
#> 4                               coarse-loamy      4
#> 5                             loamy-skeletal      5
#> 6                                      sandy      3
#> 7                             loamy-skeletal      5
#> 8                             loamy-skeletal      5
#> 9                               coarse-loamy      4
#> 10                              coarse-loamy      4
#> 11                            loamy-skeletal      5
#> 12                            sandy-skeletal      2
#> 13                              coarse-loamy      4
#> 14                              coarse-loamy      4
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment