diff --git a/CITATION.cff b/CITATION.cff index 197e64701..e83b2b488 100644 --- a/CITATION.cff +++ b/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 . 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 . + integrated into widely used tools such as SoilWeb . 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 diff --git a/DESCRIPTION b/DESCRIPTION index 4d2af2ea4..6541d1ef8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,7 +7,7 @@ Maintainer: Dylan Beaudette 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 . 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 . 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 . 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 . 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 diff --git a/NAMESPACE b/NAMESPACE index 92c3c18c9..20c7dbabc 100644 --- a/NAMESPACE +++ b/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) diff --git a/NEWS.md b/NEWS.md index bb1bb6413..3db71cddc 100644 --- a/NEWS.md +++ b/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) diff --git a/R/Class-SoilProfileCollection.R b/R/Class-SoilProfileCollection.R index 2874db2e2..a185d75ea 100644 --- a/R/Class-SoilProfileCollection.R +++ b/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 diff --git a/R/SoilProfileCollection-setters.R b/R/SoilProfileCollection-setters.R index 840fedb3e..f05bc454d 100644 --- a/R/SoilProfileCollection-setters.R +++ b/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 #' diff --git a/R/allocate.R b/R/allocate.R index d0b55e1c0..797367c7a 100644 --- a/R/allocate.R +++ b/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) +} diff --git a/R/aqp-label-placement-solvers.R b/R/aqp-label-placement-solvers.R index f2806403e..9875fe35d 100644 --- a/R/aqp-label-placement-solvers.R +++ b/R/aqp-label-placement-solvers.R @@ -101,6 +101,10 @@ overlapMetrics <- function(x, thresh) { +## TODO: incorporate ideas from N-body simulation +# https://en.wikipedia.org/wiki/N-body_simulation + + #' @title Simulation of electrostatic force between two charged particles #' @description This function computes a "force" (attraction or repulsion) between two charged "particles" (usually labels or other graphical elements), using a modification of the 1D electrostatic force equation. This function is used internally for label placement in the presence of overlap, as in [fixOverlap()]. #' @@ -134,6 +138,7 @@ overlapMetrics <- function(x, thresh) { # modified version, c/o K.C. Thompson # increase const --> dampen chaotic oscillation during simulation + # "softening" in N-body simulation res <- (Qk * Q1 * Q2 ) / (d^ex + const) return(res) diff --git a/R/data-documentation.R b/R/data-documentation.R index 13e7ae98a..baa09d78e 100644 --- a/R/data-documentation.R +++ b/R/data-documentation.R @@ -129,7 +129,7 @@ #' character vector} \item{field_ph}{a numeric vector} #' \item{hue}{a character vector} \item{value}{a numeric #' vector} \item{chroma}{a numeric vector} } -#' @references \url{http://casoilresource.lawr.ucdavis.edu/} +#' @references \url{https://casoilresource.lawr.ucdavis.edu/} #' @keywords datasets #' @examples #' @@ -190,7 +190,7 @@ NULL #' \item{b}{RGB blue component} \item{soil_color}{R-friendly #' encoding of soil color} } #' @author Dylan E. Beaudette -#' @references \url{http://casoilresource.lawr.ucdavis.edu/} +#' @references \url{https://casoilresource.lawr.ucdavis.edu/} #' @source Busacca, Alan J.; Singer, Michael J.; Verosub, Kenneth L. 1989. Late #' Cenozoic stratigraphy of the Feather and Yuba rivers area, California, with #' a section on soil development in mixed alluvium at Honcut Creek. USGS @@ -267,7 +267,7 @@ NULL #' \item{B}{color: b-coordinate, CIE-LAB colorspace (dry)} #' \item{name}{horizon name} \item{soil_color}{horizon color} } #' @keywords datasets -#' @references \url{http://casoilresource.lawr.ucdavis.edu/} +#' @references \url{https://casoilresource.lawr.ucdavis.edu/} #' @examples #' #' ## this example investigates the concept of a "median profile" diff --git a/R/plot_distance_graph.R b/R/plot_distance_graph.R index c723d2596..aecea145d 100644 --- a/R/plot_distance_graph.R +++ b/R/plot_distance_graph.R @@ -15,7 +15,7 @@ #' @return No value is returned. #' @author Dylan E Beaudette #' @seealso \code{\link{sp2}}, \code{\link{profile_compare}} -#' @references http://casoilresource.lawr.ucdavis.edu/ +#' @references https://casoilresource.lawr.ucdavis.edu/ #' @keywords hplot #' @export #' @examples diff --git a/R/segment.R b/R/segment.R index 7b27e0ce0..1134c8489 100644 --- a/R/segment.R +++ b/R/segment.R @@ -7,15 +7,16 @@ #' @param object either a `SoilProfileCollection` or `data.frame` #' @param intervals a vector of integers over which to slice the horizon data (e.g. `c(25, 100)` or `25:100`) #' @param trim logical, when `TRUE` horizons in `object` are truncated to the min/max specified in `intervals`. When `FALSE`, those horizons overlapping an interval are marked as such. Care should be taken when specifying more than one depth interval and `trim = FALSE`. -#' @param hzdepcols a character vector of length 2 specifying the names of the horizon depths (e.g. `c("hzdept", "hzdepb")`), only necessary if `object` is a `data.frame`. +#' @param depthcols a character vector of length 2 specifying the names of the horizon depths (e.g. `c("top", "bottom")`), only necessary if `object` is a +#' @param hzdepcols deprecated being replaced by depthcols. #' -#' @details `segment()` performs no aggregation or resampling of the source data, rather, labels are added to horizon records for subsequent aggregation or summary. This makes it possible to process a very large number of records outside of the constraints associated with e.g. `slice()` or `slab()`. +#' @details `hz_segment()` performs no aggregation or resampling of the source data, rather, labels are added to horizon records for subsequent aggregation or summary. This makes it possible to process a very large number of records outside of the constraints associated with e.g. `slice()` or `slab()`. #' #' @return Either a `SoilProfileCollection` or `data.frame` with the original horizon data segmented by depth intervals. There are usually more records in the resulting object, one for each time a segment interval partially overlaps with a horizon. A new column called `segment_id` identifying the depth interval is added. #' #' @author Stephen Roecker #' -#' @seealso [dice()], [glom()] +#' @seealso [dice()], [glom()], [hz_dissolve()], [hz_lag()], [hz_intersect()] #' #' @export #' @@ -28,7 +29,7 @@ #' depths(sp1) <- id ~ top + bottom #' #' # segment and trim -#' z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) +#' z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) #' #' # display segment labels #' # note that there are new horizon boundaries at segments @@ -52,7 +53,7 @@ #' #' a.slab <- slab(s, fm = ~ p1, slab.structure = c(0, 10, 20, 30), slab.fun = mean, na.rm = TRUE) #' -#' z <- segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) +#' z <- hz_segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) #' z <- horizons(z) #' z$thick <- z$bottom - z$top #' @@ -75,22 +76,22 @@ #' data(sp5) #' #' # segment by upper 25-cm -#' test1 <- segment(sp5, intervals = c(0, 100)) +#' test1 <- hz_segment(sp5, intervals = c(0, 100)) #' print(test1) #' nrow(test1) #' print(object.size(test1), units = "Mb") #' #' # segment by 1-cm increments -#' test2 <- segment(sp5, intervals = 0:100) +#' test2 <- hz_segment(sp5, intervals = 0:100) #' print(test2) #' nrow(test2) #' print(object.size(test2), units = "Mb") #' #' #' # segment and aggregate -#' test3 <- segment(horizons(sp5), +#' test3 <- hz_segment(horizons(sp5), #' intervals = c(0, 5, 15, 30, 60, 100, 200), -#' hzdepcols = c("top", "bottom") +#' depthcols = c("top", "bottom") #' ) #' test3$hzthk <- test3$bottom - test3$top #' test3_agg <- by(test3, test3$segment_id, function(x) { @@ -104,8 +105,8 @@ #' #' head(test3_agg) #' -segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { - +hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom")) { + # depth interval rules dep <- data.frame( top = intervals[-length(intervals)], @@ -119,10 +120,9 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { formatC(dep$bot, width = n, flag = 0) ) - # argument sanity check + # argument sanity check ---- test_spc <- inherits(object, 'SoilProfileCollection') test_df <- inherits(object, 'data.frame') - test_hd <- !is.null(hzdepcols) & length(hzdepcols) == 2 test_dep <- is.numeric(dep$top) & is.numeric(dep$bot) & all(dep$top < dep$bot) @@ -130,73 +130,74 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { stop("the input must be either a SoilProfileCollection or data.frame") } - if (!test_spc & (!test_df | !test_hd)) { - stop("if the input is a data.frame then hzdepcols must not be NULL and length(hzdepcols) == 2") - } + .check_depthcols_l(depthcols) if (!test_dep) { stop("intervals should be numeric and sequential (e.g. c(0, 1, 2, 3) or 0:100)") } - # standardize inputs + # standardize inputs ---- if (test_spc) { - peid <- idname(object) - hzid <- hzidname(object) - hzdepcols <- horizonDepths(object) + idcol <- idname(object) + hzidcol <- hzidname(object) + depthcols <- horizonDepths(object) h <- horizons(object) - names(h)[names(h) %in% c(peid, hzid)] <- c("peid", "hzid") + names(h)[names(h) %in% c(idcol, hzidcol)] <- c("idcol", "hzidcol") } else { h <- object } - names(h)[names(h) %in% hzdepcols] <- c("hzdept", "hzdepb") + names(h)[names(h) %in% depthcols] <- c("top", "bot") + ## TODO: consider using dice() - # filter horizons and trim + # filter horizons and trim ---- .slice <- function(h, top = NULL, bot = NULL) { - idx <- h$hzdept <= bot & h$hzdepb >= top + idx <- which(h$top < bot & h$bot > top) h <- h[idx, ] - # causing errors when FALSE + # causing errors when FALSE; fixed? if (trim == TRUE) { - h <- within(h, { - hzdept = ifelse(hzdept < top, top, hzdept) - hzdepb = ifelse(hzdepb > bot, bot, hzdepb) - }) + h$top = ifelse(h$top < top, top, h$top) + h$bot = ifelse(h$bot > bot, bot, h$bot) - h <- h[(h$hzdepb - h$hzdept) > 0, ] + # h <- h[(h$bot - h$top) > 0, ] } + # h <- h[!is.na(h$peiid), ] + return(h) } - # slice spc by intervals + # slice spc by intervals ---- # dep$df <- lapply(1:nrow(dep), function(x) h[0, ]) # pre-allocate memory - dep$df <- list(h[0, ])[rep(1, nrow(dep))] # pre-allocate memory faster + df_str <- cbind(h[0, ], segment_id = NA_character_[0]) + dep$df <- list(df_str)[rep(1, nrow(dep))] # pre-allocate memory faster h <- { split(dep, dep$id) ->.; lapply(., function(x) { - x$df[[1]] <- cbind(.slice(h, top = x$top, bot = x$bot), segment_id = x$id) + temp <- .slice(h, top = x$top, bot = x$bot) + if (nrow(temp) > 0) x$df[[1]] <- cbind(temp, segment_id = x$id) return(x$df[[1]]) }) ->.; do.call("rbind", .) ->.; } - names(h)[names(h) %in% c("hzdept", "hzdepb")] <- hzdepcols + names(h)[names(h) %in% c("top", "bot")] <- depthcols if (test_spc) { - h <- h[order(h$peid, h[[hzdepcols[1]]]), ] + h <- h[order(h$idcol, h[[depthcols[1]]]), ] # merge to re-add spc with NA - h_orig <- data.frame(peid = names(table(horizons(object)[peid])), stringsAsFactors = FALSE) - h <- merge(h_orig, h, by = "peid", all.x = TRUE, sort = FALSE) + h_orig <- data.frame(idcol = names(table(horizons(object)[idcol])), stringsAsFactors = FALSE) + h <- merge(h_orig, h, by = "idcol", all.x = TRUE, sort = FALSE) rm(h_orig) ## TODO: consider adding a flag to indicate "new" horizon records that have been added - # rebuild SPC - names(h)[names(h) == "peid"] <- peid - names(h)[names(h) == "hzid"] <- hzid + # rebuild SPC ---- + names(h)[names(h) == "idcol"] <- idcol + names(h)[names(h) == "hzidcol"] <- hzidcol h$hzID <- 1:nrow(h) replaceHorizons(object) <- h @@ -212,6 +213,14 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { } +#' @export +#' @rdname hz_segment +segment <- function(object, intervals, trim = TRUE, hzdepcols = c("top", "bottom")) { + .Deprecated("segment() will be deprecated and replaced by hz_segment()") + hz_segment(object, intervals, trim, depthcols = hzdepcols) +} + + #' @title Dissolving horizon boundaries by grouping variables #' @@ -219,19 +228,21 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { #' #' @param object a \code{data.frame} #' @param by character: column names, to be used as grouping variables, within the object. -#' @param id character: column name of the pedon ID within the object. -#' @param hztop character: column name of the horizon top depth within the object. -#' @param hzbot character: column name of the horizon bottom depth in the object. +#' @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")`). +#' @param id deprecated and replaced with idcol. +#' @param hztop deprecated and replaced by depthcols. +#' @param hzbot deprecated and replaced by depthcols. #' @param collapse logical: indicating whether to not combine grouping variables before dissolving. #' @param order logical: indicating whether or not to order the object by the id, hztop, and hzbot columns. -#' #' +#' #' @details This function assumes the profiles and horizons within the object follow the logic defined by \code{checkHzDepthLogic} (e.g. records are ordered sequentially by id, hztop, and hzbot and without gaps). If the records are not ordered, set the \code{order = TRUE}. #' -#' @return A \code{data.frame} with the original id, by grouping variables, and non-consecutive horizon depths. +#' @return A \code{data.frame} with the original idcol, by grouping variables, and non-consecutive horizon depths. #' #' @author Stephen Roecker #' -#' @seealso \code{\link{checkHzDepthLogic}} +#' @seealso [hz_lag()], [hz_intersect()], [hz_segment()] , [checkHzDepthLogic()] #' #' @export #' @@ -245,7 +256,7 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { #' spc$genhz <- generalize.hz(spc$name, c("A", "E", "B", "C"), c("A", "E", "B", "C")) #' h <- horizons(spc) #' -#' test <- dissolve_hz(h, by = c("genhz", "dep_5"), id = "id", hztop = "top", hzbot = "bottom") +#' test <- hz_dissolve(h, by = c("genhz", "dep_5"), idcol = "id", depthcols = c("top", "bottom")) #' #' vars <- c("id", "top", "bottom", "genhz", "dep_5") #' h[h$id == "92-1", vars] @@ -254,9 +265,9 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { #' #' # example 2 #' df <- data.frame( -#' peiid = 1, -#' hzdept = c(0, 5, 10, 15, 25, 50), -#' hzdepb = c(5, 10, 15, 25, 50, 100), +#' id = 1, +#' top = c(0, 5, 10, 15, 25, 50), +#' bottom = c(5, 10, 15, 25, 50, 100), #' hzname = c("A1", "A2", "E/A", "2Bt1", "2Bt2", "2C"), #' genhz = c("A", "A", "E", "2Bt", "2Bt", "2C"), #' texcl = c("sil", "sil", "sil", "sl", "sl", "s") @@ -264,15 +275,15 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { #' #' df #' -#' dissolve_hz(df, c("genhz", "texcl")) -#' dissolve_hz(df, c("genhz", "texcl"), collapse = TRUE) +#' hz_dissolve(df, c("genhz", "texcl")) +#' hz_dissolve(df, c("genhz", "texcl"), collapse = TRUE) #' -#' test <- dissolve_hz(df, "genhz") +#' test <- hz_dissolve(df, "genhz") #' subset(test, value == "2Bt") #' -dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzdepb", collapse = FALSE, order = FALSE) { +hz_dissolve <- function(object, by, idcol = "id", depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE) { # id = "peiid"; hztop = "hzdept"; hzbot = "hzdepb", collapse = FALSE, order = FALSE @@ -281,26 +292,21 @@ dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzd # test_spc <- inherits(object, 'SoilProfileCollection') # check that object & by are the right class - test_object <- inherits(object, "data.frame") - test_by <- inherits(by, "character") - - if (!any(test_object | test_by)) { - stop("the object argument must be a data.frame, and by a character", call. = FALSE) + test_object <- inherits(object, "data.frame") + if (!any(test_object)) { + stop("the object argument must be a data.frame", call. = FALSE) } - - # check that by is not NULL - if (is.null(by)) stop("the by argument must not be NULL") + # check that collapse is a logical of length 1 if (!inherits(collapse, "logical") || length(collapse) != 1) { stop("the collapse argument must be logical and a length of one", call. = FALSE) } - # check that the column names exisit within the object - var_names <- c(id = id, top = hztop, bot = hzbot, by) - if (!all(var_names %in% names(object))) { - stop("all arguments must match object names") - } + + # check that by is not NULL + if (is.null(by)) stop("the by argument must not be NULL") + # check that "by" are characters or convert if (any(!"character" %in% sapply(object[by], class))) { @@ -308,15 +314,30 @@ dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzd object[by] <- lapply(object[by], as.character) } + + # check that the column names exist within the object + .check_names(object, vars = c(idcol = idcol, top = depthcols[1], bot = depthcols[2], by)) + + + # check if previous dissolve_id exists and overwrite + nm <- names(object) + idx <- nm == "dissolve_id" + if (any(idx)) { + warning("object contains an existing column named 'dissolve_id', it will be overwritten") + object[idx] <- NULL + } + + # standardize inputs ---- - df <- object - idx_names <- sapply(var_names[1:3], function(x) which(names(df) == x)) - names(df)[idx_names] <- names(var_names)[1:3] + df_std <- .standardize_inputs(object, idcol = idcol, depthcols = depthcols) + df_conversion <- df_std$x_conversion + df <- df_std$x; rm(df_std) + # valid # vd_idx <- validate_depths(df, id = "id", hztop = "hzdept", bot = "hzdepb") if (order == TRUE) { - df <- df[order(df$id, df$top, df$bot), ] + df <- df[order(df$idcol, df$top, df$bot), ] } if (collapse == TRUE) { @@ -325,30 +346,393 @@ dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzd by <- by_co } + # var thickness ---- var_dep <- lapply(by, function(x) { - con_bot <- rle( paste(df$id, df[, x]))$length - con_top <- rle(rev(paste(df$id, df[, x])))$length + con_bot <- rle( paste(df$idcol, df[, x]))$length + con_top <- rle(rev(paste(df$idcol, df[, x])))$length bot_idx <- cumsum(con_bot) top_idx <- cumsum(con_top) vd <- data.frame( - id = df[bot_idx, "id"], - top = rev(rev(df$top)[top_idx]), - bot = df[bot_idx, "bot"], + idcol = df[bot_idx, "idcol"], + top = rev(rev(df$top)[top_idx]), + bot = df[bot_idx, "bot"], variable = x, - value = df[bot_idx, x] + value = df[bot_idx, x] ) + # vd$dissolve_id = 1:nrow(vd) return(vd) }) var_dep <- do.call("rbind", var_dep) - # undo standardization ---- - names(var_dep)[1:3] <- var_names[1:3] + # this is redundant with collapse = TRUE, and inappropriate unless the grouping by variable matches across all horizon depths, otherwise it'll generate pedons with overlapping depths + # if (direction == "wide") { + # var_dep <- reshape( + # var_dep, + # direction = "wide", + # idvar = c("id", "top", "bot"), + # timevar = "variable", + # v.names = "value" + # ) + # } + + + # append dissolve_id + n <- c( + var_dep$top, + var_dep$bot + ) |> + nchar() |> + max(na.rm = TRUE) + var_dep$dissolve_id <- paste0( + var_dep$idcol, + "_", + formatC(var_dep$top, width = n, flag = 0), + "-", + formatC(var_dep$bot, width = n, flag = 0), + "_", + var_dep$value + ) + + + # reset inputs ---- + var_dep <- .reset_inputs(var_dep, df_conversion) return(var_dep) } + + +#' @export +#' @rdname hz_dissolve + +dissolve_hz <- function(object, by, id = "idcol", hztop = "top", hzbot = "bottom", collapse = FALSE, order = FALSE) { + .Deprecated("dissolve_hz() will be deprecated and replaced by hz_dissolve()") + hz_dissolve(object, by, idcol = id, depthcols = c(hztop, hzbot), collapse, order) + } + + + +#' @title Intersecting horizon boundaries by horizon depths +#' +#' @description This function intersects two horizon tables by harmonizing their depths and merging them where they overlap. This can be useful to rejoin the results of `hz_dissolve()` to it's original horizon table, and then perform an aggregation on the dissolved variables. +#' +#' @param x a \code{data.frame} +#' @param y a \code{data.frame} +#' @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 . +#' +#' @return A \code{data.frame} with harmonized depth intervals (i.e. segment_id) and columns from both of the original \code{data.frame}. If both \code{data.frame} contain the same column names, they will both be returned (with the exception of the idcol and depthcols), and appended with either x or y to indicate which \code{data.frame} they originated from. +#' +#' @author Stephen Roecker +#' +#' @seealso [hz_dissolve()], [hz_lag()], [hz_segment()] +#' +#' @export +#' +#' @examples +#' +#' h <- data.frame( +#' id = 1, +#' top = c(0, 25, 44, 46, 50), +#' bottom = c(25, 44, 46, 50, 100), +#' by = c("Yes", "Yes", "No", "No", "Yes"), +#' clay = c(10, 12, 27, 35, 16) +#' ) +#' +#' h |> hz_dissolve("by") +#' +#' h |> hz_dissolve("by") |> hz_intersect(x = _, y = h) +#' +#' h |> +#' hz_dissolve("by") |> +#' hz_intersect(x = h, y = _) |> +#' aggregate(clay ~ dissolve_id, data = _, mean) +#' + + + +hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { + + # test inputs ---- + # argument sanity check + + # check that depthcols ---- + ## length == 2 + if (length(depthcols) != 2) stop("depthcols must length must equal 2") + + ## check for matching column names + .check_names(x, c(idcol, depthcols)) + .check_names(y, c(idcol, depthcols)) + + + # check segment_id ---- + ## if it exists, overwrite it + x_nm <- names(x) + y_nm <- names(y) + if (any(x_nm %in% "segment_id")) { + warning("x includes a column named 'segment_id', it will be overwritten") + x[x_nm == "segment_id"] <- NULL + } + + if (any(y_nm %in% "segment_id")) { + warning("y includes a column named 'segment_id', it will be overwritten") + y[y_nm == "segment_id"] <- NULL + } + + + # standardize inputs ---- + x_std <- .standardize_inputs(x, idcol = idcol, depthcols = depthcols) + x_conversion <- x_std$x_conversion + x <- x_std$x; rm(x_std) + + y <- .standardize_inputs(y, idcol = idcol, depthcols = depthcols)$x + + # intersect x & y ---- + split(x, x$idcol) ->.; + lapply(., function(x) { + xi <- x + yi <- y[which(y$idcol == xi$idcol[1]), ] + + if (nrow(yi) > 0) { + + int <- c(xi$top, xi$bot, yi$top, yi$bot) |> + sort() |> + unique() + + xi_seg <- hz_segment(xi, intervals = int, depthcols = names(x_conversion[2:3]), trim = TRUE) + yi_seg <- hz_segment(yi, intervals = int, depthcols = names(x_conversion[2:3]), trim = TRUE) + + return(list(x_seg = xi_seg, y_seg = yi_seg)) + } + }) ->.; + + + x_seg <- lapply(., function(x) x[["x_seg"]]) |> do.call("rbind", args = _) + y_seg <- lapply(., function(x) x[["y_seg"]]) |> do.call("rbind", args = _) + + + xy_int <- merge(x_seg, y_seg, by = c("segment_id", "idcol", "top", "bot"), sort = FALSE) + + + # reset inputs ---- + xy_int <- .reset_inputs(xy_int, x_conversion) + + return(xy_int) + } + + + +#' @title Find lagged horizon values +#' +#' @description This function finds adjacent values to a horizon values at lagged distances. +#' +#' @param object a \code{data.frame} +#' @param lag integer: number of horizons to lag +#' @param unit character: lag units in index or depth. +#' @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")`). +#' @param order logical: indicating whether or not to order the #' +#' @details . +#' +#' @return A \code{data.frame} with lagged values. +#' +#' @author Stephen Roecker +#' +#' @seealso [hz_dissolve()], [hz_intersect()], [hz_segment()] +#' +#' @export +#' +#' @examples +#' +#' h <- data.frame( +#' id = 1, +#' top = c(0, 25, 44, 46, 50), +#' bottom = c(25, 44, 46, 50, 100), +#' texcl = c("SL", "SL", "CL", "CL", "L"), +#' clay = c(10, 12, 27, 35, 16) +#' ) +#' +#' h |> hz_lag() +#' +#' h |> hz_lag(-1) +#' +#' h |> hz_lag(10:15, unit = "depth") +#' +#' h |> +#' hz_lag() |> +#' cbind(h, lag = _) |> +#' transform( +#' clay_dif = lag.clay_bot.1 - clay, +#' texcl_contrast = paste0(texcl, "-", lag.texcl_bot.1)) +#' + + + +hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { + + nm <- names(object) + idx_std <- which(! nm %in% c(idcol, depthcols)) + vars <- nm[idx_std] + + + # check arguments ---- + .check_depthcols_l(depthcols) + .check_names(object, vars = c(idcol, depthcols, vars)) + + + # standardize inputs ---- + x_std <- .standardize_inputs(object, idcol = idcol, depthcols = depthcols) + x_conversion <- x_std$x_conversion + x <- x_std$x; rm(x_std) + + + # check depths --- + if (unit == "depth" & max(object[[depthcols[2]]] > 1000)) { + warning("The maximum depth is greater than 1000, which is implausible and will be removed. To avoid this action either remove the offending horizon or convert the depth units to a measure which will not exceed 1000") + x <- x[x$bot < 1000, ] + } + + test <- aggregate(top ~ idcol, data = x, length)$top |> max() + if (unit == "index") { + if ((test - 1) < max(lag)) { + stop("lag can not be greater than the maximum number of horizons") + } + } + + + # order ---- + if (order) { + x <- x[order(x$idcol, x$top, x$bot), ] + } + + + # lag ---- + .lag_ind <- function(x, lag = lag) { + + nr <- nrow(x) + top <- 1:nr + if (lag >= 0) bot <- c((1 + lag):nr, rep(NA, lag)) + if (lag < 0) bot <- c(rep(NA, abs(lag)), 1:(nr + lag)) + + test_idcol <- x$idcol[top] == x$idcol[bot] + test_idcol <- ifelse(! test_idcol, NA, TRUE) + x_lag <- x[test_idcol * bot, vars] + if (lag >= 0) names(x_lag) <- paste0(vars, "_bot.", lag) + if (lag < 0) names(x_lag) <- paste0(vars, "_top.", abs(lag)) + + return(x_lag) + } + + + .lag_dep <- function(x, lag = lag) { + + n <- length(x) + x$.ID <- 1:nrow(x) + x_seg <- hz_segment(x, intervals = min(x$top):max(x$bot), trim = TRUE, depthcols = c("top", "bot")) + x_seg <- x_seg[1:(n + 1)] + + + x_seg <- lapply(lag, function(i) { + + x$bot_i <- x$bot + i + idx <- match( + paste(x$idcol, x$bot_i), + paste(x_seg$idcol, x_seg$bot) + ) + xi_seg <- x_seg[idx, ] + xi_seg <- x[xi_seg$.ID, vars, drop = FALSE] + xi_seg$.ID <- NULL + + if (i >= 0) names(xi_seg) <- paste0(names(xi_seg), "_bot.", i) + if (i < 0) names(xi_seg) <- paste0(names(xi_seg), "_top.", abs(i)) + + + return(xi_seg) + }) |> + do.call("cbind", args = _) + + return(x_seg) + } + + + if (unit == "index") { + x_lag <- lapply(lag, function(i) { + .lag_ind(x, i) + }) |> + do.call("cbind", args = _) + x_lag <- x_lag[sort(names(x_lag))] + } + if (unit == "depth") { + x_lag <- .lag_dep(x, lag) + x_lag <- x_lag[sort(names(x_lag))] + } + + + # # reset inputs ---- + x_lag <- .reset_inputs(x_lag, x_conversion) + + + return(x_lag) +} + + + +# check depthcols length +.check_depthcols_l <- function(x) { + if (length(x) != 2 & !is.null(x)) stop("depthcols must length must equal 2") +} + + +## check for matching column names +.check_names <- function(x, vars) { + + x_nm <- names(x) + + if (! all(vars %in% x_nm)) { + stop("x must contain columns with names that match the input arguments") + } +} + + +# standardize inputs +.standardize_inputs <- function(x, idcol = NULL, hzidcol = NULL, depthcols = NULL, texcl = NULL, clay = NULL, taxpartsize = NULL) { + + # set new names + var_names <- c( + idcol = idcol, + hzidcol = hzidcol, + top = depthcols[1], + bot = depthcols[2], + texcl = texcl, + clay = clay, + taxpartsize = taxpartsize + ) + + # find matches + idx_x <- sapply(var_names, function(i) which(names(x) == i)) + + # rename matching column names + names(x)[idx_x] <- names(var_names) + + return(list(x = x, x_conversion = var_names)) +} + + +.reset_inputs <- function(x, conversion) { + + # find original names + idx <- which(names(x) %in% names(conversion)) + + # reset original names + names(x)[idx] <- conversion + + return(x) +} + \ No newline at end of file diff --git a/R/soilColorIndices.R b/R/soilColorIndices.R index 0806bd8ba..11d4b8f98 100644 --- a/R/soilColorIndices.R +++ b/R/soilColorIndices.R @@ -137,7 +137,7 @@ barron.torrent.redness.LAB <- function(hue, value, chroma) { #' jacobs2000$rubif <- profileApply(jacobs2000, function(p) { #' #' # sum the melanization index over the 0-100cm interval -#' p0_100 <- segment(p, 0:100) +#' p0_100 <- hz_segment(p, 0:100) #' #' ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE) #' @@ -227,7 +227,7 @@ harden.rubification <- function(hue, chroma, hue_ref, chroma_ref) { #' jacobs2000$melan <- profileApply(jacobs2000, function(p) { #' #' # sum the melanization index over the 0-100cm interval -#' p0_100 <- segment(p, 0:100) +#' p0_100 <- hz_segment(p, 0:100) #' #' ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE) #' diff --git a/R/texture.R b/R/texture.R index 6be5cb168..25c536e5e 100644 --- a/R/texture.R +++ b/R/texture.R @@ -70,6 +70,8 @@ #' #' @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}. @@ -159,6 +161,15 @@ texcl_to_ssc <- function(texcl = NULL, clay = NULL, sample = FALSE) { 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 @@ -175,7 +186,12 @@ texcl_to_ssc <- function(texcl = NULL, clay = NULL, sample = FALSE) { 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") + 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 } } @@ -187,14 +203,6 @@ texcl_to_ssc <- function(texcl = NULL, clay = NULL, sample = FALSE) { } - # 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) - }) - - # if clay is present if (clay_not_null & sample == FALSE) { @@ -418,6 +426,8 @@ texmod_to_fragvoltot <- function(texmod = NULL, lieutex = NULL) { #' #' @return - `texture_to_taxpartsize`: a character vector containing `"taxpartsize"` classes #' +#' @seealso \code{\link{hz_to_taxpartsize}} +#' #' @rdname texture #' #' @export @@ -478,9 +488,9 @@ texture_to_taxpartsize <- function(texcl = NULL, clay = NULL, sand = NULL, fragv idx <- df$fragvoltot >= 35 if (any(idx)) { df[idx,] <- within(df[idx,], { - fpsc[texcl %in% sandytextures] = "sandy-skeletal" - fpsc[clay < 35] = "loamy-skeletal" fpsc[clay >= 35] = "clayey-skeletal" + fpsc[clay < 35] = "loamy-skeletal" + fpsc[texcl %in% sandytextures] = "sandy-skeletal" }) } @@ -507,6 +517,7 @@ texture_to_taxpartsize <- function(texcl = NULL, clay = NULL, sand = NULL, fragv } + #' Parse texmod from texture #' #' @param texmod vector of textural modifiers that conform to the USDA code @@ -1025,3 +1036,56 @@ fragvol_to_texmod <- function( return(df) } + +#' @title Ranking Systems for USDA Taxonomic Particle-Size and Substitute Classes of Mineral Soils +#' +#' @description Generate a vector of USDA Particle-Size and Substitute Classes names, sorted 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 an ordered factor +#' +#' @author Stephen Roecker +#' +#' @seealso [hz_to_taxpartsize()], [texture_to_taxpartsize()], [SoilTextureLevels()] +#' +#' @export +#' @examples +#' +#' # class codes +#' PSCS_levels() +#' + +PSCS_levels <- function() { + + fe <- c("fragmental", "pumiceous", "cindery", "sandy-skeletal", "loamy-skeletal", "gypseous-skeletal", "ashy-skeletal", "medial-skeletal", "ashy-pumiceous", "medial-pumiceous", "clay-skeletal", "sandy", "ashy", "coarse-loamy", "coarse-silty", "coarse-gypseous", "medial", "loamy", "fine-gypseous", "fine-loamy", "fine-silty", "hydrous", "fine", "clayey", "very-fine", "diatomaceous") |> rev() + + # 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) which(fe == y)) |> unlist() + # 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) + }) + + fe <- data.frame(rn = 1:length(fe), fe = fe) + sc <- data.frame(rn = unlist(idx), fe = .pscs_sc) + lu <- rbind(fe, sc) + lu <- lu[order(lu$rn), ] + lu$rn <- 1:nrow(lu) + lu$fe <- factor(lu$fe, levels = lu$fe, ordered = TRUE) + + return(lu$fe) +} + + +.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() diff --git a/R/thicknessOf.R b/R/thicknessOf.R index 9fdfbf787..3c3640568 100644 --- a/R/thicknessOf.R +++ b/R/thicknessOf.R @@ -57,8 +57,9 @@ thicknessOf <- function(x, FUN = function(x, pattern, hzdesgn, ...) grepl(pattern, x[[hzdesgn]]), na.rm = FALSE, ...) { + .internalTHK <- NULL - .interalHZM <- NULL + .internalHZM <- NULL if (is.null(hzdesgn) || !hzdesgn %in% horizonNames(x)) { stop("Horizon designation column (", hzdesgn, ") does not exist.") diff --git a/R/unroll.R b/R/unroll.R index ef6f3ebf3..6f65b03b9 100644 --- a/R/unroll.R +++ b/R/unroll.R @@ -21,7 +21,7 @@ #' defaults to FALSE #' @return a vector of "unrolled" property values #' @author Dylan E. Beaudette -#' @references http://casoilresource.lawr.ucdavis.edu/ +#' @references https://casoilresource.lawr.ucdavis.edu/ #' @keywords manip #' @export #' @examples diff --git a/R/zzz.R b/R/zzz.R index 22c5dc903..1f36d21e6 100644 --- a/R/zzz.R +++ b/R/zzz.R @@ -1,6 +1,4 @@ -# .onLoad <- function(lib, pkg) { -# packageStartupMessage("This is aqp ", utils::packageDescription("aqp", field="Version"), "\n", "see http://casoilresource.lawr.ucdavis.edu/drupal/taxonomy/term/56 for examples", appendLF = TRUE) -# } + .onAttach <- function(lib, pkg) { packageStartupMessage("This is aqp ", utils::packageDescription("aqp", field="Version"), appendLF = TRUE) diff --git a/README.Rmd b/README.Rmd index 26e2fa4e4..8c085df28 100644 --- a/README.Rmd +++ b/README.Rmd @@ -28,7 +28,7 @@ knitr::opts_chunk$set( aqp hexsticker (Paxton, Montauk, Woodbridge, Ridgebury, Whitman, Catden soil series dendogram) -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 . 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 . 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. +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 . 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 . 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. ## Installation @@ -49,7 +49,7 @@ Install suggested packages: p <- c("colorspace", "ape", "soilDB", "latticeExtra", "tactile", "compositions", "sharpshootR", "markovchain", "xtable", "testthat", "Gmedian", "farver", "Hmisc", "tibble", "RColorBrewer", "scales", "digest", -"MASS", "mpspline2", "soiltexture", "knitr", "rmarkdown") +"MASS", "mpspline2", "soiltexture", "knitr", "rmarkdown", "mvtnorm") install.packages(p) ``` @@ -94,6 +94,11 @@ plotSPC( citation("aqp") ``` +## Related Papers and Book Chapters + * Beaudette D.E., P. Roudier, and J. Skovlin. 2016. Probabilistic representation of genetic soil horizons. In Book Digital soil morphometrics. Springer. + * Maynard, J.J., S.W. Salley, D.E. Beaudette, and J.E. Herrick. 2020. Numerical soil classification supports soil identification by citizen scientists using limited, simple soil observations. Soil Science Society of America Journal 84:1675-1692. + * Beaudette, D. E., J. Skovlin, A. G. Brown, P. Roudier, and S. M. Roecker. "Algorithms for Quantitative Pedology." In Geopedology, edited by Joseph Alfred Zinck, Graciela Metternicht, Héctor Francisco del Valle, and Marcos Angelini, 201–22. Cham: Springer International Publishing, 2023. https://doi.org/10.1007/978-3-031-20667-2_11. + ## Related Packages * [soilDB](https://github.com/ncss-tech/soilDB) * [sharpshootR](https://github.com/ncss-tech/sharpshootR) @@ -102,7 +107,7 @@ citation("aqp") * [Introduction to SoilProfileCollection Objects](https://ncss-tech.github.io/aqp/articles/Introduction-to-SoilProfileCollection-Objects.html) * [Numerical Classification of Soil Profiles](https://ncss-tech.github.io/aqp/articles/NCSP.html) * [Overlapping Annotation](https://ncss-tech.github.io/aqp/articles/label-placement.html) - * [What is new in aqp 2.0?](https://ncss-tech.github.io/aqp/articles/new-in-aqp-2.html) + * [What is new in aqp 2.x?](https://ncss-tech.github.io/aqp/articles/new-in-aqp-2.html) ## Tutorials * [Soil Profile Sketches](https://ncss-tech.github.io/AQP/aqp/sketches.html) diff --git a/README.md b/README.md index 8a2e2b371..fb9598403 100644 --- a/README.md +++ b/README.md @@ -28,7 +28,7 @@ processing tasks . 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 -. Components of +. 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 @@ -117,6 +117,21 @@ citation("aqp") #> 'options(citation.bibtex.max=999)'. ``` +## Related Papers and Book Chapters + +- Beaudette D.E., P. Roudier, and J. Skovlin. 2016. Probabilistic + representation of genetic soil horizons. In Book Digital soil + morphometrics. Springer. +- Maynard, J.J., S.W. Salley, D.E. Beaudette, and J.E. Herrick. 2020. + Numerical soil classification supports soil identification by citizen + scientists using limited, simple soil observations. Soil Science + Society of America Journal 84:1675-1692. +- Beaudette, D. E., J. Skovlin, A. G. Brown, P. Roudier, and S. M. + Roecker. “Algorithms for Quantitative Pedology.” In Geopedology, + edited by Joseph Alfred Zinck, Graciela Metternicht, Héctor Francisco + del Valle, and Marcos Angelini, 201–22. Cham: Springer International + Publishing, 2023. . + ## Related Packages - [soilDB](https://github.com/ncss-tech/soilDB) @@ -131,7 +146,7 @@ citation("aqp") - [Overlapping Annotation](https://ncss-tech.github.io/aqp/articles/label-placement.html) - [What is new in aqp - 2.0?](https://ncss-tech.github.io/aqp/articles/new-in-aqp-2.html) + 2.x?](https://ncss-tech.github.io/aqp/articles/new-in-aqp-2.html) ## Tutorials diff --git a/man/PSCS_levels.Rd b/man/PSCS_levels.Rd new file mode 100644 index 000000000..dd3c70bbe --- /dev/null +++ b/man/PSCS_levels.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/texture.R +\name{PSCS_levels} +\alias{PSCS_levels} +\title{Ranking Systems for USDA Taxonomic Particle-Size and Substitute Classes of Mineral Soils} +\usage{ +PSCS_levels() +} +\value{ +an ordered factor +} +\description{ +Generate a vector of USDA Particle-Size and Substitute Classes names, sorted according to approximate particle size +} +\examples{ + +# class codes +PSCS_levels() + +} +\references{ +\href{https://nrcspad.sc.egov.usda.gov/DistributionCenter/product.aspx?ProductID=991}{Field Book for Describing and Sampling Soils, version 3.0} +} +\seealso{ +\code{\link[=hz_to_taxpartsize]{hz_to_taxpartsize()}}, \code{\link[=texture_to_taxpartsize]{texture_to_taxpartsize()}}, \code{\link[=SoilTextureLevels]{SoilTextureLevels()}} +} +\author{ +Stephen Roecker +} diff --git a/man/allocate.Rd b/man/allocate.Rd index 76b2f92e7..763816b1c 100644 --- a/man/allocate.Rd +++ b/man/allocate.Rd @@ -175,3 +175,6 @@ Richards, L.A. 1954. \href{https://www.ars.usda.gov/ARSUserFiles/20360500/hb60_p Soil Survey Staff, 2014. Keys to Soil Taxonomy, 12th ed. USDA-Natural Resources Conservation Service, Washington, D.C. } +\author{ +Stephen Roecker +} diff --git a/man/aqp-package.Rd b/man/aqp-package.Rd index 1588cf566..5f42665e8 100644 --- a/man/aqp-package.Rd +++ b/man/aqp-package.Rd @@ -3,7 +3,6 @@ \docType{package} \name{aqp-package} \alias{aqp-package} -\alias{_PACKAGE} \alias{aqp} \alias{aqp.env} \title{Algorithms for Quantitative Pedology} diff --git a/man/diagnostic_hz-set.Rd b/man/diagnostic_hz-set.Rd deleted file mode 100644 index 91d0508bf..000000000 --- a/man/diagnostic_hz-set.Rd +++ /dev/null @@ -1,40 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/SoilProfileCollection-setters.R -\docType{methods} -\name{diagnostic_hz<-} -\alias{diagnostic_hz<-} -\alias{diagnostic_hz<-,SoilProfileCollection-method} -\title{Add Data to Diagnostic Features Slot} -\usage{ -\S4method{diagnostic_hz}{SoilProfileCollection}(object) <- value -} -\arguments{ -\item{object}{A SoilProfileCollection} - -\item{value}{An object inheriting \code{data.frame}} -} -\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} -} -\examples{ - -# load test data -data(sp2) - -# promote to SPC -depths(sp2) <- id ~ top + bottom - -# assign two profiles a zone related to the mollic epipedon -newdata <- data.frame(id = c("hon-1","hon-17"), - featkind = "fixed-depth surface sample", - featdept = 0, - featdepb = 18) - -# do left join -diagnostic_hz(sp2) <- newdata - -# inspect site table: newvalue TRUE only for horizons -# with top depth equal to zero -diagnostic_hz(sp2) - -} diff --git a/man/diagnostic_hz.Rd b/man/diagnostic_hz.Rd index 074c70033..45a6089e8 100644 --- a/man/diagnostic_hz.Rd +++ b/man/diagnostic_hz.Rd @@ -1,16 +1,52 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Class-SoilProfileCollection.R +% Please edit documentation in R/Class-SoilProfileCollection.R, +% R/SoilProfileCollection-setters.R \docType{methods} \name{diagnostic_hz,SoilProfileCollection-method} \alias{diagnostic_hz,SoilProfileCollection-method} \alias{diagnostic_hz} -\title{Retrieve diagnostic data from SoilProfileCollection} +\alias{diagnostic_hz<-} +\alias{diagnostic_hz<-,SoilProfileCollection-method} +\title{Get or Set Diagnostic Horizon data in a SoilProfileCollection} \usage{ \S4method{diagnostic_hz}{SoilProfileCollection}(object) + +\S4method{diagnostic_hz}{SoilProfileCollection}(object) <- value } \arguments{ -\item{object}{a SoilProfileCollection} +\item{object}{A SoilProfileCollection} + +\item{value}{An object inheriting from \code{data.frame}} } \description{ -Get diagnostic feature data from SoilProfileCollection. Result is returned in the same \code{data.frame} class used to initially construct the SoilProfileCollection. +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. +\itemize{ +\item \code{diagnostic_hz()} (get method): Get diagnostic feature data from a SoilProfileCollection. +} + +\itemize{ +\item \verb{diagnostic_hz<-} (set method): Set diagnostic feature data for a SoilProfileCollection. The profile ID column from \code{object} (\code{idname(object)}) must be present in the replacement \code{value} object. +} +} +\examples{ + +# load test data +data(sp2) + +# promote to SPC +depths(sp2) <- id ~ top + bottom + +# assign two profiles a zone related to the mollic epipedon +newdata <- data.frame(id = c("hon-1","hon-17"), + featkind = "fixed-depth surface sample", + featdept = 0, + featdepb = 18) + +# do left join +diagnostic_hz(sp2) <- newdata + +# inspect site table: newvalue TRUE only for horizons +# with top depth equal to zero +diagnostic_hz(sp2) + } diff --git a/man/harden.melanization.Rd b/man/harden.melanization.Rd index a129727dc..0c8066b70 100644 --- a/man/harden.melanization.Rd +++ b/man/harden.melanization.Rd @@ -50,7 +50,7 @@ jacobs2000$c_horizon_color <- profileApply(jacobs2000, function(p) { jacobs2000$melan <- profileApply(jacobs2000, function(p) { # sum the melanization index over the 0-100cm interval - p0_100 <- segment(p, 0:100) + p0_100 <- hz_segment(p, 0:100) ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE) diff --git a/man/harden.rubification.Rd b/man/harden.rubification.Rd index 79a7fb792..d0c52ff25 100644 --- a/man/harden.rubification.Rd +++ b/man/harden.rubification.Rd @@ -56,7 +56,7 @@ jacobs2000$c_horizon_color <- profileApply(jacobs2000, function(p) { jacobs2000$rubif <- profileApply(jacobs2000, function(p) { # sum the melanization index over the 0-100cm interval - p0_100 <- segment(p, 0:100) + p0_100 <- hz_segment(p, 0:100) ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE) diff --git a/man/dissolve_hz.Rd b/man/hz_dissolve.Rd similarity index 61% rename from man/dissolve_hz.Rd rename to man/hz_dissolve.Rd index 44df6c0a7..742725828 100644 --- a/man/dissolve_hz.Rd +++ b/man/hz_dissolve.Rd @@ -1,15 +1,25 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/segment.R -\name{dissolve_hz} +\name{hz_dissolve} +\alias{hz_dissolve} \alias{dissolve_hz} \title{Dissolving horizon boundaries by grouping variables} \usage{ +hz_dissolve( + object, + by, + idcol = "id", + depthcols = c("top", "bottom"), + collapse = FALSE, + order = FALSE +) + dissolve_hz( object, by, - id = "peiid", - hztop = "hzdept", - hzbot = "hzdepb", + id = "idcol", + hztop = "top", + hzbot = "bottom", collapse = FALSE, order = FALSE ) @@ -19,19 +29,22 @@ dissolve_hz( \item{by}{character: column names, to be used as grouping variables, within the object.} -\item{id}{character: column name of the pedon ID within the object.} - -\item{hztop}{character: column name of the horizon top depth within the object.} +\item{idcol}{character: column name of the pedon ID within the object.} -\item{hzbot}{character: column name of the horizon bottom depth in the object.} +\item{depthcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("top", "bottom")}).} \item{collapse}{logical: indicating whether to not combine grouping variables before dissolving.} -\item{order}{logical: indicating whether or not to order the object by the id, hztop, and hzbot columns. -#'} +\item{order}{logical: indicating whether or not to order the object by the id, hztop, and hzbot columns.} + +\item{id}{deprecated and replaced with idcol.} + +\item{hztop}{deprecated and replaced by depthcols.} + +\item{hzbot}{deprecated and replaced by depthcols.} } \value{ -A \code{data.frame} with the original id, by grouping variables, and non-consecutive horizon depths. +A \code{data.frame} with the original idcol, by grouping variables, and non-consecutive horizon depths. } \description{ This function dissolves or combines horizons that have a common set of grouping variables. It only combines those horizon records that are sequential (e.g. share a horizon boundary). Thus, it can be used to identify discontinuities in the grouping variables along a profile and their unique depths. It is particularly useful for determining the depth to the top or bottom of horizons with a specific category, and should be simpler than previous methods that require aggregating over profiles. @@ -49,7 +62,7 @@ spc$dep_5 <- spc$depletion_pct >=5 spc$genhz <- generalize.hz(spc$name, c("A", "E", "B", "C"), c("A", "E", "B", "C")) h <- horizons(spc) -test <- dissolve_hz(h, by = c("genhz", "dep_5"), id = "id", hztop = "top", hzbot = "bottom") +test <- hz_dissolve(h, by = c("genhz", "dep_5"), idcol = "id", depthcols = c("top", "bottom")) vars <- c("id", "top", "bottom", "genhz", "dep_5") h[h$id == "92-1", vars] @@ -58,9 +71,9 @@ test[test$id == "92-1", ] # example 2 df <- data.frame( - peiid = 1, - hzdept = c(0, 5, 10, 15, 25, 50), - hzdepb = c(5, 10, 15, 25, 50, 100), + id = 1, + top = c(0, 5, 10, 15, 25, 50), + bottom = c(5, 10, 15, 25, 50, 100), hzname = c("A1", "A2", "E/A", "2Bt1", "2Bt2", "2C"), genhz = c("A", "A", "E", "2Bt", "2Bt", "2C"), texcl = c("sil", "sil", "sil", "sl", "sl", "s") @@ -68,15 +81,15 @@ df <- data.frame( df -dissolve_hz(df, c("genhz", "texcl")) -dissolve_hz(df, c("genhz", "texcl"), collapse = TRUE) +hz_dissolve(df, c("genhz", "texcl")) +hz_dissolve(df, c("genhz", "texcl"), collapse = TRUE) -test <- dissolve_hz(df, "genhz") +test <- hz_dissolve(df, "genhz") subset(test, value == "2Bt") } \seealso{ -\code{\link{checkHzDepthLogic}} +\code{\link[=hz_lag]{hz_lag()}}, \code{\link[=hz_intersect]{hz_intersect()}}, \code{\link[=hz_segment]{hz_segment()}} , \code{\link[=checkHzDepthLogic]{checkHzDepthLogic()}} } \author{ Stephen Roecker diff --git a/man/hz_intersect.Rd b/man/hz_intersect.Rd new file mode 100644 index 000000000..304e85a78 --- /dev/null +++ b/man/hz_intersect.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/segment.R +\name{hz_intersect} +\alias{hz_intersect} +\title{Intersecting horizon boundaries by horizon depths} +\usage{ +hz_intersect(x, y, idcol = "id", depthcols = c("top", "bottom")) +} +\arguments{ +\item{x}{a \code{data.frame}} + +\item{y}{a \code{data.frame}} + +\item{idcol}{character: column name of the pedon ID within the object.} + +\item{depthcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("top", "bottom")}).} +} +\value{ +A \code{data.frame} with harmonized depth intervals (i.e. segment_id) and columns from both of the original \code{data.frame}. If both \code{data.frame} contain the same column names, they will both be returned (with the exception of the idcol and depthcols), and appended with either x or y to indicate which \code{data.frame} they originated from. +} +\description{ +This function intersects two horizon tables by harmonizing their depths and merging them where they overlap. This can be useful to rejoin the results of \code{hz_dissolve()} to it's original horizon table, and then perform an aggregation on the dissolved variables. +} +\details{ +. +} +\examples{ + +h <- data.frame( +id = 1, +top = c(0, 25, 44, 46, 50), +bottom = c(25, 44, 46, 50, 100), +by = c("Yes", "Yes", "No", "No", "Yes"), +clay = c(10, 12, 27, 35, 16) +) + +h |> hz_dissolve("by") + +h |> hz_dissolve("by") |> hz_intersect(x = _, y = h) + +h |> +hz_dissolve("by") |> +hz_intersect(x = h, y = _) |> +aggregate(clay ~ dissolve_id, data = _, mean) + +} +\seealso{ +\code{\link[=hz_dissolve]{hz_dissolve()}}, \code{\link[=hz_lag]{hz_lag()}}, \code{\link[=hz_segment]{hz_segment()}} +} +\author{ +Stephen Roecker +} diff --git a/man/hz_lag.Rd b/man/hz_lag.Rd new file mode 100644 index 000000000..1691dc65e --- /dev/null +++ b/man/hz_lag.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/segment.R +\name{hz_lag} +\alias{hz_lag} +\title{Find lagged horizon values} +\usage{ +hz_lag( + object, + lag = 1, + unit = "index", + idcol = "id", + depthcols = c("top", "bottom"), + order = FALSE +) +} +\arguments{ +\item{object}{a \code{data.frame}} + +\item{lag}{integer: number of horizons to lag} + +\item{unit}{character: lag units in index or depth.} + +\item{idcol}{character: column name of the pedon ID within the object.} + +\item{depthcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("top", "bottom")}).} + +\item{order}{logical: indicating whether or not to order the #'} +} +\value{ +A \code{data.frame} with lagged values. +} +\description{ +This function finds adjacent values to a horizon values at lagged distances. +} +\details{ +. +} +\examples{ + +h <- data.frame( +id = 1, +top = c(0, 25, 44, 46, 50), +bottom = c(25, 44, 46, 50, 100), +texcl = c("SL", "SL", "CL", "CL", "L"), +clay = c(10, 12, 27, 35, 16) +) + +h |> hz_lag() + +h |> hz_lag(-1) + +h |> hz_lag(10:15, unit = "depth") + +h |> +hz_lag() |> +cbind(h, lag = _) |> +transform( +clay_dif = lag.clay_bot.1 - clay, +texcl_contrast = paste0(texcl, "-", lag.texcl_bot.1)) + +} +\seealso{ +\code{\link[=hz_dissolve]{hz_dissolve()}}, \code{\link[=hz_intersect]{hz_intersect()}}, \code{\link[=hz_segment]{hz_segment()}} +} +\author{ +Stephen Roecker +} diff --git a/man/segment.Rd b/man/hz_segment.Rd similarity index 74% rename from man/segment.Rd rename to man/hz_segment.Rd index bfad07f73..db41673df 100644 --- a/man/segment.Rd +++ b/man/hz_segment.Rd @@ -1,10 +1,13 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/segment.R -\name{segment} +\name{hz_segment} +\alias{hz_segment} \alias{segment} \title{Segmenting of Soil Horizon Data by Depth Interval} \usage{ -segment(object, intervals, trim = TRUE, hzdepcols = NULL) +hz_segment(object, intervals, trim = TRUE, depthcols = c("top", "bottom")) + +segment(object, intervals, trim = TRUE, hzdepcols = c("top", "bottom")) } \arguments{ \item{object}{either a \code{SoilProfileCollection} or \code{data.frame}} @@ -13,7 +16,9 @@ segment(object, intervals, trim = TRUE, hzdepcols = NULL) \item{trim}{logical, when \code{TRUE} horizons in \code{object} are truncated to the min/max specified in \code{intervals}. When \code{FALSE}, those horizons overlapping an interval are marked as such. Care should be taken when specifying more than one depth interval and \code{trim = FALSE}.} -\item{hzdepcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("hzdept", "hzdepb")}), only necessary if \code{object} is a \code{data.frame}.} +\item{depthcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("top", "bottom")}), only necessary if \code{object} is a} + +\item{hzdepcols}{deprecated being replaced by depthcols.} } \value{ Either a \code{SoilProfileCollection} or \code{data.frame} with the original horizon data segmented by depth intervals. There are usually more records in the resulting object, one for each time a segment interval partially overlaps with a horizon. A new column called \code{segment_id} identifying the depth interval is added. @@ -22,7 +27,7 @@ Either a \code{SoilProfileCollection} or \code{data.frame} with the original hor This function segments or subdivides horizon data from a \code{SoilProfileCollection} or \code{data.frame} by depth interval (e.g. \code{c(0, 10)}, \code{c(0, 50)}, or \code{25:100}). This results in horizon records being split at the specified depth intervals, which duplicates the original horizon data but also adds new horizon depths. In addition, labels (i.e. \code{"segment_id"}) are added to each horizon record that correspond with their depth interval (e.g. \code{025-100}). This function is intended to harmonize horizons to a common support (i.e. depth interval) for further aggregation or summary. See the examples. } \details{ -\code{segment()} performs no aggregation or resampling of the source data, rather, labels are added to horizon records for subsequent aggregation or summary. This makes it possible to process a very large number of records outside of the constraints associated with e.g. \code{slice()} or \code{slab()}. +\code{hz_segment()} performs no aggregation or resampling of the source data, rather, labels are added to horizon records for subsequent aggregation or summary. This makes it possible to process a very large number of records outside of the constraints associated with e.g. \code{slice()} or \code{slab()}. } \examples{ @@ -33,7 +38,7 @@ data(sp1) depths(sp1) <- id ~ top + bottom # segment and trim -z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) +z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) # display segment labels # note that there are new horizon boundaries at segments @@ -57,7 +62,7 @@ s <- combine(s) a.slab <- slab(s, fm = ~ p1, slab.structure = c(0, 10, 20, 30), slab.fun = mean, na.rm = TRUE) -z <- segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) +z <- hz_segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) z <- horizons(z) z$thick <- z$bottom - z$top @@ -80,22 +85,22 @@ res$diff < 0.001 data(sp5) # segment by upper 25-cm -test1 <- segment(sp5, intervals = c(0, 100)) +test1 <- hz_segment(sp5, intervals = c(0, 100)) print(test1) nrow(test1) print(object.size(test1), units = "Mb") # segment by 1-cm increments -test2 <- segment(sp5, intervals = 0:100) +test2 <- hz_segment(sp5, intervals = 0:100) print(test2) nrow(test2) print(object.size(test2), units = "Mb") # segment and aggregate -test3 <- segment(horizons(sp5), +test3 <- hz_segment(horizons(sp5), intervals = c(0, 5, 15, 30, 60, 100, 200), - hzdepcols = c("top", "bottom") + depthcols = c("top", "bottom") ) test3$hzthk <- test3$bottom - test3$top test3_agg <- by(test3, test3$segment_id, function(x) { @@ -111,7 +116,7 @@ head(test3_agg) } \seealso{ -\code{\link[=dice]{dice()}}, \code{\link[=glom]{glom()}} +\code{\link[=dice]{dice()}}, \code{\link[=glom]{glom()}}, \code{\link[=hz_dissolve]{hz_dissolve()}}, \code{\link[=hz_lag]{hz_lag()}}, \code{\link[=hz_intersect]{hz_intersect()}} } \author{ Stephen Roecker diff --git a/man/hz_to_taxpartsize.Rd b/man/hz_to_taxpartsize.Rd new file mode 100644 index 000000000..65e699a32 --- /dev/null +++ b/man/hz_to_taxpartsize.Rd @@ -0,0 +1,84 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/allocate.R +\name{hz_to_taxpartsize} +\alias{hz_to_taxpartsize} +\title{Allocate Particle Size Control Class for the Control Section.} +\usage{ +hz_to_taxpartsize( + x, + y, + taxpartsize = "taxpartsize", + clay = "clay", + idcol = "id", + depthcols = c("top", "bottom") +) +} +\arguments{ +\item{x}{a \code{data.frame} containing the original horizon table.} + +\item{y}{a \code{data.frame} containing the particle size control section depths for each idcol.} + +\item{taxpartsize}{\code{character} column name for taxonomic family particle size class.} + +\item{clay}{\code{character} column name for clay percent.} + +\item{idcol}{character: column name of the pedon ID within the object.} + +\item{depthcols}{a character vector of length 2 specifying the names of the horizon depths (e.g. \code{c("top", "bottom")}).} +} +\value{ +A \code{data.frame} object containing the original idcol and aggregated particle size control section allocation. +} +\description{ +This function aggregates information in the horizon table and allocates it to the particle size control section. +} +\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. +} +\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) + + +} +\seealso{ +\code{\link[=texture_to_taxpartsize]{texture_to_taxpartsize()}}, \code{\link[=PSCS_levels]{PSCS_levels()}} +} +\author{ +Stephen Roecker +} diff --git a/man/plot_distance_graph.Rd b/man/plot_distance_graph.Rd index 8ba2f8cc6..fdef71127 100644 --- a/man/plot_distance_graph.Rd +++ b/man/plot_distance_graph.Rd @@ -39,7 +39,7 @@ plot_distance_graph(d, idx=7:12) plot_distance_graph(d, idx=12:18) } \references{ -http://casoilresource.lawr.ucdavis.edu/ +https://casoilresource.lawr.ucdavis.edu/ } \seealso{ \code{\link{sp2}}, \code{\link{profile_compare}} diff --git a/man/restrictions-set.Rd b/man/restrictions-set.Rd deleted file mode 100644 index 5b4f4c83e..000000000 --- a/man/restrictions-set.Rd +++ /dev/null @@ -1,39 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/SoilProfileCollection-setters.R -\docType{methods} -\name{restrictions<-} -\alias{restrictions<-} -\alias{restrictions<-,SoilProfileCollection-method} -\title{Add Data to Restrictions Slot} -\usage{ -\S4method{restrictions}{SoilProfileCollection}(object) <- value -} -\arguments{ -\item{object}{A SoilProfileCollection} - -\item{value}{An object inheriting \code{data.frame}} -} -\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}. -} -\examples{ - -# load test data -data(sp2) - -# promote to SPC -depths(sp2) <- id ~ top + bottom - -# assign abrupt textural change to a profile -newdata <- data.frame(id = c("hon-21"), - restrkind = "abrupt textural change", - restrdep = 46) - -# do left join -restrictions(sp2) <- newdata - -# inspect site table: newvalue TRUE only for horizons -# with top depth equal to zero -restrictions(sp2) - -} diff --git a/man/restrictions.Rd b/man/restrictions.Rd index bebe38ef3..d003e07bd 100644 --- a/man/restrictions.Rd +++ b/man/restrictions.Rd @@ -1,16 +1,51 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/Class-SoilProfileCollection.R +% Please edit documentation in R/Class-SoilProfileCollection.R, +% R/SoilProfileCollection-setters.R \docType{methods} \name{restrictions,SoilProfileCollection-method} \alias{restrictions,SoilProfileCollection-method} \alias{restrictions} -\title{Retrieve restriction data from SoilProfileCollection} +\alias{restrictions<-} +\alias{restrictions<-,SoilProfileCollection-method} +\title{Get or Set Restriction data in a SoilProfileCollection} \usage{ \S4method{restrictions}{SoilProfileCollection}(object) + +\S4method{restrictions}{SoilProfileCollection}(object) <- value } \arguments{ -\item{object}{a SoilProfileCollection} +\item{object}{A SoilProfileCollection} + +\item{value}{An data.frame object containing at least a column with name \code{idname(object)}} } \description{ -Get restriction data from SoilProfileCollection. Result is returned in the same \code{data.frame} class used to initially construct the SoilProfileCollection. +Restrictions describe root-limiting features in the soil. A single profile may have multiple restrictions. +\itemize{ +\item \code{restrictions()} (get method): Get restriction data from a SoilProfileCollection. +} + +\itemize{ +\item \verb{restrictions<-} (set method): Set restriction data for a SoilProfileCollection. The profile ID column from \code{object} (\code{idname(object)}) must be present in the replacement \code{value} object. +} +} +\examples{ + +# load test data +data(sp2) + +# promote to SPC +depths(sp2) <- id ~ top + bottom + +# assign abrupt textural change to a profile +newdata <- data.frame(id = c("hon-21"), + restrkind = "abrupt textural change", + restrdep = 46) + +# do left join +restrictions(sp2) <- newdata + +# inspect site table: newvalue TRUE only for horizons +# with top depth equal to zero +restrictions(sp2) + } diff --git a/man/sp1.Rd b/man/sp1.Rd index 39a1f0b2e..4c9af583e 100644 --- a/man/sp1.Rd +++ b/man/sp1.Rd @@ -55,6 +55,6 @@ panel=panel.superpose, ylim=c(110,-5), asp=2) } \references{ -\url{http://casoilresource.lawr.ucdavis.edu/} +\url{https://casoilresource.lawr.ucdavis.edu/} } \keyword{datasets} diff --git a/man/sp2.Rd b/man/sp2.Rd index 8cc3eaddf..304c6295f 100644 --- a/man/sp2.Rd +++ b/man/sp2.Rd @@ -81,7 +81,7 @@ legend('topleft', legend=levels(sp2$surface), col=1:6, pch=15, bty='n', bg='whit } \references{ -\url{http://casoilresource.lawr.ucdavis.edu/} +\url{https://casoilresource.lawr.ucdavis.edu/} } \author{ Dylan E. Beaudette diff --git a/man/sp3.Rd b/man/sp3.Rd index 358163c19..26e0912fc 100644 --- a/man/sp3.Rd +++ b/man/sp3.Rd @@ -148,6 +148,6 @@ plotSPC( } } \references{ -\url{http://casoilresource.lawr.ucdavis.edu/} +\url{https://casoilresource.lawr.ucdavis.edu/} } \keyword{datasets} diff --git a/man/texture.Rd b/man/texture.Rd index a1cc39446..59c179810 100644 --- a/man/texture.Rd +++ b/man/texture.Rd @@ -275,6 +275,11 @@ fragvol_to_texmod(gravel = 10, cobbles = 10) \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}. } +\seealso{ +\code{\link{SoilTextureLevels}} + +\code{\link{hz_to_taxpartsize}} +} \author{ Stephen Roecker } diff --git a/man/unroll.Rd b/man/unroll.Rd index 4bad22e39..fca291a17 100644 --- a/man/unroll.Rd +++ b/man/unroll.Rd @@ -49,7 +49,7 @@ plot(x, 1:length(x), ylim=c(90,0), type='b', cex=0.5) } \references{ -http://casoilresource.lawr.ucdavis.edu/ +https://casoilresource.lawr.ucdavis.edu/ } \author{ Dylan E. Beaudette diff --git a/misc/sandbox/slice-wise-Tukey-HSD.R b/misc/sandbox/slice-wise-Tukey-HSD.R index d3f62a3c2..b095d1ace 100644 --- a/misc/sandbox/slice-wise-Tukey-HSD.R +++ b/misc/sandbox/slice-wise-Tukey-HSD.R @@ -3,11 +3,11 @@ library(aqp) library(soilDB) library(lattice) library(latticeExtra) -library(viridis) +library(tactile) library(grid) # define plotting style -tps <- list(superpose.line=list(col=c('RoyalBlue', 'DarkRed', 'DarkGreen'), lwd=2)) +tps <- tactile.theme(superpose.line=list(col=c('RoyalBlue', 'DarkRed', 'DarkGreen'), lwd=2)) # get multiple series' worth of data # TODO: add code to deal with those series that return 0 records @@ -72,7 +72,7 @@ ck <- list( # standard output from slab() -p.1 <- xyplot(top ~ p.q50 | variable, groups=taxonname, data=HSD$agg, ylab='Depth', +p.1 <- xyplot(top ~ p.q50 | variable, groups=taxonname, data=HSD$agg, ylab='Depth (cm)', xlab='median bounded by 25th and 75th percentiles', lower=HSD$agg$p.q25, upper=HSD$agg$p.q75, ylim=c(105, -5), panel=panel.depth_function, alpha=0.25, sync.colors=TRUE, @@ -94,7 +94,7 @@ p.1 <- xyplot(top ~ p.q50 | variable, groups=taxonname, data=HSD$agg, ylab='Dept # experimental HSD viz p.2 <- segplot( hzn_top ~ lwr + upr, centers = diff, level = p.adj, data = HSD$HSD, - col.regions = viridis, ylim = c(105, -5), + col.regions = hcl.colors, ylim = c(105, -5), xlab = 'HSD', at = col.at, colorkey = ck, @@ -178,7 +178,7 @@ z <- z[idx, vars] dd <- datadist(z) options(datadist="dd") -# GLS: I've used this before to parametrize correlation sturcture +# GLS: I've used this before to parametrize correlation structure (m.gls <- Gls(wmpd ~ rcs(mid) * group, data = z, correlation = corAR1(form = ~ mid | sliceID))) plot(Predict(m.gls)) diff --git a/misc/sandbox/slice-wise-entropy-example.R b/misc/sandbox/slice-wise-entropy-example.R index 6b8038f54..e62c6660c 100644 --- a/misc/sandbox/slice-wise-entropy-example.R +++ b/misc/sandbox/slice-wise-entropy-example.R @@ -1,8 +1,8 @@ library(aqp) library(latticeExtra) -library(plyr) +library(tactile) library(Hmisc) -library(reshape) +library(reshape2) data(sp3) @@ -12,7 +12,7 @@ depths(sp3) <- id ~ top + bottom # http://en.wikipedia.org/wiki/Differential_entropy # http://cran.r-project.org/web/packages/entropy/ -# calculation for continous random vars based on binning / counts +# calculation for continuous random vars based on binning / counts # http://cran.r-project.org/web/packages/entropy/entropy.pdf ## this isn't correct, and barfs when there is < 90% data available @@ -154,14 +154,16 @@ wtd.sum.qcd <- function(i){ # compute some "information" metrics -a <- slab(sp3, ~ clay + A + cec + ph, slab.fun=mean.and.sd, slab.structure=0:100) -a.1 <- slab(sp3, ~ clay + A + cec + ph, slab.fun=f.entropy, slab.structure=0:100) -a.2 <- slab(sp3, ~ clay + A + cec + ph, slab.fun=f.sig.to.noise, slab.structure=0:100) -a.3 <- slab(sp3, ~ clay + A + cec + ph, slab.fun=f.qcd, slab.structure=0:100) +a <- slab(sp3, ~ clay + A + cec + ph, slab.fun = mean.and.sd, slab.structure = 0:100) +a.1 <- slab(sp3, ~ clay + A + cec + ph, slab.fun = f.entropy, slab.structure = 0:100) +a.2 <- slab(sp3, ~ clay + A + cec + ph, slab.fun = f.sig.to.noise, slab.structure = 0:100) +a.3 <- slab(sp3, ~ clay + A + cec + ph, slab.fun = f.qcd, slab.structure = 0:100) # combine -g <- make.groups(summary=a, entropy=a.1, sig.to.noise=a.2, qcd=a.3) -g$which <- factor(g$which, labels=c('Mean +/- 1SD', 'psuedo-Entropy', 'Signal : Noise', 'QCD')) +g <- make.groups(summary = a, entropy = a.1, sig.to.noise = a.2, qcd = a.3) +g$which <- factor(g$which, labels = c('Mean +/- 1SD', 'psuedo-Entropy', 'Signal : Noise', 'QCD')) + +.tps <- tactile.theme(superpose.line=list(lwd=2, col=c('RoyalBlue', 'Orange2'))) p <- xyplot( top ~ value | which + variable, data=g, @@ -170,7 +172,7 @@ p <- xyplot( ylab='Depth (cm)', xlab='', ylim=c(100,-5), layout=c(5,3), scales=list(x=list(relation='free')), - par.settings=list(superpose.line=list(lwd=2, col=c('RoyalBlue', 'Orange2'))), + par.settings = .tps, panel=panel.depth_function, prepanel=prepanel.depth_function, auto.key=list(columns=2, lines=TRUE, points=FALSE) @@ -181,20 +183,21 @@ useOuterStrips(p, strip=strip.custom(bg=grey(0.85)), strip.left = strip.custom(h ## a "no-information QCD" must be computed from the raw data, by depth-slice -s <- slice(sp3, 0:100 ~ clay + A + cec + ph, just.the.data = TRUE) +s <- dice(sp3, fm = 0:100 ~ clay + A + cec + ph, SPC = FALSE) s.long <- melt(s, id.vars = c('top', 'bottom'), measure.vars = c('clay', 'A', 'cec', 'ph')) -qcd.ni <- ddply(s.long, c('variable'), no.information.qcd) + +qcd.ni <- lapply(split(s.long, s.long$variable), no.information.qcd) ## compute weighted mean and weighted sum QCD by variable ## note that these must be standardized by slice-wise "no-information" QCD -wm.qcd <- ddply(a.3, 'variable', .fun=wtd.mean.qcd) -ws.qcd <- ddply(a.3, 'variable', .fun=wtd.sum.qcd) - -### does this make sense? -# join QCD summaries to "no-information" baseline -ss <- join(join(wm.qcd, ws.qcd), qcd.ni) -transform(ss, mean.qcd=wt.mean.qcd / mean.ni.qcd, sum.qcd=wt.sum.qcd / sum.ni.qcd) - +# wm.qcd <- ddply(a.3, 'variable', .fun=wtd.mean.qcd) +# ws.qcd <- ddply(a.3, 'variable', .fun=wtd.sum.qcd) +# +# ### does this make sense? +# # join QCD summaries to "no-information" baseline +# ss <- join(join(wm.qcd, ws.qcd), qcd.ni) +# transform(ss, mean.qcd=wt.mean.qcd / mean.ni.qcd, sum.qcd=wt.sum.qcd / sum.ni.qcd) +# diff --git a/tests/testthat/test-segment.R b/tests/testthat/test-segment.R index 760540ec8..5e83d9391 100644 --- a/tests/testthat/test-segment.R +++ b/tests/testthat/test-segment.R @@ -1,4 +1,4 @@ -context("segment") +context("hz_segment") test_that("data.frame interface works as expected", { @@ -7,7 +7,7 @@ test_that("data.frame interface works as expected", { data(sp1) # trimming - z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE, hzdepcols = c('top', 'bottom')) + z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE, depthcols = c('top', 'bottom')) # correct object type and segment label expect_true(inherits(z, 'data.frame')) @@ -17,7 +17,7 @@ test_that("data.frame interface works as expected", { expect_true(inherits(z[['segment_id']], 'character')) # no triming - z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = FALSE, hzdepcols = c('top', 'bottom')) + z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = FALSE, depthcols = c('top', 'bottom')) # correct object type and segment label expect_true(inherits(z, 'data.frame')) @@ -35,7 +35,7 @@ test_that("SPC interface works as expected", { depths(sp1) <- id ~ top + bottom # trimming - z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) + z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE) expect_true(inherits(z, 'SoilProfileCollection')) expect_true('segment_id' %in% horizonNames(z)) @@ -44,7 +44,7 @@ test_that("SPC interface works as expected", { expect_true(inherits(z[['segment_id']], 'character')) # no trimming - z <- segment(sp1, intervals = c(0, 10, 20, 30), trim = FALSE) + z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = FALSE) expect_true(inherits(z, 'SoilProfileCollection')) expect_true('segment_id' %in% horizonNames(z)) @@ -69,8 +69,8 @@ test_that("expected outcome with NA horizon depths", { bad$top[c(1, 5)] <- NA # segment - z.bad <- segment(bad, intervals = c(0, 10, 20, 30), trim = TRUE, hzdepcols = c('top', 'bottom')) - z.good <- segment(good, intervals = c(0, 10, 20, 30), trim = TRUE, hzdepcols = c('top', 'bottom')) + z.bad <- hz_segment(bad, intervals = c(0, 10, 20, 30), trim = TRUE, depthcols = c('top', 'bottom')) + z.good <- hz_segment(good, intervals = c(0, 10, 20, 30), trim = TRUE, depthcols = c('top', 'bottom')) # label class expect_true(inherits(z.good[['segment_id']], 'character')) @@ -86,34 +86,35 @@ test_that("expected outcome with NA horizon depths", { }) -test_that("expected outcome with bogus horizon depths", { - - # init local copy of sample data - data(sp1) - - # copies - good <- sp1 - bad <- sp1 - - # add NA to horizon depths - bad$top[c(1, 5)] <- bad$bottom[c(1, 5)] - - # segment - z.bad <- segment(bad, intervals = c(0, 10, 20, 30), trim = TRUE, hzdepcols = c('top', 'bottom')) - z.good <- segment(good, intervals = c(0, 10, 20, 30), trim = TRUE, hzdepcols = c('top', 'bottom')) - - # label class - expect_true(inherits(z.good[['segment_id']], 'character')) - expect_true(inherits(z.bad[['segment_id']], 'character')) - - ## TODO: is this expected? - # row count - expect_false(nrow(z.good) == nrow(z.bad)) - - # same values - # expect_false(all(z.good$segment_id == z.bad$segment_id)) - -}) +# I think this test needs to be retired or reframed. segment previously exclude results where the thickness of the segment was zero. That feature has been removed. Upon updating segment it was removed to ensure the original data was returned regardless of the horizon errors, which should be dealt with elsewhere. +# test_that("expected outcome with bogus horizon depths", { +# +# # init local copy of sample data +# data(sp1) +# +# # copies +# good <- sp1 +# bad <- sp1 +# +# # add NA to horizon depths +# bad$top[c(1, 5)] <- bad$bottom[c(1, 5)] +# +# # segment +# z.bad <- hz_segment(bad, intervals = c(0, 10, 20, 30), trim = TRUE, depthcols = c('top', 'bottom')) +# z.good <- hz_segment(good, intervals = c(0, 10, 20, 30), trim = TRUE, depthcols = c('top', 'bottom')) +# +# # label class +# expect_true(inherits(z.good[['segment_id']], 'character')) +# expect_true(inherits(z.bad[['segment_id']], 'character')) +# +# ## TODO: is this expected? +# # row count +# expect_false(nrow(z.good) == nrow(z.bad)) +# +# # same values +# # expect_false(all(z.good$segment_id == z.bad$segment_id)) +# +# }) @@ -128,7 +129,7 @@ test_that("same results as weighted mean via slab", { a.slab <- slab(s, fm = ~ p1, slab.structure = c(0, 10, 20, 30), slab.fun = mean, na.rm = TRUE) # segment - z <- segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) + z <- hz_segment(s, intervals = c(0, 10, 20, 30), trim = TRUE) # compute horizon thickness weights z <- horizons(z) diff --git a/vignettes/Munsell-color-conversion.Rmd b/vignettes/Munsell-color-conversion.Rmd index 407d4b5af..5d608c347 100644 --- a/vignettes/Munsell-color-conversion.Rmd +++ b/vignettes/Munsell-color-conversion.Rmd @@ -61,7 +61,7 @@ axis(side = 2, las = 1) Neutral colors are commonly specified two ways in the Munsell system: `N 3/` or `N 3/0`, either format will work with `munsell2rgb()` and `parseMunsell()`. -Non-standard Munsell notation (e.g. `3.6YR 4.4 / 5.6`), possibly collected with a sensor vs. color book, can be approximated with `getClosestMunsellChip()`. A more accurate conversion can be performed with the [`munsellinterpol` package.](https://cran.r-project.org/web/packages/munsellinterpol/index.html). +Non-standard Munsell notation (e.g. `3.6YR 4.4 / 5.6`), possibly collected with a sensor vs. color book, can be approximated with `getClosestMunsellChip()`. A more accurate conversion can be performed with the [`munsellinterpol` package.](https://cran.r-project.org/package=munsellinterpol). ## Examples