Introduction

Setup R Envionment

If you have never used the aqp or soildb packages before, you will likely need to install them. This only needs to be done once.

# stable version from CRAN + dependencies
install.packages('aqp', dep=TRUE) 
install.packages('soilDB', dep=TRUE)
install.packages('sharpshootR', dep=TRUE)
install.packages('markovchain', dep=TRUE)

Once you have all of the R packages on which this document depends, it is a good idea to load them. R packages must be installed anytime you change versions of R (e.g., after an upgrade) and loaded anytime you want to access functions from within those packages.

library(markovchain)
library(aqp)
library(sharpshootR)
library(soilDB)
library(igraph)

Sample Data

While the methods outlined in this document can be applied to any collection of pedons, it is convenient to work with a standardized set of data. You can follow along with the analysis by copying code from the following blocks and running it in your R session. The sample data used in this document is based on soil profiles that have been correlated to the Loafercreek soil series from the Sierra Nevada Foothill Region of California. Note that the internal structure of the loafercreek data is identical to the structure returned by fetchNASIS() from the soilDB package. All horizon-level values are pulled from the pedon horizon table of the pedons being analyzed.

# load sample data from the soilDB package
data(loafercreek, package = 'soilDB')
pedons <- loafercreek
# plot profile sketches
par(mar=c(0,0,0,0))
plot(pedons, name='', print.id=FALSE, cex.names=0.75, axis.line.offset=-4)

Pedons correlated to the Loafercreek soil series.

Optional: Follow Along with Your Data

The following code block demonstrates how to pull data in using the fetchNASIS() function from the soilDB package.

# first load the desired data set within NASIS into your NASIS selected set
# then load data from the NASIS selected set into R
# note that the `nullFragsAreZero` argument converts NULL rock fragment class data into 0s
pedons <- fetchNASIS(nullFragsAreZero=TRUE)
# optionally subset the data by taxon name - enter your taxon name
pedons <- pedons[grep(pattern='ENTER_YOUR_TAXON_NAME', pedons$taxonname, ignore.case=TRUE), ]

Methods

sort(table(pedons$hzname), decreasing=TRUE)
## 
##    A  Bt1  Bt2   Cr  Bt3    R   BA   Oi  Crt  BCt   Bw 2Bt2 2Bt3  BAt   Bt  Bt4   A1   A2 2Bt4   2R 
##   97   93   88   49   47   40   31   24   20   16   10    8    8    8    8    8    7    7    6    6 
##    C  CBt  2Cr   BC 2Bt1 2Crt  ABt  Bw1  Bw2  2BC 2BCt 2Bt5  2CB 2CBt   AB   Ad   Ap    B  Bw3   Cd 
##    6    6    4    4    3    2    2    2    2    1    1    1    1    1    1    1    1    1    1    1 
## Cr/R   H1   Rt 
##    1    1    1

Constructing a graph of transitions from one horizon to the next (“transition probabilities”) can highlight those horizon designations that are most commonly used together. A follow-up tutorial on transition probability matrix interpretation is planned.

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

# horizon transition probabilities: row -> col transitions
(tp <- hzTransitionProbabilities(sp4, 'name'))
##     A A1 A2 AB       ABt        Bt       Bt1 Bt2 Bt3
## A   0  0  0  0 0.1111111 0.4444444 0.4444444   0   0
## A1  0  0  1  0 0.0000000 0.0000000 0.0000000   0   0
## A2  0  0  0  1 0.0000000 0.0000000 0.0000000   0   0
## AB  0  0  0  0 0.0000000 0.0000000 1.0000000   0   0
## ABt 0  0  0  0 0.0000000 0.0000000 1.0000000   0   0
## Bt  0  0  0  0 0.0000000 0.0000000 0.0000000   0   0
## Bt1 0  0  0  0 0.0000000 0.0000000 0.0000000   1   0
## Bt2 0  0  0  0 0.0000000 0.0000000 0.0000000   0   1
## Bt3 0  0  0  0 0.0000000 0.0000000 0.0000000   0   0
## attr(,"ties")
## [1] TRUE
# abuse sharpshootR code to display as graph
par(mar=c(0,0,0,0), mfcol=c(1,2))
plot(sp4)
plotSoilRelationGraph(tp, graph.mode = 'directed', edge.arrow.size=0.5, edge.scaling.factor=2, vertex.label.cex=0.75, vertex.label.family='sans')

Graph constructed from transition probabilities. Thicker lines denote more likely transitions.

More examples.

# try some other examples, seems logical
tp <- hzTransitionProbabilities(pedons, 'hzname')
# tp <- hzTransitionProbabilities(pedons, 'genhz')

# abuse sharpshootR code
# sometimes default layout doesn't look so good
par(mar=c(0,0,0,0), mfcol=c(1,2))
plotSoilRelationGraph(tp, graph.mode = 'directed', edge.arrow.size=0.5)
plotSoilRelationGraph(tp, graph.mode = 'directed', edge.arrow.size=0.5, g.layout = layout_with_lgl)