Merge branch 'master' into suggest-sp

CITATION.cff

@@ -19,7 +19,7 @@ abstract: The Algorithms for Quantitative Pedology (AQP) project was started in
   the scientist to focus on ideas rather than boilerplate data processing tasks <doi:10.1016/j.cageo.2012.10.020>.
   These functions and data structures have been extensively tested and documented,
   applied to projects involving hundreds of thousands of soil profiles, and deeply
-  integrated into widely used tools such as SoilWeb <https://casoilresource.lawr.ucdavis.edu/soilweb-apps/>.
+  integrated into widely used tools such as SoilWeb <https://casoilresource.lawr.ucdavis.edu/soilweb-apps>.
   Components of the AQP project (aqp, soilDB, sharpshootR, soilReports packages) serve
   an important role in routine data analysis within the USDA-NRCS Soil Science Division.
   The AQP suite of R packages offer a convenient platform for bridging the gap between

DESCRIPTION

@@ -7,7 +7,7 @@ Maintainer: Dylan Beaudette <[email protected]>
 Depends: R (>= 3.5.0)
 Imports: grDevices, graphics, stats, utils, methods, grid, lattice, cluster, stringr, data.table, farver
 Suggests: mvtnorm, colorspace, ape, soilDB, sp, sf, latticeExtra, tactile, compositions, sharpshootR, markovchain, xtable, testthat, Gmedian, Hmisc, tibble, RColorBrewer, scales, digest, MASS, mpspline2, soiltexture, gower, knitr, rmarkdown, plyr
-Description: The Algorithms for Quantitative Pedology (AQP) project was started in 2009 to organize a loosely-related set of concepts and source code on the topic of soil profile visualization, aggregation, and classification into this package (aqp). Over the past 8 years, the project has grown into a suite of related R packages that enhance and simplify the quantitative analysis of soil profile data. Central to the AQP project is a new vocabulary of specialized functions and data structures that can accommodate the inherent complexity of soil profile information; freeing the scientist to focus on ideas rather than boilerplate data processing tasks <doi:10.1016/j.cageo.2012.10.020>. These functions and data structures have been extensively tested and documented, applied to projects involving hundreds of thousands of soil profiles, and deeply integrated into widely used tools such as SoilWeb <https://casoilresource.lawr.ucdavis.edu/soilweb-apps/>. Components of the AQP project (aqp, soilDB, sharpshootR, soilReports packages) serve an important role in routine data analysis within the USDA-NRCS Soil Science Division. The AQP suite of R packages offer a convenient platform for bridging the gap between pedometric theory and practice.
+Description: The Algorithms for Quantitative Pedology (AQP) project was started in 2009 to organize a loosely-related set of concepts and source code on the topic of soil profile visualization, aggregation, and classification into this package (aqp). Over the past 8 years, the project has grown into a suite of related R packages that enhance and simplify the quantitative analysis of soil profile data. Central to the AQP project is a new vocabulary of specialized functions and data structures that can accommodate the inherent complexity of soil profile information; freeing the scientist to focus on ideas rather than boilerplate data processing tasks <doi:10.1016/j.cageo.2012.10.020>. These functions and data structures have been extensively tested and documented, applied to projects involving hundreds of thousands of soil profiles, and deeply integrated into widely used tools such as SoilWeb <https://casoilresource.lawr.ucdavis.edu/soilweb-apps>. Components of the AQP project (aqp, soilDB, sharpshootR, soilReports packages) serve an important role in routine data analysis within the USDA-NRCS Soil Science Division. The AQP suite of R packages offer a convenient platform for bridging the gap between pedometric theory and practice.
 License: GPL (>= 3)
 LazyLoad: yes
 Repository: CRAN

NAMESPACE

@@ -3,6 +3,7 @@
 export(HzDepthLogicSubset)
 export(L1_profiles)
 export(NCSP)
+export(PSCS_levels)
 export(ReactionClassLevels)
 export(SANN_1D)
 export(SoilProfileCollection)
@@ -95,6 +96,11 @@ export(hzOffset)
 export(hzTopographyCodeToLineType)
 export(hzTopographyCodeToOffset)
 export(hzTransitionProbabilities)
+export(hz_dissolve)
+export(hz_intersect)
+export(hz_lag)
+export(hz_segment)
+export(hz_to_taxpartsize)
 export(invertLabelColor)
 export(lunique)
 export(maxDepthOf)

NEWS.md

@@ -1,8 +1,10 @@
-# aqp 2.0.3 (2024-03-20)
+# aqp 2.0.3 (2024-04-18)
+ * CRAN release
  * `simulateColor()` gains new method `mvnorm` for simulating plausible colors
    - package mvtnorm added to SUGGESTS
  * performance improvements in `profileInformationIndex()`, `dice()`, `slab()`, `spc2mpspline()`, `fillHzGaps()`, and `flagOverlappingHz()`
  * aesthetic improvements in `huePositionCircle()`
+ * new function `thicknessOf()` used for calculating thickness of horizons within each profile of a `SoilProfileCollection` based on horizon-level logical expressions encoded in a function. Default behavior uses pattern matching on the horizon designation name.
 
 # aqp 2.0.2 (2023-11-18)
  * CRAN release
@@ -600,7 +602,7 @@ Incremental changes, should have no effect on previous code:
 
 # aqp 1.0 (2012-03-26)
    * 1.0 release, still missing condensed vignettes- should be ready soon
-   * see http://casoilresource.lawr.ucdavis.edu/drupal/taxonomy/term/56 for samples
+   * see https://casoilresource.lawr.ucdavis.edu/drupal/taxonomy/term/56 for samples
    * A small bug in profile_compare() was fixed, where slices were evaluated as 'soil' based on the bottom depth of the profile, and NOT on the presence of actual data. See ?profile_compare for details. This change will have a minor affect on profile comparisons in cases where Cr or R horizons (usually missing most characterization data) have been extended down to some arbitrary depth (usually 150 or 200 cm) AND a maximum depth of evaluation (max_d) was set beyond the actual depth of most profiles in the collection.
 
 # aqp 0.99-9.8 (2012-03-02)

R/Class-SoilProfileCollection.R

@@ -944,10 +944,12 @@ setMethod("site", signature(object = "SoilProfileCollection"),
 setGeneric("diagnostic_hz", function(object, ...)
   standardGeneric("diagnostic_hz"))
 
-#' Retrieve diagnostic data from SoilProfileCollection
-#'
-#' @description Get diagnostic feature data from SoilProfileCollection. Result is returned in the same \code{data.frame} class used to initially construct the SoilProfileCollection.
+#' Get or Set Diagnostic Horizon data in a SoilProfileCollection
 #'
+#' @description Diagnostic horizons describe features of the soil relevant to taxonomic classification. A single profile may have multiple diagnostic features or horizons, each of which may be comprised of multiple horizons.
+#' 
+#'  - `diagnostic_hz()` (get method): Get diagnostic feature data from a SoilProfileCollection. 
+#'  
 #' @param object a SoilProfileCollection
 #'
 #' @docType methods
@@ -964,9 +966,11 @@ setMethod(f = 'diagnostic_hz', signature(object = 'SoilProfileCollection'),
 setGeneric("restrictions", function(object, ...)
   standardGeneric("restrictions"))
 
-#' Retrieve restriction data from SoilProfileCollection
+#' Get or Set Restriction data in a SoilProfileCollection
 #'
-#' @description Get restriction data from SoilProfileCollection. Result is returned in the same \code{data.frame} class used to initially construct the SoilProfileCollection.
+#' @description Restrictions describe root-limiting features in the soil. A single profile may have multiple restrictions.
+#' 
+#'  - `restrictions()` (get method): Get restriction data from a SoilProfileCollection. 
 #'
 #' @param object a SoilProfileCollection
 #' @docType methods

R/SoilProfileCollection-setters.R

@@ -629,19 +629,19 @@ setReplaceMethod("horizons", signature(object = "SoilProfileCollection"),
 setGeneric('diagnostic_hz<-', function(object, value)
   standardGeneric('diagnostic_hz<-'))
 
-#' Add Data to Diagnostic Features Slot
-#'
 #' @name diagnostic_hz<-
 #'
-#' @description Diagnostic feature data in an object inheriting from \code{data.frame} can easily be added via merge (LEFT JOIN). There must be one or more same-named columns containing profile ID on the left and right hand side to facilitate the join: \code{diagnostic_hz(spc) <- newdata}
-#'
+#' @description 
+#' 
+#'  - `diagnostic_hz<-` (set method): Set diagnostic feature data for a SoilProfileCollection. The profile ID column from `object` (`idname(object)`) must be present in the replacement `value` object.
+#' 
 #' @param object A SoilProfileCollection
-#' @param value An object inheriting \code{data.frame}
+#' @param value An object inheriting from \code{data.frame}
 #'
 #' @aliases diagnostic_hz<-,SoilProfileCollection-method
 #' @docType methods
 #' @export
-#' @rdname diagnostic_hz-set
+#' @rdname diagnostic_hz
 #'
 #' @examples
 #'
@@ -716,19 +716,18 @@ setReplaceMethod("diagnostic_hz",
 setGeneric('restrictions<-', function(object, value)
   standardGeneric('restrictions<-'))
 
-#' Add Data to Restrictions Slot
-#'
 #' @name restrictions<-
 #'
-#' @description Restrictions data in an object inheriting from \code{data.frame} can easily be added via merge (LEFT JOIN). There must be one or more same-named profile ID columns on the left and right hand side to facilitate the join: \code{restrictions(spc) <- newdata}.
-#'
+#' @description 
+#' 
+#'  - `restrictions<-` (set method): Set restriction data for a SoilProfileCollection. The profile ID column from `object` (`idname(object)`) must be present in the replacement `value` object.
 #' @param object A SoilProfileCollection
-#' @param value An object inheriting \code{data.frame}
+#' @param value An data.frame object containing at least a column with name `idname(object)`
 #'
 #' @aliases restrictions<-,SoilProfileCollection-method
 #' @docType methods
 #'
-#' @rdname restrictions-set
+#' @rdname restrictions
 #' @export
 #' @examples
 #'

R/allocate.R

@@ -60,6 +60,8 @@
 #'
 #' @return A vector or \code{data.frame} object.
 #' 
+#' @author Stephen Roecker
+#' 
 #' @references 
 #' Abrol, I., Yadav, J. & Massoud, F. 1988. \href{https://www.fao.org/3/x5871e/x5871e00.htm}{Salt-affected soils and their management}. No. Bulletin 39. Rome, FAO Soils.
 #' 
@@ -413,7 +415,7 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag
 
   # combine results and subset to 0-25cm
   df_bs  <- cbind(df[vars2[1:3]], BS1 = bs1, BS2 = bs2)
-  df_bs  <- segment(df_bs, intervals = c(0, 25), hzdepcols = c("hztop", "hzbot"))
+  df_bs  <- hz_segment(df_bs, intervals = c(0, 25), depthcols = c("hztop", "hzbot"))
   df_bs  <- df_bs[df_bs$segment_id == "00-25", -6]
 
   # aggregate the horizons
@@ -635,3 +637,225 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag
   return(sp)
 }
 
+
+#' @title Allocate Particle Size Class for the Control Section.
+#' 
+#' @description This function aggregates information in the horizon table and allocates it to the particle size class for the control section.
+#'
+#' @param x a \code{data.frame} containing the original horizon table.
+#' @param y a \code{data.frame} containing the particle size control section depths for each idcol.
+#' @param taxpartsize \code{character} column name for taxonomic family particle size class.
+#' @param clay \code{character} column name for clay percent.
+# #' @param frags \code{character} column name for total rock fragments percent.
+#' @param idcol character: column name of the pedon ID within the object.
+#' @param depthcols a character vector of length 2 specifying the names of the horizon depths (e.g. `c("top", "bottom")`).
+#' 
+#' 
+#' @details
+#' This function differs from \code{\link{texture_to_taxpartsize}} in that is aggregates the results of \code{\link{texture_to_taxpartsize}}, and accounts for strongly contrasting particle size classes. 
+#' 
+#' 
+#'
+#' @return A \code{data.frame} object containing the original idcol and aggregated particle size control section allocation.
+#' 
+#' @author Stephen Roecker
+#' 
+#' @seealso [texture_to_taxpartsize()], [PSCS_levels()]
+#' 
+#' @export
+
+#' @examples
+#' 
+#' h <- data.frame(
+#' id = 1,
+#' hzname = c("A", "BA", "Bw", "BC", "C"),
+#' top    = c( 0, 10, 45, 60,  90),
+#' bottom = c(10, 45, 60, 90, 150),
+#' clay   = c(15, 16, 45, 20,  10),
+#' sand   = c(10, 35, 40, 50,  90),
+#' frags  = c( 0,  5, 10, 38,  40)
+#' )
+#'
+#' h <- cbind(
+#' h, 
+#' texcl = ssc_to_texcl(clay = h$clay, sand = h$sand)
+#' )
+#' 
+#' pscs <- data.frame(
+#' id = 1, 
+#' top = 25, 
+#' bottom = 100
+#' )
+#' 
+#' h <- cbind(
+#' h, 
+#' taxpartsize = texture_to_taxpartsize(
+#' texcl = h$texcl,
+#' clay  = h$clay, 
+#' sand  = h$sand,
+#' fragvoltot = h$frags
+#' ))
+#' 
+#' depths(h) <- id ~ top + bottom
+#' 
+#' pscs <- data.frame(id = h$id, rbind(estimatePSCS(h)))
+#' names(pscs)[2:3] <- c("top", "bottom")
+#' 
+#' hz_to_taxpartsize(horizons(h), pscs)
+#'  
+#' 
+hz_to_taxpartsize <- function(x, y, taxpartsize = "taxpartsize", clay = "clay", idcol = "id", depthcols = c("top", "bottom")) {
+  # need to incorporate fine sand for special cases of strongly contrasting classes and rock fragments (?)
+  # frags = "frags", 
+
+  # strongly contrasting
+
+  x$rn <- 1:nrow(x)
+  # xy <- hz_intersect(x, y, idcol = idcol, depthcols = depthcols)
+  # x_sub <- x[x$rn %in% xy$rn, ]
+
+
+  # check segment_id ----
+  ## if it exists, overwrite it
+  x_nm <- names(x)
+  y_nm <- names(y)
+  if (any(x_nm == "segment_id") | any(y_nm == "segment_id")) {
+    x[x_nm == "segment_id"] <- NULL
+    y[y_nm == "segment_id"] <- NULL
+  }
+
+
+  # check dissolve_id ----
+  ## if it exists, overwrite it
+  x_nm <- names(x)
+  y_nm <- names(y)
+  if (any(x_nm == "dissolve_id") | any(y_nm == "dissolve_id")) {
+    x[x_nm == "dissolve_id"] <- NULL
+    y[y_nm == "dissolve_id"] <- NULL
+  }
+
+
+  # standardize inputs
+  x_std <- .standardize_inputs(x, idcol = idcol, depthcols = depthcols, clay = clay, taxpartsize = taxpartsize)
+  x <- x_std$x; x_conv <- x_std$x_conversion
+  x_std <- NULL
+
+  y <- .standardize_inputs(y, idcol = idcol, depthcols = depthcols)$x
+
+
+  # dissolve on pscs
+  # calculate non-trimmed horizon thickness
+  x_dis <- x |>
+    hz_dissolve(by = "taxpartsize", idcol = "idcol", depthcols = c("top", "bot")) |>
+    transform(thk_o = bot - top)
+
+
+  # trim depths
+  # calculate trimmed horizon thickness
+  xy_dis <- x_dis |>
+    hz_intersect(y, idcol = "idcol", depthcols = c("top", "bot")) |>
+    transform(thk_t = bot - top)
+
+
+  # rejoin dissolved pscs to the original horizon table 
+  xy <- hz_intersect(x, xy_dis, idcol = "idcol", depthcols = c("top", "bot")) |> suppressWarnings()
+  x_dis  <- NULL
+  xy_dis <- NULL
+
+
+  # aggregate clay values within dissolved pscs
+  top <- NULL
+  bot <- NULL
+  thk_o <- NULL
+  thk_t <- NULL
+  clay_wt <- NULL
+
+  xy_agg <- data.table::as.data.table(xy)[,
+    list(
+      top     = min(top, na.rm = TRUE),
+      bot     = max(bot, na.rm = TRUE),
+      clay_wt = weighted.mean(clay, w = thk_t,  na.rm = TRUE),
+      # need to impute frags
+      # frag_wt = weighted.mean(total_frags_pct_nopf, w = thk_t), na.rm = TRUE,
+      thk_o   = sum(thk_o, na.rm = TRUE),
+      thk_t   = sum(thk_t, na.rm = TRUE)
+    ), by = c("idcol", "taxpartsize", "dissolve_id")
+    ]
+  data.table::setorder(xy_agg, idcol, top)
+  xy_agg <- as.data.frame(xy_agg)
+
+
+  # find adjacent horizons
+  xy_lag <- xy_agg |> 
+    hz_lag(idcol = "idcol", depthcols = c("top", "bot"))
+
+
+  # replace special cases of pscs with strongly contrasting classes
+  clay_wt_bot.1     <- NULL
+  taxpartsize_bot.1 <- NULL
+
+  xy_agg <- cbind(xy_agg, xy_lag) |> 
+    within({
+      clay_dif = clay_wt_bot.1 - clay_wt
+      sc       = paste0(taxpartsize, " over ", taxpartsize_bot.1)
+      sc       = gsub(" over NA$", "", sc)
+
+      sc = gsub("^fine over|^very-fine over", "clayey over", sc)
+      sc = gsub("over fine$|over very-fine$", "over clayey", sc)
+      sc = gsub("over fine over|over very-fine over", "over clayey over", sc)
+      sc = gsub("over sandy|over sandy-skeletal", "over sandy or sandy-skeletal", sc)
+
+      idx_sc     = sc %in% .pscs_sc
+      sc[idx_sc] = sc[idx_sc]
+    })
+  xy_lag <- NULL
+
+
+  # find multiple strongly contrasting ps classes within the control section
+  n_sc    <- NULL
+  n_peiid <- NULL
+
+  test <- data.table::as.data.table(xy_agg)[, list(
+    n_sc    = sum(grepl("over", sc), na.rm = TRUE),
+    n_peiid = length(idcol)
+  ),
+  by = "idcol"
+  ] |>
+    as.data.frame()
+
+
+  # pick the sc pscs with the largest contrast or pscs with the greatest thickness
+  xy_res <- xy_agg |>
+    merge(test, by = "idcol", all.x = TRUE, sort = FALSE) |>
+    transform(
+      idx_sc = grepl(" over ", sc)
+    )
+
+  xy_res <- data.table::as.data.table(xy_res)[, list(
+    pscs1 = sc[n_sc == 0 & n_peiid == 1],
+    pscs2 = sc[n_sc == 1 & n_peiid  > 1 & idx_sc],
+    pscs3 = sc[which.max(thk_t[n_sc == 0 & n_peiid == 2])],
+    pscs4 = sc[n_sc == 1 & idx_sc],
+    pscs5 = sc[which.max(abs(clay_dif[n_sc > 1 & !is.na(sc)]))]
+    ),
+    by =  "idcol"
+    ] |>
+    as.data.frame() |>
+    within({
+      # need to add fix for special case of sandy over loamy which requires fine sand percent
+      taxpartsize = paste(pscs1, pscs3, pscs4, pscs5, sep = "")
+      taxpartsize = gsub("NA", "", taxpartsize)
+      pscs1 = NULL
+      pscs2 = NULL
+      pscs3 = NULL
+      pscs4 = NULL
+      pscs5 = NULL
+    })
+
+
+  # reset inputs
+  xy_res <- .reset_inputs(xy_res, x_conv[1])
+
+
+  return(xy_res)
+}