From 33a5177b6e131d8dbc3e51e2334414dc3f17faef Mon Sep 17 00:00:00 2001 From: Beaudette Date: Fri, 5 Apr 2024 11:25:21 -0700 Subject: [PATCH 01/22] notes --- R/aqp-label-placement-solvers.R | 5 +++++ 1 file changed, 5 insertions(+) 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) From 921c84fdadb5847bdb3cbf133dd30aa245313eac Mon Sep 17 00:00:00 2001 From: Beaudette Date: Mon, 8 Apr 2024 12:48:18 -0700 Subject: [PATCH 02/22] Update slice-wise-entropy-example.R --- misc/sandbox/slice-wise-entropy-example.R | 43 ++++++++++++----------- 1 file changed, 23 insertions(+), 20 deletions(-) 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) +# From 671ce79dd0c3cc12e5fe4089a60225fd7beec63e Mon Sep 17 00:00:00 2001 From: Beaudette Date: Mon, 8 Apr 2024 12:48:30 -0700 Subject: [PATCH 03/22] update, ensuring this still works #81 --- misc/sandbox/slice-wise-Tukey-HSD.R | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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)) From 368a2d74f9c7de91ada2ced18d18ada11edf221b Mon Sep 17 00:00:00 2001 From: Beaudette Date: Mon, 8 Apr 2024 13:17:38 -0700 Subject: [PATCH 04/22] typo --- R/thicknessOf.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) 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.") From f1f1aee89bc410c7a6b81f3070c954a4cba38a27 Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Tue, 9 Apr 2024 07:34:27 -0700 Subject: [PATCH 05/22] Update NEWS.md --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index bb1bb6413..efff3c793 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,8 +1,9 @@ -# aqp 2.0.3 (2024-03-20) +# aqp 2.0.3 (2024-04-09) * `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 From 9786bf5a47aa39ebc03524f16a9c5bdbddc8c3cd Mon Sep 17 00:00:00 2001 From: Beaudette Date: Wed, 17 Apr 2024 12:06:17 -0700 Subject: [PATCH 06/22] preparing for CRAN release. --- NEWS.md | 5 +++-- README.Rmd | 9 +++++++-- README.md | 17 ++++++++++++++++- man/aqp-package.Rd | 1 - 4 files changed, 26 insertions(+), 6 deletions(-) diff --git a/NEWS.md b/NEWS.md index efff3c793..e83c62137 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,10 @@ -# aqp 2.0.3 (2024-04-09) +# aqp 2.0.3 (2024-04-17) + * 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. + * 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 diff --git a/README.Rmd b/README.Rmd index 26e2fa4e4..7ead3308c 100644 --- a/README.Rmd +++ b/README.Rmd @@ -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..ecdcc90bb 100644 --- a/README.md +++ b/README.md @@ -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/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} From 82e215cbd67c23f3d99f1d36d7a49b64af0940cf Mon Sep 17 00:00:00 2001 From: Beaudette Date: Wed, 17 Apr 2024 13:23:35 -0700 Subject: [PATCH 07/22] preparing for CRAN release --- CITATION.cff | 2 +- DESCRIPTION | 2 +- README.Rmd | 2 +- README.md | 2 +- vignettes/Munsell-color-conversion.Rmd | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) 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 6e5441358..f1c577d25 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, sp, stringr, data.table, farver Suggests: mvtnorm, colorspace, ape, soilDB, 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/README.Rmd b/README.Rmd index 7ead3308c..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 diff --git a/README.md b/README.md index ecdcc90bb..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 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 From 8bfb577d86a96cf970bfe37683426a73b32a6b1a Mon Sep 17 00:00:00 2001 From: Beaudette Date: Thu, 18 Apr 2024 09:54:27 -0700 Subject: [PATCH 08/22] Update NEWS.md --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index e83c62137..ebe6c0a91 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# aqp 2.0.3 (2024-04-17) +# aqp 2.0.3 (2024-04-18) * CRAN release * `simulateColor()` gains new method `mvnorm` for simulating plausible colors - package mvtnorm added to SUGGESTS From 7544741ed389fe2db857affec51d443061581b46 Mon Sep 17 00:00:00 2001 From: Beaudette Date: Thu, 18 Apr 2024 10:22:30 -0700 Subject: [PATCH 09/22] fixes for CRAN checks --- NEWS.md | 2 +- R/data-documentation.R | 6 +++--- R/plot_distance_graph.R | 2 +- R/unroll.R | 2 +- R/zzz.R | 4 +--- man/plot_distance_graph.Rd | 2 +- man/sp1.Rd | 2 +- man/sp2.Rd | 2 +- man/sp3.Rd | 2 +- man/unroll.Rd | 2 +- 10 files changed, 12 insertions(+), 14 deletions(-) diff --git a/NEWS.md b/NEWS.md index ebe6c0a91..3db71cddc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -602,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/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/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/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/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/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 From 2db9f707641d663c887027074724736bb57bd2ee Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 21 May 2024 12:25:30 -0500 Subject: [PATCH 10/22] adding hz_intersect() and renaming segment() and dissolve_hz() also renaming several segment() and dissolve_hz() arguments to be consistent with SoilProfileCollection() --- NAMESPACE | 3 + R/allocate.R | 2 +- R/segment.R | 330 ++++++++++++++++++++----- R/soilColorIndices.R | 4 +- man/harden.melanization.Rd | 2 +- man/harden.rubification.Rd | 2 +- man/{dissolve_hz.Rd => hz_dissolve.Rd} | 52 ++-- man/hz_intersect.Rd | 48 ++++ man/{segment.Rd => hz_segment.Rd} | 37 ++- tests/testthat/test-segment.R | 73 +++--- 10 files changed, 417 insertions(+), 136 deletions(-) rename man/{dissolve_hz.Rd => hz_dissolve.Rd} (64%) create mode 100644 man/hz_intersect.Rd rename man/{segment.Rd => hz_segment.Rd} (76%) diff --git a/NAMESPACE b/NAMESPACE index ea7f94412..f6c233992 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -95,6 +95,9 @@ export(hzOffset) export(hzTopographyCodeToLineType) export(hzTopographyCodeToOffset) export(hzTransitionProbabilities) +export(hz_dissolve) +export(hz_intersect) +export(hz_segment) export(invertLabelColor) export(lunique) export(maxDepthOf) diff --git a/R/allocate.R b/R/allocate.R index d0b55e1c0..3b75611bd 100644 --- a/R/allocate.R +++ b/R/allocate.R @@ -413,7 +413,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 diff --git a/R/segment.R b/R/segment.R index 7b27e0ce0..7f43e55ba 100644 --- a/R/segment.R +++ b/R/segment.R @@ -7,9 +7,10 @@ #' @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. #' @@ -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,7 +105,12 @@ #' #' head(test3_agg) #' -segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { +hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom"), hzdepcols = NULL) { + + if (!is.null(hzdepcols) & is.null(depthcols)) { + .Deprecated("hzdepcols is being replaced with depthcols to be consistent with SoilProfileCollection()") + depthcols <- hzdepcols + } # depth interval rules dep <- data.frame( @@ -119,10 +125,10 @@ 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_hd <- !is.null(depthcols) & length(depthcols) == 2 test_dep <- is.numeric(dep$top) & is.numeric(dep$bot) & all(dep$top < dep$bot) @@ -131,7 +137,7 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { } 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") + stop("if the input is a data.frame then depthcols must not be NULL and length(depthcols) == 2") } if (!test_dep) { @@ -139,64 +145,67 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { } - # 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 +221,14 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = NULL) { } +#' @export +#' @rdname hz_segment +segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom"), hzdepcols = NULL) { + .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,15 +236,17 @@ 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 #' @@ -245,7 +264,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 +273,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,43 +283,47 @@ 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", id = NULL, hztop = NULL, hzbot = NULL, depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE) { # id = "peiid"; hztop = "hzdept"; hzbot = "hzdepb", collapse = FALSE, order = FALSE # test inputs ---- # argument sanity check + # check idcol and hzdepcols + if (!is.null(hztop) & !is.null(hzbot)) { + .Deprecated("hztop and hztop parameters will be deprecated and replaced by hzdepcols in order to be consistent with SoilProfileCollection()") + depthcols <- c(hztop, hzbot) + } + if (!is.null(id) & is.null(idcol)) { + .Deprecated("id will be deprecated and replaced by idcol in order to be consistent with SoilProfileCollection()") + depthcols <- c(hztop, hzbot) + } + # test_spc <- inherits(object, 'SoilProfileCollection') # check that object & by are the right class - test_object <- inherits(object, "data.frame") - test_by <- inherits(by, "character") + test_object <- inherits(object, "data.frame") + # test_by <- inherits(by, "character") - if (!any(test_object | test_by)) { + if (!any(test_object)) { stop("the object argument must be a data.frame, and by a character", 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,6 +331,21 @@ 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 + var_names <- c(idcol = idcol, top = depthcols[1], bot = depthcols[2], by) + if (!all(var_names %in% names(object))) { + stop("all arguments must match object names") + } + + # 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)) @@ -316,7 +354,7 @@ dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzd # 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) { @@ -328,27 +366,183 @@ dissolve_hz <- function(object, by, id = "peiid", hztop = "hzdept", hzbot = "hzd # 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) + + # 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" + # ) + # } + + # undo standardization ---- names(var_dep)[1:3] <- var_names[1:3] + # append dissolve_id + n <- c( + var_dep[["top"]], + var_dep[["bottom"]] + ) |> + nchar() |> + max() + var_dep$dissolve_id <- paste0( + var_dep[[idcol]], + "_", + formatC(var_dep[["top"]], width = n, flag = 0), + "-", + formatC(var_dep[["bottom"]], width = n, flag = 0), + "_", + var_dep[["variable"]] + ) + + return(var_dep) } + + +#' @export +#' @rdname hz_dissolve + +dissolve_hz <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot = NULL, depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE) { + .Deprecated("dissolve_hz() will be deprecated and replaced by hz_dissolve()") + hz_dissolve(object, by, idcol, id, hztop, hzbot, depthcols, 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 +#' +#' +#' @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 + x_deps <- names(x[depthcols]) + y_deps <- names(y[depthcols]) + + if (!all.equal( + names(x[depthcols]), + names(y[depthcols])) + ) stop("both x and y must contain columns with names matching 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 ---- + var_names <- c(idcol = idcol, top = depthcols[1], bot = depthcols[2]) + idx_x <- sapply(var_names[1:3], function(i) which(names(x) == i)) + idx_y <- sapply(var_names[1:3], function(i) which(names(y) == i)) + names(x)[idx_x] <- names(var_names)[1:3] + names(y)[idx_y] <- names(var_names)[1:3] + + + # intersect x & y ---- + split(x, x$idcol) ->.; + lapply(., function(x) { + xi <- x + yi <- y[which(y$idcol == xi$idcol[1]), ] + + int <- c(xi$top, xi$bot, yi$top, yi$bot) |> + sort() |> + unique() + + xi_seg <- hz_segment(xi, intervals = int, depthcols = names(var_names[2:3]), trim = TRUE) + yi_seg <- hz_segment(yi, intervals = int, depthcols = names(var_names[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")) + + + # undo standardization ---- + names(xy_int)[names(xy_int) %in% names(var_names)] <- var_names + + + return(xy_int) + } + + + \ 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/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 64% rename from man/dissolve_hz.Rd rename to man/hz_dissolve.Rd index 44df6c0a7..9fc244679 100644 --- a/man/dissolve_hz.Rd +++ b/man/hz_dissolve.Rd @@ -1,15 +1,30 @@ % 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", + id = NULL, + hztop = NULL, + hzbot = NULL, + depthcols = c("top", "bottom"), + collapse = FALSE, + order = FALSE +) + dissolve_hz( object, by, - id = "peiid", - hztop = "hzdept", - hzbot = "hzdepb", + idcol = "id", + id = NULL, + hztop = NULL, + hzbot = NULL, + depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE ) @@ -19,19 +34,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{idcol}{character: column name of the pedon ID within the object.} + +\item{id}{deprecated and replaced with idcol.} + +\item{hztop}{deprecated and replaced by depthcols.} -\item{hztop}{character: column name of the horizon top depth within the object.} +\item{hzbot}{deprecated and replaced by depthcols.} -\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.} } \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 +67,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 +76,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,10 +86,10 @@ 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") } diff --git a/man/hz_intersect.Rd b/man/hz_intersect.Rd new file mode 100644 index 000000000..b30948f1c --- /dev/null +++ b/man/hz_intersect.Rd @@ -0,0 +1,48 @@ +% 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) + +} +\author{ +Stephen Roecker +} diff --git a/man/segment.Rd b/man/hz_segment.Rd similarity index 76% rename from man/segment.Rd rename to man/hz_segment.Rd index bfad07f73..02f16afea 100644 --- a/man/segment.Rd +++ b/man/hz_segment.Rd @@ -1,10 +1,25 @@ % 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"), + hzdepcols = NULL +) + +segment( + object, + intervals, + trim = TRUE, + depthcols = c("top", "bottom"), + hzdepcols = NULL +) } \arguments{ \item{object}{either a \code{SoilProfileCollection} or \code{data.frame}} @@ -13,7 +28,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 +39,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 +50,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 +74,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 +97,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) { 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) From 4caa286597ae2a2feff8e63bd1256e8cce45f5ce Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 21 May 2024 19:08:58 -0500 Subject: [PATCH 11/22] adding hz_lag() --- NAMESPACE | 1 + R/segment.R | 221 ++++++++++++++++++++++++++++++++++---------- man/hz_dissolve.Rd | 23 ++--- man/hz_intersect.Rd | 7 +- man/hz_lag.Rd | 60 ++++++++++++ man/hz_segment.Rd | 16 +--- 6 files changed, 247 insertions(+), 81 deletions(-) create mode 100644 man/hz_lag.Rd diff --git a/NAMESPACE b/NAMESPACE index f6c233992..78adeb991 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -97,6 +97,7 @@ export(hzTopographyCodeToOffset) export(hzTransitionProbabilities) export(hz_dissolve) export(hz_intersect) +export(hz_lag) export(hz_segment) export(invertLabelColor) export(lunique) diff --git a/R/segment.R b/R/segment.R index 7f43e55ba..c482ee8ab 100644 --- a/R/segment.R +++ b/R/segment.R @@ -105,13 +105,8 @@ #' #' head(test3_agg) #' -hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom"), hzdepcols = NULL) { +hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom")) { - if (!is.null(hzdepcols) & is.null(depthcols)) { - .Deprecated("hzdepcols is being replaced with depthcols to be consistent with SoilProfileCollection()") - depthcols <- hzdepcols - } - # depth interval rules dep <- data.frame( top = intervals[-length(intervals)], @@ -128,7 +123,6 @@ hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bot # argument sanity check ---- test_spc <- inherits(object, 'SoilProfileCollection') test_df <- inherits(object, 'data.frame') - test_hd <- !is.null(depthcols) & length(depthcols) == 2 test_dep <- is.numeric(dep$top) & is.numeric(dep$bot) & all(dep$top < dep$bot) @@ -136,9 +130,7 @@ hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bot 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 depthcols must not be NULL and length(depthcols) == 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)") @@ -223,9 +215,9 @@ hz_segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bot #' @export #' @rdname hz_segment -segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom"), hzdepcols = NULL) { +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) + hz_segment(object, intervals, trim, depthcols = hzdepcols) } @@ -291,51 +283,41 @@ segment <- function(object, intervals, trim = TRUE, depthcols = c("top", "bottom #' -hz_dissolve <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot = NULL, depthcols = c("top", "bottom"), 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 # test inputs ---- # argument sanity check - # check idcol and hzdepcols - if (!is.null(hztop) & !is.null(hzbot)) { - .Deprecated("hztop and hztop parameters will be deprecated and replaced by hzdepcols in order to be consistent with SoilProfileCollection()") - depthcols <- c(hztop, hzbot) - } - if (!is.null(id) & is.null(idcol)) { - .Deprecated("id will be deprecated and replaced by idcol in order to be consistent with SoilProfileCollection()") - depthcols <- c(hztop, hzbot) - } - # 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)) { - stop("the object argument must be a data.frame, and by a character", call. = FALSE) + stop("the object argument must be a data.frame", call. = FALSE) } + # 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 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))) { message("non-character grouping variables are being converted to characters") object[by] <- lapply(object[by], as.character) } + # check that the column names exist within the object - var_names <- c(idcol = idcol, top = depthcols[1], bot = depthcols[2], by) - if (!all(var_names %in% names(object))) { - stop("all arguments must match object names") - } + .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) @@ -347,9 +329,10 @@ hz_dissolve <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot # 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") @@ -363,6 +346,7 @@ hz_dissolve <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot by <- by_co } + # var thickness ---- var_dep <- lapply(by, function(x) { @@ -398,8 +382,8 @@ hz_dissolve <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot # } - # undo standardization ---- - names(var_dep)[1:3] <- var_names[1:3] + # reset inputs ---- + var_dep <- .reset_inputs(var_dep, df_conversion) # append dissolve_id @@ -427,9 +411,9 @@ hz_dissolve <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot #' @export #' @rdname hz_dissolve -dissolve_hz <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot = NULL, depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE) { +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, hztop, hzbot, depthcols, collapse, order) + hz_dissolve(object, by, idcol = id, depthcols = c(hztop, hzbot), collapse, order) } @@ -466,9 +450,10 @@ dissolve_hz <- function(object, by, idcol = "id", id = NULL, hztop = NULL, hzbot #' #' h |> hz_dissolve("by") |> hz_intersect(x = _, y = h) #' -#' h |> hz_dissolve("by") |> -#' hz_intersect(x = h, y = _) |> -#' aggregate(clay ~ dissolve_id, data = _, mean) +#' h |> +#' hz_dissolve("by") |> +#' hz_intersect(x = h, y = _) |> +#' aggregate(clay ~ dissolve_id, data = _, mean) #' @@ -508,12 +493,11 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { # standardize inputs ---- - var_names <- c(idcol = idcol, top = depthcols[1], bot = depthcols[2]) - idx_x <- sapply(var_names[1:3], function(i) which(names(x) == i)) - idx_y <- sapply(var_names[1:3], function(i) which(names(y) == i)) - names(x)[idx_x] <- names(var_names)[1:3] - names(y)[idx_y] <- names(var_names)[1:3] + 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) ->.; @@ -525,8 +509,8 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { sort() |> unique() - xi_seg <- hz_segment(xi, intervals = int, depthcols = names(var_names[2:3]), trim = TRUE) - yi_seg <- hz_segment(yi, intervals = int, depthcols = names(var_names[2:3]), trim = TRUE) + 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)) }) ->.; @@ -537,12 +521,149 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { xy_int <- merge(x_seg, y_seg, by = c("segment_id", "idcol", "top", "bot")) - # undo standardization ---- - names(xy_int)[names(xy_int) %in% names(var_names)] <- var_names - + # 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 vars character: column names, to lag. +#' @param lag integer: number of horizons to lag +#' @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 +#' +#' +#' @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(c("texcl", "clay"), 1) +#' +#' h |> +#' hz_lag(c("texcl", "clay"), 1) |> +#' cbind(h) |> +#' transform( +#' clay_dif = clay - clay_lag.1, +#' texcl_contrast = paste0(texcl, "-", texcl_lag.1)) +#' + + + +hz_lag <- function(object, vars, lag = 1, idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { + + # 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) + + + # order ---- + if (order) { + x[order(x$idcol, x$top, x$bot), ] + } + + + # lag ---- + .lag <- function(x, lag = lag, var = NULL) { + + nr <- nrow(x) + top <- 1:nr + bot <- c((1 + lag):nr, rep(NA, lag)) + + test_idcol <- x$idcol[top] == x$idcol[bot] + # test_deps <- x$bot[top] == x$top[bot] + lag_var <- ifelse( + test_idcol & !is.na(test_idcol), + x[bot, var], + NA) + + return(lag_var) + } + + x_lag <- lapply(vars, function(y) .lag(x, lag = lag, y)) |> + do.call("data.frame", args = _) + names(x_lag) <- paste0(vars, "_lag.", lag) + + # # reset inputs ---- + # x <- .reset_inputs(cbind(x, 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) { + + # set new names + var_names <- c( + idcol = idcol, + hzidcol = hzidcol, + top = depthcols[1], + bot = depthcols[2] + ) + + # 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/man/hz_dissolve.Rd b/man/hz_dissolve.Rd index 9fc244679..19fc9f1f2 100644 --- a/man/hz_dissolve.Rd +++ b/man/hz_dissolve.Rd @@ -9,9 +9,6 @@ hz_dissolve( object, by, idcol = "id", - id = NULL, - hztop = NULL, - hzbot = NULL, depthcols = c("top", "bottom"), collapse = FALSE, order = FALSE @@ -20,11 +17,9 @@ hz_dissolve( dissolve_hz( object, by, - idcol = "id", - id = NULL, - hztop = NULL, - hzbot = NULL, - depthcols = c("top", "bottom"), + id = "idcol", + hztop = "top", + hzbot = "bottom", collapse = FALSE, order = FALSE ) @@ -36,17 +31,17 @@ dissolve_hz( \item{idcol}{character: column name of the pedon ID within the object.} -\item{id}{deprecated and replaced with idcol.} - -\item{hztop}{deprecated and replaced by depthcols.} - -\item{hzbot}{deprecated and replaced by depthcols.} - \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{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 idcol, by grouping variables, and non-consecutive horizon depths. diff --git a/man/hz_intersect.Rd b/man/hz_intersect.Rd index b30948f1c..23503244d 100644 --- a/man/hz_intersect.Rd +++ b/man/hz_intersect.Rd @@ -38,9 +38,10 @@ 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) +h |> +hz_dissolve("by") |> +hz_intersect(x = h, y = _) |> +aggregate(clay ~ dissolve_id, data = _, mean) } \author{ diff --git a/man/hz_lag.Rd b/man/hz_lag.Rd new file mode 100644 index 000000000..bdd73f10c --- /dev/null +++ b/man/hz_lag.Rd @@ -0,0 +1,60 @@ +% 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, + vars, + lag = 1, + idcol = "id", + depthcols = c("top", "bottom"), + order = FALSE +) +} +\arguments{ +\item{object}{a \code{data.frame}} + +\item{vars}{character: column names, to lag.} + +\item{lag}{integer: number of horizons to lag} + +\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(c("texcl", "clay"), 1) + +h |> +hz_lag(c("texcl", "clay"), 1) |> +cbind(h) |> +transform( +clay_dif = clay - clay_lag.1, +texcl_contrast = paste0(texcl, "-", texcl_lag.1)) + +} +\author{ +Stephen Roecker +} diff --git a/man/hz_segment.Rd b/man/hz_segment.Rd index 02f16afea..d2071bf15 100644 --- a/man/hz_segment.Rd +++ b/man/hz_segment.Rd @@ -5,21 +5,9 @@ \alias{segment} \title{Segmenting of Soil Horizon Data by Depth Interval} \usage{ -hz_segment( - object, - intervals, - trim = TRUE, - depthcols = c("top", "bottom"), - hzdepcols = NULL -) +hz_segment(object, intervals, trim = TRUE, depthcols = c("top", "bottom")) -segment( - object, - intervals, - trim = TRUE, - depthcols = c("top", "bottom"), - hzdepcols = NULL -) +segment(object, intervals, trim = TRUE, hzdepcols = c("top", "bottom")) } \arguments{ \item{object}{either a \code{SoilProfileCollection} or \code{data.frame}} From f6ee56fa672a4b5e0416cbd392088482d4c28bd1 Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 21 May 2024 22:55:45 -0500 Subject: [PATCH 12/22] Update hz_lag() --- R/segment.R | 32 ++++++++++++++++++++------------ man/hz_lag.Rd | 10 +++++----- 2 files changed, 25 insertions(+), 17 deletions(-) diff --git a/R/segment.R b/R/segment.R index c482ee8ab..b70530de8 100644 --- a/R/segment.R +++ b/R/segment.R @@ -558,10 +558,10 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { #' clay = c(10, 12, 27, 35, 16) #' ) #' -#' h |> hz_lag(c("texcl", "clay"), 1) +#' h |> hz_lag() #' #' h |> -#' hz_lag(c("texcl", "clay"), 1) |> +#' hz_lag() |> #' cbind(h) |> #' transform( #' clay_dif = clay - clay_lag.1, @@ -570,7 +570,12 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { -hz_lag <- function(object, vars, lag = 1, idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { +hz_lag <- function(object, lag = 1, vars = NULL, idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { + + nm <- names(object) + idx <- which(! nm %in% c(idcol, depthcols)) + if (is.null(vars)) vars <- nm[idx] + # check arguments ---- .check_depthcols_l(depthcols) @@ -590,7 +595,7 @@ hz_lag <- function(object, vars, lag = 1, idcol = "id", depthcols = c("top", "bo # lag ---- - .lag <- function(x, lag = lag, var = NULL) { + .lag <- function(x, lag = lag, vars = NULL) { nr <- nrow(x) top <- 1:nr @@ -598,17 +603,20 @@ hz_lag <- function(object, vars, lag = 1, idcol = "id", depthcols = c("top", "bo test_idcol <- x$idcol[top] == x$idcol[bot] # test_deps <- x$bot[top] == x$top[bot] - lag_var <- ifelse( - test_idcol & !is.na(test_idcol), - x[bot, var], - NA) + # lag_vars <- ifelse( + # test_idcol & !is.na(test_idcol), + # x[bot, var], + # NA) + lag_vars <- x[test_idcol * bot, vars] + names(lag_vars) <- paste0(vars, "_lag.", lag) - return(lag_var) + return(lag_vars) } - x_lag <- lapply(vars, function(y) .lag(x, lag = lag, y)) |> - do.call("data.frame", args = _) - names(x_lag) <- paste0(vars, "_lag.", lag) + # x_lag <- lapply(vars, function(y) .lag(x, lag = lag, y)) |> + # do.call("data.frame", args = _) + # names(x_lag) <- paste0(vars, "_lag.", lag) + x_lag <- .lag(x, lag, vars) # # reset inputs ---- # x <- .reset_inputs(cbind(x, x_lag), x_conversion) diff --git a/man/hz_lag.Rd b/man/hz_lag.Rd index bdd73f10c..f234da021 100644 --- a/man/hz_lag.Rd +++ b/man/hz_lag.Rd @@ -6,8 +6,8 @@ \usage{ hz_lag( object, - vars, lag = 1, + vars = NULL, idcol = "id", depthcols = c("top", "bottom"), order = FALSE @@ -16,10 +16,10 @@ hz_lag( \arguments{ \item{object}{a \code{data.frame}} -\item{vars}{character: column names, to lag.} - \item{lag}{integer: number of horizons to lag} +\item{vars}{character: column names, to lag.} + \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")}).} @@ -45,10 +45,10 @@ texcl = c("SL", "SL", "CL", "CL", "L"), clay = c(10, 12, 27, 35, 16) ) -h |> hz_lag(c("texcl", "clay"), 1) +h |> hz_lag() h |> -hz_lag(c("texcl", "clay"), 1) |> +hz_lag() |> cbind(h) |> transform( clay_dif = clay - clay_lag.1, From c4d37e9280e5131d7bd755b7f36b3c6649cf9f32 Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Wed, 22 May 2024 16:54:34 -0500 Subject: [PATCH 13/22] update hz_lag() --- R/segment.R | 98 ++++++++++++++++++++++++++++++++++++++++----------- man/hz_lag.Rd | 14 +++++--- 2 files changed, 86 insertions(+), 26 deletions(-) diff --git a/R/segment.R b/R/segment.R index b70530de8..4dc7a0965 100644 --- a/R/segment.R +++ b/R/segment.R @@ -534,8 +534,8 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { #' @description This function finds adjacent values to a horizon values at lagged distances. #' #' @param object a \code{data.frame} -#' @param vars character: column names, to lag. #' @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 #' @@ -560,21 +560,25 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { #' #' h |> hz_lag() #' +#' h |> hz_lag(-1) +#' +#' h |> hz_lag(10:15, unit = "depth") +#' #' h |> #' hz_lag() |> -#' cbind(h) |> +#' cbind(h, lag = _) |> #' transform( -#' clay_dif = clay - clay_lag.1, -#' texcl_contrast = paste0(texcl, "-", texcl_lag.1)) +#' clay_dif = lag.clay_bot.1 - clay, +#' texcl_contrast = paste0(texcl, "-", lag.texcl_bot.1)) #' -hz_lag <- function(object, lag = 1, vars = NULL, idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { +hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c("top", "bottom"), order = FALSE) { nm <- names(object) - idx <- which(! nm %in% c(idcol, depthcols)) - if (is.null(vars)) vars <- nm[idx] + idx_std <- which(! nm %in% c(idcol, depthcols)) + vars <- nm[idx_std] # check arguments ---- @@ -586,6 +590,20 @@ hz_lag <- function(object, lag = 1, vars = NULL, idcol = "id", depthcols = c("to 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 ---- @@ -595,31 +613,69 @@ hz_lag <- function(object, lag = 1, vars = NULL, idcol = "id", depthcols = c("to # lag ---- - .lag <- function(x, lag = lag, vars = NULL) { + .lag_ind <- function(x, lag = lag) { nr <- nrow(x) top <- 1:nr - bot <- c((1 + lag):nr, rep(NA, lag)) + 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_deps <- x$bot[top] == x$top[bot] - # lag_vars <- ifelse( - # test_idcol & !is.na(test_idcol), - # x[bot, var], - # NA) lag_vars <- x[test_idcol * bot, vars] - names(lag_vars) <- paste0(vars, "_lag.", lag) - + if (lag >= 0) names(lag_vars) <- paste0(vars, "_bot.", lag) + if (lag < 0) names(lag_vars) <- paste0(vars, "_top.", abs(lag)) + return(lag_vars) } - # x_lag <- lapply(vars, function(y) .lag(x, lag = lag, y)) |> - # do.call("data.frame", args = _) - # names(x_lag) <- paste0(vars, "_lag.", lag) - x_lag <- .lag(x, lag, vars) + + .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 <- .reset_inputs(cbind(x, x_lag), x_conversion) + x <- .reset_inputs(cbind(x, x_lag), x_conversion) + return(x_lag) } diff --git a/man/hz_lag.Rd b/man/hz_lag.Rd index f234da021..2a5f32df0 100644 --- a/man/hz_lag.Rd +++ b/man/hz_lag.Rd @@ -7,7 +7,7 @@ hz_lag( object, lag = 1, - vars = NULL, + unit = "index", idcol = "id", depthcols = c("top", "bottom"), order = FALSE @@ -18,7 +18,7 @@ hz_lag( \item{lag}{integer: number of horizons to lag} -\item{vars}{character: column names, to lag.} +\item{unit}{character: lag units in index or depth.} \item{idcol}{character: column name of the pedon ID within the object.} @@ -47,12 +47,16 @@ 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) |> +cbind(h, lag = _) |> transform( -clay_dif = clay - clay_lag.1, -texcl_contrast = paste0(texcl, "-", texcl_lag.1)) +clay_dif = lag.clay_bot.1 - clay, +texcl_contrast = paste0(texcl, "-", lag.texcl_bot.1)) } \author{ From 63a61048548afe72f71ddcd9efbfcc21904578ac Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Wed, 22 May 2024 17:10:43 -0500 Subject: [PATCH 14/22] Update segment.R --- R/segment.R | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/R/segment.R b/R/segment.R index 4dc7a0965..9fab4a6fc 100644 --- a/R/segment.R +++ b/R/segment.R @@ -382,28 +382,27 @@ hz_dissolve <- function(object, by, idcol = "id", depthcols = c("top", "bottom") # } - # reset inputs ---- - var_dep <- .reset_inputs(var_dep, df_conversion) - - # append dissolve_id n <- c( - var_dep[["top"]], - var_dep[["bottom"]] + var_dep$top, + var_dep$bot ) |> nchar() |> - max() + max(na.rm = TRUE) var_dep$dissolve_id <- paste0( - var_dep[[idcol]], + var_dep$idcol, "_", - formatC(var_dep[["top"]], width = n, flag = 0), + formatC(var_dep$top, width = n, flag = 0), "-", - formatC(var_dep[["bottom"]], width = n, flag = 0), + formatC(var_dep$bot, width = n, flag = 0), "_", - var_dep[["variable"]] + var_dep$variable ) + # reset inputs ---- + var_dep <- .reset_inputs(var_dep, df_conversion) + return(var_dep) } From 7e61dc5fb9a629f7ee71411fdb3e90e710995f1d Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Wed, 22 May 2024 17:38:13 -0500 Subject: [PATCH 15/22] Update segment.R --- R/segment.R | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/R/segment.R b/R/segment.R index 9fab4a6fc..da52ff50c 100644 --- a/R/segment.R +++ b/R/segment.R @@ -466,14 +466,9 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { ## length == 2 if (length(depthcols) != 2) stop("depthcols must length must equal 2") - ## check for matching column names - x_deps <- names(x[depthcols]) - y_deps <- names(y[depthcols]) - - if (!all.equal( - names(x[depthcols]), - names(y[depthcols])) - ) stop("both x and y must contain columns with names matching depthcols") + ## check for matching column names + .check_names(x, c(idcol, depthcols)) + .check_names(y, c(idcol, depthcols)) # check segment_id ---- @@ -504,19 +499,24 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { xi <- x yi <- y[which(y$idcol == xi$idcol[1]), ] - int <- c(xi$top, xi$bot, yi$top, yi$bot) |> - sort() |> + 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)) + + 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")) From 6344ca6ab266d0fcd7380b88b3af829af32cf2cc Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Thu, 23 May 2024 12:35:17 -0500 Subject: [PATCH 16/22] Update segment.R --- R/segment.R | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/R/segment.R b/R/segment.R index da52ff50c..96d3b4f1d 100644 --- a/R/segment.R +++ b/R/segment.R @@ -517,7 +517,7 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { 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")) + xy_int <- merge(x_seg, y_seg, by = c("segment_id", "idcol", "top", "bot"), sort = FALSE) # reset inputs ---- @@ -607,7 +607,7 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( # order ---- if (order) { - x[order(x$idcol, x$top, x$bot), ] + x <- x[order(x$idcol, x$top, x$bot), ] } @@ -619,12 +619,13 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( 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] - lag_vars <- x[test_idcol * bot, vars] - if (lag >= 0) names(lag_vars) <- paste0(vars, "_bot.", lag) - if (lag < 0) names(lag_vars) <- paste0(vars, "_top.", abs(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(lag_vars) + return(x_lag) } @@ -673,7 +674,7 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( # # reset inputs ---- - x <- .reset_inputs(cbind(x, x_lag), x_conversion) + x_lag <- .reset_inputs(x_lag, x_conversion) return(x_lag) From 8ea05879cf5336b28a1957af15db66b332aac20d Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 28 May 2024 14:37:05 -0500 Subject: [PATCH 17/22] adding aloc_taxpartsize() --- NAMESPACE | 1 + R/allocate.R | 202 ++++++++++++++++++++++++++++++++++++++++ R/segment.R | 8 +- R/texture.R | 29 +++--- man/aloc_taxpartsize.Rd | 81 ++++++++++++++++ 5 files changed, 307 insertions(+), 14 deletions(-) create mode 100644 man/aloc_taxpartsize.Rd diff --git a/NAMESPACE b/NAMESPACE index 78adeb991..1a141c7fe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,6 +15,7 @@ export(aggregateColor) export(aggregateSoilDepth) export(alignTransect) export(allocate) +export(aloc_taxpartsize) export(aqp.env) export(argillic.clay.increase.depth) export(barron.torrent.redness.LAB) diff --git a/R/allocate.R b/R/allocate.R index 3b75611bd..790a3e8d3 100644 --- a/R/allocate.R +++ b/R/allocate.R @@ -635,3 +635,205 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag return(sp) } + +#' @title Allocate Particle Size Control Class for the Control Section. +#' +#' @description This function aggregates information in the horizon table and allocates it to the particle size 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. +#' +#' @seealso \code{\link{texture_to_taxpartsize}} +#' +#' @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") +#' +#' aloc_taxpartsize(horizons(h), pscs) +#' +#' +aloc_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 + pscs_sc <- NULL + 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() + + x$rn <- 1:nrow(x) + # xy <- hz_intersect(x, y, idcol = idcol, depthcols = depthcols) + # x_sub <- x[x$rn %in% xy$rn, ] + + + # standardize inputs + x_std <- .standardize_inputs(x, idcol = idcol, depthcols = depthcols, clay = clay) + 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/segment.R b/R/segment.R index 96d3b4f1d..b416b8a34 100644 --- a/R/segment.R +++ b/R/segment.R @@ -396,7 +396,7 @@ hz_dissolve <- function(object, by, idcol = "id", depthcols = c("top", "bottom") "-", formatC(var_dep$bot, width = n, flag = 0), "_", - var_dep$variable + var_dep$value ) @@ -700,14 +700,16 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( # standardize inputs -.standardize_inputs <- function(x, idcol = NULL, hzidcol = NULL, depthcols = NULL) { +.standardize_inputs <- function(x, idcol = NULL, hzidcol = NULL, depthcols = NULL, texcl = NULL, clay = NULL) { # set new names var_names <- c( idcol = idcol, hzidcol = hzidcol, top = depthcols[1], - bot = depthcols[2] + bot = depthcols[2], + texcl = texcl, + clay = clay ) # find matches diff --git a/R/texture.R b/R/texture.R index 6be5cb168..f1c8e90bb 100644 --- a/R/texture.R +++ b/R/texture.R @@ -159,6 +159,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 +184,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 +201,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) { @@ -478,9 +484,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 +513,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 diff --git a/man/aloc_taxpartsize.Rd b/man/aloc_taxpartsize.Rd new file mode 100644 index 000000000..d560eb9ae --- /dev/null +++ b/man/aloc_taxpartsize.Rd @@ -0,0 +1,81 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/allocate.R +\name{aloc_taxpartsize} +\alias{aloc_taxpartsize} +\title{Allocate Particle Size Control Class for the Control Section.} +\usage{ +aloc_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") + +aloc_taxpartsize(horizons(h), pscs) + + +} +\seealso{ +\code{\link{texture_to_taxpartsize}} +} From 6aef52934a23d038fa7eb3a5e6ea4ae96c6dc838 Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 28 May 2024 15:07:39 -0500 Subject: [PATCH 18/22] fix when dissolve_id already present --- R/allocate.R | 25 +++++++++++++++++++++++-- R/segment.R | 5 +++-- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/R/allocate.R b/R/allocate.R index 790a3e8d3..115628a8d 100644 --- a/R/allocate.R +++ b/R/allocate.R @@ -714,8 +714,28 @@ aloc_taxpartsize <- function(x, y, taxpartsize = "taxpartsize", clay = "clay", i # 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) + 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 @@ -750,7 +770,8 @@ aloc_taxpartsize <- function(x, y, taxpartsize = "taxpartsize", clay = "clay", i clay_wt <- NULL xy_agg <- data.table::as.data.table(xy)[, - list(top = min(top, na.rm = TRUE), + 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 diff --git a/R/segment.R b/R/segment.R index b416b8a34..43d10c415 100644 --- a/R/segment.R +++ b/R/segment.R @@ -700,7 +700,7 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( # standardize inputs -.standardize_inputs <- function(x, idcol = NULL, hzidcol = NULL, depthcols = NULL, texcl = NULL, clay = NULL) { +.standardize_inputs <- function(x, idcol = NULL, hzidcol = NULL, depthcols = NULL, texcl = NULL, clay = NULL, taxpartsize = NULL) { # set new names var_names <- c( @@ -709,7 +709,8 @@ hz_lag <- function(object, lag = 1, unit = "index", idcol = "id", depthcols = c( top = depthcols[1], bot = depthcols[2], texcl = texcl, - clay = clay + clay = clay, + taxpartsize = taxpartsize ) # find matches From 81f3cbb3e6f86b441bd748b8320b590fbd32b932 Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 28 May 2024 17:40:35 -0500 Subject: [PATCH 19/22] renaming aloc_taxpartsize() & adding PSCS_levels() --- NAMESPACE | 3 +- R/allocate.R | 15 +++++----- R/segment.R | 6 ++-- man/PSCS_levels.Rd | 29 +++++++++++++++++++ man/allocate.Rd | 3 ++ man/hz_dissolve.Rd | 2 +- man/hz_intersect.Rd | 3 ++ man/hz_lag.Rd | 3 ++ man/hz_segment.Rd | 2 +- ...oc_taxpartsize.Rd => hz_to_taxpartsize.Rd} | 13 +++++---- man/texture.Rd | 5 ++++ 11 files changed, 67 insertions(+), 17 deletions(-) create mode 100644 man/PSCS_levels.Rd rename man/{aloc_taxpartsize.Rd => hz_to_taxpartsize.Rd} (88%) diff --git a/NAMESPACE b/NAMESPACE index 1a141c7fe..1e845e096 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) @@ -15,7 +16,6 @@ export(aggregateColor) export(aggregateSoilDepth) export(alignTransect) export(allocate) -export(aloc_taxpartsize) export(aqp.env) export(argillic.clay.increase.depth) export(barron.torrent.redness.LAB) @@ -100,6 +100,7 @@ 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/R/allocate.R b/R/allocate.R index 115628a8d..2f95505f1 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. #' @@ -656,7 +658,9 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag #' #' @return A \code{data.frame} object containing the original idcol and aggregated particle size control section allocation. #' -#' @seealso \code{\link{texture_to_taxpartsize}} +#' @author Stephen Roecker +#' +#' @seealso [texture_to_taxpartsize()], [PSCS_levels()] #' #' @export @@ -697,17 +701,14 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag #' pscs <- data.frame(id = h$id, rbind(estimatePSCS(h))) #' names(pscs)[2:3] <- c("top", "bottom") #' -#' aloc_taxpartsize(horizons(h), pscs) +#' hz_to_taxpartsize(horizons(h), pscs) #' #' -aloc_taxpartsize <- function(x, y, taxpartsize = "taxpartsize", clay = "clay", idcol = "id", depthcols = c("top", "bottom")) { +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 - pscs_sc <- NULL - 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() x$rn <- 1:nrow(x) # xy <- hz_intersect(x, y, idcol = idcol, depthcols = depthcols) @@ -804,7 +805,7 @@ aloc_taxpartsize <- function(x, y, taxpartsize = "taxpartsize", clay = "clay", i 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 + idx_sc = sc %in% .pscs_sc sc[idx_sc] = sc[idx_sc] }) xy_lag <- NULL diff --git a/R/segment.R b/R/segment.R index 43d10c415..1134c8489 100644 --- a/R/segment.R +++ b/R/segment.R @@ -16,7 +16,7 @@ #' #' @author Stephen Roecker #' -#' @seealso [dice()], [glom()] +#' @seealso [dice()], [glom()], [hz_dissolve()], [hz_lag()], [hz_intersect()] #' #' @export #' @@ -242,7 +242,7 @@ segment <- function(object, intervals, trim = TRUE, hzdepcols = c("top", "bottom #' #' @author Stephen Roecker #' -#' @seealso \code{\link{checkHzDepthLogic}} +#' @seealso [hz_lag()], [hz_intersect()], [hz_segment()] , [checkHzDepthLogic()] #' #' @export #' @@ -432,6 +432,7 @@ dissolve_hz <- function(object, by, id = "idcol", hztop = "top", hzbot = "bottom #' #' @author Stephen Roecker #' +#' @seealso [hz_dissolve()], [hz_lag()], [hz_segment()] #' #' @export #' @@ -544,6 +545,7 @@ hz_intersect <- function(x, y, idcol = "id", depthcols = c("top", "bottom")) { #' #' @author Stephen Roecker #' +#' @seealso [hz_dissolve()], [hz_intersect()], [hz_segment()] #' #' @export #' 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/hz_dissolve.Rd b/man/hz_dissolve.Rd index 19fc9f1f2..742725828 100644 --- a/man/hz_dissolve.Rd +++ b/man/hz_dissolve.Rd @@ -89,7 +89,7 @@ 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 index 23503244d..304e85a78 100644 --- a/man/hz_intersect.Rd +++ b/man/hz_intersect.Rd @@ -43,6 +43,9 @@ 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 index 2a5f32df0..1691dc65e 100644 --- a/man/hz_lag.Rd +++ b/man/hz_lag.Rd @@ -58,6 +58,9 @@ 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/hz_segment.Rd b/man/hz_segment.Rd index d2071bf15..db41673df 100644 --- a/man/hz_segment.Rd +++ b/man/hz_segment.Rd @@ -116,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/aloc_taxpartsize.Rd b/man/hz_to_taxpartsize.Rd similarity index 88% rename from man/aloc_taxpartsize.Rd rename to man/hz_to_taxpartsize.Rd index d560eb9ae..65e699a32 100644 --- a/man/aloc_taxpartsize.Rd +++ b/man/hz_to_taxpartsize.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/allocate.R -\name{aloc_taxpartsize} -\alias{aloc_taxpartsize} +\name{hz_to_taxpartsize} +\alias{hz_to_taxpartsize} \title{Allocate Particle Size Control Class for the Control Section.} \usage{ -aloc_taxpartsize( +hz_to_taxpartsize( x, y, taxpartsize = "taxpartsize", @@ -72,10 +72,13 @@ depths(h) <- id ~ top + bottom pscs <- data.frame(id = h$id, rbind(estimatePSCS(h))) names(pscs)[2:3] <- c("top", "bottom") -aloc_taxpartsize(horizons(h), pscs) +hz_to_taxpartsize(horizons(h), pscs) } \seealso{ -\code{\link{texture_to_taxpartsize}} +\code{\link[=texture_to_taxpartsize]{texture_to_taxpartsize()}}, \code{\link[=PSCS_levels]{PSCS_levels()}} +} +\author{ +Stephen Roecker } 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 } From fd2ce0b774563b91ef4a15d35d347926b34f6afb Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Tue, 28 May 2024 17:41:03 -0500 Subject: [PATCH 20/22] adding PSCS_levels() --- R/texture.R | 54 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/R/texture.R b/R/texture.R index f1c8e90bb..df1a84ea6 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}. @@ -424,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 @@ -1032,3 +1036,53 @@ 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") + + # 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 <- stats::lag(idx) + idx <- ifelse(idx < l & !is.na(l), idx * -1, idx * 1) + n <- length(idx) + idx = sum(idx * c(1, 0.1, 0.01)[1:n]) + }) + + 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() From 1b040b87f41423a806314d7eea582a0bd6af4beb Mon Sep 17 00:00:00 2001 From: Andrew Gene Brown Date: Tue, 28 May 2024 16:08:59 -0700 Subject: [PATCH 21/22] Consolidate/update diagnostic_hz/restrictions documentation --- R/Class-SoilProfileCollection.R | 14 ++++++---- R/SoilProfileCollection-setters.R | 23 ++++++++-------- man/diagnostic_hz-set.Rd | 40 ---------------------------- man/diagnostic_hz.Rd | 44 ++++++++++++++++++++++++++++--- man/restrictions-set.Rd | 39 --------------------------- man/restrictions.Rd | 43 +++++++++++++++++++++++++++--- 6 files changed, 99 insertions(+), 104 deletions(-) delete mode 100644 man/diagnostic_hz-set.Rd delete mode 100644 man/restrictions-set.Rd diff --git a/R/Class-SoilProfileCollection.R b/R/Class-SoilProfileCollection.R index 6ed74f8ca..4ea830da4 100644 --- a/R/Class-SoilProfileCollection.R +++ b/R/Class-SoilProfileCollection.R @@ -948,10 +948,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 @@ -968,9 +970,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/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/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) + } From 32839473b7da79f7f7c94753be5df0a508b055f2 Mon Sep 17 00:00:00 2001 From: Stephen Roecker Date: Wed, 29 May 2024 11:17:24 -0500 Subject: [PATCH 22/22] fix sorting --- R/allocate.R | 4 ++-- R/texture.R | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/R/allocate.R b/R/allocate.R index 2f95505f1..797367c7a 100644 --- a/R/allocate.R +++ b/R/allocate.R @@ -638,9 +638,9 @@ allocate <- function(..., to = c("FAO Salt Severity", "FAO Black Soil", "ST Diag } -#' @title Allocate Particle Size Control Class for the Control Section. +#' @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 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. diff --git a/R/texture.R b/R/texture.R index df1a84ea6..25c536e5e 100644 --- a/R/texture.R +++ b/R/texture.R @@ -1058,7 +1058,7 @@ fragvol_to_texmod <- function( 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") + 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") @@ -1067,10 +1067,13 @@ PSCS_levels <- function() { idx <- lapply(test, function(x) { idx <- sapply(x, function(y) which(fe == y)) |> unlist() - l <- stats::lag(idx) + # 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)[1:n]) + 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)