Algorithms for Quantitative Pedology

Description

The aqp (Algorithms for Quantitative Pedology) package for R was developed to address some of the difficulties associated with processing soils information, specifically related to visualization, aggregation, and classification of soil profile data. This package is based on a mix of S3/S4 functions and classes, and most functions use basic dataframes as input, where rows represent soil horizons and columns define properties of those horizons. Common to most functions are the requirements that horizon boundaries are defined as depth from 0, and that profiles are uniquely defined by an id column. The aqp package defines an S4 class, "SoilProfileCollection", for storage of profile-level metadata, as well as summary, print, and plotting methods that have been customized for common tasks related to soils data.

Details

Demos: demo(aqp)

Project homepage

Author(s)

Dylan E. Beaudette debeaudette@ucdavis.edu, Pierre Roudier, Andrew G. Brown

See Also

⁠depths<-()⁠, SoilProfileCollection(), sp1, sp2, sp3, sp4, sp5, sp6

Get data from column of horizon or site data in a SoilProfileCollection

Description

Get the data from a column accessed by name x$name. Column names other than profile ID are not shared between site and horizons.

Usage

## S4 method for signature 'SoilProfileCollection'
x$name

Arguments

Examples

data(sp1)

depths(sp1) <- id ~ top + bottom

# get data from a column by name (prop)
sp1$prop

Set data in column of horizon or site data in a SoilProfileCollection

Description

Set the data in a column accessed by name spc$name. Column names other than profile ID are not shared between site and horizons.

When using $<-, the length of input and output matching either the number of sites or number of horizons is used to determine which slot new columns are assigned to. Use site(x)$name <- value or horizons(x)$name <- value to be explicit about which slot is being accessed.

Usage

## S4 replacement method for signature 'SoilProfileCollection'
x$name <- value

Arguments

Get or Set Generalized Horizon Label (GHL) Column Name

Description

GHL(): Get column name containing generalized horizon labels

⁠GHL<-⁠: Set generalized horizon label column name

Usage

## S4 method for signature 'SoilProfileCollection'
GHL(object, required = FALSE)

## S4 replacement method for signature 'SoilProfileCollection'
GHL(object, required = FALSE) <- value

Arguments

Details

Store the column name containing generalized horizon labels in the metadata slot of the SoilProfileCollection.

Examples

data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# set horizon designation column
GHL(sp1) <- "name"

# get horizon designation column
GHL(sp1)

Subset SoilProfileCollection Objects or Horizons via checkHzDepthLogic

Description

This function removes profiles or horizons from a SoilProfileCollection that are flagged as having invalid horizon depth logic by checkHzDepthLogic. Invalid profiles may be created when setting byhz = TRUE; use caution as some functions may not work properly in the presence of gaps. Consider using fillHzGaps to fill these gaps.

Usage

HzDepthLogicSubset(x, byhz = FALSE)

Arguments

Value

a SoilProfileCollection object

Note

This function cannot identify (and remove) overlapping horizons when byhz = TRUE.

Create Representative Soil Profiles via L1 Estimator

Description

The L1 estimator, or geometric median, is a multivariate generalization of the (univariate) median concept. This function performs a multivariate aggregation (via L1 estimator) according to a suite of ratio-scale soil properties. The L1 estimator is applied to soil profile data that have been sliced to a 1-depth-unit basis. Data should be well stratified by groups defined in fm, otherwise the L1 median may not make any sense.

See the L1 Profiles Tutorial for additional examples.

Usage

L1_profiles(
  x,
  fm,
  basis = 1,
  method = c("regex", "simple", "constant"),
  maxDepthRule = c("max", "min"),
  maxDepthConstant = NULL
)

Arguments

Details

See this related tutorial for additional examples. The method, maxDepthRule, and maxDepthConstant arguments set the maximum depth (over the entire collection) of analysis used to build "L1 profiles". The following rules are available:

• method = 'regex' uses pattern matching on horizon designations (note that hzdesgnname metadata must be set with hzdesgnname(x) <- 'columnname')

• method = 'simple' uses min or max as applied to x, no accounting for non-soil horizons (e.g. Cr or R)

• method = 'constant' uses a fixed depth value supplied by maxDepthConstant

The maxDepthRule argument sets depth calculation constraint, applied to soil depths computed according to method (min or max).

Value

a SoilProfileCollection object

Note

This function requires the Gmedian package.

References

Cardot, H., Cenac, P. and Zitt, P-A. (2013). Efficient and fast estimation of the geometric median in Hilbert spaces with an averaged stochastic gradient algorithm. Bernoulli, 19, 18-43.

Numerical Classification of Soil Profiles

Description

Replacement for profile_compare().

Performs a numerical comparison of soil profiles using named properties, based on a weighted, summed, depth-segment-aligned dissimilarity calculation.

Variability in soil depth can interfere significantly with the calculation of between-profile dissimilarity– what is the numerical “distance” (or dissimilarity) between a slice of soil from profile A and the corresponding, but missing, slice from a shallower profile B? Gower's distance metric would yield a NULL distance, despite the fact that intuition suggests otherwise: shallower soils should be more dissimilar from deeper soils. For example, when a 25 cm deep profile is compared with a 50 cm deep profile, numerical distances are only accumulated for the first 25 cm of soil (distances from 26 - 50 cm are NULL). When summed, the total distance between these profiles will generally be less than the distance between two profiles of equal depth. Our algorithm has an option (setting replace_na=TRUE) to replace NULL distances with the maximum distance between any pair of profiles for the current depth slice. In this way, the numerical distance between a slice of soil and a corresponding slice of non-soil reflects the fact that these two materials should be treated very differently (i.e. maximum dissimilarity).

This alternative calculation of dissimilarities between soil and non-soil slices solves the problem of comparing shallow profiles with deeper profiles. However, it can result in a new problem: distances calculated between two shallow profiles will be erroneously inflated beyond the extent of either profile's depth. Our algorithm has an additional option (setting add_soil_flag=TRUE) that will preserve NULL distances between slices when both slices represent non-soil material. With this option enabled, shallow profiles will only accumulate mutual dissimilarity to the depth of the deeper profile.

Slices are classified as 'soil' down to the maximum depth to which at least one of variables used in the dissimilarity calculation is not NA. This will cause problems when profiles within a collection contain all NAs within the columns used to determine dissimilarity. An approach for identifying and removing these kind of profiles is presented in the examples section below.

A notice is issued if there are any NA values within the matrix used for distance calculations, as these values are optionally replaced by the max dissimilarity.

Our approach builds on the work of (Moore, 1972) and the previously mentioned depth-slicing algorithm.

Usage

NCSP(
  x,
  vars,
  fm = NULL,
  weights = rep(1, times = length(vars)),
  maxDepth = max(x),
  k = 0,
  isColor = FALSE,
  rescaleResult = FALSE,
  progress = TRUE,
  verbose = TRUE,
  returnDepthDistances = FALSE
)

Arguments

Note

NCSP() will overwrite the removed.profiles metadata from x.

Author(s)

Dylan E. Beaudette and Jon Maynard

References

• J.J Maynard, S.W. Salley, D.E. Beaudette, J.E Herrick. Numerical soil classification supports soil identification by citizen scientists using limited, simple soil observations. Soil Sci. Soc. Am. J. 2020; 84: 1675-1692. doi:10.1002/saj2.20119.

• D.E. Beaudette, P. Roudier, A.T. O'Geen, Algorithms for quantitative pedology: A toolkit for soil scientists, Computers & Geosciences, Volume 52, 2013, Pages 258-268, ISSN 0098-3004, doi:10.1016/j.cageo.2012.10.020.

• Moore, A.; Russell, J. & Ward, W. Numerical analysis of soils: A comparison of three soil profile models with field classification. Journal of Soil Science, 1972, 23, 194-209.

See Also

dice(), cluster::daisy(), compareSites()

Average Hydraulic Parameters from the ROSETTA Model by USDA Soil Texture Class

Description

Average soil hydraulic parameters generated by the first stage predictions of the ROSETTA model by USDA soil texture class. These data were extracted from ROSETTA documentation and re-formatted for ease of use.

Usage

data(ROSETTA.centroids)

Format

A data frame:

texture

soil texture class, ordered from low to high available water holding capacity

theta_r

average saturated water content

theta_s

average residual water content

alpha

average value, related to the inverse of the air entry suction, log10-transformed values with units of cm

npar

average value, index of pore size distribution, log10-transformed values with units of 1/cm

theta_r_sd

1 standard deviation of theta_r

theta_s_sd

1 standard deviation of theta_s

alpha_sd

1 standard deviation of alpha

npar_sd

1 standard deviation of npar

sat

approximate volumetric water content at which soil material is saturated

fc

approximate volumetric water content at which matrix potential = -33kPa

pwp

approximate volumetric water content at which matrix potential = -1500kPa

awc

approximate available water holding capacity: VWC(-33kPa)

• VWC(-1500kPa)

Details

Theoretical water-retention parameters for uniform soil material of each texture class have been estimated via van Genuchten model.

See the related tutorial

Source

ROSETTA Class Average Hydraulic Parameters

References

van Genuchten, M.Th. (1980). "A closed-form equation for predicting the hydraulic conductivity of unsaturated soils". Soil Science Society of America Journal. 44 (5): 892-898.

Schaap, M.G., F.J. Leij, and M.Th. van Genuchten. 2001. ROSETTA: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 251(3–4): 163-176.

Examples

## Not run: 

library(aqp)
library(soilDB)
library(latticeExtra)

data("ROSETTA.centroids")

# iterate over horizons and generate VG model curve
res <- lapply(1:nrow(ROSETTA.centroids), function(i) {
  m <- KSSL_VG_model(VG_params = ROSETTA.centroids[i, ], phi_min = 10^-3, phi_max=10^6)$VG_curve
  # copy generalized hz label
  m$hz <- ROSETTA.centroids$hz[i]
  # copy ID
  m$texture_class <- ROSETTA.centroids$texture[i]
  return(m)
})

# copy over lab sample number as ID
res <- do.call('rbind', res)

# check: OK
str(res)


# visual check: OK
xyplot(
  phi ~ theta | texture_class, data=res,
  type=c('l', 'g'),
  scales=list(alternating=3, x=list(tick.number=10), y=list(log=10, tick.number=10)),
  yscale.components=yscale.components.logpower,
  ylab=expression(Suction~~(kPa)),
  xlab=expression(Volumetric~Water~Content~~(cm^3/cm^3)),
  par.settings = list(superpose.line=list(col='RoyalBlue', lwd=2)),
  strip=strip.custom(bg=grey(0.85)),
  as.table=TRUE
)



## End(Not run)

Fix Overlap within a Sequence via Simulated Annealing

Description

This function makes small adjustments to elements of x until overlap defined by thresh is removed, or until maxIter is reached. Rank order and boundary conditions (defined by min.x and max.x) are preserved. The underlying algorithm is based on simulated annealing. The "cooling schedule" parameters T0 and k can be used to tune the algorithm for specific applications.

Usage

SANN_1D(
  x,
  thresh = 0.6,
  adj = thresh * 2/3,
  min.x = min(x) - 0.2,
  max.x = max(x) + 0.2,
  maxIter = 1000,
  trace = FALSE,
  tiny = 1e-04,
  T0 = 500,
  k = 10,
  ...
)

Arguments

Details

Ideas for solving difficult overlap scenarios:

• widen the boundary conditions by adjusting min.x and max.x beyond the original scale of x

• reduce the allowable overlap threshold thresh

• reduce the magnitude of perturbations (adj) and increase maxIter

• increase k

Value

When trace = FALSE, a vector of the same length as x, preserving rank-ordering and boundary conditions. When trace = TRUE a list containing the new sequence along with information about objective functions and decisions made during iteration.

Author(s)

D.E. Beaudette and K.C. Thompson

See Also

electroStatics_1D(), fixOverlap()

Examples

x <- c(1, 2, 3, 3.4, 3.5, 5, 6, 10)

# easy
z <- fixOverlap(x, thresh = 0.2, trace = TRUE)

# harder
z <- fixOverlap(x, thresh = 0.6, trace = TRUE)

# much harder
z <- fixOverlap(x, thresh = 0.9, trace = TRUE)


# interpret `trace` output

# relatively challenging
x <- c(1, 2, 3.4, 3.4, 3.4, 3.4, 6, 8, 10, 12, 13, 13, 15, 15.5)

# fix overlap, return debugging information
set.seed(10101)
z <- fixOverlap(x, thresh = 0.8, trace = TRUE)

# setup plot device
par(mar = c(4, 4, 1, 1))
layout(matrix(c(1,2,3)), widths = 1, heights = c(1,1,2))

# objective function = overlap + SSD
plot(
  seq_along(z$stats), z$stats, 
  type = 'h', las = 1,
  xlab = 'Iteration', ylab = 'Overlap',
  cex.axis = 0.8
)

# SSD: deviation from original configuration 
plot(
  seq_along(z$ssd), z$ssd, 
  type = 'h', las = 1,
  xlab = 'Iteration', ylab = 'Deviation',
  cex.axis = 0.8
)
# adjustments at each iteration
matplot(
  z$states, type = 'l', 
  lty = 1, las = 1, 
  xlab = 'Iteration', ylab = 'x-position'
)

# trace log
# B: boundary condition violation
# O: rank (order) violation
# +: accepted perturbation
# -: rejected perturbation
table(z$log)

Example SoilProfileCollection with Overlapping Horizons

Description

A SoilProfileCollection with overlapping horizons, derived from a Dynamic Soil Properties project.

Usage

data(SPC.with.overlap)

Format

A SoilProfileCollection

An S4 object representation of a group of soil profiles.

Description

In general, one should use depths() to initiate a SoilProfileCollection object from data. However, sometimes there are instances where either an empty, or very specific, object is needed. If that is the case, the general constructor SoilProfileCollection is available.

Usage

SoilProfileCollection(
  idcol = "id",
  hzidcol = "hzID",
  depthcols = c("top", "bottom"),
  metadata = list(aqp_df_class = "data.frame", aqp_group_by = "", aqp_hzdesgn = "",
    aqp_hztexcl = "", stringsAsFactors = FALSE),
  horizons = data.frame(id = character(0), hzID = character(0), top = numeric(0), bottom
    = numeric(0), stringsAsFactors = FALSE),
  site = data.frame(id = character(0), stringsAsFactors = FALSE),
  diagnostic = data.frame(stringsAsFactors = FALSE),
  restrictions = data.frame(stringsAsFactors = FALSE)
)

Arguments

Details

After aqp 2.0.2, the ⁠@sp⁠ slot was removed from the SoilProfileCollection object. If you run into errors related to old object definitions, use rebuildSPC() on the problematic object.

Slots

idcol

character.

hzidcol

character.

depthcols

character.

metadata

list.

horizons

data.frame.

site

data.frame.

diagnostic

data.frame.

restrictions

data.frame.

Author(s)

Pierre Roudier, Dylan E. Beaudette, Andrew G. Brown

Examples

## structure of default, empty SoilProfileCollection
str(SoilProfileCollection())


## use the depths() formula interface to specify
## profile ID, top and bottom depth and set up
## a SPC that is topologically correct and complete

d <- do.call('rbind', lapply(1:10, random_profile))

# promote to SoilProfileCollection and plot
depths(d) <- id ~ top + bottom
plot(d)

# split into new SoilProfileCollection objects by index
d.1 <- d[1, ]
d.2 <- d[2, ]
d.345 <- d[3:5, ]

# combine profile collections
# note that profiles are sorted according to ID
d.new <- c(d.345, d.1, d.2)
plot(d.new)

data(sp1)

## init SoilProfileCollection objects from data.frame
depths(sp1) <- id ~ top + bottom

## depth units
du <- depth_units(sp1)
depth_units(sp1) <- 'in'
depth_units(sp1) <- du

## horizon designation column
hzdesgnname(sp1) <- "name"
hzdesgnname(sp1)

## all designations in an SPC (useful for single profile SPC)
hzDesgn(sp1)

## horizon texture class column
hztexclname(sp1) <- "texture"
hztexclname(sp1)

## get/set metadata on SoilProfileCollection objects
# this is a 1-row data.frame
m <- metadata(sp1)
m$sampler <- 'Dylan'
metadata(sp1) <- m

## extract horizon data from SoilProfileCollection objects as data.frame
h <- horizons(sp1)

# also merge (left-join) of new columns and
# replacement of existing columns via horizons<-
horizons(sp1) <- h

# get number of horizons
nrow(sp1)


## getting site-level data
site(sp1)

## setting site-level data
# site-level data from horizon-level data (stored in @horizons)
site(sp1) <- ~ group


# make some fake site data, and append from data.frame
# a matching ID column must be present in both @site and new data
# note that IDs should all be character class
d <- data.frame(id=profile_id(sp1), p=runif(n=length(sp1)), stringsAsFactors=FALSE)
site(sp1) <- d

# edit horizon depths
horizonDepths(sp1) <- c('t', 'b')
horizonDepths(sp1)

# edit profile IDs
p <- sprintf("%s-new", profile_id(sp1))
profile_id(sp1) <- p
profile_id(sp1)

Ranking Systems for USDA Soil Texture Classes

Description

Generate a vector of USDA soil texture codes or class names, sorted according to approximate particle size

Usage

SoilTextureLevels(which = "codes", simplify = FALSE)

Arguments

Value

an ordered factor

References

Field Book for Describing and Sampling Soils, version 3.0

Examples

# class codes
SoilTextureLevels()

# class names
SoilTextureLevels(which = 'names')

# simpler class names
SoilTextureLevels(which = 'names', simplify = TRUE)
 

Matrix/data.frame-like access to profiles and horizons in a SoilProfileCollection

Description

You can access the contents of a SoilProfileCollection by profile and horizon "index", i and j, respectively: spc[i, j, ...]. Subset operations are propagated to other slots (such as diagnostics or spatial) when they result in removal of sites from a collection.

• i refers to the profile position within the collection. By default the order is based on the C SORT order of the variable that you specified as your unique profile ID at time of object construction. Note that if your ID variable was numeric, then it has been sorted as a character.

• j refers to the horizon or "slice" index. This index is most useful when either a) working with slice'd SoilProfileCollection or b) working with single-profile collections. j returns the layer in the specified index positions for all profiles in a collection.

• ... is an area to specify an expression that is evaluated in the subset. Currently supported

  • .LAST (last horizon in each profile): return the last horizon from each profile. This uses i but ignores the regular j index.

  • .FIRST (first horizon in each profile): return the last horizon from each profile. This uses i but ignores the regular j index.

  • .HZID (horizon index): return the horizon indices corresponding to i+j+... ("k") constraints

  • .NHZ (number of horizons): return the number of horizons in the profiles resulting from i+j+... ("k") constraints

Usage

## S4 method for signature 'SoilProfileCollection'
x[i, j, ..., drop = TRUE]

Arguments

Get column of horizon or site data in a SoilProfileCollection

Description

Get the data from a column accessed by name. Column names other than profile ID are not shared between site and horizons. Bonus: [[ gives access to all site and horizon level variables in tab complete for RStudio using the magrittr pipe operator!

Usage

## S4 method for signature 'SoilProfileCollection'
x[[i, j]]

Arguments

Examples

data(sp2)
depths(sp2) <- id ~ top + bottom
site(sp2) <- ~ surface

# get with [[
sp2[['surface']]

# get using "unknown" expression:
#  "2nd + 3rd horizon column names"
for(i in horizonNames(sp2)[2:3])
 print(sp2[[i]])

data(sp5)

# some column names to work with
rgb.columns <- c("R25","G25","B25")

res <- lapply(rgb.columns, function(x) {

  # [[ allows you to access column names in a loop
  round(sp5[[x]] * 255)

})

# rename scaled results
names(res) <- paste0(rgb.columns,"_scl")

# add horizon ID to results
result <- data.frame(hzID = hzID(sp5), do.call('cbind', res))
head(result)

# join result back into horizons
horizons(sp5) <- result

Add or change column of horizon or site data in a SoilProfileCollection

Description

Add or change the data from a column accessed by name. Column names other than profile ID are not shared between site and horizons. The benefit of using double bracket setter over $ is that name can be calculated, whereas with $, it must be known a priori and hard coded.

When using the double bracket setter the length of input and output matching either the number of sites or number of horizons is used to determine which slot new columns are assigned to.

Usage

## S4 replacement method for signature 'SoilProfileCollection'
x[[i]] <- value

Arguments

Accumulate horizon depths, and reflect reversed depths, relative to new datum

Description

Fix old-style organic horizon depths or depths with a non-standard datum by the "depth accumulation" method.

Usage

accumulateDepths(
  x,
  id = NULL,
  hzdepths = NULL,
  hzname = NULL,
  hzdatum = 0,
  seqnum = NULL,
  pattern = "O",
  fix = TRUE
)

Arguments

Details

The "depth accumulation" method calculates thicknesses of individual horizons and then cumulative sums them after putting them in id + top depth order. The routine tries to determine context based on hzname and pattern. The main transformation is if a top depth is deeper than the bottom depth, the depths are reflected on the Z-axis (made negative). The data are then id + top depth sorted again, the thickness calculated and accumulated to replace the old depths.

This function uses several heuristics to adjust data before transformation and thickness calculation:

Regex matching of horizon designation patterns and similar

• matches of pattern where both top and bottom depth NA -> ⁠[0,1]⁠ ⁠[top,bottom]⁠ depth

• REMOVE horizons that do not match pattern where both top and bottom depths NA

Over-ride hzname handling with the sequence column argument seqnum

• if seqnum column specified "first record with NA hzname" is considered a pattern match if seqnum == 1

Trigger "fixing" with the fix argument:

• Add 1 cm to bottom-most horizons with NA bottom depth

• Add 1 cm thickness to horizons with top and bottom depth equal

• Add 1 cm thickness to horizons with NA top depth and bottom depth 0

Value

A horizon-level data.frame, suitable for promoting to SPC with ⁠depths<-⁠, or a SoilProfileCollection, depending on the class of x.

Examples

# example using hzdatum argument
data(sp4)
depths(sp4) <- id ~ top + bottom
hz <- accumulateDepths(sp4,
                       id = "id",
                       hzdepths = c("top", "bottom"),
                       hzname = "name",
                       hzdatum = 5 * 1:length(sp4))
plot(hz)
  
# example using old-style O horizons
hz <- read.table(text = "peiidref hzdept hzdepb hzname seqnum phiid
                    1        11      0      5      A      2   295
                    2        11      1      0     Oe      1   294
                    3        11      5     13     C1      3   296
                    4        11     13     58     C2      4   297
                    5        11     58    152     C3      5   298
                    6        13      0      5      A      2   303
                    7        13      1      0     Oe      1   302
                    8        13      5     25     Bw      3   304
                    9        13     25     61      C      4   305
                    10       13     61     NA      R      5   306
                    11      136      0     13     A1      3   695
                    12      136      1      0     Oe      2   694
                    13      136      2      1     Oi      1   693
                    14      136     13     61     C1      4   696
                    15      136     61     76     C2      5   697")
                    
depths(hz) <- peiidref ~ hzdept + hzdepb

hz_fixed <- accumulateDepths(hz,
                             id = "peiidref",
                             hzdepths = c("hzdept", "hzdepb"),
                             hzname = "hzname")
is_valid <- checkHzDepthLogic(hz_fixed)$valid

test0 <- subset(hz_fixed, !is_valid)
test1 <- subset(hz_fixed, is_valid)

origO <- subset(hz, grepl("O", hzname))
fixedO <- subset(hz_fixed, grepl("O", hzname))

par(mfrow = c(2, 1), mar = c(0, 0, 3, 2))

plotSPC(origO, max.depth = 25)
plotSPC(fixedO, max.depth = 25)

Add Depth Brackets

Description

Add depth brackets to soil profile sketches.

Usage

addBracket(
  x,
  label.cex = 0.75,
  tick.length = 0.05,
  arrow.length = 0.05,
  offset = -0.3,
  missing.bottom.depth = NULL,
  ...
)

Arguments

Details

x may contain multiple records per profile. Additional examples can be found in this tutorial.

Note

This is a low-level plotting function: you must first plot a SoilProfileCollection object before using this function.

Author(s)

D.E. Beaudette

See Also

addDiagnosticBracket, plotSPC

Examples

# sample data
data(sp1)

# add color vector
sp1$soil_color <- with(sp1, munsell2rgb(hue, value, chroma))

# promote to SoilProfileCollection
depths(sp1) <- id ~ top + bottom

# plot profiles
par(mar = c(0, 0, 0, 1))
plotSPC(sp1, width = 0.3)

# extract min--max depths associated with all A horizons
# result is a single-row data.frame / profile
combinedBracket <- function(i) {
  h <- horizons(i)
  idn <- idname(i)
  this.id <- h[[idn]][1]
  
  idx <- grep('^A', h$name)
  
  res <- data.frame(
    id = this.id,
    top = min(h$top[idx]), 
    bottom = max(h$bottom[idx], na.rm=TRUE)
  )
  names(res)[1] <- idn
  
  return(res)
}

# return matching horizon top / bottom depths for A or C horizons
# result is a 0 or more row data.frame / profile
individualBrackets <- function(i) {
  h <- horizons(i)
  idn <- idname(i)
  this.id <- h[[idn]][1]
  
  idx <- grep('^A|^C', h$name)
  
  res <- data.frame(
    id = this.id,
    top = h$top[idx], 
    bottom = h$bottom[idx]
  )
  names(res)[1] <- idn
  
  return(res)
}

# combined brackets
b1 <- profileApply(sp1, combinedBracket, frameify = TRUE)

# individual brackets
b2 <- profileApply(sp1, individualBrackets, frameify = TRUE)

# plot in reverse order
plotSPC(sp1, plot.order = rev(1:length(sp1)), width = 0.25)

# note that plotting order is derived from the call to `plotSPC(sp1)`
addBracket(b1, col='red', offset = -0.35)

# plot in reverse order
plotSPC(sp1, plot.order = rev(1:length(sp1)), width = 0.25)

# note that plotting order is derived from the call to `plotSPC(sp1)`
addBracket(b2, col='red', offset = -0.35)

Annotate Diagnostic Features

Description

Annotate diagnostic features within a sketch of soil profiles.

Usage

addDiagnosticBracket(
  s,
  kind,
  feature = "featkind",
  top = "featdept",
  bottom = "featdepb",
  ...
)

Arguments

Details

Additional examples can be found in this tutorial.

Note

This is a low-level plotting function: you must first plot a SoilProfileCollection object before using this function.

Author(s)

D.E. Beaudette

See Also

addBracket(), plotSPC()

Symbolize Volume Fraction within a Soil Profile Collection Plot

Description

Symbolize volume fraction on an existing soil profile collection plot.

Usage

addVolumeFraction(
  x,
  colname,
  res = 10,
  cex.min = 0.1,
  cex.max = 0.5,
  pch = 1,
  col = "black"
)

Arguments

Details

This function can only be called after plotting a SoilProfileCollection object. Details associated with a call to plotSPC() are automatically accounted for within this function: e.g. plot.order, width, etc..

Note

It may be necessary to adjust both res, cex.min, and cex.max for optimal legibility.

Author(s)

D.E. Beaudette

See Also

plotSPC()

Summarize Soil Colors

Description

Summarize soil color data, weighted by occurrence and horizon thickness.

Usage

aggregateColor(
  x,
  groups = "genhz",
  col = "soil_color",
  k = NULL,
  profile_wt = NULL,
  mixingMethod = c("estimate", "exact")
)

Arguments

Details

Weights are computed by: w_i = sqrt(sum(thickness_i)) * n_i where w_i is the weight associated with color i, thickness_i is the total thickness of all horizons associated with the color i, and n_i is the number of horizons associated with color i. Weights are computed within groups specified by groups.

See the related tutorial for additional examples.

Value

A list with the following components:

• scaled.data: a list of colors and associated weights, one item for each generalized horizon label with at least one color specified in the source data

• aggregate.data: a data.frame of weighted-mean colors, one row for each generalized horizon label with at least one color specified in the source data

Author(s)

D.E. Beaudette

See Also

generalize.hz(), aggregateColorPlot()

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# load some example data
data(sp1, package = 'aqp')

# upgrade to SoilProfileCollection and convert Munsell colors
sp1$soil_color <- with(sp1, munsell2rgb(hue, value, chroma))
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# generalize horizon names
n <- c('O', 'A', 'B', 'C')
p <- c('O', 'A', 'B', 'C')
sp1$genhz <- generalize.hz(sp1$name, n, p)

# aggregate colors over horizon-level attribute: 'genhz'
a <- aggregateColor(sp1, groups = 'genhz', col = 'soil_color')

# check results
str(a)

# simple visualization
aggregateColorPlot(a)

Plot aggregate soil color data

Description

Generate a plot from summaries generated by aggregateColor().

Usage

aggregateColorPlot(
  x,
  print.label = TRUE,
  label.font = 1,
  label.cex = 0.65,
  label.orientation = c("v", "h"),
  buffer.pct = 0.02,
  print.n.hz = FALSE,
  rect.border = "black",
  horizontal.borders = FALSE,
  horizontal.border.lwd = 2,
  x.axis = TRUE,
  y.axis = TRUE,
  ...
)

Arguments

Details

See the related tutorial for additional examples.

Value

no data are returned, function is called for graphical output

Author(s)

D.E. Beaudette

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# load some example data
data(sp1, package = 'aqp')

# upgrade to SoilProfileCollection and convert Munsell colors
sp1$soil_color <- with(sp1, munsell2rgb(hue, value, chroma))
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# generalize horizon names
n <- c('O', 'A', 'B', 'C')
p <- c('O', 'A', 'B', 'C')
sp1$genhz <- generalize.hz(sp1$name, n, p)

# aggregate colors over horizon-level attribute: 'genhz'
a <- aggregateColor(sp1, groups = 'genhz', col = 'soil_color')

# check results
str(a)

# simple visualization
aggregateColorPlot(a)

Probabilistic Estimation of Soil Depth within Groups

Description

Estimate the most-likely depth to contact within a collection of soil profiles. Consider getSoilDepthClass followed by group-wise percentile estimation as a faster alternative.

Usage

aggregateSoilDepth(
  x,
  groups,
  crit.prob = 0.9,
  name = hzdesgnname(x),
  p = "Cr|R|Cd",
  ...
)

Arguments

Details

This function computes a probability-based estimate of soil depth by group. If no grouping variable exists, a dummy value can be used to compute a single estimate. The crit.prob argument sets the critical probability (e.g. 0.9) at which soil depth within a group of profiles is determined. For example, a crit.prob of 0.95 might result in an estimated soil depth (e.g. 120cm) where 95% of the profiles (by group) had depths that were less than or equal to 120cm.

Value

A data.frame is returned, with as many rows as there are unique group labels, as specified in groups.

Author(s)

D.E. Beaudette

See Also

estimateSoilDepth() slab()

Examples

data(sp1)
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# set horizon designation in SPC
hzdesgnname(sp1) <- 'name'

aggregateSoilDepth(sp1, 'group', crit.prob = 0.9)

Calculate Relative Positions from Transect Data

Description

This function is used to support relative positioning of soil profiles by plotSPC, based on transect or gradient values typically associated with a site level attribute (e.g. elevation). Gradient values specified in x are translated to the range used by plotSPC (usually ⁠1, length(SPC)⁠) specified in x.min and x.max.

Usage

alignTransect(x, x.min, x.max, fix = TRUE, ...)

Arguments

Details

See the Pair-Wise Distances by Generalized Horizon Labels tutorial for additional examples.

Value

list containing:

• grad: values of x in ascending order

• order: ordering vector of x

• relative.pos: elements of x translated to the new relative scale defined by x.min and x.max

Examples

data("sierraTransect")

# split transects
g <- subset(sierraTransect, transect == 'Granite')
a <- subset(sierraTransect, transect == 'Andesite')

g.p <- alignTransect(g$elev, x.min = 1, x.max = length(g), fix = FALSE)
a.p <- alignTransect(a$elev, x.min = 1, x.max = length(a), fix = FALSE)

op <- par(mar=c(2,0,0,2), mfrow=c(2,1))

plotSPC(g, width=0.25, name.style='center-center', 
        cex.names=0.75, 
        relative.pos = g.p$relative.pos, plot.order = g.p$order)

axis(1, at = g.p$relative.pos, labels = g.p$grad, line = -1.5)

plotSPC(a, width=0.25, name.style='center-center', 
        cex.names=0.75, 
        relative.pos = a.p$relative.pos, plot.order = a.p$order)

axis(1, at = a.p$relative.pos, labels = a.p$grad, line = -1.5)


par(op)

Allocate soil properties within various classification systems.

Description

Generic function to allocate soil properties to different classification schemes.

Usage

allocate(
  ...,
  to = c("FAO Salt Severity", "FAO Black Soil", "ST Diagnostic Features"),
  droplevels = FALSE
)

Arguments

Details

This function is intended to allocate a set of soil properties to an established soil classification scheme, such as Salt Severity or Black Soil. Allocation is semantically different from classification. While classification is the 'act' of developing a grouping scheme, allocation is the assignment or identification of measurements to a established class (Powell, 2008).

Usage Details

Each classification scheme (to argument) uses a different set of arguments.

• ⁠FAO Salt Severity⁠

  • EC: electrical conductivity column name, dS/m

  • pH: pH column name, saturated paste extract

  • ESP: exchangeable sodium percentage column name, percent

• ⁠FAO Black Soils⁠

  • object: a data.frame or SoilProfileCollection

  • pedonid: pedon ID column name, required when object is a data.frame

  • hztop: horizon top depth column name, required when object is a data.frame

  • hzbot: horizon bottom depth column name, required when object is a data.frame

  • OC: organic carbon column name, percent

  • m_chroma: moist Munsell chroma column name

  • m_value: moist Munsell value column name

  • d_value: dry Munsell value column name

  • CEC: cation exchange capacity column name (NH4OAc at pH 7), units of cmol(+)/kg soil

  • BS: base saturation column name (NH4OAc at pH 7), percent

  • tropical: logical, data are associated with "tropical soils"

• ⁠ST Diagnostic Features⁠

  • object: a data.frame or SoilProfileCollection

  • pedonid: pedon ID column name, required when object is a data.frame

  • hzname: horizon name column, required when object is a data.frame

  • hztop: horizon top depth column name, required when object is a data.frame

  • hzbot: horizon bottom depth column name, required when object is a data.frame

  • texcl: soil texture class (USDA) column name

  • rupresblkcem: rupture resistance column name

  • m_value: moist Munsell value column name

  • m_chroma: moist Munsell chroma column name

  • d_value: dry Munsell value column name

  • BS: base saturation column name (method ??), percent

  • OC: organic carbon column name, percent

  • n_value: ??

  • featkind: ??

Value

A vector or data.frame object.

Note

The results returned by allocate(to = "ST Diagnostic Features") currently return a limited set of diagnostic features that are easily defined. Also, the logic implemented for some features does not include all the criteria defined in the Keys to Soil Taxonomy.

Author(s)

Stephen Roecker

References

Abrol, I., Yadav, J. & Massoud, F. 1988. Salt-affected soils and their management. No. Bulletin 39. Rome, FAO Soils.

FAO. 2006. Guidelines for soil description. Rome, Food and Agriculture Organization of the United Nations.

FAO. 2020. DEFINITION | What is a black soil? (online). (Cited 28 December 2020). http://www.fao.org/global-soil-partnership/intergovernmental-technical-panel-soils/gsoc17-implementation/internationalnetworkblacksoils/more-on-black-soils/definition-what-is-a-black-soil/es/

Powell, B., 2008. Classifying soil and land, in: McKenzie, N.J., Grundy, M.J., Webster, R., Ringrose-Voase, A.J. (Eds.), Guidelines for Survey Soil and Land Resources, Australian Soil and Land Survey Handbook Series. CSIRO, Melbourne, p. 572.

Richards, L.A. 1954. Diagnosis and Improvement of Saline and Alkali Soils. U. S. Government Printing Office. 166 pp.

Soil Survey Staff, 2014. Keys to Soil Taxonomy, 12th ed. USDA-Natural Resources Conservation Service, Washington, D.C.

Examples

# Salt Severity
test <- expand.grid(
  EC  = sort(sapply(c(0, 0.75, 2, 4, 8, 15, 30), function(x) x + c(0, -0.05, 0.05))),
  pH  = c(8.1, 8.2, 8.3, 8.4, 8.5, 8.6),
  ESP = sort(sapply(c(0, 15, 30, 50, 70, 100), function(x) x + c(0, 0.1, -0.1)))
)
test$ss      <- with(test, allocate(EC = EC, pH = pH, ESP = ESP, to = "FAO Salt Severity"))
table(test$ss)

# Black Soil Category 1 (BS1)
test <- expand.grid(
  dept = seq(0, 50, 10),
  OC   = sort(sapply(c(0, 0.6, 1.2, 20, 40), function(x) x + c(0, -0.05, 0.05))),
  chroma_moist  = 2:4,
  value_moist   = 2:4,
  value_dry     = 4:6,
  thickness     = 24:26,
  CEC           = 24:26,
  BS            = 49:51,
  tropical      = c(TRUE, FALSE)
)
test$pedon_id <- rep(1:21870, each = 6)
test$depb     <- test$dept + 10

bs1 <- allocate(test, pedonid = "pedon_id", hztop = "dept", hzbot = "depb", 
                OC = "OC", m_chroma = "chroma_moist", m_value = "value_moist", 
                d_value = "value_dry", CEC = "CEC", BS = "BS", 
                to = "FAO Black Soil"
)

table(BS1 = bs1$BS1, BS2 = bs1$BS2)


# SoilProfileCollection interface

data(sp3)
depths(sp3) <- id ~ top + bottom
hzdesgnname(sp3) <- 'name'

# fake base saturation
horizons(sp3)$bs <- 75

plotSPC(sp3)

allocate(
  sp3, 
  to = 'FAO Black Soil', 
  OC = 'tc', 
  m_chroma = 'chroma', 
  m_value = 'value', 
  d_value = 'value',
  CEC = 'cec',
  BS = 'bs'
)

# make a copy and edit horizon values
x <- sp3
x$value <- 2
x$chroma <- 2
x$cec <- 26
x$tc <- 2

x$soil_color <- munsell2rgb(x$hue, x$value, x$chroma)

plotSPC(x)

allocate(
  x, 
  to = 'FAO Black Soil', 
  OC = 'tc', 
  m_chroma = 'chroma', 
  m_value = 'value', 
  d_value = 'value',
  CEC = 'cec',
  BS = 'bs'
)


# Soil Taxonomy Diagnostic Features
data(sp1)
sp1$texcl = gsub("gr|grv|cbv", "", sp1$texture)
df <- allocate(object = sp1, pedonid = "id", hzname = "name", 
               hzdept = "top", hzdepb = "bottom", texcl = "texcl", 
               to = "ST Diagnostic Features"
)
aggregate(featdept ~ id, data = df, summary)

Get aqp_df_class entry from metadata or return a safe value.

Description

This is an accessor and replacement method for the aqp_df_class entry in the metadata slot. This entry is used internally by methods that interact with data.frame objects and slots to ensure that the same class used to promote to the SoilProfileCollection initially is used throughout the process.

Usage

## S4 method for signature 'SoilProfileCollection'
aqp_df_class(object)

## S4 replacement method for signature 'SoilProfileCollection'
aqp_df_class(object) <- value

Arguments

Return upper boundary of argillic horizon

Description

Returns the top depth of the argillic horizon as a numeric vector.

Usage

argillic.clay.increase.depth(p, clay.attr = "clay")

Arguments

Details

Uses crit.clay.argillic to determine threshold clay increase, and get.increase.matrix to determine where increase is met within a vertical distance of 30 cm.

Value

A numeric vector containing top depth of argillic horizon, if present, or NA.

Author(s)

Andrew Gene Brown

See Also

getArgillicBounds, get.increase.matrix, crit.clay.argillic

Examples

data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

p <- sp1[1]
attr <- 'prop' # clay contents 
foo <- argillic.clay.increase.depth(p, clay.attr = attr)
foo

Coerce SoilProfileCollection with as()

Description

SoilProfileCollections can be coerced to other R object types using as(spc, 'type').

Possible endpoints include: list, data.frame, SpatialPointsDataFrame and SpatialPoints.

Usage

## S4 method for signature 'SoilProfileCollection'
as.data.frame(x)

Arguments

Value

list

data.frame

tbl_df

data.table

SpatialPointsDataFrame

sf

SpatialPoints

Examples

# load example data stored as SoilProfileCollection
data(sp5)

# sp5
str(sp5)

# list output
str(as(sp5, 'list'))

# data.frame output
str(as(sp5, 'data.frame'))

# Spatial Objects
# make some random coordinate data for each profile
sp5$x <- sp5$y <- rnorm(length(sp5))
initSpatial(sp5, crs = "OGC:CRS84") <- ~ x + y

if (requireNamespace("sf")) {

  # sf output
  str(as(sp5, 'sf'))

  # SpatialPointsDataFrame output
  str(as(sp5, 'SpatialPointsDataFrame'))

  # SpatialPoints output
  str(as(sp5, 'SpatialPoints'))
  
}

Barron & Torrent (1986) Redness Index in LAB color space

Description

Calculate Redness Index after Barron & Torrent (1986) "Use of the Kubelka—Munk Theory to Study the Influence of Iron Oxides on Soil Colour" using Munsell colors converted to LAB. DOI: 10.1111/j.1365-2389.1986.tb00382.x. Accepts vectorized inputs for hue, value and chroma, produces vector output.

Usage

barron.torrent.redness.LAB(hue, value, chroma)

Arguments

Value

A numeric vector of horizon redness index (higher values = redder).

Author(s)

Andrew G. Brown

References

Barron, V. and Torrent, J. (1986), Use of the Kubelka-Munk theory to study the influence of iron oxides on soil colour. Journal of Soil Science, 37: 499-510. doi:10.1111/j.1365-2389.1986.tb00382.x

Bootstrap Soil Texture Data

Description

Simulate realistic sand/silt/clay values (a composition) using multivariate Normal distribution or Dirichlet distribution. Simulations from the multivariate Normal distribution are based on the compositional mean and variance-covariance matrix. Simulations from the Dirichlet distribution are based on maximum likelihood estimation of alpha parameters.

Usage

bootstrapSoilTexture(ssc, method = c("dirichlet", "normal"), n = 100)

Arguments

Details

Simulations from the multivariate normal distribution will more closely track the marginal distributions of sand, silt, and clay–possibly a better fit for "squished" compositions (TODO elaborate). However, these simulations can result in extreme (unlikely) estimates.

Simulations from the Dirichlet distribution will usually be a better fit (fewer extreme estimates) but require a fairly large number of records in ssc (n >= 30?) for a reliable fit.

Additional examples will be added to this tutorial.

Value

a list containing:

• samples - data.frame of simulated sand, silt, clay values

• mean - compositional mean

• var - compositional variance-covariance matrix

• D.alpha - (fitted) alpha parameters of the Dirichlet distribution, NULL when method = 'normal'

Author(s)

D.E. Beaudette

References

Aitchison, J. (1986) The Statistical Analysis of Compositional Data Monographs on Statistics and Applied Probability. Chapman & Hall Ltd., London (UK). 416p.

Aitchison, J, C. Barcel'o-Vidal, J.J. Egozcue, V. Pawlowsky-Glahn (2002) A concise guide to the algebraic geometric structure of the simplex, the sample space for compositional data analysis, Terra Nostra, Schriften der Alfred Wegener-Stiftung, 03/2003

Malone Brendan, Searle Ross (2021) Updating the Australian digital soil texture mapping (Part 1*): re-calibration of field soil texture class centroids and description of a field soil texture conversion algorithm. Soil Research. https://www.publish.csiro.au/SR/SR20283

Malone Brendan, Searle Ross (2021) Updating the Australian digital soil texture mapping (Part 2*): spatial modelling of merged field and lab measurements. Soil Research. https://doi.org/10.1071/SR20284

Examples

if(
requireNamespace("compositions") &
  requireNamespace("soiltexture")
) {
  
  # sample data, data.frame
  data('sp4')
  
  # filter just Bt horizon data
  ssc <- sp4[grep('^Bt', sp4$name), c('sand', 'silt', 'clay')]
  names(ssc) <- toupper(names(ssc))
  
  # simulate 100 samples
  s <- bootstrapSoilTexture(ssc, n = 100)
  s <- s$samples
  
  # empty soil texture triangle
  TT <- soiltexture::TT.plot(
    class.sys= "USDA-NCSS.TT",
    main= "",
    tri.sum.tst=FALSE,
    cex.lab=0.75,
    cex.axis=0.75,
    frame.bg.col='white',
    class.lab.col='black',
    lwd.axis=1.5,
    arrows.show=TRUE,
    new.mar = c(3, 0, 0, 0)
  )
  
  # add original data points
  soiltexture::TT.points(
    tri.data = s, geo = TT, col='firebrick', 
    pch = 3, cex = 0.5, lwd = 1, 
    tri.sum.tst = FALSE
  )
  
  # add simulated points
  soiltexture::TT.points(
    tri.data = ssc, geo = TT, bg='royalblue', 
    pch = 22, cex = 1, lwd = 1, 
    tri.sum.tst = FALSE
  )
  
  # simple legend
  legend('top', 
         legend = c('Source', 'Simulated'), 
         pch = c(22, 3), 
         col = c('black', 'firebrick'), 
         pt.bg = c('royalblue', NA), 
         horiz = TRUE, bty = 'n'
  )
  
  
}

Multinominal Brier Score

Description

Compute a multinominal Brier score from predicted class probabilities and observed class label. Lower values are associated with a more accurate classifier.

Usage

brierScore(x, classLabels, actual = "actual")

Arguments

Value

a single Brier score, representative of data in x

Author(s)

D.E. Beaudette

References

Brier, Glenn W. 1950. "Verification of Forecasts Expressed in Terms of Probability." Monthly Weather Review 78 (1): 1-3. doi:10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2.

Examples

# columns 'a', 'b', 'c' contain predicted probabilities
# column 'actual' contains observed class label

# a good classifier
d.good <- data.frame(
  a = c(0.05, 0.05, 0.10),
  b = c(0.90, 0.85, 0.75),
  c = c(0.05, 0.10, 0.15),
  actual = c('b', 'b', 'b'),
  stringsAsFactors = FALSE
)

# a rather bad classifier
d.bad <- data.frame(
  a = c(0.05, 0.05, 0.10),
  b = c(0.90, 0.85, 0.75),
  c = c(0.05, 0.10, 0.15),
  actual = c('c', 'c', 'c'),
  stringsAsFactors = FALSE
)

# class labels are factors
d.factors <- data.frame(
  a = c(0.05, 0.05, 0.10),
  b = c(0.90, 0.85, 0.75),
  c = c(0.05, 0.10, 0.15),
  actual = c('b', 'b', 'b'),
  stringsAsFactors = TRUE
)

# relatively low value = accurate
brierScore(x = d.good, classLabels = c('a', 'b', 'c'), actual = 'actual')

# high values = not accuate
brierScore(x = d.bad, classLabels = c('a', 'b', 'c'), actual = 'actual')

# message related to conversion of factor -> character
brierScore(x = d.factors, classLabels = c('a', 'b', 'c'), actual = 'actual')

Buntley-Westin (1965) Index

Description

Calculate "Color Development Equivalent" by the method of Buntley & Westin (1965) "A Comparative Study of Developmental Color in a Chestnut-Chernozem-Brunizem Soil Climosequence" DOI: 10.2136/sssaj1965.03615995002900050029x. Originally developed for Mollisols, the Buntley-Westin index has been used as a tool to separate soils based on depth to particular colors.

Usage

buntley.westin.index(hue, chroma)

Arguments

Value

A numeric vector reflecting horizon color development.

Author(s)

Andrew G. Brown

References

Buntley, G.J. and Westin, F.C. (1965), A Comparative Study of Developmental Color in a Chestnut-Chernozem-Brunizem Soil Climosequence. Soil Science Society of America Journal, 29: 579-582. doi:10.2136/sssaj1965.03615995002900050029x

Combine SoilProfileCollection objects

Description

Combine SoilProfileCollection objects or lists of SoilProfileCollection objects. This method provides ... expansion for the pbindlist method.

Usage

## S4 method for signature 'SoilProfileCollection'
c(x, ...)

## S4 method for signature 'SoilProfileCollection'
combine(...)

## S4 method for signature 'list'
combine(...)

Arguments

Value

A SoilProfileCollection

Examples

# example data
spc1 <- random_profile(1, SPC = TRUE)
spc2 <- random_profile(2, SPC = TRUE)
spc3 <- random_profile('A', SPC = TRUE)

# combine into a single SPC, ... interface
spc <- combine(spc1, spc2, spc3)

# combine into a single SPC, list interface
spc <- combine(list(spc1, spc2, spc3))

# input are combined into a single SPC
spc <- c(spc1, spc2, spc3)

# result is a list when a mixture of objects are provided
spc <- c(spc1, bar=spc2, baz="foo")

Soil Data from the Central Sierra Nevada Region of California

Description

Site and laboratory data from soils sampled in the central Sierra Nevada Region of California.

Usage

data(ca630)

Format

List containing:

$site : A data frame containing site information.

user_site_id

national user site id

mlra

the MLRA

county

the county

ssa

soil survey area

lon

longitude, WGS84

lat

latitude, WGS84

pedon_key

national soil profile id

user_pedon_id

local soil profile id

cntrl_depth_to_top

control section top depth (cm)

cntrl_depth_to_bot

control section bottom depth (cm)

sampled_taxon_name

soil series name

$lab : A data frame containing horizon information.

pedon_key

national soil profile id

layer_key

national horizon id

layer_sequence

horizon sequence number

hzn_top

horizon top (cm)

hzn_bot

horizon bottom (cm)

hzn_desgn

horizon name

texture_description

USDA soil texture

nh4_sum_bases

sum of bases extracted by ammonium acetate (pH 7)

ex_acid

exchangeable acidity [method ?]

CEC8.2

cation exchange capacity by sum of cations method (pH 8.2)

CEC7

cation exchange capacity by ammonium acetate (pH 7)

bs_8.2

base saturation by sum of cations method (pH 8.2)

bs_7

base saturation by ammonium acetate (pH 7)

Details

These data were extracted from the NSSL database. ca630 is a list composed of site and lab data, each stored as data.frame objects. These data are modeled by a 1:many (site:lab) relation, with the pedon_id acting as the primary key in the site table and as the foreign key in the lab table.

Note

These data are out of date. Pending some new data + documentation. Use with caution

Source

https://ncsslabdatamart.sc.egov.usda.gov/

Examples

## Not run: 
library(tactile)
library(lattice)
library(Hmisc)
library(sp)

# check the data out:
data(ca630)
str(ca630)

# note that pedon_key is the link between the two tables

# make a copy of the horizon data
ca <- ca630$lab

# promote to a SoilProfileCollection class object
depths(ca) <- pedon_key ~ hzn_top + hzn_bot

# add site data, based on pedon_key
site(ca) <- ca630$site

# ID data missing coordinates: '|' is a logical OR
(missing.coords.idx <- which(is.na(ca$lat) | is.na(ca$lon)))

# remove missing coordinates by safely subsetting
if(length(missing.coords.idx) > 0)
	ca <- ca[-missing.coords.idx, ]

# register spatial data
initSpatial(ca) <- ~ lon + lat

# assign a coordinate reference system
prj(ca) <- 'EPSG:4269'

# check the result
print(ca)

# aggregate %BS 7 for all profiles into 1 cm slices
a <- slab(ca, fm= ~ bs_7)

# plot median & IQR by 1 cm slice
xyplot(
top ~ p.q50, 
data = a, 
lower=a$p.q25, 
upper=a$p.q75,
alpha=0.5,
ylim=c(160,-5), 
scales = list(alternating = 1, y = list(tick.num = 7)),
panel = panel.depth_function, 
prepanel = prepanel.depth_function,
ylab='Depth (cm)', xlab='Base Saturation at pH 7',
par.settings = tactile.theme(superpose.line = list(col = 'black', lwd = 2))
)

# aggregate %BS at pH 8.2 for all profiles by MLRA, along 1 cm slices
# note that mlra is stored in @site
a <- slab(ca, mlra ~ bs_8.2)

# keep only MLRA 18 and 22
a <- subset(a, subset=mlra %in% c('18', '22'))

# plot median & IQR by 1 cm slice, using different colors for each MLRA
xyplot(
top ~ p.q50, 
groups = factor(mlra), 
data = a,
lower=a$p.q25, 
upper=a$p.q75,
alpha=0.25,
sync.colors = TRUE,
ylim=c(160,-5), 
scales = list(alternating = 1, y = list(tick.num = 7)),
panel = panel.depth_function, 
prepanel = prepanel.depth_function,
ylab='Depth (cm)', xlab='Base Saturation at pH 7',
par.settings = tactile.theme(superpose.line = list(lwd = 2)),
auto.key = list(lines = TRUE, points = FALSE, columns = 2)
)



# Extract the 2nd horizon from all profiles as SPDF
ca.2 <- ca[, 2]

# subset profiles 1 through 10
ca.1.to.10 <- ca[1:10, ]

# basic plot method: profile plot
par(mar = c(0, 0, 3, 1))
plotSPC(ca.1.to.10, name='hzn_desgn', color = 'CEC7')

## End(Not run)

Check a SoilProfileCollection object for errors in horizon depths.

Description

This function inspects a SoilProfileCollection object, looking for four common errors in horizon depths:

1. bottom depth shallower than top depth

2. equal top and bottom depth

3. missing top or bottom depth (e.g. NA)

4. gap or overlap between adjacent horizons (only if byhz = FALSE)

Usage

checkHzDepthLogic(
  x,
  hzdepths = NULL,
  idname = NULL,
  fast = FALSE,
  byhz = FALSE
)

Arguments

Value

A data.frame containing profile IDs, validity boolean (valid) and test results if fast = FALSE.

The data.frame will have as many rows as profiles in x (length(x)).

• id : Profile IDs, named according to idname(x)

• valid : boolean, profile passes all of the following tests

• depthLogic : boolean, errors related to depth logic

• sameDepth : boolean, errors related to same top/bottom depths

• missingDepth : boolean, NA in top / bottom depths

• overlapOrGap : boolean, gaps or overlap in adjacent horizons (NA when byhz = TRUE)

Author(s)

D.E. Beaudette, A.G. Brown, S.M. Roecker

Examples

## sample data

data(sp3)
depths(sp3) <- id ~ top + bottom

# these data should be clean
res <- checkHzDepthLogic(sp3)

head(res)

# less memory if only concerned about net validity
res <- checkHzDepthLogic(sp3, fast = TRUE)

head(res)

Test for a valid SoilProfileCollection

Description

Test for a valid SoilProfileCollection

Usage

checkSPC(x)

Arguments

Details

Test for valid SoilProfileCollection by checking for slots defined in the class prototype. Likely only used between major versions of aqp where internal structure of SoilProfileCollection has changed. Use checkHzDepthLogic to check for common errors in horizon depths.

Value

TRUE or FALSE. Consider using rebuildSPC() if FALSE.

Author(s)

D.E. Beaudette

See Also

rebuildSPC, checkHzDepthLogic

Convert colors into Munsell Notation

Description

Lookup the n closest Munsell chips from the munsell lookup table from various color notations. This function replaces rgb2munsell().

Usage

col2Munsell(col, space = c("sRGB", "CIELAB"), nClosest = 1)

Arguments

Value

an (NA-padded) data.frame containing hue, value, chroma, and CIE delta-E 2000 color contrast metric between source and nearest matching color(s).

Note

This function is fully vectorized and will pad output with NA-records when NA are present in color.

Author(s)

D.E. Beaudette

References

http://ncss-tech.github.io/AQP/ http://www.brucelindbloom.com/index.html?ColorCalcHelp.html https://www.munsellcolourscienceforpainters.com/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html http://www.cis.rit.edu/mcsl/online/munsell.php

Examples

# vector of named R colors
col2Munsell(c('red', 'green', 'blue'))

# sRGB matrix in the range of 0-255
col2Munsell(cbind(255, 0, 0))

# sRGB matrix in the range of 0-1
col2Munsell(cbind(1, 0, 0))

# 10YR 5/6 in CIELAB
col2Munsell(
  cbind(51.4337, 9.917916, 38.6889), 
  space = 'CIELAB'
)

# 2.5YR 6/8 in hex notation
col2Munsell("#D18158FF")

# 7.5YR 8/1 in sRGB {0, 1}
col2Munsell(
  cbind(0.8240707, 0.7856834, 0.7541048)
)

# 7.5YR 8/1 in sRGB {0, 255}
col2Munsell(
  cbind(0.8240707, 0.7856834, 0.7541048) * 255
)

# multple colors in CIELAB
col2Munsell(
  parseMunsell(c('10BG 6/6', '2.5YR 4/6'), returnLAB = TRUE),
  space = 'CIELAB'
)

# data.frame input
col2Munsell(
  data.frame(r = 1, g = 0, b = 0),
  space = 'sRGB'
)

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# Munsell notation to sRGB triplets {0, 1} 
color <- munsell2rgb(
  the_hue = c('10YR', '2.5YR', '5YR'), 
  the_value = c(3, 5, 2.5), 
  the_chroma = c(5, 6, 2), 
  return_triplets = TRUE
)

# result is a data.frame of sRGB {0, 1}
color

# back-transform sRGB -> closest Munsell color
# sigma is the dE00 color contrast metric
col2Munsell(color, space = 'sRGB')

Collapse Horizons within Profiles Based on Pattern Matching

Description

Combines layers and aggregates data by grouping adjacent horizons which match pattern in hzdesgn or, alternately, share a common value in by argument. Numeric properties are combined using the weighted average, and other properties are derived from the dominant condition based on thickness of layers and values in each group.

Usage

collapseHz(
  x,
  pattern = NULL,
  by = NULL,
  hzdesgn = hzdesgnname(x, required = TRUE),
  FUN = function(x, pattern, hzdesgn, ...) grepl(pattern, x[[hzdesgn]], ignore.case =
    FALSE),
  ...,
  AGGFUN = NULL,
  ignore_numerics = NULL,
  na.rm = FALSE
)

Arguments

Details

If a custom matching function (FUN) is used, it should accept arbitrary additional arguments via an ellipsis (...). It is not necessary to do anything with arguments, but the result should match the number of horizons found in the input SoilProfileCollection x.

Custom aggregation functions defined in the AGGFUN argument should either return a single vector value for each group*column combination, or should return a data.frame object with named columns. If the input column name is used as a column name in the result data.frame, then the values of that column name in the result SoilProfileCollection will be replaced by the output of the aggregation function. See examples.

Value

A SoilProfileCollection

Author(s)

Andrew G. Brown

See Also

hz_dissolve()

Examples

data(jacobs2000)

# calculate a new SPC with genhz column based on patterns
new_labels <- c("A", "E", "Bt", "Bh", "C")
patterns <- c("A", "E", "B.*t", "B.*h", "C")
jacobs2000_gen <- generalizeHz(jacobs2000, new = new_labels, pattern = patterns)

# use existing generalized horizon labels
i <- collapseHz(jacobs2000_gen, by = "genhz")

profile_id(i) <- paste0(profile_id(i), "_collapse")

plot(
  c(i, jacobs2000),
  color = "genhz",
  name = "name",
  name.style = "center-center",
  cex.names = 1
)
 
# custom pattern argument  
j <- collapseHz(jacobs2000,
                c(
                  `A` = "^A",
                  `E` = "E",
                  `Bt` = "[ABC]+t",
                  `C` = "^C",
                  `foo` = "bar"
                ))
profile_id(j) <- paste0(profile_id(j), "_collapse")
plot(c(j, jacobs2000), color = "clay")

# custom aggregation function for matrix_color_munsell
k <- collapseHz(jacobs2000,
                pattern = c(
                  `A` = "^A",
                  `E` = "E",
                  `Bt` = "[ABC]+t",
                  `C` = "^C",
                  `foo` = "bar"
                ),
                AGGFUN = list(
                  matrix_color_munsell = function(x, top, bottom) {
                    thk <- bottom - top
                    if (length(x) > 1) {
                      xord <- order(thk, decreasing = TRUE)
                      paste0(paste0(x[xord], " (t=", thk[xord], ")"), collapse = ", ")
                    } else
                      x
                  }
                )
              )
profile_id(k) <- paste0(profile_id(k), "_collapse_custom")

unique(k$matrix_color_munsell)

# custom aggregation function for matrix_color_munsell (returns data.frame)
m <- collapseHz(jacobs2000,
                pattern = c(
                  `A` = "^A",
                  `E` = "E",
                  `Bt` = "[ABC]+t",
                  `C` = "^C",
                  `foo` = "bar"
                ),
                AGGFUN = list(
                  matrix_color_munsell = function(x, top, bottom) {
                    thk <- bottom - top
                    if (length(x) > 1) {
                      xord <- order(thk, decreasing = TRUE)
                      data.frame(matrix_color_munsell = paste0(x, collapse = ";"),
                                 n_matrix_color = length(x))
                    } else {
                      data.frame(matrix_color_munsell = x, 
                                 n_matrix_color = length(x))
                    }
                  }
                )
              )
profile_id(m) <- paste0(profile_id(m), "_collapse_custom")

m$matrix_color_munsell.n_matrix_color

Visualize soil colors in Munsell notation according to within-group frequency.

Description

Visualize soil colors in Munsell notation according to within-group frequency.

Usage

colorChart(
  m,
  g = factor("All"),
  size = TRUE,
  annotate = FALSE,
  chip.cex = 3,
  chip.cex.min = 0.1,
  chip.cex.max = 1.5,
  chip.border.col = "black",
  annotate.cex = chip.cex * 0.25,
  annotate.type = c("count", "percentage"),
  threshold = NULL
)

Arguments

Value

a trellis object

Examples

library(lattice)

# two hue pages
ric <- expand.grid(
  hue = c('5YR', '7.5YR'),
  value = 2:8,
  chroma = 2:8
)

# combine hue, value, chroma into standard Munsell notation
ric <- sprintf("%s %s/%s", ric$hue, ric$value, ric$chroma)

# note that chip frequency-based size is disabled
# because all chips have equal frequency
colorChart(ric, chip.cex = 4, size = TRUE)

# annotation of frequency
colorChart(ric, chip.cex = 4, annotate = TRUE)

# bootstrap to larger size
ric.big <- sample(ric, size = 100, replace = TRUE)

# frequency can be encoded in size
colorChart(ric.big, chip.cex = 3)
colorChart(ric.big, chip.cex = 5, annotate = TRUE)

# constant size
colorChart(ric.big, chip.cex = 3, size = FALSE)
colorChart(ric.big, chip.cex = 3, size = FALSE, chip.border.col = 'NA')

# simulate colors based dE00 thresholding
p <- list(
  list(m = '10YR 4/4', thresh = 10, hues = c('10YR', '7.5YR'))
)

# perform 500 simulations
s <- simulateColor(method = 'dE00', n = 500, parameters = p)

# result is a list, use the first element
colorChart(s[[1]], chip.cex = 4)

# increase the possible range of color chip sizes
colorChart(s[[1]], chip.cex = 4, chip.cex.min = 0.01, chip.cex.max = 2)

# slightly funky support for neutral hues
N <- sprintf('N %s/', 2:8)
cols <- c(rep(N, times = 5), ric.big)

# note special panel used to show neutral hues
colorChart(cols, size = FALSE, annotate = TRUE)

# filter proportions below given threshold
colorChart(cols, size = FALSE, annotate = TRUE, threshold = 0.01,
           chip.cex = 4, annotate.type = 'percentage')

Metrics of Contrast Suitable for Comparing Soil Colors

Description

Pair-wise comparisons of Munsell color specifications, based on the NCSS color contrast classes (Soil Survey Technical Note 2) and CIE delta-E 2000 metric.

Usage

colorContrast(m1, m2)

Arguments

Details

This function is fully vectorized but expects input to be of the same length. Use expand.grid() to generate suitable input from 1:many or many:1 type comparisons. See this tutorial for an expanded discussion and more examples. Neutral colors are not mentioned in SSTN2: in this function any comparison to a neutral color (e.g. 'N 3/') are assigned a delta-hue of 1. Since SSTN2 expects hues to be counted clock wise from 5R, it possible to get very large delta-hue values for otherwise adjacent colors: '5R' vs. '2.5R'. This will be addressed in an update to the standards.

The most meaningful representation of color contrast is the CIE2000 (dE00) metric.

Value

data.frame with the following columns:

• m1: Munsell color 1

• m2: Munsell color 2

• dH: delta-hue, as computed by huePosition

• dV: delta-value, absolute value of difference in Munsell value (m1 vs. m2)

• dc: delta-chroma, absolute value of difference in Munsell chroma (m1 vs. m2)

• dE00: delta-E00, e.g. the CIE delta-E as refined in 2000

• cc: soil color contrast class, as specified in Soil Survey Technical Note 2.

Note

delta-E00 is computed by the farver package.

Author(s)

D.E. Beaudette

References

1. https://en.wikipedia.org/wiki/Color_difference

See Also

colorContrastPlot, huePosition, huePositionCircle

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# two sets of colors to compare
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3')

# contrast metrics
colorContrast(m1, m2)

# adjacent chips
colorContrast('10YR 3/3', '10YR 3/4')
colorContrast('10YR 3/3', '7.5YR 3/3')

# highly contrasting colors
# http://colour.granjow.net/fabercastell-polychromos.html
colorContrastPlot('10B 4/13', '10YR 10/15',
labels = c('helioblue-reddish', 'light cadmium yellow')
)


## Note: neutral hues aren't defined in TN2
# approximation / extension of the concept
colorContrast(m1 = 'N 3/', m2 = 'N 6/')

#
colorContrast(m1 = '10YR 3/3', m2 = 'N 3/')

m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', 'N 3/')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3')
colorContrast(m1, m2)

Color Contrast Plot

Description

A simple display of two sets of colors, NCSS color contrast class and CIE delta-E00.

Usage

colorContrastPlot(
  m1,
  m2,
  col.cex = 1,
  col.font = 2,
  d.cex = 1,
  cc.font = 3,
  dE00.font = 1,
  labels = c("m1", "m2"),
  label.cex = 1,
  label.font = 1,
  printMetrics = TRUE,
  ...
)

Arguments

Note

This function requires the farver package for calculation of CIE delta-E00.

Author(s)

D.E. Beaudette

See Also

colorContrast()

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# two sets of colors to compare
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3')

# contrast metrics
colorContrast(m1, m2)

# graphical display
colorContrastPlot(m1, m2)

Soil Color Range via Quantiles

Description

Estimate central tendency and spread of soil color using marginal quantiles and L1 median of CIELAB coordinates.

Usage

colorQuantiles(soilColors, p = c(0.05, 0.5, 0.95))

Arguments

Details

Colors are converted from sRGB to CIELAB (D65 illuminant), marginal quantiles of (L,A,B) coordinates are estimated, and L1 median (L,A,B) is estimates. The closest Munsell chips (via Munsell/CIELAB lookup table provided by munsell) and R colors are determined by locating chips closest to the marginal quantiles and L1 median.

The results can be conveniently inspected using plotColorQuantiles().

Value

A List containing the following elements:

• marginal: data.frame containing marginal quantiles in CIELAB (D65), closest Munsell chips, and dE00

• L1: L1 median CIELAB (D65) values, closest Munsell chip, and dE00

Author(s)

D.E. Beaudette

Examples

## Not run: 
# example data, see manual page for details
data(sp5)

# slice top 25 cm
# 24-25cm is the last slice
s <- dice(sp5, 0:24 ~ .)

# check some of the data
par(mar=c(0,0,0,0))
plotSPC(sample(s, 25), divide.hz = FALSE, name = '', print.id = FALSE, width = 0.5)

# colors
previewColors(unique(s$soil_color))

# compute marginal quantiles and L1 median
cq <- colorQuantiles(s$soil_color)

# simple graphical display of results
plotColorQuantiles(cq)

## End(Not run)

Compare Site Level Attributes of a SoilProfileCollection

Description

Compare site level attributes of a SoilProfileCollection object, returning a distance matrix conformal with the output from NCSP(). Values are within the range of 0-1.

Usage

compareSites(x, vars, weights = rep(1, times = length(vars)), ...)

Arguments

Details

This function is typically used in conjunction with the output from NCSP().

Value

dissimilarity / dist class object containing pair-wise distances, row/column names derived from profile_id(x)

See Also

NCSP() cluster::daisy()

Return a list representation of site and horizon level data

Description

compositeSPC() is a convenience function that returns a named list representation of the columns from the @site and @horizons slots.

Usage

compositeSPC(object)

Arguments

Value

A list.

Author(s)

Andrew G. Brown.

Confusion Index

Description

Calculate the confusion index of Burrough et al., 1997.

Usage

confusionIndex(x)

Arguments

Value

A single numeric value.

Author(s)

D.E. Beaudette

References

Burrough, P.A., P.F.M. van Gaans, and R. Hootsmans. 1997. "Continuous Classification in Soil Survey: Spatial Correlation, Confusion and Boundaries." Geoderma 77: 115-35. doi:10.1016/S0016-7061(97)00018-9.

Examples

# a very simple example
p <- c(0.25, 0.25, 0.4, 0.05, 0.05)
confusionIndex(p)

# for comparison
shannonEntropy(p)

Color Contrast Chart

Description

Compare one or more pages from a simulated Munsell book of soil colors to a reference color.

Usage

contrastChart(
  m,
  hues,
  ccAbbreviate = 1,
  style = "hue",
  gridLines = FALSE,
  de00.cex = 0.6,
  cc.cex = 0.6,
  thresh = NULL,
  returnData = FALSE
)

Arguments

Details

A simulated Munsell color book page or pages are used to demonstrate color contrast between all chips and the reference color m (highlighted in red). NCSS color contrast class and CIE delta-E00 values are printed below all other color chips. Munsell color chips for chroma 5 and 7 are omitted, but axis labels are retained as a reminder of this fact.

Setting style='hue' emphasizes the contrast classes and CIE delta-E00 of chips adjacent to m. Setting style='CC' emphasizes adjacent chips according to respective contrast class via lattice panels.

Two-way panels are used when multiple hues are provided and style='CC'. The default output can be greatly enhanced via:

latticeExtra::useOuterStrips(..., strip = strip.custom(bg=grey(0.85)), strip.left = strip.custom(bg=grey(0.85)) )

Author(s)

D.E. Beaudette

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# single hue page
contrastChart(m = '10YR 3/3', hues = '10YR')

# multiple hue pages
contrastChart(m = '10YR 3/3', hues = c('10YR', '2.5Y'))

# contrast class, single hue
contrastChart(m = '10YR 3/3', hues = '10YR', style='CC')

# contrast class, multiple hues
# consider latticeExtra::useOuterStrips()
contrastChart(m = '10YR 5/6', hues = c('10YR', '2.5Y'), style='CC')

Soil Color Contrast

Description

Determine soil color contrast class according to methods outlined in the Soil Survey Manual. This function is typically called from colorContrast() which is simpler to use and provides more information.

Usage

contrastClass(v1, c1, v2, c2, dH, dV, dC, verbose = FALSE)

Arguments

Details

This function is fully vectorized but expects all inputs have the same length.

Value

A vector of color contrast classes (ordered factor). A list when verbose is TRUE.

Author(s)

D.E. Beaudette

References

• Soil Survey Technical Note 2 wayback machine URL

See Also

colorContrast

Examples

## standard use, result is an ordered factor
# 10YR 6/3 vs 5YR 3/4
contrastClass(v1=6, c1=3, v2=3, c2=4, dH=2, dV=3, dC=1)

## verbose output, useful for testing rules/cases
# 10YR 6/3 vs 5YR 3/4
contrastClass(v1=6, c1=3, v2=3, c2=4, dH=2, dV=3, dC=1, verbose = TRUE)

Apply rock fragment or salt correction to available water content

Description

Apply rock fragment or salt correction to available water content

Usage

correctAWC(
  awc,
  total_rf = numeric(length(awc)),
  gravel = NULL,
  ec = NULL,
  nullFragsAreZero = TRUE
)

Arguments

Value

A numeric vector (double) containing estimated available water capacities corrected for rock fragments and salts

Examples

# medium organic matter, loam texture 
base.awc <- 0.18 # estimateAWC(texcl = "l", omcl = 2, na.rm = TRUE)

# medium organic matter, loam texture w/ 23% rock fragments by volume 
corrected.awc <- correctAWC(base.awc, total_rf = 23)
corrected.awc

# medium organic matter, loam texture w/ 0% frags by volume and 8 mmhos/cm salts
salty.awc <- correctAWC(base.awc, total_rf = 0, ec = 8)
salty.awc

Determines threshold (minimum) clay content for argillic upper bound

Description

Given a vector or matrix of "eluvial" horizon clay contents (\ crit.clay.argillic() returns a vector or matrix of minimum clay contents (thresholds) that must be met for an argillic horizon clay increase.

Usage

crit.clay.argillic(eluvial_clay_content)

Arguments

Details

Uses the standard equations for clay contents less than 15 \ and 40 \ the definition of the argillic horizon from 12th Edition Keys to Soil Taxonomy (Soil Survey Staff, 2014).

Value

A vector or matrix (input-dependent) containing minimum "illuvial" horizon clay contents (thresholds) to be met for argillic horizon clay increase.

Note

This function is intended for identifying clay content threshold required for an argillic horizon. These thresholds may not apply depending on the specifics of your soil. E.g. if the upper part of argillic has been plowed (has Ap immediately over upper boundary) the clay increase requirement can be waived (Soil Survey Staff, 2014).

Author(s)

Andrew Gene Brown

References

Soil Survey Staff. 2014. Keys to Soil Taxonomy, 12th ed. USDA-Natural Resources Conservation Service, Washington, DC.

See Also

getArgillicBounds, get.increase.matrix

Examples

# crit.clay.argillic uses different equations for clay content
# less than 15 %, between 15 and 40 %, and >40 %

crit.clay.argillic(eluvial_clay_content=c(5, 20, 45))

Create a (redundant) horizon-level attribute from a site-level attribute

Description

Create a (redundant) horizon-level attribute from a site-level attribute. Specify a SoilProfileCollection and a site-level attribute from that SPC (by name) to receive a vector of length equal to the number of horizons containing the site-level values. This vector is directly usable with the SoilProfileCollection horizon setter.

denormalize is the inverse operation for the formula interface that "normalizes" a horizon level variable to site level:

site(object) <- ~ horizonvar

Usage

denormalize(object, attr)

Arguments

Details

"Denormalization" is the process of trying to improve the read performance of a database, at the expense of losing some write performance, by adding redundant copies of data or by grouping data. Sometimes it is beneficial to have site-level attributes denormalized for grouping of horizon-level data in analyses. denormalize achieves this result for SoilProfileCollections.

Value

A vector of values of equal length to the number of rows in the horizon table of the input SPC.

Author(s)

Andrew G. Brown, Dylan Beaudette

Examples

data(sp1)

# create a SoilProfileCollection from horizon data
depths(sp1) <- id ~ top + bottom

# create random site-level attribute `sitevar` with a binary (0/1) outcome
sp1$sitevar <- round(runif(length(sp1)))

# use denormalize() to create a mirror of sitevar in the horizon table
# name the attribute something different (e.g. `hz.sitevar`) to 
# prevent collision with the site attribute
# the attributes can have the same name but you will then need 
# site() or horizons() to access explicitly
sp1$hz.sitevar <- denormalize(sp1, 'sitevar')

# compare number of profiles to number of sitevar assignments
length(sp1)
table(sp1$sitevar)

# compare number of horizons to number of horizon-level copies of sitevar `hz.'sitevar`
nrow(sp1)
table(sp1$hz.sitevar)

Get top or bottom depths of horizons matching a regular expression pattern

Description

The depthOf family of functions calculate depth of occurrence of a horizon designation pattern, or any other value that can be coerced to character and matched with a regular expression.

If you need all depths of occurrence for a particular pattern, depthOf is what you are looking for. minDepthOf and maxDepthOf are wrappers around depthOf that return the minimum and maximum depth. They are all set up to handle missing values and missing "contacts" with the target pattern.

Usage

depthOf(
  p,
  pattern,
  FUN = NULL,
  top = TRUE,
  hzdesgn = hzdesgnname(p, required = TRUE),
  no.contact.depth = NULL,
  no.contact.assigned = NA_real_,
  na.rm = TRUE,
  simplify = TRUE
)

maxDepthOf(
  p,
  pattern,
  top = TRUE,
  hzdesgn = hzdesgnname(p, required = TRUE),
  no.contact.depth = NULL,
  no.contact.assigned = NA,
  na.rm = TRUE,
  simplify = TRUE
)

minDepthOf(
  p,
  pattern,
  top = TRUE,
  hzdesgn = hzdesgnname(p, required = TRUE),
  no.contact.depth = NULL,
  no.contact.assigned = NA,
  na.rm = TRUE,
  simplify = TRUE
)

Arguments

Value

a numeric vector containing specified depth(s) of horizons matching a pattern. If length(p) > 1 then a data.frame containing profile ID, horizon ID, top or bottom depths, horizon designation and pattern.

Author(s)

Andrew G. Brown

Examples

# construct a fake profile
spc <- data.frame(id=1, taxsubgrp = "Lithic Haploxerepts",
                  hzname   = c("A","AB","Bw","BC","R"),
                  hzdept   = c(0,  20, 32, 42,  49),
                  hzdepb   = c(20, 32, 42, 49, 200),
                  clay     = c(19, 22, 22, 21,  NA),
                  texcl    = c("l","l","l", "l","br"),
                  d_value  = c(5,   5,  5,  6,  NA),
                  m_value  = c(2.5, 3,  3,  4,  NA),
                  m_chroma = c(2,   3,  4,  4,  NA))

# promote to SoilProfileCollection
depths(spc) <- id ~ hzdept + hzdepb
hzdesgnname(spc) <- 'hzname'
hztexclname(spc) <- 'texcl'

# multiple horizons contain B
depthOf(spc, "B")

# deepest top depth of horizon containing B
maxDepthOf(spc, "B")

# shallowest top depth
minDepthOf(spc, "B")

# deepest bottom depth
maxDepthOf(spc, "B", top = FALSE)

# deepest bottom depth above 35cm
maxDepthOf(spc, "B", top = FALSE, no.contact.depth = 35)

# assign infinity (Inf) if B horizon does not start within 10cm
minDepthOf(spc, "B", no.contact.depth = 10, no.contact.assigned = Inf)

Return a vector of contributing fractions over a depth interval

Description

depthWeights() calculates the contributing fraction for each pair of horizon top and bottom depths, given an upper and lower boundary.

Usage

depthWeights(top, bottom, upper, lower)

Arguments

Value

A named list.

Author(s)

Andrew G. Brown.

Get depth units from metadata

Description

Get units of depth measurement from metadata. Default value is centimeters.

Usage

## S4 method for signature 'SoilProfileCollection'
depth_units(object)

## S4 replacement method for signature 'SoilProfileCollection'
depth_units(object) <- value

Arguments

Examples

data(sp5)

## get depth units
du <- depth_units(sp5)

# set alternate units; e.g. inches
depth_units(sp5) <- 'in'

# replace original value (cm)
depth_units(sp5) <- du

Initialize a SoilProfileCollection from data.frame

Description

⁠depths(<data.frame>) <- <formula>⁠: Initialize SoilProfileCollection

⁠depths(<SoilProfileCollection>)⁠: Extract profile ID and horizon depths from SoilProfileCollection

Usage

## S4 method for signature 'SoilProfileCollection'
depths(x, hzID = FALSE, ...)

## S4 replacement method for signature 'SoilProfileCollection'
depths(object) <- value

## S4 replacement method for signature 'data.frame'
depths(object) <- value

Arguments

Details

The input horizon data, and the resulting profile order, is sorted based on unique profile ID and top depth. ID columns are converted to character, depth columns are converted to integer. If NA values exist in all of the top depths, a prototype with 1 horizon per profile ID is returned, with NA in all non-essential columns. If the input object has 0 rows, a prototype with 0 horizons and 0 rows, but same column names as object, is returned.

Value

a data.frame containing profile ID, top depth, and bottom depth

See Also

horizons() idname() hzidname() horizonDepths()

Examples

# load a SoilProfileCollection
data(jacobs2000, package = "aqp")

depths(jacobs2000)
## init SoilProfileCollection objects from data.frame of horizon data

# load demo data
data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# plot
plot(sp1)

# number of profiles
length(sp1)

# number of horizons
nrow(sp1)

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.

• ⁠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.

Usage

## S4 method for signature 'SoilProfileCollection'
diagnostic_hz(object)

## S4 replacement method for signature 'SoilProfileCollection'
diagnostic_hz(object) <- value

Arguments

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)

Efficient Slicing of SoilProfileCollection Objects

Description

Cut ("dice") soil horizons into 1-unit thick slices. This function replaces aqp::slice(), which will be deprecated in aqp 2.0.

Usage

## S4 method for signature 'SoilProfileCollection'
dice(
  x,
  fm = NULL,
  SPC = TRUE,
  pctMissing = FALSE,
  fill = FALSE,
  strict = TRUE,
  byhz = TRUE,
  verbose = FALSE
)

Arguments

Details

For large and potentially messy collections that may include missing horizon depth logic errors, consider using repairMissingHzDepths() before dice(). Consider using accumulateDepths() before invoking dice() on collections that may contain old-style O horizon notation (e.g. 5-0cm).

Value

a SoilProfileCollection object, or data.frame when SPC = FALSE

Author(s)

D.E. Beaudette and A.G. Brown

See Also

repairMissingHzDepths(), accumulateDepths(), fillHzGaps()

Duplicate Profiles of a SoilProfileCollection

Description

A simple function to duplicate the contents of a SoilProfileCollection object. Old profile IDs are saved as a site-level attribute (oldID) and new IDs are generated using a numeric serial number.

Usage

duplicate(x, times = 3, oldID = ".oldID")

Arguments

Value

a SoilProfileCollection object

Author(s)

D.E. Beaudette

Examples

# sample data
data('sp4')

# promote to SPC
depths(sp4) <- id ~ top + bottom

# duplicate each profile 2 times
d <- duplicate(sp4, times = 2)

# graphical check
par(mar = c(0, 0, 3, 1))
plotSPC(d, color = 'Ca', width = 0.25)

Label placement based on a simulation of electrostatic forces

Description

This function attempts to move labels along a 1D coordinate system such that overlap (as specified by threshold) is minimized. An electrostatic simulation applies forces of repulsion between labels that are within thresh (e.g. overlapping) and forces of attraction to a uniformly spaced sequence to iteratively perturb affected labels until either no overlap is reported, or a maximum number of iterations (maxIter) has been reached.

Usage

electroStatics_1D(
  x,
  thresh,
  q = 1,
  chargeDecayRate = 0.01,
  QkA_GrowthRate = 0.05,
  maxIter = 100,
  tiny = 1e-04,
  const = 0.001,
  trace = FALSE,
  ...
)

Arguments

Details

Difficult overlap problems can be addressed by reducing thresh and increasing q. Large values of q can lead to chaotic results.

This function will generate unpredictable output when x contains duplicate values.

This function requires input to be pre-sorted, although interesting "artistic" simulations will often result from unsorted x.

Value

When trace = TRUE a list, otherwise numeric vector with converged attribute.

Author(s)

D.E. Beaudette and K.C. Thompson

See Also

fixOverlap(), SANN_1D()

Examples

# vector of object locations, with potential overlap
x <- c(1, 2, 3, 3.3, 3.8, 5, 6, 7, 8, 9, 10)

# full diagnostic output
z <- electroStatics_1D(x, thresh = 0.65, trace = TRUE, q = 1)
txt <- sprintf("Converged %s (%s iterations)", z$converged, length(z$cost))

plot(
seq_along(z$cost),
z$cost, 
las = 1, 
xlab = 'Iteration', 
ylab = 'Overlap Cost', 
type = 'b', 
main = txt
)

abline(h = 0, lty = 2, col = 2)

# final configuration only
xnew <- electroStatics_1D(x, thresh = 0.65, q = 1)

# check for convergence
attr(xnew, 'converged')

# compare original vs. modified
data.frame(orig = x, new = round(xnew, 2))

Identify "equivalent" (whole number value/chroma) Munsell chips

Description

Uses a pre-calculated lookup list (equivalent_munsell) based on pair-wise CIE2000 contrast (dE00) of LAB color with D65 illuminant for all whole value/chroma "chips" in the aqp::munsell data set.

The intention is to identify Munsell chips that may be "functionally equivalent" to some other given whole value/chroma chip elsewhere in the Munsell color space – as discretized in the aqp::munsell data table. This basic assumption needs to be validated against your end goal: probably by visual inspection of some or all of the resulting sets. See colorContrast and colorContrastPlot.

"Equivalent" chips table are based (fairly arbitrarily) on the 0.001 probability level of dE00 (default Type 7 quantile) within the upper triangle of the 8467x8467 contrast matrix. This corresponds to a dE00 contrast threshold of approximately 2.16.

Usage

equivalentMunsellChips(hue = NULL, value = NULL, chroma = NULL)

Arguments

Value

A named list; Each list element contains a data.frame with one or more rows of "equivalent" Munsell, RGB and LAB color coordinates from munsell data set.

References

Gaurav Sharma, Wencheng Wu, Edul N. Dalal. (2005). The CIEDE2000 Color-Difference Formula: Implementation Notes, Supplementary Test Data, and Mathematical Observations. COLOR research and application. 30(1):21-30. http://www2.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf

Thomas Lin Pedersen, Berendea Nicolae and Romain François (2020). farver: High Performance Colour Space Manipulation. R package version 2.0.3. https://CRAN.R-project.org/package=farver

Dong, C.E., Webb, J.B., Bottrell, M.C., Saginor, I., Lee, B.D. and Stern, L.A. (2020). Strengths, Limitations, and Recommendations for Instrumental Color Measurement in Forensic Soil Characterization. J Forensic Sci, 65: 438-449. https://doi.org/10.1111/1556-4029.14193

See Also

colorContrast colorContrastPlot equivalent_munsell

Examples

# 7.5YR 4/4 (the one and only)

equivalentMunsellChips("7.5YR", 4, 4)
#>
#> $`7.5YR 4/4`
#>        hue value chroma         r        g         b        L       A       B
#> 8330 7.5YR     4      4 0.4923909 0.352334 0.2313328 41.26403 10.8689 23.5914

# 7.5YR 1/1 (two chips are equivalent; 3 row result)

equivalentMunsellChips("7.5YR", 1, 1)
#>
#> $`7.5YR 1/1`
#>        hue value chroma         r         g          b        L        A        B
#> 1983  10YR     1      1 0.1345633 0.1087014 0.07606787 10.64787 1.621323 6.847629
#> 6189   5YR     1      1 0.1330994 0.1076359 0.09450179 10.63901 2.489012 3.515146
#> 8303 7.5YR     1      1 0.1329483 0.1082380 0.08862581 10.64210 2.065514 4.623922

# 10YR 6/8 (two chips are equivalent; 3 row result)

equivalentMunsellChips("10YR", 6, 8)
#>
#> $`10YR 6/8`
#>       hue value chroma         r         g         b        L        A        B
#> 2039 10YR     6      7 0.7382230 0.5512957 0.2680260 61.76795 10.50886 44.78574
#> 2040 10YR     6      8 0.7519872 0.5472116 0.2157209 61.77496 11.83215 51.15496
#> 2041 10YR     6      9 0.7642826 0.5433189 0.1559069 61.78085 13.09599 57.49773

# compare visually a very red color

veryred <- equivalentMunsellChips("10R", 6, 28)[[1]]

par(mar=c(0,0,1,1))

pie(rep(1, nrow(veryred)), col = with(veryred, munsell2rgb(hue, value, chroma)),
    label = with(veryred, sprintf("%s %s/%s", hue, value, chroma)))

table(veryred$hue) # 2 hues
#> 
#>  10R 7.5R 
#>    8   17

table(veryred$value) # 2 values
#> 
#>  5  6 
#> 11 14

table(veryred$chroma) # 10 chromas
#> 
#> 21 22 23 24 25 26 27 28 29 30 
#>  1  2  2  3  3  4  3  3  2  2

Indices of "equivalent" Munsell chips in the munsell data set

Description

A pre-calculated lookup list (made with farver::compare_colour) based on pair-wise color contrast (CIE2000 or dE00) evaluated over all "chips" in the aqp::munsell data set.

The intention is to identify Munsell chips that may be "functionally equivalent" to some other given whole chip elsewhere in the Munsell color space – as discretized in the aqp::munsell lookup table.

"Equivalent" chips are based (fairly arbitrarily) on the 0.001 probability level of dE00 (default Type 7 quantile) within the upper triangle of the 8467x8467 contrast matrix. This corresponds to a dE00 threshold of approximately 2.15.

This is a naive (to the subtleties of human color perception, and overall magnitude of contrast between some of the "chips") but computationally consistent approach. Using the lookup list, as opposed to manual contrast via e.g. farver::compare_colour may have some benefits for efficiency in certain applications where the exact contrast value is not as important as the concept of having some threshold that is non-zero, but very small.

Usage

data(equivalent_munsell)

Format

A named list with 8467 elements, each containing a numeric vector of indices corresponding to the munsell data set, which has 8467 rows (unique, whole-number chips). Names have the format HUE VALUE/CHROMA, e.g. "7.5YR 4/4"

References

Gaurav Sharma, Wencheng Wu, Edul N. Dalal. (2005). The CIEDE2000 Color-Difference Formula: Implementation Notes, Supplementary Test Data, and Mathematical Observations. COLOR research and application. 30(1):21-30. http://www2.ece.rochester.edu/~gsharma/ciede2000/ciede2000noteCRNA.pdf

Thomas Lin Pedersen, Berendea Nicolae and Romain Francois (2020). farver: High Performance Colour Space Manipulation. R package version 2.0.3. https://CRAN.R-project.org/package=farver

Dong, C.E., Webb, J.B., Bottrell, M.C., Saginor, I., Lee, B.D. and Stern, L.A. (2020). Strengths, Limitations, and Recommendations for Instrumental Color Measurement in Forensic Soil Characterization. J Forensic Sci, 65: 438-449. https://doi.org/10.1111/1556-4029.14193

See Also

equivalentMunsellChips

Examples

data(equivalent_munsell)

Estimate available water capacity for fine-earth fraction

Description

Estimate available water capacity for fine-earth fraction

Usage

estimateAWC(texcl, omcl, precision = 2, FUN = mean, ...)

Arguments

Value

A numeric vector double containing estimated available water capacities for fine-earth fraction.

Examples

# organic matter, loam texture, low medium and high OM
base.awc <- estimateAWC(c("l","l","l"), c(1, 2, 3), na.rm = TRUE)
base.awc

Estimate boundaries of the U.S Soil Taxonomy Particle Size Control Section

Description

Estimates the upper and lower boundary of the particle size control section for Mineral or Organic soilsby applying a programmatic version of the particle size control section key from the Keys to Soil Taxonomy (13th edition). See details for assumptions and required data elements.

Usage

estimatePSCS(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  clay.attr = hzmetaname(p, "clay", required = TRUE),
  texcl.attr = hztexclname(p, required = TRUE),
  lieutex = hzmetaname(p, "lieutex"),
  tax_order_field = "tax_order",
  bottom.pattern = "Cr|R|Cd|m",
  simplify = TRUE,
  ...
)

Arguments

Details

Requires information to identify argillic horizons (clay contents, horizon designations) with getArgillicBounds() as well as the presence of plow layers and surface organic soil material. Any getArgillicBounds() arguments may be passed to estimatePSCS.

Also, requires information on taxonomic order to handle Andisols.

In aqp 2.1, particle size control sections of organic soils Histosols and Histels are supported. This requires setting the "lieutex" horizon metadata column using ⁠hzmetaname<-()⁠ Horizon designations containing "f" or "W" are recognized as permafrost and water layers, respectively, for application of the organic soils key to control sections. In lieu textures "SPM" and "PEAT" are used to identify low density organic materials used for surface tier thickness. To avoid using the 160 cm surface tier, set the "lieutex" column to any column that does not contain "SPM" or "PEAT" values.

WARNING: Soils in arenic or grossarenic subgroups, with fragipans, or with strongly contrasting PSCs may not be classified correctly. The author would welcome a dataset to develop this functionality for.

Value

A numeric vector (when simplify=TRUE) containing the top and bottom depth of the particle size control section. First value is top, second value is bottom. If p contains more than one profile, the result is a data.frame with profile ID plus PSCS top and bottom depths.

Author(s)

Andrew Gene Brown

References

Soil Survey Staff. 2014. Keys to Soil Taxonomy, 12th ed. USDA-Natural Resources Conservation Service, Washington, DC.

See Also

getArgillicBounds(), getSurfaceHorizonDepth()

Examples

data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# set required metadata
hzdesgnname(sp1) <- 'name'
hztexclname(sp1) <- 'texture'
hzmetaname(sp1, 'clay') <- 'prop'

x <- estimatePSCS(sp1)
x

# change horizon texture and set inlieu texture column to turn
# first profile into an organic soil
sp1$name[1:6] <- c("Oi1", "Oi2", "Oi3", "Oaf", "Cf1", "Cf2")
sp1$texture <- as.character(sp1$texture)
sp1$texture[1:6] <- c("PEAT", "PEAT", "PEAT", "MUCK", "GRVLS", "GRVLS")
sp1$bottom[6] <- 200
hzmetaname(sp1, 'lieutex') <- 'texture'

y <- estimatePSCS(sp1[1, ], simplify = FALSE)

# 74cm lower boundary is 25cm past the upper boundary of permafrost (49cm)
# but minimum depth is 100cm unless there is a root-limiting layer
y

Estimate dry soil colors from moist soil colors and vice versa.

Description

Soil color is typically described at dry and moist conditions. This function attempts to estimate soil color at dry or moist condition when one is missing. Estimation proceeds as:

• convert Munsell notation to CIELAB color coordinates via munsell2rgb()

• apply scaling, rotation, and translation parameters in CIELAB color space

• locate closest Munsell chip to CIELAB coordinates via col2munsell()

Estimation of dry from moist soil color state is not guaranteed to be symmetric with estimation of moist from dry.

Usage

estimateSoilColor(hue, value, chroma, sourceMoistureState = c("dry", "moist"))

Arguments

Details

Scaling, rotation, and translation parameters for shifting between dry <–> moist CIELAB coordinates was determined using vegan::procrustes(), from those official series descriptions (OSD) where moist and dry soil colors were available.

Estimates for colors having a (dry or moist) Munsell value of 10 are not likely correct.

This is still a work in progress.

Value

data.frame of estimated colors in Munsell notation. The sigma column contains CIE2000 color contrast metric values describing the perceptual distance between estimated color in CIELAB coordinates and closest Munsell chip.

Author(s)

D.E. Beaudette

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

estimateSoilColor(hue = '10YR', value = 3, chroma = 3, sourceMoistureState = 'moist')

# note that estimation is not symmetric
estimateSoilColor(hue = '10YR', value = 5, chroma = 3, sourceMoistureState = 'dry')

Estimate Soil Depth

Description

Estimate the soil depth of a single profile within a SoilProfileCollection object. This function would typically be called by profileApply().

Usage

estimateSoilDepth(
  f,
  name = hzdesgnname(f),
  p = "Cr|R|Cd",
  selection = min,
  no.contact.depth = NULL,
  no.contact.assigned = NULL
)

Arguments

Details

The choice of a selection function usually follows:

min: the top of the first matching horizon, max: the top bot the last matching horizon, or possibly mean: somewhere in-between.

Value

single value representing the depth to contact or no.contact.assigned

Author(s)

D.E. Beaudette and J.M. Skovlin

See Also

getSoilDepthClass, profileApply

Examples

## consider a situation where there were multiple candidate
## "contacts": 2 Cd horizons over an R

# init hypothetical profile
d <- data.frame(
  id = '1',
  top = c(0, 10, 20, 30, 40, 50, 60),
  bottom = c(10, 20, 30, 40, 50, 60, 80),
  name = c('A', 'Bt1', 'Bt2', 'BC', 'Cd1', 'Cd2', 'R'),
  stringsAsFactors = FALSE
)

# upgrade to SPC
depths(d) <- id ~ top + bottom

# init horizon designation
hzdesgnname(d) <- 'name'

# visual check
par(mar = c(0, 0, 0, 2))
plotSPC(d, hz.depths = TRUE, name.style = 'center-center', cex.names = 1, width = 0.1)

# top of the first Cd
estimateSoilDepth(d, name = 'name')

# top of the first Cd
estimateSoilDepth(d, name = 'name', selection = min)

# top of the R
estimateSoilDepth(d, name = 'name', selection = max)

# top of the second Cd
estimateSoilDepth(d, name = 'name', selection = max, p = 'Cd')


## another example

data(sp1)
depths(sp1) <- id ~ top + bottom

# init horizon designation
hzdesgnname(d) <- 'name'

# apply to each profile in a collection, and save as site-level attribute
sp1$depth <- profileApply(sp1, estimateSoilDepth, name='name')

# this function can be used to "find" depth to any feature 
# that can be defined via REGEX pattern matching on the horizon name
# for example, locate the depth to the top "Bt" horizon
# returning NA when there is no match
sp1$top_Bt <- profileApply(
  sp1, estimateSoilDepth, 
  name='name', 
  p='Bt', 
  no.contact.depth=0, 
  no.contact.assigned=NA
)

# reduced margins
par(mar=c(1,1,1,2))
# adjust default y-offset and depth scaling for following examples
plotSPC(sp1, y.offset=10, scaling.factor=0.5)

# get plotting parameters for profile widths and depth scaling factors
lsp <- get("last_spc_plot", envir = aqp.env)

# positions on x-axis, same for both depth and top "Bt" horizon
x.positions <- (1:length(sp1)) - lsp$width

# annotate contact with unicode right-arrow
# y-position is adjusted based on plot y-offset and scaling factor
y.positions <- lsp$y.offset + (sp1$depth * lsp$scaling.factor)
text(x.positions, y.positions, '\u2192', col='red', adj=1, cex=1.25, lwd=2)

# annotate top "Bt" depth with unicode right-arrow
# y-position is adjusted based on plot y-offset and scaling factor
y.positions <- lsp$y.offset + (sp1$top_Bt * lsp$scaling.factor)
text(x.positions, y.positions, '\u2192', col='blue', adj=1, cex=1.25, lwd=2)


## Not run: 
  # sample data
  data(gopheridge, package='soilDB')
  
  # run on a single profile
  estimateSoilDepth(gopheridge[1, ], name = 'hzname')
  
  # apply to an entire collection
  profileApply(gopheridge, estimateSoilDepth, name = 'hzname')

## End(Not run)

Evaluate Generalized Horizon Labels

Description

Data-driven evaluation of generalized horizon labels using nMDS and silhouette width.

Usage

evalGenHZ(
  obj,
  genhz = GHL(obj, required = TRUE),
  vars,
  non.matching.code = "not-used",
  stand = TRUE,
  metric = "euclidean"
)

Arguments

Details

Classic multidimensional scaling is performed via stats::cmdscale(). The input distance matrix is generated by cluster::daisy() using (complete cases of) horizon-level attributes from obj as named in vars.

Silhouette widths are computed via cluster::silhouette(). The input distance matrix is generated by cluster::daisy() using (complete cases of) horizon-level attributes from obj as named in vars. Note that observations with genhz labels specified in non.matching.code are removed filtered before calculation of the distance matrix.

Value

a list is returned containing:

• horizons: ⁠c('mds.1', mds.2', 'sil.width', 'neighbor')⁠

• stats: mean and standard deviation vars, computed by generalized horizon label

• dist: the distance matrix as passed to stats::cmdscale()

Author(s)

D.E. Beaudette

See Also

get.ml.hz()

Evaluate Missing Data within a SoilProfileCollection

Description

Evaluate missing data within a SoilProfileCollection object

Data completeness is evaluated by profile or by horizon. Profile-level evaluation is based on the thickness of horizons (method = absolute) with complete horizon-level attributes (vars), optionally divided by the total thickness (method = relative). The REGEX pattern (p) is used to filter non-soil horizons from the calculation.

Usage

evalMissingData(
  x,
  vars,
  name = hzdesgnname(x),
  p = "Cr|R|Cd",
  method = c("relative", "absolute", "horizon")
)

Arguments

Value

A vector values ranging from 0 to 1 (method = 'relative') or 0 to maximum depth in specified depth units (method = 'absolute') representing the quantity of non-missing data (as specified in vars) for each profile. When method = 'horizon' a non-missing data fraction is returned for each horizon.

Author(s)

D.E. Beaudette

Examples

# example data
data("jacobs2000")

# fully populated
plotSPC(jacobs2000, name.style = 'center-center', 
        cex.names = 0.8, color = 'time_saturated')

# missing some data
plotSPC(jacobs2000, name.style = 'center-center', 
        cex.names = 0.8, color = 'concentration_color')

# very nearly complete
plotSPC(jacobs2000, name.style = 'center-center', 
        cex.names = 0.8, color = 'matrix_color')


# variables to consider
v <- c('time_saturated', 'concentration_color', 'matrix_color')

# compute data completeness by profile
# ignore 2C horizons
jacobs2000$data.complete <- evalMissingData(
  jacobs2000, 
  vars = v, 
  method = 'relative',
  p = '2C'
)

jacobs2000$data.complete.abs <- evalMissingData(
  jacobs2000, 
  vars = v, 
  method = 'absolute',
  p = '2C'
)

# compute data completeness by horizon
# ignore 2C horizons
jacobs2000$hz.data.complete <- evalMissingData(
  jacobs2000, 
  vars = v, 
  method = 'horizon',
  p = '2C'
)


# "fraction complete" by horizon
plotSPC(
  jacobs2000, name.style = 'center-center', 
  cex.names = 0.8, color = 'hz.data.complete'
)


# rank on profile completeness
new.order <- order(jacobs2000$data.complete)

# plot along data completeness ranking
plotSPC(
  jacobs2000, name.style = 'center-center', 
  cex.names = 0.8, color = 'concentration_color', 
  plot.order = new.order
)

# add relative completeness axis
# note re-ordering of axis labels
axis(
  side = 1, at = 1:length(jacobs2000), 
  labels = round(jacobs2000$data.complete[new.order], 2),
  line = 0, cex.axis = 0.75
)

# add absolute completeness (cm)
axis(
  side = 1, at = 1:length(jacobs2000), 
  labels = jacobs2000$data.complete.abs[new.order],
  line = 2.5, cex.axis=0.75
)

Visual Explanation for plotSPC

Description

Create a visual explanation for the many arguments to plotSPC. Call this function instead of plotSPC, all objects after x are passed on to plotSPC. Nearly all of the figures in the Introduction to SoilProfileCollection Objects tutorial are created with this function.

Usage

explainPlotSPC(x, ...)

Arguments

Value

a list of internally-used ordering vectors and graphical offsets / scaling factors

Author(s)

D.E. Beaudette

See Also

plotSPC

Examples

# sample data
data(sp4)
depths(sp4) <- id ~ top + bottom

# proposed vector of relative positions, overlap likely
pos <- c(1, 1.1, 3, 4, 5, 5.2, 7, 8, 9, 10)

# try it
explainPlotSPC(sp4, name = 'name', relative.pos=pos)

# attempt to fix using an integer sequence, short-circut will prevent adjustments
explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(1:10))

# attempt to adjust using defaults
explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos))

# attempt to adjust and tinker with defaults
explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos, adj = 0.2))

# enforce larger space between
explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos, thresh = 0.7))

# more complex adjustments required
pos <- c(1, 2, 3, 3.3, 5, 5.1, 5.5, 8, 9, 10)

# tinker
explainPlotSPC(sp4, name = 'name', relative.pos = pos)
explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos))

explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos, 
thresh = 0.7))

explainPlotSPC(sp4, name = 'name', relative.pos = fixOverlap(pos, 
thresh = 0.7, adj = 0.2))

# SANN: solution requires many iterations, and will not always converge
explainPlotSPC(sp4, name = 'name', 
relative.pos = fixOverlap(pos, thresh = 0.85, adj = 0.2)
)

# electrostatics: solution requires larger charge (q)
explainPlotSPC(sp4, name = 'name', 
relative.pos = fixOverlap(pos, thresh = 0.85, method = 'E', q = 2)
)

Find and Fill Horizon Gaps

Description

This function attempts to find "gaps" in the horizon records of a SoilProfileCollection object and fill with placeholder horizons (profile ID, horizon ID, to/bottom depths, all else NA). Missing horizon records between the top of each profile and to_top, or the bottom of each profile and to_bottom are treated as gaps when those arguments are not NULL. You can use this function to prepare a potentially messy SoilProfileCollection for subsequent analyses that are sensitive to horizon sequence inconsistencies or require a conformal "rectangle" of data spanning known depths.

Gaps are defined as:

• within each profile, for horizons i to n_hz:

• ⁠bottom_i != top_i+1 (but only to i = 1:(n_hz - 1)⁠

Usage

fillHzGaps(x, flag = TRUE, to_top = 0, to_bottom = max(x))

Arguments

Value

a possibly modified SoilProfileCollection object

Author(s)

A.G. Brown and D.E. Beaudette

Examples

data(sp4)
depths(sp4) <- id ~ top + bottom

# introduce depth logic errors
idx <- c(2, 6:7, 8, 12)
sp4$top[idx] <- NA

# check
horizons(sp4)[idx, ]

# create gaps by removing logic errors
x <- HzDepthLogicSubset(sp4, byhz = TRUE)

# check on removed horizons (hzID values)
metadata(x)$removed.horizons

# inspect
par(mar = c(0, 0, 0, 2))
plotSPC(x, width = 0.3, default.color = 'royalblue', 
name = 'hzID', name.style = 'center-center', cex.names = 0.8,
cex.id = 0.66)

# fill gaps left by HzDepthLogicSubset()
z <- fillHzGaps(x, flag = TRUE)

# graphical check
plotSPC(z, width = 0.3, color = '.filledGap', name = 'hzID', 
show.legend = FALSE, name.style = 'center-center', cex.names = 0.8,
cex.id = 0.66)

# fill top to 0 cm
z2 <- fillHzGaps(x, flag = TRUE, to_top = 0)
plotSPC(z2, width = 0.3, color = '.filledGap', name = 'hzID', show.legend = FALSE)

# fill bottom to max(SPC)
z3 <- fillHzGaps(x, flag = TRUE, to_top = 0, to_bottom = max(x))
plotSPC(z3, width = 0.3, color = '.filledGap', name = 'hzID', show.legend = FALSE)

## another example
data(sp4)
depths(sp4) <- id ~ top + bottom

# remove 1st horizons from profiles 1:4
idx <- sp4[,, .FIRST, .HZID]
replaceHorizons(sp4) <- horizons(sp4)[-idx[1:4], ]

# prepare for dice()
z <- fillHzGaps(sp4, to_top = 0, to_bottom = 50, flag = TRUE) 

# empty-horizon padding is in place for formula interface to dice()
d <- dice(z, fm = 0:50 ~ .)
plotSPC(d, color = 'Ca', show.legend = FALSE)
plotSPC(d, color = '.filledGap', show.legend = FALSE)

Find Overlap within a Sequence

Description

Establish which elements within a vector of horizontal positions overlap beyond a given threshold

Desc.

Usage

findOverlap(x, thresh)

overlapMetrics(x, thresh)

Arguments

Value

unique index to affected (overlapping) elements in x

Examples

x <- c(1, 2, 3, 3.4, 3.5, 5, 6, 10)

findOverlap(x, thresh = 0.5)


x <- c(1, 2, 3, 3.4, 3.5, 5, 6, 10)

overlapMetrics(x, thresh = 0.5)

Fix Overlap within a Sequence

Description

Fix Overlap within a Sequence

Usage

fixOverlap(x, thresh = 0.6, method = c("S", "E"), trace = FALSE, ...)

Arguments

Value

When trace = FALSE, a vector of the same length as x, preserving rank-ordering and boundary conditions. When trace = TRUE a list containing the new sequence along with information about objective functions and decisions made during adjustment of x.

See Also

electroStatics_1D(), SANN_1D()

Examples

s <- c(1, 2, 2.3, 4, 5, 5, 7)

# simulated annealing, solution is non-deterministic
fixOverlap(s, thresh = 0.6, method = 'S')

# electrostatics-inspired simulation of particles
# solution is deterministic
fixOverlap(s, thresh = 0.6, method = 'E')


# create a very busy profile with lots of possible overlapping
# depth annotation
x <- quickSPC(
  "SPC:AAA|BBB|CCC|D|EEEEE|FF|GG|HH|I|I|JJ|KK|LL|M|N|O|P|QQQQ|RR|S|TTTTTT|U", 
  interval = 1
)

# convert horizon ID to numeric
x$z <- as.numeric(x$hzID)

# plotSPC arguments
.a <- list(
  width = 0.2, 
  hz.depths = TRUE, 
  name.style = 'center-center', 
  cex.names = 1.5, 
  depth.axis = FALSE, 
  name = NA,
  color = 'z',
  show.legend = FALSE,
  print.id = FALSE,
  col.palette = hcl.colors(n = 25, palette = 'Spectral', rev = TRUE)
)

# set plotSPC default arguments
options(.aqp.plotSPC.args = .a)

# wrapper function to test label collision solutions
testIt <- function(x, ...) {
  
  plotSPC(x, ...)
  
  # a normalized index of label adjustment
  .txt <- sprintf(
    "LAI: %0.3f", 
    get('last_spc_plot', envir = aqp.env)$hz.depth.LAI
  )
  mtext(.txt, side = 1, at = 1, line = -2, cex = 0.8)
  
}


# compare and contrast
op <- par(mar = c(0, 0, 0, 0), mfcol = c(1, 6))

testIt(x)
title('ES (defaults)', line = -3)

testIt(x, fixOverlapArgs = list(method = 'S'))
title('SANN (defaults)', line = -3)

testIt(x, fixOverlapArgs = list(method = 'E', q = 1.5))
title('ES (q = 1.5)', line = -3)

testIt(x, fixOverlapArgs = list(method = 'E', q = 1))
title('ES (q = 1)', line = -3)

testIt(x, fixOverlapArgs = list(method = 'E', q = 0.5))
title('ES (q = 0.5)', line = -3)

testIt(x, fixOverlapArgs = list(method = 'E', q = 0.1))
title('ES (q = 0.1)', line = -3)

par(op)

Flag perfectly overlapping horizons within a SoilProfileCollection

Description

Flag perfectly overlapping horizons within a SoilProfileCollection

Usage

flagOverlappingHz(x)

Arguments

Details

Horizons with NA depths can be flagged as overlapping. Consider finding these horizons with checkHzDepthLogic(byhz=TRUE) and removing or fixing them.

Value

logical vector with length (and order) matching the horizons of x

Author(s)

D.E. Beaudette, A.G. Brown

See Also

checkHzDepthLogic() fillHzGaps()

Examples

# two overlapping horizons
z <- data.frame(
  id = 'SPC',
  top = c(0, 25, 25, 50, 75, 100, 100),
  bottom = c(25, 50, 50, 75, 100, 125, 125)
)

# init SPC
depths(z) <- id ~ top + bottom

# flag perfectly overlapping horizons
z$.overlapFlag <- flagOverlappingHz(z)

# thematic sketches
plotSPC(z, color = '.overlapFlag', hz.depths = TRUE, 
depth.axis = FALSE, cex.names = 0.85)

Coarse Fragment Class Labels and Diameter

Description

This is a convenience function for accessing coarse fragment class labels and associated diameter (mm), as defined in various classification systems such as USDA, Unified, and AASHTO.

Usage

fragmentClasses(
  sys = c("usda_simplified", "usda", "international", "unified", "aashto",
    "mod.wentworth"),
  flat = FALSE,
  rounded = FALSE
)

Arguments

Value

named vector of fragment diameter in mm

References

Schoeneberger, P.J., D.A. Wysocki, E.C. Benham, and Soil Survey Staff. 2012. Field book for describing and sampling soils, Version 3.0. Natural Resources Conservation Service, National Soil Survey Center, Lincoln, NE.

See Also

fragmentSieve()

Examples

# use default system: "usda_simplified"
fragmentClasses()
fragmentClasses(flat = TRUE)

fragmentClasses(sys = 'usda')
fragmentClasses(sys = 'USDA', flat = TRUE)

fragmentClasses(sys = 'international')

fragmentClasses(sys = 'unified')

fragmentClasses(sys = 'aashto')
fragmentClasses(sys = 'aashto', rounded = TRUE)

fragmentClasses(sys = 'mod.wentworth')

Sieve the Coarse Fraction of Soil

Description

Sieve applies thresholds to a numeric vector of fragment diameter values, returning fragment size classes. Particle diameter thresholds are evaluated as d < threshold.

Usage

fragmentSieve(
  diameter,
  sieves = NULL,
  ordered = FALSE,
  prefix = "",
  new_names = NULL,
  ...
)

Arguments

Value

character. Size class labels based on names of sieves, new_names, and prefix (if specified).

References

Soil Science Division Staff. 2017. Soil survey manual. C. Ditzler, K. Scheffe, and H.C. Monger (eds.). USDA Handbook 18. Government Printing Office, Washington, D.C.

See Also

fragmentClasses()

Examples

# use a simplified version of the USDA system
# common within NRCS/SPSD and NCSS
fragmentSieve(c(30, 125, 180, 500, 1000))

# pararock fragments
fragmentSieve(c(30, 125, 180, 500, 1000), prefix = 'para')

# result as an ordered factor
fragmentSieve(c(30, 125, 180, 500, 1000), ordered = TRUE)

# USDA system, flat size classes
fragmentSieve(c(30, 125, 180, 500, 1000), flat = TRUE)

# alternative classification systems
fragmentSieve(c(30, 125, 180, 500, 1000), sys = 'usda')
fragmentSieve(c(30, 125, 180, 500, 1000), sys = 'international')
fragmentSieve(c(30, 125, 180, 500, 1000), sys = 'unified')
fragmentSieve(c(30, 125, 180, 500, 1000), sys = 'aashto')
fragmentSieve(c(30, 125, 180, 500, 1000), sys = 'mod.wentworth')

# custom fragment labels / diameter
fragmentSieve(
  c(30, 125, 180, 500, 1000),
  sieves = c(clumps = 50, chunks = 300, blocks = 100000)
)

# unnamed sieves, generic labels used
fragmentSieve(c(10, 50), sieves = c(30, 70))

fragmentSieve(c(10, 50), sieves = c(30, 70), ordered = TRUE)
 

Generate Labels for Slabs

Description

This method is used by slab() for generating labels that assign IDs to layers in a SoilProfileCollection

Usage

genSlabLabels(
  slab.structure = 1,
  max.d = NULL,
  n.profiles = NULL,
  spc = NULL,
  diced = NULL,
  ...
)

Arguments

Details

The new routine used in aqp 2.0 requires that, at a minimum, the spc and slab.structure arguments be specified.

Value

factor. slab IDs, labels are the segment top and bottom depth separated by "-"

See Also

slab()

Generalize Horizon Names

Description

Generalize a vector of horizon names, based on new classes, and REGEX patterns. Or create a new column ghl in a SoilProfileCollection (requires a horizon designation name to be defined for the collection, see details)

Usage

generalize.hz(
  x,
  new,
  pattern,
  non.matching.code = "not-used",
  hzdepm = NULL,
  ordered = !missing(hzdepm),
  na.rm = TRUE,
  ...
)

## S4 method for signature 'character'
generalizeHz(
  x,
  new,
  pattern,
  non.matching.code = "not-used",
  hzdepm = NULL,
  ordered = !missing(hzdepm),
  ...
)

## S4 method for signature 'SoilProfileCollection'
generalizeHz(
  x,
  new,
  pattern,
  non.matching.code = "not-used",
  hzdepm = NULL,
  ordered = !missing(hzdepm),
  ghl = "genhz",
  ...
)

Arguments

Details

When x is a SoilProfileCollection the ghl column will be updated with the factor results. This requires that the "horizon designation name" metadata be defined for the collection to set the column for input designations.

Value

factor (an ordered factor when ordered=TRUE) of the same length as x (if character) or as number of horizons in x (if SoilProfileCollection)

Author(s)

D.E. Beaudette

References

Beaudette, D.E., Roudier, P., Skovlin, J. (2016). Probabilistic Representation of Genetic Soil Horizons. In: Hartemink, A., Minasny, B. (eds) Digital Soil Morphometrics. Progress in Soil Science. Springer, Cham. https://doi.org/10.1007/978-3-319-28295-4_18

See Also

hzdesgnname()

Examples

data(sp1)

# check original distribution of hz designations
table(sp1$name)

# generalized horizon labels
# character vector input
sp1$genhz <- generalizeHz(
  sp1$name,
  new = c('O','A','B','C','R'),
  pattern = c('O', '^A','^B','C','R'),
  ordered = TRUE
)

# see how we did / what we missed
table(sp1$genhz, sp1$name)


## a more advanced example, requries `perl = TRUE`
# example data
x <- c('A', 'AC', 'Bt1', '^AC', 'C', 'BC', 'CB')

# new labels
n <- c('A', '^AC', 'C')

# patterns:
# "A anywhere in the name"
# "literal '^A' anywhere in the name"
# "C anywhere in name, but without preceding A"
p <- c('A', '^A', '(?<!A)C')

# note additional argument
res <- generalizeHz(
  x, 
  new = n, 
  pattern = p, 
  perl = TRUE
)

# double-check: OK
table(res, x)

## apply to a SoilProfileCollection
data(sp1)
depths(sp1) <- id ~ top + bottom

# must set horizon designation metadata
hzdesgnname(sp1) <- 'name'

# result is a SoilProfileCollection
x <- generalizeHz(
  sp1,
  new = c('O','A','B','C','R'),
  pattern = c('O', '^A','^B','C','R'),
  ordered = TRUE
)

# GHL stored in 'genhz' column
x$genhz

# GHL metadata is set
GHL(x)

Convert cross-tabulation to adjacency matrix.

Description

Convert a cross-tabulation of e.g. original horizon designations vs. generalized horizon labels to adjacency matrix form.

Usage

genhzTableToAdjMat(tab)

Arguments

Value

matrix of numeric values

Author(s)

D.E. Beaudette

Compute Pair-wise Distances of Soil Properties over Depth

Description

Computes pair-wise distance matrix to determine where an attribute increases within a specified vertical distance threshold.

get.increase.depths performs the conversion of the square matrix output of get.increase.matrix back to horizon top depth for where criteria were met.

Usage

get.increase.matrix(p, attr, threshold.fun, vertical.distance)

get.increase.depths(p, attr, threshold.fun, vertical.distance)

Arguments

Details

Uses matrix outer product to determine all pair-wise differences in attr for the horizons of p. Supplies attr to threshold.fun to determine the minimum value criterion to return TRUE in output matrix for an "increase". Also, computes all pair-wise distances in depth dimension to determine whether the vertical distance criteria have been met simultaneously with attr increase.

This function assumes that the threshold.fun supplied by the user returns either a constant or a vector of equal length to its input.

Note that the threshold.fun result is allowed to contain NA, but that will result in no output for affected cells.

get.increase.depths() performs the conversion of the square matrix output of get.increase.matrix back to horizon top depth for where criteria were met.

Note that the threshold.fun result is allowed to contain NA, but that will result in no output for affected cells.

Value

Returns a square logical matrix reflecting where the increase criteria were met.

get.increase.depths converts to horizon top depth by using above matrix output to determine depths where increase is met.

Returns a numeric vector of depths where the increase requirement is met. For the argillic, the first is the one of interest.

get.increase.depths() converts to horizon top depth by using above matrix output to determine depths where increase is met.

Author(s)

Andrew Gene Brown

See Also

getArgillicBounds(), crit.clay.argillic()

getArgillicBounds() crit.clay.argillic()

Examples

data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

p <- sp1[1]
attr <- 'prop' # clay contents
foo <- get.increase.matrix(p, threshold.fun = crit.clay.argillic,
                           attr = attr, vertical.distance = 30)
foo


data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

p <- sp1[1]
attr <- 'prop' # clay contents
foo <- get.increase.depths(p, threshold.fun = crit.clay.argillic,
                           attr = attr, vertical.distance = 30)
foo

Determine ML Horizon Boundaries

Description

This function accepts input from slab() (a data.frame) along with a vector of horizon names, and returns a data.frame of the most likely horizon boundaries.

This function expects that x is a data.frame generated by slab(). If x was not generated by slab(), then o.names is required.

Usage

get.ml.hz(x, o.names = attr(x, which = "original.levels"))

Arguments

Value

A data.frame with the following columns:

• hz: horizon names

• top: horizon top depth

• bottom: horizon bottom depth

• confidence: integrated probability over thickness of each ML horizon, rounded to the nearest integer

• pseudo.brier: A "pseudo"" Brier Score for a multi-class prediction, where the most-likely horizon label is treated as the "correct" outcome. Details on the calculation for traditional Brier Scores here: https://en.wikipedia.org/wiki/Brier_score. Lower values suggest better agreement between ML horizon label and class-wise probabilities.

• mean.H: mean Shannon entropy (bits), derived from probabilities within each most-likely horizon. Larger values suggest more confusion within each ML.

Author(s)

D.E. Beaudette

References

Beaudette, D. E., Roudier, P., & Skovlin, J. (2016). Probabilistic representation of genetic soil horizons. Digital soil morphometrics, 281-293.

See Also

slab()

Examples

# init SPC
data(sp1)
depths(sp1) <- id ~ top + bottom

# set horizon designation metadata
hzdesgnname(sp1) <- 'name'

# generalize horizon designations from character vector
# result is an ordered factor
sp1$genhz <- generalizeHz(
  sp1$name,
  new = c('O','A','B','C'),
  pat = c('O', '^A','^B','C'),
  ordered = TRUE
)

# compute slice-wise GHL probability
# so that it sums to contributing fraction
# from 0-150cm
a <- slab(sp1, fm = ~ genhz, cpm = 1, slab.structure = 0:150)

# note original GHL names are set by slab()
attr(a, 'original.levels')

# generate table of ML horizonation
get.ml.hz(a)

Estimate upper and lower boundary of argillic diagnostic subsurface horizon

Description

getArgillicBounds estimates the upper and lower boundary of argillic diagnostic subsurface horizon for a profile in a single-profile SoilProfileCollection object (p).

The upper boundary is where the clay increase threshold is met. The function uses crit.clay.argillic as the threshold function for determining whether a clay increase occurs and get.increase.matrix to determine whether the increase is met, whether vertical distance of increase is sufficiently small, and in which horizon.

Usage

getArgillicBounds(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  clay.attr = hzmetaname(p, "clay", required = TRUE),
  texcl.attr = hztexclname(p, required = TRUE),
  require_t = TRUE,
  bottom.pattern = "Cr|R|Cd",
  lower.grad.pattern = "^[2-9]*B*CB*[^rtd]*[1-9]*$",
  sandy.texture.pattern = "-S$|^S$|COS$|L[^V]FS$|[^L]VFS$|LS$|LFS$",
  vertical.distance = 30,
  simplify = TRUE,
  verbose = FALSE
)

Arguments

Details

The lower boundary is first approximated as the depth to a lithic/paralithic/densic contact, or some other horizon matchable by a custom regular expression pattern. Subsequently, that boundary is extended upwards to the end of "evidence of illuviation."

The depth to contact is estimated using 'bottom.pattern' "Cr|R|Cd" by default. It matches anything containing Cr, R or Cd.

The lower gradational horizon regular expression ‘lower.grad.pattern' default is ⁠^[2-9]*B*CB*[^rtd]*[1-9]*$}⁠. It matches anything that starts with a lithologic discontinuity (or none) and a C master horizon designation. May contain B as second horizon designation in transitional horizon. May not contain 'r' or 't' subscript.

The minimum thickness of the argillic horizon is dependent on whether all subhorizons are "sandy" or not. The sandy.texture.pattern default ⁠-S$|^S$|COS$|L[^V]FS$|[^L]VFS$|LS$|LFS$⁠ captures USDA textural class fine earth fractions that meet "sandy" particle size class criteria.

There also is an option ‘require_t' to omit the requirement for evidence of eluviation in form of 't' subscript in 'hzdesgn'. Even if "t" subscript is not required for positive identification, the presence of lower gradational C horizons lacking 't' will still be used to modify the lower boundary upward from a detected contact, if needed. If this behavior is not desired, just set 'lower.grad.pattern' to something that will not match any horizons in your data.

Value

Returns a numeric vector; first value is top depth, second value is bottom depth. If as.list is TRUE, returns a list with top depth named "ubound" and bottom depth named "lbound". If p has more than one profile or if simplify = FALSE the result is a data.frame containing profile ID, upper and lower boundary columns.

Author(s)

Andrew G. Brown

Examples

data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# set required metadata
hzdesgnname(sp1) <- 'name'
hztexclname(sp1) <- 'texture'
hzmetaname(sp1, 'clay') <- 'prop'

x <- getArgillicBounds(sp1)
x

Find all intervals that are potentially part of a Cambic horizon

Description

Find all intervals that are potentially part of a Cambic horizon excluding those that are part of an argillic horizon (defined either by depth interval or getArgillicBounds()).

There may be multiple cambic horizons (indexes) in a profile. Each cambic index has a top and bottom depth associated: cambic_top and cambic_bottom. This result is designed to be used for single profiles, or with profileApply(..., frameify = TRUE)

Usage

getCambicBounds(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  texcl.attr = hztexclname(p, required = TRUE),
  clay.attr = hzmetaname(p, "clay", required = TRUE),
  argi_bounds = NULL,
  d_value = "d_value",
  m_value = "m_value",
  m_chroma = "m_chroma",
  sandy.texture.pattern = "-S$|^S$|COS$|L[^V]FS$|[^L]VFS$|LS$|LFS$",
  ...
)

Arguments

Value

A data.frame containing profile, cambic indexes, along with top and bottom depths.

Author(s)

Andrew G. Brown

Examples

# construct a fake profile
spc <- data.frame(id=1, taxsubgrp = "Lithic Haploxerepts",
                  hzname   = c("A","AB","Bw","BC","R"),
                  hzdept   = c(0,  20, 32, 42,  49),
                  hzdepb   = c(20, 32, 42, 49, 200),
                  clay     = c(19, 22, 22, 21,  NA),
                  texcl    = c("l","l","l", "l","br"),
                  d_value  = c(5,   5,  5,  6,  NA),
                  m_value  = c(2.5, 3,  3,  4,  NA),
                  m_chroma = c(2,   3,  4,  4,  NA))

# promote to SoilProfileCollection
depths(spc) <- id ~ hzdept + hzdepb

# set required metadata
hzdesgnname(spc) <- 'hzname'
hztexclname(spc) <- 'texcl'
hzmetaname(spc, 'clay') <- 'clay'

# print results in table
getCambicBounds(spc)

Get Approximate Munsell Chip

Description

Non-standard Munsell notation ('7.9YR 2.7/2.0') can be matched (nearest-neighbor, no interpolation) to the closest color within the munsell sRGB/CIELAB look-up table via getClosestMunsellChip(). A more accurate estimate of sRGB values from non-standard notation can be achieved with the munsellinterpol package. For example, conversion from Munsell to CIELAB, assuming a D65 illuminant via: MunsellToLab('0.1Y 3.3/4.4', white='D65', adapt='Bradford').

Usage

getClosestMunsellChip(munsellColor, convertColors = TRUE, ...)

Arguments

Value

a data.frame when convertColors=TRUE, otherwise character vector

Examples

# convert a non-standard color to closest "chip" in `munsell` look-up table
getClosestMunsellChip('7.9YR 2.7/2.0', convertColors = FALSE)

# convert directly to R color
getClosestMunsellChip('7.9YR 2.7/2.0')

# special case for 2.5 value -> no rounding, we have these records in the conversion LUT
getClosestMunsellChip('7.5YR 2.5/2', convertColors = FALSE)


getClosestMunsellChip('7.5YR 6.8/4.4', convertColors = FALSE)

Get IDs of Deepest Horizons by Profile

Description

Return horizon IDs of the deepest horizon within each profile of a SoilProfileCollection. IDs are returned in the same order as profile_id(x). Horizon top depths are used because there are cases where bottom depths may be missing.

Usage

getLastHorizonID(x)

Arguments

Generate Soil Depth Class Matrix

Description

Generate a boolean matrix of soil depth classes, actual soil depth class, and estimate of soil depth from a SoilProfileCollection object. Soil depths are estimated using pattern matching applied to horizon designations, by estimateSoilDepth(). The default REGEX pattern (p = 'Cr|R|Cd') will match most "contacts" described using the USDA / Soil Taxonomy horizon designation conventions.

Usage

getSoilDepthClass(
  f,
  depth.classes = c(very.shallow = 25, shallow = 50, mod.deep = 100, deep = 150,
    very.deep = 10000),
  ...
)

Arguments

Value

a data.frame containing soil depth and depth class for each profile, see examples

Author(s)

D.E. Beaudette and J.M. Skovlin

See Also

estimateSoilDepth

Examples

data(sp1)
depths(sp1) <- id ~ top + bottom

# generate depth-class matrix
sdc <- getSoilDepthClass(sp1, name = 'name')

# inspect
head(sdc)

# join back into sp1 as site-level data
site(sp1) <- sdc

## Not run: 
# sample data
data(gopheridge, package='soilDB')

getSoilDepthClass(gopheridge, name = 'hzname')

## End(Not run)

Determine thickness of horizons (continuous from surface) matching a pattern

Description

Find the thickness of horizon designations, or any other character patterns, that are continuous from the soil surface (depth = 0 or shallowest depth in profile).

Usage

getSurfaceHorizonDepth(
  p,
  pattern,
  hzdesgn = hzdesgnname(p, required = TRUE),
  simplify = TRUE
)

getMineralSoilSurfaceDepth(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  pattern = "O",
  simplify = TRUE
)

getPlowLayerDepth(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  pattern = "^Ap[^b]*",
  simplify = TRUE
)

Arguments

Details

The horizon designation to match is specified with the regular expression pattern 'pattern'. All horizons matching that pattern, that are continuous from the soil surface, count towards the depth / thickness value that is ultimately returned. For instance: horizon designations: A1-A2-A3-C-Ab , would return A3 bottom depth given pattern = "^A[1-9]*$".

getSurfaceHorizonDepth is used by getPlowLayerDepth for matching Ap horizons; and, it is used by getMineralSoilSurfaceDepth to find the thickness of O horizons in lieu of lab data.

Value

a numeric value corresponding to the bottom depth of the last horizon matching 'pattern' that is contiguous with other matching horizons up to the soil surface. If length(p) > 1 then a data.frame containing profile ID, horizon ID, top or bottom depths, horizon designation and pattern.

Author(s)

Andrew G. Brown

Examples

library(aqp)
data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

p <- sp1[1]
q <- sp1[2]

# look at horizon designations in p and q
p$name
q$name

# thickness of all surface horizons containing A
getSurfaceHorizonDepth(p, pattern = 'A', hzdesgn = 'name')

# thickness of all surface horizons that start with A
getSurfaceHorizonDepth(p, pattern = '^A', hzdesgn = 'name')

# thickness of all surface horizons that start with A, and the A is not followed by B
getSurfaceHorizonDepth(p, pattern = '^A[^B]*', hzdesgn = 'name')

# thickness of all surface horizons that start with A
#  followed by a number from _2_ to 9 (returns ZERO)
getSurfaceHorizonDepth(p, pattern = '^A[2-9]*', hzdesgn = 'name')

# getPlowLayerDepth matches first two horizons in fake Ap horizon data with "buried Ap"
p$aphorizons <- c("Ap1","Ap2","AB", rep('C', nrow(p) - 4), "Apb")
getPlowLayerDepth(p, hzdesgn = 'aphorizons')

# getMineralSoilSurfaceDepthmatches first 3 horizons in fake O horizon data
p$ohorizons <- c("Oi1","Oi2","Oe", rep('C', nrow(p) - 4), "2C")
getMineralSoilSurfaceDepth(p, hzdesgn='ohorizons')

# matches first Oi horizon with original horizon designations of pedon 2
getMineralSoilSurfaceDepth(q, hzdesgn='name')

Subset soil horizon data using a depth or depth interval

Description

Make a "clod" of horizons from a SoilProfileCollection given a point or a depth interval to intersect. The interval ⁠[z1,z2]⁠ may be profile-specific (equal in length to p), or may be recycled over all profiles (if boundaries are length 1). For "point" intersection, z2 may be left as the default value NULL.

trunc() is a wrapper method (using S4 generic) for glom() where truncate=TRUE

Usage

## S4 method for signature 'SoilProfileCollection'
glom(
  p,
  z1,
  z2 = NULL,
  ids = FALSE,
  df = FALSE,
  truncate = FALSE,
  invert = FALSE,
  fill = FALSE,
  modality = "all",
  drop = !fill,
  ...
)

## S4 method for signature 'SoilProfileCollection'
trunc(x, z1, z2, ...)

Arguments

Details

"To glom" is "to steal" or to "become stuck or attached to". The word is related to the compound "glomalin", which is a glycoprotein produced by mycorrhizal fungi in soil.

The full depth range of horizons included within the interval are returned (a "ragged" SoilProfileCollection) unless the truncate argument is set as TRUE. Horizon intersection is based on unique ID hzidname(spc) and depth range of interest. Profiles that lack data in the range of interest will be dropped from the resulting SoilProfileCollection.

If inverting results with invert, it is possible that thick horizons (whose boundaries span wider than the specified interval) will be split into two horizons, where previously they were one. This may make the results from ids = TRUE different from what you expect, as they will be based on a profile with an "extra" horizon and re-calculated unique horizon ID (hzidname(spc)) "hzID".

Value

a SoilProfileCollection, data.frame, or a vector of horizon IDs. NULL if no result.

Author(s)

Andrew G. Brown

See Also

glomApply trunc

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

data(sp1, package = 'aqp')
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

p <- glom(sp1, 25, 150)

# 28 horizons
nrow(p) 

# inspect graphically
par(mar = c(1,1,3,1))
plot(p, color = "prop", max.depth = 200)
abline(h = c(25, 100), lty = 2)

## glom(..., truncate = TRUE)

p2 <- glom(sp1, 25, 150, truncate = TRUE)

# 28 horizons
nrow(p2) 

# inspect graphically
par(mar = c(1,1,3,1))
plot(p2, color = "prop", max.depth = 200)
abline(h = c(25, 100), lty = 2)

## glom(..., truncate = TRUE, invert = TRUE)

p3 <- glom(sp1, 25, 150, truncate = TRUE, invert = TRUE)

# 45 horizons
nrow(p3) 

# inspect graphically
par(mar = c(1,1,3,1))
plot(p3, color = "prop", max.depth = 200)
abline(h = c(25, 100), lty = 2)

## profile-specific interval, using expressions evaluated within sp1@site

# calculate some new site-level variables containing target interval
sp1$glom_top <- (1:9) * 10
sp1$glom_bottom <- 10 + sp1$glom_top

# glom evaluates non-standard expressions using siteNames(sp1) column names
p4 <- glom(sp1, glom_top / 2, glom_bottom * 1.2, truncate = TRUE)

# inspect graphically
par(mar = c(1,1,3,1))
plot(p4, color = "prop", max.depth = 200)


# load sample data
data("sp3")

# promote to SPC
depths(sp3) <- id ~ top + bottom

### TRUNCATE all profiles in sp3 to [0,25]

# set up plot parameters
par(mfrow=c(2,1), mar=c(0,0,0,0))

# full profiles
plot(sp3)

# trunc'd profiles
plot(trunc(sp3, 0, 25))

Subset an SPC by applying glom to each profile

Description

glomApply() is a function used for subsetting SoilProfileCollection objects by depth. It is a wrapper around glom which is intended to subset single-profile SPCs based on depth intervals/intersection.

glomApply works by accepting a function .fun as argument. This function is used on each profile to process a multi-profile SPC for input to glom (via profileApply). For each profile, .fun returns a 2-length numeric vector of top and bottom boundaries glom arguments: z1, z2.

glomApply provides the option to generate profile-specific glom depths for a large SPC and handles iteration and rebuilding of a subset SPC object. Optional arguments include: truncate to cut the boundaries to specified [z1, z2]; invert to the portion outside [z1, z2], modality to either "all" horizons or "thickest" horizon in the glom interval. ... are various expressions you can run on the individual profiles using NSE, similar to mutate.

Usage

glomApply(
  object,
  .fun = NULL,
  truncate = FALSE,
  invert = FALSE,
  modality = "all",
  ...,
  chunk.size = 100
)

Arguments

Value

A SoilProfileCollection.

Author(s)

Andrew G. Brown.

See Also

glom trunc

glom glomApply

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

data(sp3)
depths(sp3) <- id ~ top + bottom

# init horizon designation column in metadata, used by estimateSoilDepth
hzdesgnname(sp3) <- 'name'

# constant depths, whole horizon returns by default
plot(glomApply(sp3, function(p) c(25,100)))

# constant depths, truncated
#(see aqp::trunc for helper function)
plot(glomApply(sp3, function(p) c(25,30), truncate = TRUE))

# constant depths, inverted
plot(glomApply(sp3, function(p) c(25,100), invert = TRUE))

# constant depths, inverted + truncated (same as above)
plot(glomApply(sp3, function(p) c(25,30), invert = TRUE, truncate=TRUE))

# random boundaries in each profile
plot(glomApply(sp3, function(p) round(sort(runif(2, 0, max(sp3))))))

# random boundaries in each profile (truncated)
plot(glomApply(sp3, function(p) round(sort(runif(2, 0, max(sp3)))), truncate = TRUE))

# calculate some boundaries as site level attribtes
sp3$glom_top <- profileApply(sp3, getMineralSoilSurfaceDepth)
sp3$glom_bottom <- profileApply(sp3, estimateSoilDepth)

# use site level attributes for glom intervals for each profile
plot(glomApply(sp3, function(p) return(c(p$glom_top, p$glom_bottom))))

Subset SPC with pattern-matching for text-based attributes

Description

grepSPC() is a shorthand function for subsetting SoilProfileCollection objects. For example, by filter(grepl(spc, ...)) or filter(stringr::str_detect(spc, ...)). It provides pattern matching for a single text-based site or horizon level attribute.

Usage

grepSPC(object, attr, pattern, ...)

Arguments

Value

A SoilProfileCollection.

Author(s)

Andrew G. Brown.

Store groupings within a profile collection.

Description

Store groupings within a profile collection.

Usage

groupSPC(object, ...)

Arguments

Grouped Soil Profile Plot

Description

Plot a collection of soil profiles, sorted by group.

The left-right ordering of groups can be adjusted by converting groups into a factor and explicitly setting factor levels. Alpha-numeric ordering is used for all other types.

Usage

groupedProfilePlot(
  x,
  groups,
  group.name.offset = -5,
  group.name.cex = 0.75,
  group.line.col = "RoyalBlue",
  group.line.lwd = 2,
  group.line.lty = 2,
  break.style = c("line", "arrow", "both"),
  break.offset = 0.5,
  arrow.offset = group.name.offset + 5,
  arrow.length = 0.1,
  ...
)

Arguments

Author(s)

D.E. Beaudette

See Also

plotSPC

Examples

# sample data
data(sp1)
# convert colors from Munsell to hex-encoded RGB
sp1$soil_color <- with(sp1, munsell2rgb(hue, value, chroma))

# promote to SoilProfileCollection
depths(sp1) <- id ~ top + bottom
site(sp1) <- ~ group

# add a groups
sp1$group.2 <- sprintf("%s-%s", rev(LETTERS[1:3]), sp1$group)

# convert fake groupt to factor with new levels
sp1$group.3 <- factor(sp1$group.2, levels=c('C-2', 'B-2', 'A-2', 'C-1', 'B-1', 'A-1'))

# plot profiles, sorted and annotated by 'group' (integers)
par(mar=c(1,1,1,1))
groupedProfilePlot(sp1, groups='group', max.depth=150, group.name.offset = -5, id.style='side')

# plot profiles, sorted and annotated by 'group.2' (characters)
par(mar=c(1,1,1,1))
groupedProfilePlot(sp1, groups='group.2', max.depth=150, group.name.offset = -5, id.style='side')

# plot profiles, sorted and annotated by 'group.3' (characters)
par(mar=c(1,1,1,1))
groupedProfilePlot(sp1, groups='group.3', max.depth=150, group.name.offset = -5, id.style='side')


# make fake site-level attribute and adjust levels
sp1$new.group <- sample(letters[1:3], size=length(sp1), replace=TRUE)

# tabulate pedons / group
tab <- table(sp1$new.group)

# sort large -> small
tab <- sort(tab, decreasing = TRUE)

# set levels based on sorted tabulation
# assign custom labels
sp1$new.group <- factor(sp1$new.group, levels=names(tab),
labels=paste0(names(tab), ' (', tab, ')'))

groupedProfilePlot(sp1, groups='new.group', max.depth=150,
group.name.offset = -10, id.style='side')

# offsets can be set using a vector of values, recycled as needed
groupedProfilePlot(sp1, groups='new.group', max.depth=150,
group.name.offset=c(-10, -5), id.style='side')

# annotate with arrows instead of vertical lines
groupedProfilePlot(sp1, groups='new.group', max.depth=150,
group.name.offset=c(-10, -12), break.style='arrow', arrow.offset=-3,
group.line.lty = 1, group.line.lwd = 1, id.style='side')


## Not run: 
# more complete example using data from soilDB package
data(loafercreek, package='soilDB')
par(mar=c(1,1,1,1))
# lines
groupedProfilePlot(loafercreek, groups='hillslopeprof', group.name.cex = 0.5,
group.name.offset = -10)

# arrows
groupedProfilePlot(loafercreek, groups='hillslopeprof', group.name.cex = 0.5,
group.name.offset = -10, break.style ='arrow', group.line.lty = 1,
group.line.lwd = 1)

# both
groupedProfilePlot(loafercreek, groups='hillslopeprof', group.name.cex = 0.5,
group.name.offset = -10, break.style ='both', group.line.lty = 1,
group.line.lwd = 1)

## End(Not run)

Guess Appropriate Ordering for Generalized Horizon Labels

Description

This function makes an (educated) guess at an appropriate set of levels for generalized horizon labels using the median of horizon depth mid-points.

Usage

guessGenHzLevels(x, hz = GHL(x, required = TRUE))

Arguments

Details

This function is useful when groups of horizons have been generalized via some method other than generalize.hz. For example, it may be useful to generalize horizons using labels derived from slice depths. The default sorting of these labels will not follow a logical depth-wise sorting when converted to a factor. guessGenHzLevels does a good job of "guessing" the proper ordering of these labels based on median horizon depth mid-point.

Value

a list:

Author(s)

D.E. Beaudette

See Also

generalize.hz

Examples

# load some example data
data(sp1, package='aqp')

# upgrade to SoilProfileCollection
depths(sp1) <- id ~ top + bottom

# generalize horizon names
n <- c('O', 'A', 'B', 'C')
p <- c('O', 'A', 'B', 'C')
sp1$genhz <- generalize.hz(sp1$name, n, p)

# note: levels are in the order in which originally defined:
levels(sp1$genhz)

# generalize horizons by depth slice
s <- dice(sp1, c(5, 10, 15, 25, 50, 100, 150) ~ .)
s$slice <- paste0(s$top, ' cm')
# not a factor
levels(s$slice)

# the proper ordering of these new labels can be guessed from horizon depths
guessGenHzLevels(s, 'slice')

# convert to factor, and set proper order
s$slice <- factor(s$slice, levels=guessGenHzLevels(s, 'slice')$levels)

# that is better
levels(s$slice)

Guess Horizon Slot Column Names

Description

guessHzAttrName(): Guess the horizon column name where possible/preferred formative elements are known. There is a preference for records where more optional requirements are met to handle cases where there will be many matches. For example, working with soil data one might have "low, RV and high" total clay, as well as clay fractions. One could distinguish between these different measurements using standard formative elements for column names from the database of interest. Result is the first match in horizonNames(x) with the most required plus optional patterns matched.

e.g. guessHzAttrName(x, attr="clay", optional=c("total", "_r")) matches (claytotal_r == totalclay_r) over (clay_r == claytotal == totalclay) over clay.

guessHzDesgnName(): DEPRECATED This follows the historic convention used by aqp::plotSPC() looking for "hzname" or other column names containing the regular expression "name". If the pattern "name" is not found, the pattern "desgn" is searched as a fallback, as "hzdesgn" or "hz_desgn" are other common column naming schemes for horizon designation name.

guessHzTexClName(): DEPRECATED This function is used to provide a texture class attribute column name to functions. It will use regular expressions to match "texcl" which is typically the texture of the fine earth fraction, without modifiers or in-lieu textures. Alternately, it will match "texture" for cases where "texcl" is absent (e.g. in NASIS Component Horizon).

Usage

guessHzAttrName(x, attr, optional = NULL, verbose = TRUE, required = FALSE)

guessHzDesgnName(x, required = FALSE)

guessHzTexClName(x, required = FALSE)

Arguments

Value

Character containing horizon attribute column name. Result is the first match in horizonNames(x) with the most required plus optional patterns matched.

Author(s)

Andrew G. Brown

Examples

# a has the required attr pattern, but none of the optional
a <- data.frame(id = 1, top = c(0,10), bottom=c(10,40),
                clay=c(18,19))
depths(a) <- id ~ top + bottom

guessHzAttrName(a, attr="clay", optional=c("total", "_r"))

# b has requried attr pattern, and one of the opional patterns
#   notice that it also contains "clay" but preferentially matches more optional patterns
b <- data.frame(id = 1, top = c(0,10), bottom=c(10,40),
                clay=c(0.18,0.19), clay_r=c(18,19))
depths(b) <- id ~ top + bottom

guessHzAttrName(b, attr="clay", optional=c("total", "_r"))

# c has total and _r (both optional) on either side of clay
# having all of the optional patterns plus required is best evidence, and first
# column containing that combination will be returned
c <- data.frame(id = 1, top = c(0,10), bottom=c(10,40),
                totalclay_r=c(18,19), claytotal_r=c(0.18,0.19))
depths(c) <- id ~ top + bottom

guessHzAttrName(c, attr="clay", optional=c("total", "_r"))

Harden (1982) Melanization

Description

Calculate "melanization" component of "Profile Development Index" after Harden (1982) "A quantitative index of soil development from field descriptions: Examples from a chronosequence in central California". Accepts vectorized inputs for value and reference value to produce vector output. A convenient use case would be to apply this on a profile-specific basis, where the value_ref has a single value, and value is a vector of length equal to the number of horizons within the upper 100 cm.

Usage

harden.melanization(value, value_ref)

Arguments

Details

In Harden (1982), "melanization" is calculated relative to a reference parent material for all horizons within 100cm of the soil surface. In addition, several other non-color components are normalized relative to a maximum value and summed to obtain the overall Profile Development Index.

Value

A numeric vector reflecting horizon darkening relative to a reference (e.g. parent) material.

Author(s)

Andrew G. Brown

References

Harden, J.W. (1982) A quantitative index of soil development from field descriptions: Examples from a chronosequence in central California. Geoderma. 28(1) 1-28. doi: 10.1016/0016-7061(82)90037-4

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

library(aqp)
data("jacobs2000", package="aqp")

# LEFT JOIN hue, value, chroma matrix color columns
horizons(jacobs2000) <- cbind(horizons(jacobs2000)[,c(idname(jacobs2000), hzidname(jacobs2000))],
                              parseMunsell(jacobs2000$matrix_color_munsell, convertColors = FALSE))

# calculate a mixed 150-200cm color ~"parent material"

jacobs2000$c_horizon_color <- profileApply(jacobs2000, function(p) {

  # and derive the parent material from the 150-200cm interval
  p150_200 <- glom(p, 150, 200, truncate = TRUE)
  p150_200$thickness <- p150_200$bottom - p150_200$top

  # mix colors
  clrs <- na.omit(horizons(p150_200)[,c('matrix_color_munsell','thickness')])
  mixMunsell(clrs$matrix_color_munsell, w = clrs$thickness)$munsell

})

# segment profile into 1cm slices (for proper depth weighting)
jacobs2000$melan <- profileApply(jacobs2000, function(p) {

  # sum the melanization index over the 0-100cm interval
  p0_100 <- hz_segment(p, 0:100)

  ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE)

  sum(harden.melanization(
    value = as.numeric(p0_100$value),
    value_ref = as.numeric(ccol$value)), na.rm = TRUE)

})

jacobs2000$melanorder <- order(jacobs2000$melan)

# Plot in order of increasing Melanization index

plotSPC(jacobs2000, 
        color = "matrix_color",
        label = "melan",
        plot.order = jacobs2000$melanorder,
        max.depth = 250
        )

segments(
  x0 = 0.5, 
  x1 = length(jacobs2000) + 0.5, 
  y0 = c(0,100,150,200), 
  y1 = c(0,100,150,200), 
  lty = 2
)

# Add [estimated] parent material color swatches
lapply(seq_along(jacobs2000$c_horizon_color), function(i) {
  rect(i - 0.15, 250, i + 0.15, 225,
       col = parseMunsell(jacobs2000$c_horizon_color[jacobs2000$melanorder[i]]))
})

Harden (1982) Rubification

Description

Calculate "rubification" component of "Profile Development Index" after Harden (1982) "A quantitative index of soil development from field descriptions: Examples from a chronosequence in central California". Accepts vectorized inputs for hue and chroma to produce vector output.

In Harden (1982) "rubification" is calculated relative to a reference parent material. Several other non-color components are normalized relative to a maximum value and summed to obtain the overall Profile Development Index.

Usage

harden.rubification(hue, chroma, hue_ref, chroma_ref)

Arguments

Value

A numeric vector reflecting horizon redness increase relative to a reference (e.g. parent) material.

Author(s)

Andrew G. Brown

References

Harden, J.W. (1982) A quantitative index of soil development from field descriptions: Examples from a chronosequence in central California. Geoderma. 28(1) 1-28. doi: 10.1016/0016-7061(82)90037-4

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

library(aqp)
data("jacobs2000", package="aqp")

# LEFT JOIN hue, value, chroma matrix color columns
horizons(jacobs2000) <- cbind(horizons(jacobs2000)[,c(idname(jacobs2000), hzidname(jacobs2000))],
                              parseMunsell(jacobs2000$matrix_color_munsell, convertColors = FALSE))

#' # calculate a mixed 150-200cm color ~"parent material"
jacobs2000$c_horizon_color <- profileApply(jacobs2000, function(p) {

  # and derive the parent material from the 150-200cm interval
  p150_200 <- glom(p, 150, 200, truncate = TRUE)
  p150_200$thickness <- p150_200$bottom - p150_200$top
  
  # subset colors and thickness
  clrs <- na.omit(horizons(p150_200)[,c('matrix_color_munsell','thickness')])
  
  # simulate a subtractive mixture using thickness as weight
  mixMunsell(
  clrs$matrix_color_munsell, 
  w = clrs$thickness, 
  mixingMethod = 'exact')$munsell
})

# segment profile into 1cm slices (for proper depth weighting)
jacobs2000$rubif <- profileApply(jacobs2000, function(p) {

  # sum the melanization index over the 0-100cm interval
  p0_100 <- hz_segment(p, 0:100)

  ccol <- parseMunsell(p$c_horizon_color, convertColors = FALSE)

  sum(harden.rubification(
    hue = p0_100$hue,
    chroma = as.numeric(p0_100$chroma),
    hue_ref = ccol$hue,
    chroma_ref = as.numeric(ccol$chroma)
  ), na.rm = TRUE)

})

jacobs2000$rubiforder <- order(jacobs2000$rubif)

# Plot in order of increasing Rubification index

plotSPC(jacobs2000,
color = "matrix_color",
label = "rubif",
plot.order = jacobs2000$rubiforder,
max.depth = 250
)

segments(
  x0 = 0.5, 
  x1 = length(jacobs2000) + 0.5, 
  y0 = c(0,100,150,200), 
  y1 = c(0,100,150,200), 
  lty = 2
)

# Add [estimated] parent material color swatches
trash <- sapply(seq_along(jacobs2000$c_horizon_color), function(i) {
  rect(i - 0.15, 250, i + 0.15, 225,
       col = parseMunsell(jacobs2000$c_horizon_color[jacobs2000$rubiforder[i]]))
})

Harmonize a property by profile-level denormalization for convenient visualization or analysis of ranges

Description

It is sometimes convenient to be able to "denormalize" to a SoilProfileCollection with fewer attributes but more profiles. This is helpful wherever calculations are made on a profile basis and ranges or repeated measures are depicted with multiple attributes per soil horizon.

harmonize is most commonly used for creating "comparison" soil profile sketches with plotSPC–where the thematic attribute is derived from multiple data sources or summary statistics (such as quantiles of a property for Low-RV-High). However, the method more generally applies wherever one wants to alias between multiple columns containing "similar" data as input to an algorithm.

Data are "harmonized" to a common attribute names specified by the names of list elements in x.names. Profiles are essentially duplicated. In order to satisfy uniqueness constraints of the SoilProfileCollection, the label from the sub-elements of x.names are used to disambiguate profiles. A new column in the site table is calculated to reflect these groupings and facilitate filtering. See examples below.

Usage

## S4 method for signature 'SoilProfileCollection'
harmonize(x, x.names, keep.cols = NULL, grp.name = "hgroup")

Arguments

Details

If attributes reflecting the same or similar property within a soil layer have different names (e.g. socQ05, socQ50, socQ95) it is sometimes inconvenient to work with them as multiple attributes within the same profile. These similar attributes may need to be analyzed together, or in sequence by profile, displayed using the same name or using a common scale. It is also useful to be able to alias different data sources that have the same attributes with different names.

Each list element in x.names specifies a single "harmonization," which is comprised of one or more mappings from new to old. Each named "sub-element" of x.names specifies the name and attribute to use for updating the profile ID and site table of the duplicated profiles.

Value

A (redundant) SoilProfileCollection, with one profile for each set of harmonizations specified by x.names.

Author(s)

Andrew G. Brown

Examples

### single source "harmonization" of single-profile with range -> single attribute, multi-profile

# make some test data
spc <- combine(lapply(1:10, random_profile, SPC = TRUE))

# assume that p1, p2 and p3 are the low RV and high quantiles for a hypothetical property "foo"
h1 <- harmonize(spc, x.names = list(foo = c(q05 = "p1", q50 = "p2", q95 = "p3")))

# inspect result
plotSPC(h1, color = "foo")

# filter with calculated "harmonized group" to get just RV profiles
plotSPC(subset(h1, hgroup == "q50"), color="foo")

### single source, two properties at once; with common labels: "method1" "method2"

# assume that p1, p2 are measurements by two (!=) methods for a hypothetical property "foo"
#             p3, p4 are measurements by same two methods for a hypothetical property "bar"
h3 <- harmonize(spc, x.names = list(foo = c(method1 = "p1", method2 = "p2"),
                                    bar = c(method1 = "p3", method2 = "p4")))
plotSPC(h3, color = "foo")
plotSPC(h3, color = "bar")
head(horizons(h3))

# a slight modification, "method 1" onlyused for "foo" and "method 3" for "bar"
h3 <- harmonize(spc, x.names = list(foo = c(method1 = "p1", method2 = "p2"),
                                    bar = c(method2 = "p3", method3 = "p4")))
plotSPC(h3, color = "foo") # note the pattern of values missing for foo (*_method3)
plotSPC(h3, color = "bar") #  likewise for bar (*_method1)

#' the new labels need not match across harmonizations -- not sure how useful this is but it works
h3 <- harmonize(spc, x.names = list(foo = c(method1 = "p1", method2 = "p2"),
                                    bar = c(method3 = "p3", method4 = "p4")))
plotSPC(h3, color = "foo") # note the pattern of values missing for foo (*_method 3 + 4)
plotSPC(h3, color = "bar") #  likewise for bar (*_method 1 + 2)

### two-source harmonization

# make test data
spc1 <- combine(lapply(LETTERS[1:5], random_profile, SPC = TRUE))
spc2 <- combine(lapply(letters[1:5], random_profile, SPC = TRUE))

h4 <- combine(list(harmonize(spc1, list(foo = c(transect1 = "p4"))),   # foo is p4 in dataset 1
                      harmonize(spc2, list(foo = c(transect2 = "p2")))))  # foo is p2 in dataset 2

# same property with different name in two different datasets
plotSPC(h4, color = "foo")

### many source harmonization

# make test datasets (n=10); highly redundant IDs (1:3 repeated)
spcs <- lapply(1:10, function(x) pbindlist(lapply(1:3, random_profile, SPC = TRUE)))

# randomly varying column name for demo (in each dataset, foo could could be p1 thru p5)
rcolname <- paste0("p", round(runif(10, 0.5, 5.5)))

# iterate over data sources
bigspc <- combine(lapply(1:length(spcs), function(i) {

  # assume each data source has a unique name for the property "foo"
  xn <- rcolname[i]

  # set names attribute to be equal to index i [creating unique profile IDs]
  #   i.e. 2_10 will be profile ID 2 from 10th dataset
  names(xn) <- i

  # harmonize each data source, using unique column name and target name "foo"
  harmonize(spcs[[i]], x.names = list(foo = xn))
}))

# inspect a subset
plotSPC(bigspc[1:30,], color = "foo")

Find horizons with colors darker than a Munsell hue, value, chroma threshold

Description

hasDarkColors returns a boolean value by horizon representing whether darkness thresholds are met. The code is fully vectorized and deals with missing data and optional thresholds.

Default arguments are set up for "5-3-3 colors" – the basic criteria for Mollic/Umbric epipedon/mineral soil darkness. Any of the thresholds or column names can be altered. Any thresholds that are set equal to NA will be ignored.

Usage

hasDarkColors(
  p,
  d_hue = NA,
  m_hue = NA,
  d_value = 5,
  d_chroma = NA,
  m_value = 3,
  m_chroma = 3,
  dhuenm = "d_hue",
  dvalnm = "d_value",
  dchrnm = "d_chroma",
  mhuenm = "m_hue",
  mvalnm = "m_value",
  mchrnm = "m_chroma"
)

Arguments

Value

Boolean value (for each horizon in p) reflecting whether "darkness" criteria are met.

Author(s)

Andrew G. Brown

Examples

# construct a fake profile
spc <- data.frame(id=1, taxsubgrp = "Lithic Haploxeralfs",
                  hzdesgn  = c("A","AB","Bt","BCt","R"),
                  hzdept   = c(0, 20, 32, 42,  49),
                  hzdepb   = c(20, 32, 42, 49, 200),
                  d_value  = c(5,   5,  5,  6,  NA),
                  m_value  = c(2.5, 3,  3,  4,  NA),
                  m_chroma = c(2,   3,  4,  4,  NA))

# promote to SoilProfileCollection
depths(spc) <- id ~ hzdept + hzdepb

# print results in table
data.frame(id = spc[[idname(spc)]],
           hz_desgn = spc$hzdesgn,
           has_dark_colors = hasDarkColors(spc))

Horizon Color Indices

Description

Calculate basic horizon-level color indices for a SoilProfileCollection. Basic indices do not require aggregation over the whole profile or comparison to a "reference" (e.g. parent material) color. Includes Hurst (1977) Redness Index, Barron-Torrent Redness Index (1986) and Buntley-Westin Index (1965). This is a wrapper method around several horizon-level indices. See the individual functions for more details.

Usage

horizonColorIndices(p, hue = "m_hue", value = "m_value", chroma = "m_chroma")

Arguments

Value

A data.frame containing unique pedon and horizon IDs and horizon-level color indices.

Author(s)

Andrew G. Brown

See Also

hurst.redness barron.torrent.redness.LAB buntley.westin.index

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

data(sp1)

# promote sp1 data to SoilProfileCollection
depths(sp1) <- id ~ top + bottom

# move site data
site(sp1) <- ~ group

# use Munsell color notation as horizon name
sp1$m <- sprintf("%s %s/%s", sp1$hue, sp1$value, sp1$chroma)

# compute indices
# merged into `sp1` with left-join on hzidname(sp1)
horizons(sp1) <- horizonColorIndices(sp1, hue="hue", value="value", chroma="chroma")

# visualize
par(mar = c(0, 1, 3, 1))
plotSPC(sp1, color='hurst_redness', name = 'm')
plotSPC(sp1, color='barron_torrent_redness', name = 'm')
plotSPC(sp1, color='buntley_westin', name = 'm')

Set horizon depth column names

Description

Set column name containing horizon ID

Get column names containing horizon depths

Usage

## S4 replacement method for signature 'SoilProfileCollection'
horizonDepths(object) <- value

## S4 method for signature 'SoilProfileCollection'
horizonDepths(object)

Arguments

Set horizon column names

Description

Set horizon column names

Get names of columns in horizon table.

Usage

## S4 replacement method for signature 'SoilProfileCollection'
horizonNames(object) <- value

## S4 method for signature 'SoilProfileCollection'
horizonNames(object)

Arguments

Retrieve horizon data from SoilProfileCollection

Description

Get horizon data from SoilProfileCollection. Result is returned in the same data.frame class used to initially construct the SoilProfileCollection.

Horizon data in an object inheriting from data.frame can easily be added via merge (LEFT JOIN). There must be one or more same-named columns (with at least some matching data) on the left and right hand side to facilitate the join: horizons(spc) <- newdata

Usage

## S4 method for signature 'SoilProfileCollection'
horizons(object)

## S4 replacement method for signature 'SoilProfileCollection'
horizons(object) <- value

Arguments

Examples

# load test data
data(sp2)

# promote to SPC
depths(sp2) <- id ~ top + bottom

# assign true to surface horizon
newdata <- data.frame(top = 0,
                      newvalue = TRUE)

# do left join
horizons(sp2) <- newdata

# inspect site table: newvalue TRUE only for horizons
#  with top depth equal to zero
horizons(sp2)

Munsell Hue Reference and Position Searching

Description

The 40 Munsell hues are typically arranged from 5R to 2.5R moving clock wise on the unit circle. This function matches a vector of hues to positions on that circle, with options for setting a custom origin or search direction.

This function is fully vectorized.

Usage

huePosition(
  x,
  returnHues = FALSE,
  includeNeutral = FALSE,
  origin = "5R",
  direction = c("cw", "ccw")
)

Arguments

Value

A vector of integer hue positions is returned, of the same length and order as x. If returnHues = TRUE, then all hue names and ordering are returned and x is ignored.

Author(s)

D.E. Beaudette

References

• Soil Survey Technical Note 2 wayback machine URL

• Munsell book of color. 1976. Macbeth, a Division of Kollmorgen Corp., Baltimore, MD.

See Also

colorContrast, huePositionCircle

Examples

# get hue ordering for setting levels of a factor
huePosition(returnHues = TRUE)

# get hue ordering including N (neutral)
huePosition(returnHues = TRUE, includeNeutral = TRUE)

# get position of the '10YR' hue, relative to standard origin of '5R'
# should be 7
huePosition(x = '10YR')

# get position of the '10YR' hue, relative to standard origin of '5YR'
# should be 3
huePosition(x = '10YR', origin = '5YR')

# visualize
op <- par(mar = c(0, 0, 0, 0), fg = 'white', bg = 'black')

huePositionCircle(huePosition(returnHues = TRUE, origin = '5YR'))

par(op)

Visual Description of Munsell Hue Ordering

Description

Munsell hues are arranged on the unit circle with "neutral" at the center.

Usage

huePositionCircle(
  hues = huePosition(returnHues = TRUE),
  value = 6,
  chroma = 10,
  chip.cex = 5.5,
  label.cex = 0.66,
  seg.adj = 0.8,
  seg.col = grey(0.4),
  plot = TRUE,
  simulateCVD = NULL,
  CVDseverity = 1
)

Arguments

Value

an invisible data.frame of data used to create the figure

Note

The best results are obtained when setting margins to zero, and inverting foreground / background colors. For example: par(mar = c(0, 0, 0, 0), fg = 'white', bg = 'black').

References

Munsell book of color. 1976. Macbeth, a Division of Kollmorgen Corp., Baltimore, MD.

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))


# better graphics defaults
op <- par(
mar = c(0, 0, 0, 0), 
fg = 'white', 
bg = 'black', 
xpd = NA
)

# full set of hues, as generated by huePosition(returnHues = TRUE)
huePositionCircle()

# just a few hues
huePositionCircle(hues = c('5R', '5Y', '5G', '5B', '5P'))

# adjust Munsell value and chroma
huePositionCircle(value = 3, chroma = 6)

# reset graphics state
par(op)

Hurst (1977) Redness Index

Description

Calculate Redness Index after Hurst (1977) "Visual estimation of iron in saprolite" DOI: 10.1130/0016-7606(1977)88<174:VEOIIS>2.0.CO;2. Accepts vectorized inputs for hue, value and chroma, produces vector output.

Usage

hurst.redness(hue, value, chroma)

Arguments

Value

A numeric vector of horizon redness index (lower values = redder).

Author(s)

Andrew G. Brown

References

Hurst, V.J. (1977) Visual estimation of iron in saprolite. GSA Bulletin. 88(2): 174–176. doi: https://doi.org/10.1130/0016-7606(1977)88<174:VEOIIS>2.0.CO;2

Horizons Above or Below

Description

Horizons Above or Below

Usage

hzAbove(x, ..., offset = 1, SPC = TRUE, simplify = SPC)

hzBelow(x, ..., offset = 1, SPC = TRUE, simplify = SPC)

hzOffset(x, hzid, offset, SPC = FALSE, simplify = TRUE)

Arguments

Details

To minimize likelihood of issues with non-standard evaluation context, especially when using hzAbove()/hzBelow() inside another function, all expressions used in ... should be in terms of variables that are in the horizon data frame.

Value

A SoilProfileCollection (when SPC = TRUE) or a vector of horizon row indices (when SPC = FALSE and simplify = TRUE) or a list (when SPC = FALSE and simplify = FALSE))

Examples

data(sp4)
depths(sp4) <- id ~ top + bottom

# get the horizon above the last horizon (j-index of bottom horizon minus 1)
hzAbove(sp4, hzID(sp4) %in% getLastHorizonID(sp4))

# get horizons below the last horizon (none; j-index of bottom horizon plus 1)
hzBelow(sp4, hzID(sp4) %in% getLastHorizonID(sp4))

Tests of horizon depth logic

Description

Function used internally by checkHzDepthLogic(), glom() and various other functions that operate on horizon data from single soil profiles and require a priori depth logic checks. Checks for bottom depths less than top depth / bad top depth order ("depthLogic"), bottom depths equal to top depth ("sameDepth"), overlaps/gaps ("overlapOrGap") and missing depths ("missingDepth"). Use names(res)[res] on result res of hzDepthTest() to to determine type of logic error(s) found – see examples below.

Usage

hzDepthTests(top, bottom = NULL)

Arguments

Value

A named logical vector containing TRUE for each type of horizon logic error found in the given data.

Author(s)

Andrew G. Brown & Dylan E. Beaudette

Examples

# no logic errors
res <- hzDepthTests(top = c(0,10,20,30), bottom = c(10,20,30,50))
names(res)[res]

# bottom < top
hzDepthTests(top = c(10,20,30,50), bottom = c(0,10,20,30))
names(res)[res]

# bottom == top
hzDepthTests(top = c(10,20,30,50), bottom = c(0,20,20,30))
names(res)[res]

# overlap
hzDepthTests(top = c(0,5,20,30), bottom = c(10,20,30,50))
names(res)[res]

# gap
hzDepthTests(top = c(0,15,20,30), bottom = c(10,20,30,50))
names(res)[res]

# missing
hzDepthTests(c(0,15,NA,30),c(10,NA,30,50))
names(res)[res]

Get horizon designation column name

Description

Get horizon designation names

Usage

## S4 method for signature 'SoilProfileCollection'
hzDesgn(object)

Arguments

Convert Horizon Boundary Distinctness to Vertical Offset

Description

This function will convert USDA-NCSS horizon boundary distinctness codes into vertical (+/-) offsets in cm, based on the Field Book for Describing and Sampling Soils, version 3.0.

Usage

hzDistinctnessCodeToOffset(
  x,
  codes = c("very abrupt", "abrupt", "clear", "gradual", "diffuse"),
  offset = c(0.5, 2, 5, 15, 20)/2
)

Arguments

Details

The default offsets are based on the high-end of ranges presented in "transitional zone thickness criteria" from the Field Book version 3.0 (page 2-6). Offsets are returned as 1/2 of the transitional zone thickness so that horizon boundaries can be adjusted up/down from horizon depths. See plotSPC, specifically the hz.distinctness.offset argument for visualization ideas. Missing data in x (NA) or codes that are not defined in codes are returned as 0 offsets.

Either format (or mixture) are accepted, case insensitive:

• terms: c('very abrupt', 'abrupt', 'clear', 'gradual', 'diffuse')

• coded values: ⁠c('v', 'a', 'c', 'g', d')⁠

Additional examples are available in the Visualization of Horizon Boundaries tutorial.

Value

vector of offsets with same length as x

Author(s)

D.E. Beaudette

References

Field Book for Describing and Sampling Soils, version 3.0

See Also

plotSPC

Examples

# example data
data(sp1)

# compute 1/2 transitional zone thickness from distinctness codes
sp1$hzdo <- hzDistinctnessCodeToOffset(sp1$bound_distinct)

# convert colors from Munsell to hex-encoded RGB
sp1$soil_color <- with(sp1, munsell2rgb(hue, value, chroma))

# promote to SoilProfileCollection
depths(sp1) <- id ~ top + bottom
hzdesgnname(sp1) <- 'name'

# adjust margins
op <- par(mar=c(0,0,0,1.5))

# sketches, adjust width, adjust text size, include coded hz distinctness offsets
plotSPC(sp1, width=0.3, cex.names=0.75, hz.distinctness.offset = 'hzdo')

# clean-up
par(op)

Set horizon IDs

Description

Set vector containing horizon IDs

Get vector containing horizon IDs

Usage

## S4 replacement method for signature 'SoilProfileCollection'
hzID(object) <- value

## S4 method for signature 'SoilProfileCollection'
hzID(object)

Arguments

Get horizon-level metadata

Description

Get idname(object) and hzidname(object), with hzdesgnname(object), hztexclname(object) (if defined)

Usage

## S4 method for signature 'SoilProfileCollection'
hzMetadata(object)

Arguments

Convert Horizon Boundary Topography to Line Type

Description

This function will convert USDA-NCSS horizon boundary topography codes into line types, based on the Field Book for Describing and Sampling Soils, version 3.0.

Usage

hzTopographyCodeToLineType(
  x,
  codes = c("smooth", "wavy", "irregular", "broken"),
  lty = c(1, 2, 3, 4)
)

Arguments

Details

Visualization of horizon boundary topography can be difficult, line type offers an additional visual cue. See hzTopographyCodeToOffset for an offset-based approach. Additional examples are available in the Visualization of Horizon Boundaries tutorial. Missing data in x (NA) or codes that are not defined in codes are returned as line type 1.

Either format (or mixture) are accepted, case insensitive:

• terms: c('smooth', 'wavy', 'irregular', 'broken')

• coded values: c('s', 'w', 'i', 'b')

Value

vector of line types with same length as x

Author(s)

D.E. Beaudette

References

Field Book for Describing and Sampling Soils, version 3.0

See Also

plotSPC, hzTopographyCodeToOffset

Convert Horizon Boundary Topography to Vertical Offset

Description

This function will convert USDA-NCSS horizon boundary topography codes into a vertical offset, suitable for use in plotSPC. Default values are reasonable starting points for encoding smooth, wavy, irregular, or broken style horizon boundary topography as defined in Field Book for Describing and Sampling Soils, version 3.0.

Usage

hzTopographyCodeToOffset(
  x,
  codes = c("smooth", "wavy", "irregular", "broken"),
  offset = c(0, 4, 8, 12)
)

Arguments

Details

Additional examples are available in the Visualization of Horizon Boundaries tutorial. Missing data in x (NA) or codes that are not defined in codes are returned with an offset of 0.

Either format (or mixture) are accepted, case insensitive:

• terms: c('smooth', 'wavy', 'irregular', 'broken')

• coded values: c('s', 'w', 'i', 'b')

Value

vector of vertical offsets with same length as x

Author(s)

D.E. Beaudette

References

Field Book for Describing and Sampling Soils, version 3.0

See Also

plotSPC

Horizon Transition Probabilities

Description

Functions for creating and working with horizon (sequence) transition probability matrices.

See the following tutorials for some ideas:

• horizon designation TP

• soil color TP

Usage

hzTransitionProbabilities(
  x,
  name = GHL(x, required = TRUE),
  loopTerminalStates = FALSE
)

mostLikelyHzSequence(mc, t0, maxIterations = 10)

Arguments

Value

A square matrix of transition probabilities. See examples.

The function genhzTableToAdjMat() returns a square adjacency matrix. See examples.

The function mostLikelyHzSequence() returns the most likely sequence of horizons, given a markovchain object initialized from horizon transition probabilities and an initial state, t0. See examples.

Note

These functions are still experimental and subject to change.

Author(s)

D.E. Beaudette

See Also

generalize.hz()

Examples

data(sp4)
depths(sp4) <- id ~ top + bottom

# horizon transition probabilities: row -> col transitions
(tp <- hzTransitionProbabilities(sp4, 'name'))

Dissolving horizon boundaries by grouping variables

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.

Usage

hz_dissolve(
  object,
  by,
  idcol = "id",
  depthcols = c("top", "bottom"),
  collapse = FALSE,
  order = FALSE
)

dissolve_hz(
  object,
  by,
  id = "idcol",
  hztop = "top",
  hzbot = "bottom",
  collapse = FALSE,
  order = FALSE
)

Arguments

Details

This function assumes the profiles and horizons within the object follow the logic defined by checkHzDepthLogic (e.g. records are ordered sequentially by id, hztop, and hzbot and without gaps). If the records are not ordered, set the order = TRUE.

Value

A data.frame with the original idcol, by grouping variables, and non-consecutive horizon depths.

Author(s)

Stephen Roecker

See Also

hz_lag(), hz_intersect(), hz_segment() , checkHzDepthLogic()

Examples

# example 1
data(jacobs2000)
spc <- jacobs2000

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 <- 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]
test[test$id == "92-1", ]


# example 2
df <- data.frame(
    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")
    )

df

hz_dissolve(df, c("genhz", "texcl"))
hz_dissolve(df, c("genhz", "texcl"), collapse = TRUE)

test <- hz_dissolve(df, "genhz")
subset(test, value == "2Bt")

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.

Usage

hz_intersect(x, y, idcol = "id", depthcols = c("top", "bottom"))

Arguments

Details

.

Value

A data.frame with harmonized depth intervals (i.e. segment_id) and columns from both of the original data.frame. If both 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 data.frame they originated from.

Author(s)

Stephen Roecker

See Also

hz_dissolve(), hz_lag(), hz_segment()

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

hz_dissolve(h, "by")

hz_intersect(x = hz_dissolve(h, "by"), y = h)

hi <- hz_intersect(x = h, y = hz_dissolve(h, "by"))
aggregate(clay ~ dissolve_id, data = hi, mean)

Find lagged horizon values

Description

This function finds adjacent values to a horizon values at lagged distances.

Usage

hz_lag(
  object,
  lag = 1,
  unit = "index",
  idcol = "id",
  depthcols = c("top", "bottom"),
  order = FALSE
)

Arguments

Details

.

Value

A data.frame with lagged values.

Author(s)

Stephen Roecker

See Also

hz_dissolve(), hz_intersect(), hz_segment()

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

hz_lag(h)

hz_lag(h, -1)

hz_lag(h, 10:15, unit = "depth")

transform(cbind(h, lag = hz_lag(h)), 
  clay_dif = lag.clay_bot.1 - clay,
  texcl_contrast = paste0(texcl, "-", lag.texcl_bot.1)
)

Segmenting of Soil Horizon Data by Depth Interval

Description

This function segments or subdivides horizon data from a SoilProfileCollection or data.frame by depth interval (e.g. c(0, 10), c(0, 50), or 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. "segment_id") are added to each horizon record that correspond with their depth interval (e.g. 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.

Usage

hz_segment(object, intervals, trim = TRUE, depthcols = c("top", "bottom"))

segment(object, intervals, trim = TRUE, hzdepcols = c("top", "bottom"))

Arguments

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().

Value

Either a SoilProfileCollection or data.frame with the original horizon data segmented by depth intervals. There are usually more records in the resulting object, one for each time a segment interval partially overlaps with a horizon. A new column called segment_id identifying the depth interval is added.

Author(s)

Stephen Roecker

See Also

dice(), glom(), hz_dissolve(), hz_lag(), hz_intersect()

Examples

# example data
data(sp1)

# upgrade to SPC
depths(sp1) <- id ~ top + bottom

# segment and trim
z <- hz_segment(sp1, intervals = c(0, 10, 20, 30), trim = TRUE)

# display segment labels
# note that there are new horizon boundaries at segments
par(mar = c(0, 0, 3, 1))
plotSPC(z, color = 'segment_id', width = 0.3)

# highlight new horizon records
par(mar = c(0, 0, 2, 1))
plotSPC(z, color = NA, default.color = NA, width = 0.3, lwd = 1)
plotSPC(sp1, color = NA, default.color = NA, 
width = 0.3, lwd = 3, add = TRUE, name = NA, print.id = FALSE)
legend('top', horiz = TRUE, 
legend = c('original', 'segmented'), 
lwd = c(1, 3), cex = 0.85, bty = 'n')


# same results as slab()
# 10 random profiles
s <- lapply(1:10, random_profile, n_prop = 1, SPC = TRUE, method = 'random_walk')
s <- combine(s)

a.slab <- slab(s, fm = ~ p1, slab.structure = c(0, 10, 20, 30), slab.fun = mean, na.rm = TRUE)

z <- hz_segment(s, intervals = c(0, 10, 20, 30), trim = TRUE)
z <- horizons(z)
z$thick <- z$bottom - z$top

a.segment <- sapply(split(z, z$segment_id), function(i) {
  weighted.mean(i$p1, i$thick)
})


res <- data.frame(
  slab = a.slab$value,
  segment = a.segment,
  diff = a.slab$value - a.segment
)

print(res)
res$diff < 0.001



data(sp5)

# segment by upper 25-cm
test1 <- hz_segment(sp5, intervals = c(0, 100))
print(test1)
nrow(test1)
print(object.size(test1), units = "Mb")

# segment by 1-cm increments
test2 <- hz_segment(sp5, intervals = 0:100)
print(test2)
nrow(test2)
print(object.size(test2), units = "Mb")


# segment and aggregate
test3 <- hz_segment(horizons(sp5), 
                 intervals = c(0, 5, 15, 30, 60, 100, 200), 
                 depthcols = c("top", "bottom")
)
test3$hzthk <- test3$bottom - test3$top
test3_agg <- by(test3, test3$segment_id, function(x) {
  data.frame(
    hzID = x$hzID[1],
    segment_id = x$segment_id[1],
    average = weighted.mean(x$clay, w = x$hzthk)
  )
})
test3_agg <- do.call("rbind", test3_agg)

head(test3_agg)

Allocate Particle Size Class for the Control Section.

Description

This function aggregates information in the horizon table and allocates it to the particle size class for the control section.

Usage

hz_to_taxpartsize(
  x,
  y,
  taxpartsize = "taxpartsize",
  clay = "clay",
  idcol = "id",
  depthcols = c("top", "bottom")
)

Arguments

Details

This function differs from texture_to_taxpartsize in that is aggregates the results of texture_to_taxpartsize, and accounts for strongly contrasting particle size classes.

Value

A data.frame object containing the original idcol, the aggregated particle size control section allocation, and an aniso column to indicate more than one contrasting class.

Author(s)

Stephen Roecker

See Also

texture_to_taxpartsize(), lookup_taxpartsize()

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

# set required metadata for estimatePSCS()
hzdesgnname(h) <- "hzname"
hztexclname(h) <- "texcl"
hzmetaname(h, "clay") <- "clay"

pscs <- data.frame(id = h$id, rbind(estimatePSCS(h)))
names(pscs)[2:3] <- c("top", "bottom")

hz_to_taxpartsize(horizons(h), pscs)
 

Get or Set Horizon Designation Column Name

Description

hzdesgnname(): Get column name containing horizon designations

⁠hzdesgnname<-⁠: Set horizon designation column name

Usage

## S4 method for signature 'SoilProfileCollection'
hzdesgnname(object, required = FALSE)

## S4 replacement method for signature 'SoilProfileCollection'
hzdesgnname(object, required = FALSE) <- value

Arguments

Details

Store the column name containing horizon designations or other identifiers in the metadata slot of the SoilProfileCollection.

See Also

hzDesgn()

Examples

data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# set horizon designation column
hzdesgnname(sp1) <- "name"

# get horizon designation column
hzdesgnname(sp1)

Set horizon ID column name

Description

Set unique horizon ID column name

Get column name containing unique horizon ID

Usage

## S4 replacement method for signature 'SoilProfileCollection'
hzidname(object) <- value

## S4 method for signature 'SoilProfileCollection'
hzidname(object)

Arguments

Examples

data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# create new horizon ID
sp1$hzIDrev <- rev(sp1$hzID)

# set horizon designation column
hzidname(sp1) <- "hzIDrev"

# get horizon designation column
hzidname(sp1)

Get or Set Horizon Metadata Column Name

Description

hzmetaname(): Get column name containing horizon data of interest

⁠hzmetaname<-⁠: Set horizon designation column name

Usage

## S4 method for signature 'SoilProfileCollection'
hzmetaname(object, attr, required = FALSE)

## S4 replacement method for signature 'SoilProfileCollection'
hzmetaname(object, attr, required = FALSE) <- value

Arguments

Details

Store the column name containing a specific type of horizon data in the metadata slot of the SoilProfileCollection.

See Also

guessHzAttrName()

Examples

data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# set important metadata columns
hzdesgnname(sp1) <- "name"
hztexclname(sp1) <- "texture"

# set custom horizon property (clay content) column
hzmetaname(sp1, "clay") <- "prop"

# inspect metadata list
metadata(sp1)

# get horizon clay content column
hzmetaname(sp1, "clay")

# uses hzdesgname(), hztexclname(), hzmetaname(attr="clay") in function definition
estimatePSCS(sp1)

Get or Set Horizon Texture Class Column Name

Description

hztexclname(): Get column name containing horizon designation name

⁠hztexclname<-⁠: Set horizon texture class column name for a SoilProfileCollection

Usage

## S4 method for signature 'SoilProfileCollection'
hztexclname(object, required = FALSE)

## S4 replacement method for signature 'SoilProfileCollection'
hztexclname(object, required = FALSE) <- value

Arguments

Details

Store the column name containing horizon texture classes or other identifiers in the metadata slot of the SoilProfileCollection.

Examples

data(sp1)

# promote to SPC
depths(sp1) <- id ~ top + bottom

# set horizon texture class column
hztexclname(sp1) <- "texture"

# get horizon texture class column
hztexclname(sp1)

Get profile ID column name

Description

Get column name containing unique profile IDs

Usage

## S4 method for signature 'SoilProfileCollection'
idname(object)

Arguments

Initialize Spatial Data in a SoilProfileCollection

Description

⁠initSpatial()<-⁠: Set the column names containing spatial data and the corresponding coordinate reference system for a SoilProfileCollection.

getSpatial(): Get spatial data associated with a SoilProfileCollection

Usage

## S4 replacement method for signature 'SoilProfileCollection,ANY,ANY'
initSpatial(object, crs = NULL) <- value

## S4 replacement method for signature 'SoilProfileCollection,ANY,character'
initSpatial(object, crs = NULL) <- value

## S4 method for signature 'SoilProfileCollection'
getSpatial(object)

## S4 method for signature 'SoilProfileCollection'
coordinates(obj)

## S4 replacement method for signature 'SoilProfileCollection,ANY'
coordinates(object) <- value

## S4 replacement method for signature 'SoilProfileCollection,character'
coordinates(object) <- value

Arguments

See Also

prj()

Examples

data(sp5)

# coordinates are stored in x and y column of site
sp5$x <- rnorm(length(sp5))
sp5$y <- rnorm(length(sp5))

# coordinates takes a formula object as input
initSpatial(sp5) <- ~ x + y

# optionally specify Coordinate Reference System (crs) on left-hand side
initSpatial(sp5, crs = "OGC:CRS84") <- ~ x + y

Make High Contrast Label Colors

Description

Generate a vector of white or black label colors conditioned on a vector of colors to maximize label contrast.

Usage

invertLabelColor(colors, threshold = 0.65)

Arguments

Value

vector of label colors

Author(s)

D.E. Beaudette

Examples

# test with shades of grey
s <- seq(0, 1, by = 0.05)
cols <- grey(s)
soilPalette(cols, lab = as.character(s))

# test with 10YR x/3
m <- sprintf('10YR %s/3', 1:8)
cols <- parseMunsell(m)
soilPalette(cols, lab = m)

Check for "empty" profiles in a SoilProfileCollection

Description

"Empty" profiles are used as placeholders for positions in a SoilProfileCollection These profiles result from operations that remove or extract portions of horizons from source profiles.

Usage

## S4 method for signature 'SoilProfileCollection'
isEmpty(object, ...)

Arguments

Details

In a SoilProfileCollection an empty profile occurs when it has one horizon, with NA top and bottom depths. Generally all non-profile ID site and horizon-level values are all also NA, but only the depths are checked by isEmpty().

Value

logical. Vector of length equal to number of profiles in object. Returns TRUE when a profile has one horizon with NA top and bottom depths

Soil Morphologic Data from Jacobs et al. 2002.

Description

Select soil morphologic data from "Redoximorphic Features as Indicators of Seasonal Saturation, Lowndes County, Georgia". This is a useful sample dataset for testing the analysis and visualization of redoximorphic features.

Usage

data(jacobs2000)

Format

A SoilProfileCollection object.

References

Jacobs, P. M., L. T. West, and J. N. Shaw. 2002. Redoximorphic Features as Indicators of Seasonal Saturation, Lowndes County, Georgia. Soil Sci. Soc. Am. J. 66:315-323. doi:doi:10.2136/sssaj2002.3150

Examples

# keep examples from using more than 2 cores
data.table::setDTthreads(Sys.getenv("OMP_THREAD_LIMIT", unset = 2))

# load
data(jacobs2000)

# basic plot
par(mar = c(0, 1, 3, 1.5))
plotSPC(jacobs2000, name='name', color='matrix_color', width=0.3)

# add concentrations
addVolumeFraction(jacobs2000, 'concentration_pct',
                  col = jacobs2000$concentration_color, pch = 16, cex.max = 0.5)

# add depletions
plotSPC(jacobs2000, name='name', color='matrix_color', width=0.3)
addVolumeFraction(jacobs2000, 'depletion_pct',
                  col = jacobs2000$depletion_color, pch = 16, cex.max = 0.5)

# time saturated
plotSPC(jacobs2000, color='time_saturated', cex.names=0.8, col.label = 'Time Saturated')


# color contrast: matrix vs. concentrations
cc <- colorContrast(jacobs2000$matrix_color_munsell, jacobs2000$concentration_munsell)
cc <- na.omit(cc)

cc <- cc[order(cc$dE00), ]
cc <- unique(cc)

par(bg = 'black', fg = 'white')
colorContrastPlot(cc$m1[1:10], cc$m2[1:10], labels = c('matrix', 'concentration'))
colorContrastPlot(cc$m1[11:21], cc$m2[11:21], labels = c('matrix', 'concentration'))


# color contrast: depletion vs. concentrations
cc <- colorContrast(jacobs2000$depletion_munsell, jacobs2000$concentration_munsell)
cc <- na.omit(cc)

cc <- cc[order(cc$dE00), ]
cc <- unique(cc)

par(bg = 'black', fg = 'white')
colorContrastPlot(cc$m1, cc$m2, labels = c('depletion', 'concentration'))

Get the number of profiles in a SoilProfileCollection

Description

Get the number of profiles in a SoilProfileCollection

Usage

## S4 method for signature 'SoilProfileCollection'
length(x)

Arguments

Ranking Systems for USDA Taxonomic Particle-Size and Substitute Classes of Mineral Soils

Description

Generate a lookup table of USDA Particle-Size and Substitute Classes names, ranked according to approximate particle size

Usage

lookup_taxpartsize()

Value

A data.frame with a rank column, taxonomic family particle size class, and a flag for contrasting.

Author(s)

Stephen Roecker

References

Field Book for Describing and Sampling Soils, version 3.0

See Also

hz_to_taxpartsize(), texture_to_taxpartsize(), SoilTextureLevels()

Examples

# class codes
lu <- lookup_taxpartsize()

idx <- lu$contrasting == FALSE

lu$taxpartsize[idx]

lu$rank[as.integer(lu$taxpartsize)[idx]]

Eliminate duplicate instances of profile IDs in a list of SoilProfileCollections

Description

@description Experimental function to "clean" list input where duplicates exist (that would otherwise prevent pbindlist). Useful for queries that may have overlapping instances of the same data, for instance a list of SoilProfileCollections where each list element contains profiles gathered from a set of (potentially overlapping) extents.

Usage

lunique(l)

Arguments

Value

A list of SoilProfileCollections, with duplicate profile IDs removed.

Author(s)

Andrew G. Brown

Examples

data(sp5)

# EXAMPLE #1 -- resolving overlap

# 6 profiles in four sets, and 5,6,7 are missing
input <- lapply(list(c(1,3,4), c(2,2,3), NA, c(8,9,1)), function(idx) {
      if(!all(is.na(idx)))
       sp5[idx,]
})

output <- lunique(input)

# 6 profiles are in final SPC; 5,6,7 are missing
match(profile_id(pbindlist(output)), profile_id(sp5))

# EXAMPLE #2 -- exact duplicates

# deliberately duplicate an SPC
sp5_2 <- sp5
res <- lunique(list(sp5, sp5_2))

# the number of profiles in first element is equal to number in sp5
length(res[[1]]) == length(sp5)

# second list element contains NA b/c all uniques are in #1
res[[2]]

Get the maximum bottom depth in a SoilProfileCollection

Description

Get the deepest depth of description out of all profiles in a SoilProfileCollection. Data missing one or more of: bottom depth, profile ID, or any optional attribute are omitted using complete.cases.

Usage

## S4 method for signature 'SoilProfileCollection'
max(x, v = NULL, na.rm = TRUE)

Arguments

Retrieve metadata from SoilProfileCollection

Description

Get metadata from SoilProfileCollection. Result is a list. Two entries (aqp_df_class, depth_units) should not be edited in the metadata list directly. There are methods that facilitate changing them – and propagating their changes throughout the collection. Otherwise, metadata list is a free-form slot used to store arbitrary information about the data, how it was collected, citations, etc.

Usage

## S4 method for signature 'SoilProfileCollection'
metadata(object)

## S4 replacement method for signature 'SoilProfileCollection'
metadata(object) <- value

Arguments

Examples

data(sp5)

# replace default metadata with itself
metadata(sp5) <- metadata(sp5)

# set new metadata attribute value
metadata(sp5)$newvalue <- 'foo'

# get metadata attribute
metadata(sp5)$newvalue

Get the minimum bottom depth in a SoilProfileCollection

Description

Get the shallowest depth of description out of all profiles in a SoilProfileCollection. Data missing one or more of: bottom depth, profile ID, or any optional attribute are omitted using complete.cases.

Usage

## S4 method for signature 'SoilProfileCollection'
min(x, v = NULL, na.rm = TRUE)

Arguments

Missing Data Grid

Description

Generate a levelplot of missing data from a SoilProfileCollection object.

Usage

missingDataGrid(
  s,
  max_depth,
  vars,
  filter.column = NULL,
  filter.regex = NULL,
  cols = NULL,
  ...
)

Arguments

Details

This function evaluates a ⁠missing data fraction⁠ based on slice-wise evaluation of named variables in a SoilProfileCollection object.

Value

A data.frame describing the percentage of missing data by variable.

Note

A lattice graphic is printed to the active output device.

Author(s)

D.E. Beaudette

See Also

slice

Examples

# 10 random profiles
set.seed(10101)
s <- lapply(as.character(1:10), random_profile)
s <- do.call('rbind', s)

# randomly sprinkle some missing data
s[sample(nrow(s), 5), 'p1'] <- NA
s[sample(nrow(s), 5), 'p2'] <- NA
s[sample(nrow(s), 5), 'p3'] <- NA

# set all p4 and p5 attributes of `soil 1' to NA
s[which(s$id == '1'), 'p5'] <- NA
s[which(s$id == '1'), 'p4'] <- NA

# upgrade to SPC
depths(s) <- id ~ top + bottom

# plot missing data via slicing + levelplot
missingDataGrid(
  s,
  max_depth = 100,
  vars = c('p1', 'p2', 'p3', 'p4', 'p5'),
  main='Missing Data Fraction'
)

Mix Munsell Colors via Spectral Library

Description

Simulate mixing of colors in Munsell notation, similar to the way in which mixtures of pigments operate.

Usage

mixMunsell(
  x,
  w = rep(1, times = length(x))/length(x),
  mixingMethod = c("exact", "reference", "estimate", "adaptive"),
  n = 1,
  keepMixedSpec = FALSE,
  distThreshold = 0.025,
  ...
)

Arguments

Details

See the expanded tutorial for examples.

An accurate simulation of pigment mixtures ("subtractive" color mixtures) is incredibly complex due to factors that aren't easily measured or controlled: pigment solubility, pigment particle size distribution, water content, substrate composition, and physical obstruction to name a few. That said, it is possible to simulate reasonable, subtractive color mixtures given a reference spectra library (350-800nm) and some assumptions about pigment qualities and lighting. For the purposes of estimating a mixture of soil colors (these are pigments after all) we can relax these assumptions and assume a standard light source. The only missing piece is the spectral library for all Munsell chips in our color books.

Thankfully, Scott Burns has outlined the entire process, and Paul Centore has provided a Munsell color chip reflectance spectra library. The estimation of a subtractive mixture of soil colors can proceed as follows:

1. look up the associated spectra for each color in x

2. compute the weighted (w argument) geometric mean of the spectra

3. convert the spectral mixture to the closest Munsell color via:

• search for the closest n matching spectra in the reference library (mixtureMethod = 'reference')

• direct conversion of spectra to closest Munsell color via spec2Munsell() (mixtureMethod = 'exact')

1. suggest resulting Munsell chip(s) as the best candidate for a simulated mixture

Key assumptions include:

• similar particle size distribution

• similar mineralogy (i.e. pigmentation qualities)

• similar water content.

For the purposes of estimating (for example) a "mixed soil color within the top 18cm of soil" these assumptions are usually valid. Again, these are estimates that are ultimately "snapped" to the nearest chip and not do not need to approach the accuracy of paint-matching systems.

A message is printed when scaledDistance is larger than 1.

Value

A data.frame with the closest matching Munsell color(s):

• munsell: Munsell notation of the n-closest spectra

• distance: spectral (Gower) distance to the n-closest spectra

• scaledDistance: spectral distance scaled by distThreshold

• mixingMethod: method used for each mixture

When keepMixedSpec = TRUE then a list:

• mixed: a data.frame containing the same elements as above

• spec: spectra for the 1st closest match

Author(s)

D.E. Beaudette

References

Marcus, R.T. (1998). The Measurement of Color. In K. Nassau (Ed.), Color for Science, Art, and Technology (pp. 32-96). North-Holland.

• inspiration / calculations based on the work of Scott Burns

• related discussion on Stack Overflow

• spectral library source

See Also

munsell.spectra

Calculate the minimum thickness requirement for Mollic epipedon

Description

Utilize horizon depths, designations and textures in a profile to estimate the thickness requirement for the Mollic or Umbric epipedon, per criterion 6 in the U.S. Keys to Soil Taxonomy (12th Edition).

Usage

mollic.thickness.requirement(
  p,
  hzdesgn = hzdesgnname(p, required = TRUE),
  texcl.attr = hztexclname(p, required = TRUE),
  clay.attr = hzmetaname(p, "clay", required = TRUE),
  truncate = TRUE
)

Arguments

Value

A unit length numeric vector containing Mollic or Umbric epipedon minimum thickness requirement.

Author(s)

Andrew G. Brown

Examples

# construct a fake profile
spc <- data.frame(id=1, taxsubgrp = "Lithic Haploxeralfs",
                  hzname   = c("A","AB","Bt","BCt","R"),
                  hzdept   = c(0,  20, 32, 42,  49),
                  hzdepb   = c(20, 32, 42, 49, 200),
                  prop     = c(18, 22, 28, 24,  NA),
                  texcl    = c("l","l","cl", "l","br"),
                  d_value  = c(5,   5,  5,  6,  NA),
                  m_value  = c(2.5, 3,  3,  4,  NA),
                  m_chroma = c(2,   3,  4,  4,  NA))

# promote to SoilProfileCollection
depths(spc) <- id ~ hzdept + hzdepb
hzdesgnname(spc) <- 'hzname'
hztexclname(spc) <- 'texcl'

# print results in table
data.frame(id = spc[[idname(spc)]],
           thickness_req = mollic.thickness.requirement(spc, clay.attr='prop'),
           thickness_req_nobound = mollic.thickness.requirement(spc,
                                        clay.attr='prop', truncate=FALSE))

Munsell to sRGB Lookup Table for Common Soil Colors

Description

A lookup table of interpolated Munsell color chips for common soil colors.

Usage

data(munsell)

Format

A data.frame with 8825 rows.

• hue: Munsell Hue, upper case

• value: Munsell Value

• chroma: Munsell Chroma

• r: sRGB "red" value (0-1)

• g: sRGB "green" value (0-1)

• b: sRGB "blue" value (0-1)

• L: CIELAB "L" coordinate

• A: CIELAB "A" coordinate

• B: CIELAB "B" coordinate

Details

See munsell2rgb for conversion examples. Values are referenced to the D65 standard illuminant.

Source

Color chip XYZ values: https://www.rit.edu/science/munsell-color-science-lab-educational-resources#munsell-renotation-data

References

• Color conversion equations

  • http://www.brucelindbloom.com/index.html?ColorCalcHelp.html

• Methods used to generate this table

  • http://dx.doi.org/10.1016/j.cageo.2012.10.020

Examples

data(munsell)

Spectral Library of Munsell Colors

Description

The original database "SpectralReflectancesOf2007MunsellBookOfColorGlossy.txt" was provided by Paul Centore and downloaded July, 2020. Reflectance values for odd chroma have been interpolated from adjacent chips. See aqp/misc/utils/Munsell/ for the entire set of processing steps.

Munsell value typically ranges from 2-9, and chroma from 1-12. Ranges vary by hue. Run aqp:::.summarizeMunsellSpectraRanges() for a detailed listing by hue.

The original database contains the following description:

This file contains spectral reflectance measurements of X-Rite's 2007 Munsell Book of Color (Glossy Finish). The measurements were made in 2012 with a ColorMunki spectrophotometer. The first column is the Munsell name. The remaining columns give reflectance values for 380 nm to 730 nm, in steps of 10 nm. The reflectance is a value between 0 (indicating that no light at that wavelength is reflected) and 1 (indicating that all the light at that wavelength is reflected). Occasionally an entry is slightly greater than 1. The likely cause is random variability, and those entries can be adjusted to 1 with negligible loss. In all, 1485 colour samples were measured. Researchers are invited to analyze the data in this file.

Usage

data(munsell.spectra)

Format

A data frame with 89496 rows and 10 variables:

munsell

munsell color

hue

hue component

value

value component

chroma

chroma component

wavelength

wavelength (nm)

reflectance

reflectance

References

Centore, Paul. Colour Tools for Painters. https://www.munsellcolourscienceforpainters.com/.

Convert Munsell Color Notation to other Color Space Coordinates (sRGB and CIELAB)

Description

Color conversion based on a look-up table of common soil colors.

Usage

munsell2rgb(
  the_hue,
  the_value,
  the_chroma,
  alpha = 1,
  maxColorValue = 1,
  return_triplets = FALSE,
  returnLAB = FALSE
)

Arguments

Details

This function is vectorized without recycling: i.e. the length of each argument must be the same. Both functions will pad output with NA if there are any NA present in the inputs.

Neutral hues are approximated by greyscale shades ranging from 20\

Gley soil colors that are missing a chroma will not be correctly interpreted. Consider using a chroma of 1. Non-standard Munsell notation (e.g. '7.9YR 2.7/2.0') can be matched (nearest-neighbor, no interpolation) to the closest color within the munsell sRGB/CIELAB look-up table via getClosestMunsellChip(). A more accurate estimate of sRGB values from non-standard notation can be achieved with the munsellinterpol package.

See examples below.

Value

A vector of R colors is returned that is the same length as the input data. When return_triplets = TRUE and/or returnLAB = TRUE, then a data.frame (of sample length as input) is returned.

Note

Care should be taken when using the resulting sRGB values; they are close to their Munsell counterparts, but will vary based on your monitor and ambient lighting conditions. Also, the value used for maxColorValue will affect the brightness of the colors. Th default value (1) will usually give acceptable results, but can be adjusted to force the colors closer to what the user thinks they should look like.

Author(s)

D.E. Beaudette

References

• http://www.brucelindbloom.com/index.html?ColorCalcHelp.html

• https://www.munsellcolourscienceforpainters.com/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html

• https://www.rit.edu/science/munsell-color-lab

Examples

# neutral hues (N) can be defined with chroma of 0 or NA 
g <- expand.grid(hue = 'N', value = 2:8, chroma = 0, stringsAsFactors = FALSE)
(m <- munsell2rgb(g$hue, g$value, g$chroma))
soilPalette(m)

# back-transform to Munsell notation
col2Munsell(t(col2rgb(m)) / 255)


# basic example
d <- expand.grid(hue = '10YR', value = 2:8, chroma = 1:8, stringsAsFactors = FALSE)
d$color <- with(d, munsell2rgb(hue, value, chroma))

# similar to the 10YR color book page
plot(value ~ chroma, data = d, col = d$color, pch = 15, cex = 3, las = 1)

# multiple pages of hue:
hues <- c('2.5YR', '5YR', '7.5YR', '10YR')
d <- expand.grid(
  hue = hues, 
  value = c(2, 2.5, 3:8), 
  chroma = seq(2, 8, by = 2), 
  stringsAsFactors = FALSE
)
# convert Munsell -> sRGB
d$color <- with(d, munsell2rgb(hue, value, chroma))

# extract CIELAB coordinates
with(d, munsell2rgb(hue, value, chroma, returnLAB = TRUE))

# plot: note that we are setting panel order from red --> yellow
library(lattice)

xyplot(
  value ~ factor(chroma) | factor(hue, levels = hues),
  main = "Common Soil Colors", layout = c(4, 1), scales = list(alternating = 1),
  strip = strip.custom(bg = grey(0.85)),
  data = d, as.table = TRUE, subscripts = TRUE,
  xlab = 'Chroma', ylab = 'Value',
  panel = function(x, y, subscripts, ...) {
    panel.xyplot(x, y, pch = 15, cex = 4, col = d$color[subscripts])
  }
)


# convert a non-standard color to closest "chip" in `munsell` look-up table
getClosestMunsellChip('7.9YR 2.7/2.0', convertColors = FALSE)

# convert directly to hex notation of sRGB
getClosestMunsellChip('7.9YR 2.7/2.0')

Merge Munsell Hue, Value, Chroma converted to sRGB & CIELAB into a SoilProfileCollection

Description

Convert Munsell hue, value and chroma into sRGB (rgb_R, rgb_G, rgb_B) and CIELAB (lab_L, lab_A, lab_B) color coordinates using munsell2rgb. The converted values are stored in the horizons() slot unless as.spc is FALSE, in which case the results are combined with profile and horizon ID columns and returned as the data.frame subclass used by the SPC.

Usage

## S4 method for signature 'SoilProfileCollection'
munsell2spc(
  object,
  hue = "hue",
  value = "value",
  chroma = "chroma",
  .data = NULL,
  as.spc = TRUE
)

Arguments

Value

A SoilProfileCollection or data.frame-like object

See Also

parseMunsell rgb2munsell munsell2rgb

Examples

data(sp3)
depths(sp3) <- id ~ top + bottom

# inspect input data
horizons(sp3)[,c("hue","value","chroma")]

# do color conversions to sRGB and LAB, join into horizon data
sp3 <- munsell2spc(sp3)

# plot rgb "R" coordinate by horizon
plot(sp3, color = "rgb_R")

# plot lab "A" coordinate by horizon
plot(sp3, color = "lab_A")

# note that `lab_A` values do not exactly match the original `A` values
# this is because `lab_A` was computed from the (field determined) Munsell color notation,
# while `A` was directly measured in the lab by colorimeter
plot(sp3$A, sp3$lab_A, xlab = 'Measured', ylab = 'Converted from Field Observed Munsell')

Munsell Hue Position Reference

Description

Position data for the 40 standard Munsell hues (and neutral). Data include angular positions (compass-style, origin at ⁠[x = 0, y = 1]⁠, CW rotation) and Cartesian coordinates on the unit circle.

Usage

data(munsellHuePosition)

Format

An object of class data.frame with 41 rows and 4 columns.

References

Munsell book of color. 1976. Macbeth, a Division of Kollmorgen Corp., Baltimore, MD.

Transform a SPC (by profile) with a set of expressions

Description

mutate_profile() is a function used for transforming SoilProfileCollections. Each expression is applied to site or horizon level attributes of individual profiles. This distinguishes this function from transform, which is applied to all values in a collection, regardless of which profile they came from.

Usage

mutate_profile(object, ..., col_names = NULL, horizon_level = NULL)

mutate_profile_raw(object, expr, col_names = NULL, horizon_level = NULL)

Arguments

Details

If the length an expression's result matches the number of horizons, the result is stored as a horizon-level variable. If the result has length 1, it is stored as a site-level variable. In the ambiguous case where the first and last profile have only one horizon, the results are stored in the horizon slot by default. To force results into site slot use horizon_level = FALSE.

Value

A SoilProfileCollection.

Author(s)

Andrew G. Brown.

Examples

data(sp4)
depths(sp4) <- id ~ top + bottom

mutate_profile(sp4, clay_wtd_average = weighted.mean(clay, bottom - top))

data(jacobs2000)

set.seed(123)

## col_names allows for column names to be calculated
x <- mutate_profile(jacobs2000, bottom - top / sum(bottom - top),
                    col_names = paste0("relthk", floor(runif(1, 0, 100))))
x$relthk28

# mutate_profile_raw allows for lists of expressions to be evaluated
master_desgn <- c("O", "A", "E", "B", "C", "R", "L", "M")
thk_names <- paste0("thk_", master_desgn)

# calculate thickness for each horizon
x$thk <- x$bottom - x$top

## construct an arbitrary number of expressions using variable inputs
ops <- lapply(master_desgn, function(x) {
  substitute(sum(thk[grepl(PATTERN, name)], na.rm = TRUE), list(PATTERN = x))
})
names(ops) <- thk_names

# do mutation
y <- mutate_profile_raw(x, ops)

site(y)[c(idname(y), thk_names)]

Get names of columns in site and horizons table

Description

Get names of columns in site and horizons table of a SoilProfileCollection.

Usage

## S4 method for signature 'SoilProfileCollection'
names(x)

Arguments

Get the number of horizons in a SoilProfileCollection

Description

Get the number of horizons in a SoilProfileCollection

Usage

## S4 method for signature 'SoilProfileCollection'
nrow(x)

Arguments

Example Output from soilDB::fetchOSD()

Description

An example SoilProfileCollection object created by soilDB::fetchOSD(), derived from the Cecil, Appling, and Bonneau Official Series Descriptions.

Usage

data(osd)

Format

A SoilProfileCollection

Lattice Panel Function for Soil Profiles

Description

Panel function for plotting grouped soil property data, along with upper and lower estimates of uncertainty.

This function can be used to replace panel.superpose when plotting depth function data. When requested, contributing fraction data are printed using colors the same color as corresponding depth function lines unless a single color value is given via cf.col.

This function is not able to apply transformations (typically log = 10) applied in the scales argument to xyplot to upper/lower bounds. These will have to be manually applied. See examples.

Usage

panel.depth_function(
  x,
  y,
  id,
  upper = NA,
  lower = NA,
  subscripts = NULL,
  groups = NULL,
  sync.colors = FALSE,
  cf = NA,
  cf.col = NA,
  cf.interval = 20,
  ...
)

prepanel.depth_function(
  x,
  y,
  upper = NA,
  lower = NA,
  subscripts,
  groups = NULL,
  ...
)

Arguments

Author(s)

D.E. Beaudette

See Also

sp1, slice, slab

Examples

library(lattice)
data(sp1)

# 1. plotting mechanism for step-functions derived from soil profile data
xyplot(
  cbind(top, bottom) ~ prop,
  data = sp1,
  id = sp1$id,
  panel = panel.depth_function,
  ylim = c(250, -10),
  scales = list(y = list(tick.number = 10)),
  xlab = 'Property',
  ylab = 'Depth (cm)',
  main = 'panel.depth_function() demo'
)

# 1.1 include groups argument to leverage lattice styling framework
sp1$group <- factor(sp1$group, labels = c('Group 1', 'Group2'))

xyplot(
  cbind(top, bottom) ~ prop,
  groups = group,
  data = sp1,
  id = sp1$id,
  panel = panel.depth_function,
  ylim = c(250, -10),
  scales = list(y = list(tick.number = 10)),
  xlab = 'Property',
  ylab = 'Depth (cm)',
  main = 'panel.depth_function() demo',
  auto.key = list(
    columns = 2,
    points = FALSE,
    lines = TRUE
  ),
  par.settings = list(superpose.line = list(col = c(
    'Orange', 'RoyalBlue'
  )))
)


# more complex examples, using step functions with grouped data
# better looking figures with less customization via tactile package
if(requireNamespace('tactile')) {
  
  library(data.table)
  library(lattice)
  library(tactile)
  
  # example data
  data(sp6)
  
  # a single profile
  x <- sp6[1:5, ]
  
  # wide -> long format
  x.long <- data.table::melt(
    data.table::data.table(x), 
    id.vars = c('id', 'top', 'bottom'), 
    measure.vars = c('sand', 'silt', 'clay')
  )
  
  # (optional) convert back to data.frame
  x.long <- as.data.frame(x.long)
  
  # three variables sharing a common axis
  # factor levels set by melt()
  xyplot(
    cbind(top, bottom) ~ value | id,
    groups = variable,
    data = x.long,
    id = x.long$id,
    ylim = c(200, -5), xlim = c(10, 60),
    scales = list(alternating = 1, y = list(tick.number = 10)),
    par.settings = tactile.theme(superpose.line = list(lwd = 2)),
    xlab = 'Sand, Silt, Clay (%)',
    ylab = 'Depth (cm)',
    panel = panel.depth_function,
    auto.key = list(columns = 3, lines = TRUE, points = FALSE),
    asp = 1.5
  )
  
  
  # all profiles
  x <- sp6
  
  # wide -> long format
  x.long <- data.table::melt(
    data.table::data.table(x), 
    id.vars = c('id', 'top', 'bottom'), 
    measure.vars = c('sand', 'silt', 'clay')
  )
  
  # (optional) convert back to data.frame
  x.long <- as.data.frame(x.long)
  
  # three variables sharing a common axis
  # factor levels set by melt()
  xyplot(
    cbind(top, bottom) ~ value | id,
    groups = variable,
    data = x.long,
    id = x.long$id,
    ylim = c(200, -5), xlim = c(0, 70),
    scales = list(alternating = 1, y = list(tick.number = 10)),
    par.settings = tactile.theme(superpose.line = list(lwd = 2)),
    xlab = 'Sand, Silt, Clay (%)',
    ylab = 'Depth (cm)',
    panel = panel.depth_function,
    auto.key = list(columns = 3, lines = TRUE, points = FALSE),
    as.table = TRUE
  )
  
  xyplot(
    cbind(top, bottom) ~ value,
    groups = variable,
    data = x.long,
    id = x.long$id,
    ylim = c(200, -5), xlim = c(0, 70),
    scales = list(alternating = 1, y = list(tick.number = 10)),
    par.settings = tactile.theme(superpose.line = list(lwd = 2)),
    xlab = 'Sand, Silt, Clay (%)',
    ylab = 'Depth (cm)',
    panel = panel.depth_function,
    auto.key = list(columns = 3, lines = TRUE, points = FALSE),
    as.table = TRUE
  )
  
  xyplot(
    cbind(top, bottom) ~ value | variable,
    groups = variable,
    data = x.long,
    id = x.long$id,
    ylim = c(200, -5), xlim = c(0, 70),
    scales = list(alternating = 1, y = list(tick.number = 10)),
    par.settings = tactile.theme(superpose.line = list(lwd = 2)),
    xlab = 'Sand, Silt, Clay (%)',
    ylab = 'Depth (cm)',
    panel = panel.depth_function,
    auto.key = list(columns = 3, lines = TRUE, points = FALSE),
    as.table = TRUE
  )
  
  
  xyplot(
    cbind(top, bottom) ~ value | variable,
    data = x.long,
    id = x.long$id,
    ylim = c(200, -5), xlim = c(0, 70),
    scales = list(alternating = 1, y = list(tick.number = 10)),
    par.settings = tactile.theme(superpose.line = list(lwd = 2)),
    xlab = 'Sand, Silt, Clay (%)',
    ylab = 'Depth (cm)',
    panel = panel.depth_function,
    auto.key = list(columns = 3, lines = TRUE, points = FALSE),
    as.table = TRUE
  )
  
  
}

Parse Munsell Color Notation

Description

Split Munsell color notation into "hue", "value", and "chroma", with optional conversion to sRGB hex notation, sRGB coordinates, and CIELAB coordinates. Conversion is performed by munsell2rgb().

Usage

parseMunsell(munsellColor, convertColors = TRUE, delim = NA, ...)

Arguments

Value

a data.frame object

Author(s)

P. Roudier and D.E. Beaudette

Examples

# just sRGB
parseMunsell("10YR 3/5", return_triplets = TRUE)

# sRGB + CIELAB (D65 illuminant)
parseMunsell("10YR 3/5", return_triplets = TRUE, returnLAB = TRUE)

# CIELAB only
parseMunsell("10YR 3/5", return_triplets = FALSE, returnLAB = TRUE)

# neutral hue
# note chroma encoded as '0'
parseMunsell('N 3/', convertColors = FALSE)