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)
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)
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), ]
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')
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)