Introduction to SoilProfileCollection Objects

This is a very basic introduction to the SoilProfileCollection class object defined in the aqp package for R. The SoilProfileCollection class was designed to simplifiy the process of working with the collection of data associated with soil profiles: site-level data, horizon-level data, spatial data, diagnostic horizon data, metadata, etc. Examples listed below are meant to be copied/pasted from this document and interactively run within R. Comments (green text) briefly describe what the code in each line does. This document assumes a basic level of proficiency with R which can be gained by reviewing some of the material in tutorials like this. Further documentation on objects and functions from the aqp package can be accessed by typing help(aqp) (or more generally, ?function_name) at the R console.

Object Creation

SoilProfileCollection objects are typically created by “promoting” data.frame objects (rectangular tables of data) that contain at least three essential columns:

The data.frame should be pre-sorted according to the profile ID and horizon top boundary. Formula notation is used to define the columns used to promote a data.frame object:

idcolumn ~ hz_top_column + hz_bottom_column

In this tutorial we will use some sample data included with the aqp package, based on characterization data from 10 soils sampled on serpentinitic parent material as described in McGahan et al, 2009.

# load required packages, you may have to install these if missing:
# install.packages('aqp', dep=TRUE)
library(aqp)
library(Hmisc)
library(lattice)
library(MASS)

# load sample data set, a data.frame object with horizon-level data from 10 profiles
data(sp4)
str(sp4)
## 'data.frame':    30 obs. of  13 variables:
##  $ id         : chr  "colusa" "colusa" "colusa" "colusa" ...
##  $ name       : chr  "A" "ABt" "Bt1" "Bt2" ...
##  $ top        : int  0 3 8 30 0 9 0 4 13 0 ...
##  $ bottom     : int  3 8 30 42 9 34 4 13 40 6 ...
##  $ K          : num  0.3 0.2 0.1 0.1 0.2 0.3 0.2 0.6 0.8 0.4 ...
##  $ Mg         : num  25.7 23.7 23.2 44.3 21.9 18.9 12.1 12.1 17.7 16.4 ...
##  $ Ca         : num  9 5.6 1.9 0.3 4.4 4.5 1.4 7 4.4 24.1 ...
##  $ CEC_7      : num  23 21.4 23.7 43 18.8 27.5 23.7 18 20 31.1 ...
##  $ ex_Ca_to_Mg: num  0.35 0.23 0.08 0.01 0.2 0.2 0.58 0.51 0.25 1.47 ...
##  $ sand       : int  46 42 40 27 54 49 43 36 27 43 ...
##  $ silt       : int  33 31 28 18 20 18 55 49 45 42 ...
##  $ clay       : int  21 27 32 55 25 34 3 15 27 15 ...
##  $ CF         : num  0.12 0.27 0.27 0.16 0.55 0.84 0.5 0.75 0.67 0.02 ...

# optionally read about it...
# ?sp4

# upgrade to SoilProfileCollection
# 'id' is the name of the column containing the profile ID
# 'top' is the name of the column containing horizon upper boundaries
# 'bottom' is the name of the column containing horizon lower boundaries
depths(sp4) <- id ~ top + bottom

# check it out:
class(sp4) # class name
## [1] "SoilProfileCollection"
## attr(,"package")
## [1] "aqp"
print(sp4)
## Object of class SoilProfileCollection
## Number of profiles: 10
## Depth range: 16-49 cm
## 
## Horizon attributes:
##       id name top bottom   K   Mg  Ca CEC_7 ex_Ca_to_Mg sand silt clay   CF
## 1 colusa    A   0      3 0.3 25.7 9.0  23.0        0.35   46   33   21 0.12
## 2 colusa  ABt   3      8 0.2 23.7 5.6  21.4        0.23   42   31   27 0.27
## 3 colusa  Bt1   8     30 0.1 23.2 1.9  23.7        0.08   40   28   32 0.27
## 4 colusa  Bt2  30     42 0.1 44.3 0.3  43.0        0.01   27   18   55 0.16
## 5  glenn    A   0      9 0.2 21.9 4.4  18.8        0.20   54   20   25 0.55
## 6  glenn   Bt   9     34 0.3 18.9 4.5  27.5        0.20   49   18   34 0.84
## 
## Sampling site attributes:
##          id
## 1    colusa
## 2     glenn
## 3     kings
## 4  mariposa
## 5 mendocino
## 6      napa

Accessing, Setting, and Replacing Data

“Accessor” functions are used to extract specfic components from within SoilProfileCollection objects.

# methods for object inspection
idname(sp4) # self-explanitory
horizonDepths(sp4) # self-explanitory
depth_units(sp4) # defaults to 'cm'
metadata(sp4) # not much to start with
profile_id(sp4) # vector of profile IDs

# further inspection with common function overloads
length(sp4) # number of profiles in the collection
nrow(sp4) # number of horizons in the collection
names(sp4) # column names from site and horizon data
min(sp4) # shallowest profile depth in collection
max(sp4) # deepest profile depth in collection

Several accessor functions can also be used to define, modify, or replace components within SoilProfileCollection objects. For example,

horizons(sp4)

will extract horizon data and

# careful!
horizons(sp4) <- x

will replace horizon data. Further examples include:

# extraction of soil profile components
site(sp4)  # get or set site data
diagnostic_hz(sp4) # get or set diagnostic horizons
proj4string(sp4)
coordinates(sp4)

Horizon and Site Data

Typically, horizon and site data are the most important components of SoilProfileCollection objects. Both are internally stored as data.frame objects; with one or more rows (per profile ID) in the horizon table and one row (per profile ID) in the site table. Columns from either table can be accessed with the $ style notation of data.frames. New data can be assigned to either table in the same manner, as long as the length of the new data is either:

# assignment of new data to existing or new attributes
sp4$elevation <- rnorm(n=length(sp4), mean=1000, sd=150) # site-level, based on length of assigned data
sp4$thickness <- sp4$bottom - sp4$top # horizon-level

# extraction of specific attributes by name
sp4$clay # vector of clay content (horizon data)
##  [1] 21 27 32 55 25 34  3 15 27 32 25 31 33 13 21 23 15 17 12 19 14 14 22 25 40 51 67 24 25 32
sp4$elevation # vector of simulated elevation (site data)
##  [1] 1070.1 1135.0 1137.8  972.6 1125.4 1206.6 1153.0 1163.1  866.5 1295.1

# assign a single single value into horizon-level attributes
sp4$constant <- rep(1, times=nrow(sp4))

# promote horizon-level data to site-level data (when it makes sense to do so)
# note that this _moves_ the named column from horizon to site
site(sp4) <- ~ constant 

Horizon and site data can also be modified via extraction to data.frame followed by replacement (horizon data) or join (site data). Note that while this approach gives the most flexibility, it is also the most dangerous– replacement of horizon data with new data that don't exactly conform to the original sorting may corrupt your SoilProfileCollection.

# extract horizon data to data.frame
h <- horizons(sp4)

# add a new column and save back to original object
h$random.numbers <- rnorm(n=nrow(h), mean=0, sd=1)

# _replace_ original horizon data with modified version
# ! row-order should not be altered !
horizons(sp4) <- h

# extract site data to data.frame
s <- site(sp4)

# add a fake group to the site data
s$group <- rep(c('A', 'B'), length.out=nrow(s))

# join new site data with previous data: old data are _not_ replaced
site(sp4) <- s

# check:
sp4
## Object of class SoilProfileCollection
## Number of profiles: 10
## Depth range: 16-49 cm
## 
## Horizon attributes:
##       id name top bottom   K   Mg  Ca CEC_7 ex_Ca_to_Mg sand silt clay   CF thickness
## 1 colusa    A   0      3 0.3 25.7 9.0  23.0        0.35   46   33   21 0.12         3
## 2 colusa  ABt   3      8 0.2 23.7 5.6  21.4        0.23   42   31   27 0.27         5
## 3 colusa  Bt1   8     30 0.1 23.2 1.9  23.7        0.08   40   28   32 0.27        22
## 4 colusa  Bt2  30     42 0.1 44.3 0.3  43.0        0.01   27   18   55 0.16        12
## 5  glenn    A   0      9 0.2 21.9 4.4  18.8        0.20   54   20   25 0.55         9
## 6  glenn   Bt   9     34 0.3 18.9 4.5  27.5        0.20   49   18   34 0.84        25
##   random.numbers
## 1         0.6372
## 2        -0.3140
## 3         1.7968
## 4         0.5328
## 5        -0.7985
## 6         0.3723
## 
## Sampling site attributes:
##          id elevation constant group
## 1    colusa    1070.1        1     A
## 2     glenn    1135.0        1     B
## 3     kings    1137.8        1     A
## 4  mariposa     972.6        1     B
## 5 mendocino    1125.4        1     A
## 6      napa    1206.6        1     B

Diagnostic Horizons

Diagnostic horizons typically span several genetic horizons and may or may not be present in all profiles. To accomodate the wide range of possibilities, diagnostic horizon data are stored as a data.frame in “long format”: each row corresponds to a diagnostic horizon, identified with a column matching the ID column used to initialize the SoilProfileCollection object.

# manually create some diagnostic horizon data
# there is no restrictions on data format, as long as each row has an ID that exists within the collection
# be sure to use the ID column name that was used to initialize the SoilProfileCollection object
# check via: idname(sp4)
dh <- data.frame(id='colusa', kind='argillic', top=8, bottom=42, stringsAsFactors=FALSE)

# overwrite any existing diagnostic horizon data
diagnostic_hz(sp4) <- dh

# append to diagnostic horizon data
dh <- diagnostic_hz(sp4)
dh.new <- data.frame(id='napa', kind='argillic', top=6, bottom=20, stringsAsFactors=FALSE)

# overwrite existing diagnostic horizon data with appended data
diagnostic_hz(sp4) <- rbind(dh, dh.new)

Spatial Data

Spatial data can be explicitly stored within a SoilProfileCollection object and accessed with methods imported from the sp package. The use of sp classes (in this case SpatialPoints and SpatialPointsDataFrame objects) simplifies operations such as plotting spatial data, coordinate system transformations, and spatial queries.

# generate some fake coordinates as site level attributes
sp4$x <- rnorm(n=length(sp4), mean=354000, sd=100)
sp4$y <- rnorm(n=length(sp4), mean=4109533, sd=100)

# initialize spatial coordinates
coordinates(sp4) <- ~ x + y

# extract coordinates as matrix
coordinates(sp4)
##         x       y
## 1  354277 4109635
## 2  353990 4109565
## 3  354287 4109379
## 4  353921 4109633
## 5  354152 4109521
## 6  354050 4109669
## 7  353955 4109681
## 8  353963 4109519
## 9  354142 4109561
## 10 353939 4109451

# get/set spatial reference system using PROJ4 syntax
proj4string(sp4) <- '+proj=utm +zone=11 +datum=NAD83'
proj4string(sp4)
## [1] "+proj=utm +zone=11 +datum=NAD83"

# extract spatial data + site level attribtutes
# see ?SpatialPointsDataFrame for details
sp4.sp <- as(sp4, 'SpatialPointsDataFrame')

# plot the fake coordinates
plot(sp4.sp)
box()

plot of chunk spatial

Subsetting SoilProfileCollection Objects

SoilProfileCollection objects can be subset using the familiar [-style notiation used by matrix and data.frame objects, such that: spc[i, j] will return profiles identified by the integer vector i, and horizons identified by the integer vector j. Ommitting either index will result in all profiles (i ommitted) or all horizons (j ommitted). Typically, site-level attributes will be used as the subsetting criteria. Functions that return an index to matches (such as grep() or which()) provide the link between attributes and an index to matching profiles. Some examples:

# explicit string matching
idx <- which(sp4$group == 'A')

# numerical expressions
idx <- which(sp4$elevation < 1000)

# regular expression, matches any profile ID containing 'shasta'
idx <- grep('shasta', profile_id(sp4), ignore.case=TRUE)

# perform subset based on index
sp4[idx, ]

The subsetProfiles() function can be used to subset collections based on criteria from both horizon and site level attributes. Several examples are given below.

# matrix-style subsetting
sp4[1:2, ] # profiles 1--2 from the collection
sp4[, 1:2] # horizons 1--2 from each profile

# in the presence of spatial data, subsetting can either result in 
# 1. a new SoilProfileCollection (> 1 horizon / profile)
# 2. a new SpatialPointsDataFrame (1 horizon / profile)
class(sp4[1, ]) # profile 1 from the collection
class(sp4[, 1]) # horizon 1 from each profile

# filtering criteria from horizon and site level data
subsetProfiles(sp4, h='clay < 20', s='group == "A"')

Concatenation of SoilProfileCollection Objects

If the internal structures of several SoilProfileCollection objects are identical (e.g. same ID name, horizon, site, diagnostic horizon, and spatial attributes) they can be safely combined into a new object. This is typically only possible when re-combining subsets of the same original object. Note that concatenation will re-order all input SoilProfileCollection objects such that the resulting IDs are in alpha-numeric order.

# subset data into chunks
s1 <- sp4[1:2, ]
s2 <- sp4[4, ]
s3 <- sp4[c(6, 8, 9), ]

# re-combine: think "row-bind"
s <- rbind(s1, s2, s3)

Plotting SoilProfileCollection Objects

The plot() method for SoilProfileCollection objects generates sketches of profiles within the collection based on horizon boundaries, vertically aligned to an integer sequence from 1 to the number of profiles. Horizon names are automatically extracted from a horizon-level attribute name (if present), or via an alternate attributed given as an argument: name='column.name'. Horizon colors are automatically generated from the horizon-level attribute soil_color, or any other attribute of R-compatible color description given as an argument: color='column.name'. This function is highly customizable, therefore, it is prudent to consult help(plotSPC) from time to time. Soil colors in Munsell notation can be converted to R-compatible colors via munsell2rgb().

# plot uses base graphics, so it is simple to add low-level annotation layers
plot(sp4, name='name')

# inspect plotting area, very simple to overlay graphical elements
abline(v=1:length(sp4), lty=3, col='blue')

# profiles are centered at integers, from 1 to length(obj)
axis(1, line=0, at=1:10, cex.axis=0.75, font=4, col='blue', col.axis='blue', lwd=2)
mtext('profile axis', side=1, line=2, font=4, col='blue')

# y-axis is based on profile depths
axis(2, line=-1, at=pretty(1:max(sp4)), cex.axis=0.75, font=4, las=1, col='blue', col.axis='blue', lwd=2)
mtext('depth axis', side=2, line=1, font=4, col='blue')

plot of chunk plotting

Horizon-level attributes can be symbolized with color, in this case using the horizon-level attribute “clay”:

# plot again, this time using new colors
par(mar=c(0,0,3,0)) # tighter figure margins
plot(sp4, name='name', color='clay')

plot of chunk plotting-with-color

Horizon-level attributes that represent a volume fraction (e.g. coarse-fragment percentage) can be added to an existing figure:

par(mar=c(0,0,3,0)) # tighter figure margins
sp4$frag_pct <- sp4$CF * 100 # convert fraction -> percent
plot(sp4, name='name', col='frag_pct')
addVolumeFraction(sp4, 'frag_pct') # symbolize volume fraction data

plot of chunk plotting-vol-fraction

Annotation of depth-intervals can be accomplished using:

# generate some fake depth brackets:
fake.tops <- rep(0, times=length(sp4))
fake.bottoms <- runif(n=length(sp4), min=0, max=20)
par(mar=c(0,0,3,0)) # tighter figure margins
plot(sp4, name='name')
addBracket(fake.tops, fake.bottoms, col='red') # add depth brackets

plot of chunk plotting-depth-interval

Iterating Over Profiles in a Collection

The profileApply() function is an extension of the familiar apply() family of functions that operate on vectors (sapply and tapply), matrices (apply), and lists (lapply)– extended to SoilProfileCollection objects. The function named in the FUN argument is evaluated once for each profile in the collection, typicially returning a single value per profile. In this case, the ordering of the results would match the ordering of values in the site level attribute table.

# max() returns the depth of a soil profile
sp4$soil.depth <- profileApply(sp4, FUN=max)
# max() with additional argument give max depth to non-missing 'clay'
sp4$soil.depth.clay <- profileApply(sp4, FUN=max, v='clay')
# nrow() returns the number of horizons
sp4$n.hz <- profileApply(sp4, FUN=nrow)
# compute the mean clay content by profile using an inline function
sp4$mean.clay <- profileApply(sp4, FUN=function(i) mean(i$clay))
# estimate soil depth based on horizon designation
sp4$soil.depth <- profileApply(sp4, estimateSoilDepth, name='name', top='top', bottom='bottom')

When FUN returns a vector of the same length as the number of horizons in a profile, profileApply() can be used to create new horizon-level attributes. For example, the change in clay content with depth could be calculated via:

# save as horizon-level attribute
sp4$delta.clay <- profileApply(sp4, FUN=function(i) c(NA, diff(i$clay)))

# check results:
horizons(sp4)[1:6, c('id', 'top', 'bottom', 'clay', 'delta.clay')]
##       id top bottom clay delta.clay
## 1 colusa   0      3   21         NA
## 2 colusa   3      8   27          6
## 3 colusa   8     30   32          5
## 4 colusa  30     42   55         23
## 5  glenn   0      9   25         NA
## 6  glenn   9     34   34          9

More complex summaries can be generated by writing a custom function that is then called by profileApply(). Note that each profile is passed into this function and accessed via a temporary variable (i), which is a SoilProfileCollection object containing a single profile. A list of SoilProfileCollection objects returned from a custom function can be re-constituted into a single SoilProfileCollection object via rbind. See help(profileApply) for details.

# compute hz-thickness weighted mean exchangeable-Ca:Mg
wt.mean.ca.mg <- function(i) {
    # use horizon thickness as a weight
    thick <- i$bottom - i$top
    # compute the weighted mean, accounting for the possibility of missing data
    m <- wtd.mean(i$ex_Ca_to_Mg, weights=thick, na.rm=TRUE)
    return(m)
    }

# apply our custom function and save results as a site-level attribute
sp4$wt.mean.ca.to.mg <- profileApply(sp4, wt.mean.ca.mg)

We can now use our some of our new site-level attributes to order the profiles when plotting. In this case profiles are ordered based on the horizon-thickness weighted mean, exchageable Ca:Mg values. Horizons are colored by exchangeable Ca:Mg values.

# the result is an index of rank
new.order <- order(sp4$wt.mean.ca.to.mg)
# plot the data using our new order based on Ca:Mg
par(mar=c(4,0,3,0)) # tighten figure margins
plot(sp4, name='name', color='ex_Ca_to_Mg', plot.order=new.order)
# add an axis labeled with the sorting criteria
axis(1, at=1:length(sp4), labels=round(sp4$wt.mean.ca.to.mg, 3), cex.axis=0.75)
mtext(1, line=2.25, text='Horizon Thickness Weighted Mean Ex. Ca:Mg', cex=0.75)

plot of chunk profileApply-4

Slicing Soil Profile Collections

Collections of soil profiles can be sliced (or re-sampled) into regular depth-intervals with the slice() function. The slicing structure and variables of interest are defined via formula notation:

# slice select horizon-level attributes
seq ~ var.1 + var.2 + var.3 + ...
# slice all horizon-level attributes
seq ~ .

where seq is a sequence of integers (e.g. 0:15, c(5,10,15,20), etc.) and var.1 + var.2 + var.3 + ... are horizon-level attributes to slice. Both continuous and categorical variables can be named on the right-hand-side of the formula. The results returned by slice() is either a SoilProfileCollection, or data.frame when called with the optional argument just.the.data=TRUE. For example, to slice our sample data set into 1-cm intervals, from 0–15 cm:

# slice data
sliced <- slice(sp4, fm= 0:15 ~ sand + silt + clay + name + ex_Ca_to_Mg)
# check the result
class(sliced)
## [1] "SoilProfileCollection"
## attr(,"package")
## [1] "aqp"
# plot sliced data
par(mar=c(0,0,3,0)) # tighten figure margins
plot(sliced, name='name', color='ex_Ca_to_Mg')

plot of chunk slice

Once soil profile data have been sliced, it is simple to extract “chunks” of data by depth interval via matrix-style subsetting:

# slice from 0 to max depth in the collection
sliced <- slice(sp4, fm= 0:max(sp4) ~ sand + silt + clay + name + ex_Ca_to_Mg)
# extract all data over the range of 5--10 cm:
plot(sliced[, 5:10])
# extract all data over the range of 25--50 cm:
plot(sliced[, 25:50])
# extract all data over the range of 10--20 and 40--50 cm:
plot(sliced[, c(10:20, 40:50)])

Aggregating Soil Profile Collections Along Regular “Slabs”

Depth-wise summary of horizon-level attributes is performed with the slab() function. Profile grouping criterial and horizon attribute selection is parametrized via formula: either groups ~ var1 + var2 + var3 where named variables are aggregated within groups OR where named variables are aggregated across the entire collection ~ var1 + var2 + var3. The default summary function (slab.fun) computes the 5th, 25th, 50th, 75th, and 95th percentiles via Harrell-Davis quantile estimator.

The depth structure (“slabs”) over which summaries are computed is defined with the slab.structure argument using:

# aggregate a couple of the horizon-level attributes, 
# across the entire collection, 
# from 0--10 cm
# computing the mean value ignoring missing data
slab(sp4, fm= ~ sand + silt + clay, slab.structure=c(0,10), slab.fun=mean, na.rm=TRUE)
##   variable all.profiles value top bottom contributing_fraction
## 1     sand            1 47.63   0     10                     1
## 2     silt            1 31.15   0     10                     1
## 3     clay            1 21.11   0     10                     1

# again, this time within groups defined by a site-level attribute:
slab(sp4, fm= group ~ sand + silt + clay, slab.structure=c(0,10), slab.fun=mean, na.rm=TRUE)
##   variable group value top bottom contributing_fraction
## 1     sand     A 48.26   0     10                     1
## 2     silt     A 31.52   0     10                     1
## 3     clay     A 20.30   0     10                     1
## 4     sand     B 47.00   0     10                     1
## 5     silt     B 30.78   0     10                     1
## 6     clay     B 21.92   0     10                     1

# again, this time over several depth ranges
slab(sp4, fm= ~ sand + silt + clay, slab.structure=c(0,10,25,40), slab.fun=mean, na.rm=TRUE)
##   variable all.profiles value top bottom contributing_fraction
## 1     sand            1 47.63   0     10                1.0000
## 2     sand            1 42.39  10     25                0.8733
## 3     sand            1 32.15  25     40                0.5933
## 4     silt            1 31.15   0     10                1.0000
## 5     silt            1 29.41  10     25                0.8733
## 6     silt            1 31.35  25     40                0.5933
## 7     clay            1 21.11   0     10                1.0000
## 8     clay            1 28.11  10     25                0.8733
## 9     clay            1 36.27  25     40                0.5933

# again, this time along 1-cm slices, computing quantiles
agg <- slab(sp4, fm= ~ Mg + Ca + ex_Ca_to_Mg + CEC_7 + clay)

# see ?slab for details on the default aggregate function
head(agg)
##   variable all.profiles  p.q5 p.q25 p.q50 p.q75 p.q95 top bottom contributing_fraction
## 1       Mg            1 4.169 10.79 15.38 22.22 27.71   0      1                     1
## 2       Mg            1 4.169 10.79 15.38 22.22 27.71   1      2                     1
## 3       Mg            1 4.174 11.33 19.02 25.60 28.05   2      3                     1
## 4       Mg            1 4.255 12.53 20.31 25.90 32.56   3      4                     1
## 5       Mg            1 4.255 12.53 20.31 25.90 32.56   4      5                     1
## 6       Mg            1 4.255 12.53 20.46 27.00 32.88   5      6                     1

# plot median +/i bounds defined by the 25th and 75th percentiles
# this is lattice graphics, syntax is a little rough
xyplot(top ~ p.q50 | variable, data=agg, ylab='Depth',
             xlab='median bounded by 25th and 75th percentiles',
             lower=agg$p.q25, upper=agg$p.q75, ylim=c(42,-2),
             panel=panel.depth_function,
             alpha=0.25, sync.colors=TRUE,
             par.settings=list(superpose.line=list(col='RoyalBlue', lwd=2)),
             prepanel=prepanel.depth_function,
             cf=agg$contributing_fraction, cf.col='black', cf.interval=5, 
             layout=c(5,1), strip=strip.custom(bg=grey(0.8)),
             scales=list(x=list(tick.number=4, alternating=3, relation='free'))
             )

plot of chunk slab

Depth-wise aggregation can be useful for visual evaluation of multivariate similarity among groups of profiles.

# processing the "napa" and tehama profiles
idx <- which(profile_id(sp4) %in% c('napa', 'tehama'))
napa.and.tehama <- slab(sp4[idx, ], fm= ~ Mg + Ca + ex_Ca_to_Mg + CEC_7 + clay)

# combine with the collection-wide aggregate data
g <- make.groups(collection=agg, napa.and.tehama=napa.and.tehama)

# compare graphically:
xyplot(top ~ p.q50 | variable, groups=which, data=g, ylab='Depth',
             xlab='median bounded by 25th and 75th percentiles',
             lower=g$p.q25, upper=g$p.q75, ylim=c(42,-2),
             panel=panel.depth_function,
             alpha=0.25, sync.colors=TRUE, cf=g$contributing_fraction, cf.interval=10, 
             par.settings=list(superpose.line=list(col=c('RoyalBlue', 'Red4'), lwd=2, lty=c(1,2))),
             prepanel=prepanel.depth_function,
             layout=c(5,1), strip=strip.custom(bg=grey(0.8)),
             scales=list(x=list(tick.number=4, alternating=3, relation='free')),
             auto.key=list(columns=2, lines=TRUE, points=FALSE)
             )

plot of chunk slab-2

Pair-Wise Dissimilarity

Calculation of between-profile dissimilarity is performed using the profile_compare() function. Dissimilarity values depend on attributes selection (e.g. clay, CEC, pH , etc.), optional depth-weighting parameter (k), and a maximum depth of evaluation (max_d). This topic is further documented in a related tutorial, the function manual page (type help(profile_compare) in the R console), and this paper.

# eval dissimilarity:
# using Ex-Ca:Mg and CEC at pH 7
# with no depth-weighting (k=0)
# to a maximum depth of 40 cm
d <- profile_compare(sp4, vars=c('ex_Ca_to_Mg', 'CEC_7'), k=0, max_d=40)

# check distance matrix:
round(d, 1)
## Dissimilarities :
##                colusa glenn kings mariposa mendocino napa san benito shasta shasta-trinity
## glenn            13.5                                                                     
## kings            16.0  12.7                                                               
## mariposa          8.4  11.3  16.5                                                         
## mendocino        11.5   8.0  16.4     15.0                                                
## napa             30.4  24.1  29.4     29.2      21.6                                      
## san benito       25.7  20.6  26.3     28.2      15.8 18.0                                 
## shasta           17.2  13.3   8.7     17.6      17.1 33.7       22.2                      
## shasta-trinity    6.4  16.6  22.3      9.6      16.5 29.8       27.2   23.3               
## tehama           28.7  22.9  27.9     27.3      20.0  8.8       15.1   31.4           27.9
## 
## Metric :  mixed ;  Types = I, I 
## Number of objects : 10

# vizualize dissimilarity matrix via hierarchical clustering
h <- hclust(d)
plot(as.dendrogram(h))

plot of chunk dissimilarity

Object Metadata

# print metadata:
metadata(sp4)

# alter the depth unit metadata attribute
depth_units(sp4) <- 'inches' # units are really 'cm'

# more generic interface for adjusting metadata
md <- metadata(sp4) # save original metadata

# add columns
md$describer <- 'DGM'
md$date <- as.Date('2009-01-01')
md$citation <- 'McGahan, D.G., Southard, R.J, Claassen, V.P. 2009. Plant-Available Calcium Varies Widely in Soils on Serpentinite Landscapes. Soil Sci. Soc. Am. J. 73: 2087-2095.'

# re-assign user defined metadata to original object
metadata(sp4) <- md

# check:
metadata(sp4)

# fix depth units, set back to 'cm'
depth_units(sp4) <- 'cm'

Object Structure and Coercion

# check our work by viewing the internal structure
str(sp4)

# deconstruct SoilProfileCollection into a data.frame, with horizon+site data
as(sp4, 'data.frame')

# extraction of site + spatial data as SpatialPointsDataFrame
as(sp4, 'SpatialPointsDataFrame')

This document is based on aqp version 1.7.