Introduction

This document demonstrates how to use the soilDB package to download KSSL data from SoilWeb. These data are from the September 2018 snapshot, and will be updated as future snapshots are released. Comparisons are made via graphical summaries of key soil properties with depth, using data structures and functions from the aqp package.

Installation

With a recent version of R (>= 2.15), it is possible to get all of the packages that this tutorial depends on via:

# stable packages from CRAN
install.packages('maps', dep=TRUE)
install.packages('plyr', dep=TRUE)
install.packages('reshape2', dep=TRUE)
install.packages('soilDB', dep=TRUE)

# latest versions from GitHub
devtools::install_github("ncss-tech/aqp", dependencies=FALSE, upgrade_dependencies=FALSE)
devtools::install_github("ncss-tech/soilDB", dependencies=FALSE, upgrade_dependencies=FALSE)

Quick Example: getting lab characterization (KSSL) and basic morphology (NASIS)

KSSL and NASIS data are from the September 2018 snapshot. Details pending.

library(soilDB)
library(plyr)
library(reshape2)

# get lab and morphologic data
s <- fetchKSSL(series='auburn', returnMorphologicData = TRUE)

# the result is a list, check it out
str(s, 2)
## List of 2
##  $ SPC  :Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##  $ morph:List of 4
##   ..$ phcolor    :'data.frame':  139 obs. of  7 variables:
##   ..$ phfrags    :'data.frame':  40 obs. of  10 variables:
##   ..$ phpores    :'data.frame':  95 obs. of  5 variables:
##   ..$ phstructure:'data.frame':  57 obs. of  7 variables:
# check out the "raw" morphologic data:
lapply(s$morph, head)
## $phcolor
##    phiid labsampnum colorpct colorhue colorvalue colorchroma colormoistst
## 1 820051   05N02182       NA     10YR          4           4          Dry
## 2 820051   05N02182       NA     10YR          3           3        Moist
## 3 820052   05N02183       NA    7.5YR          3           3        Moist
## 4 820052   05N02183       NA    7.5YR          4           6          Dry
## 5 820053   05N02184       NA      5YR          4           4          Dry
## 6 820053   05N02184       NA      5YR          3           3        Moist
## 
## $phfrags
##     phiid labsampnum fragvol                   fragkind fragsize_l fragsize_r fragsize_h fragshp
## 1  820051   05N02182      14 Metamorphic rock fragments          2          5         76      NA
## 2  820052   05N02183      10 Metamorphic rock fragments          2         20         76      NA
## 3  820053   05N02184      18 Metamorphic rock fragments          2          5         76      NA
## 4 1842112   10N04396       2       Greenstone fragments          2          3          5      NA
## 5 1842114   10N04397       4       Greenstone fragments          2         10         20      NA
## 6 1842114   10N04397       2       Greenstone fragments          2         10         20      NA
##    fraground          fraghard
## 1 Subrounded         Indurated
## 2 Subrounded         Indurated
## 3 Subangular         Indurated
## 4       <NA>   Weakly cemented
## 5       <NA> Strongly cemented
## 6       <NA>   Weakly cemented
## 
## $phpores
##     phiid labsampnum poreqty  poresize           poreshp
## 1  820051   05N02182     5.0      Fine      Interstitial
## 2  820051   05N02182     7.0    Medium           Tubular
## 3  820052   05N02183     5.0 Very fine Dendritic tubular
## 4  820052   05N02183     6.0    Medium           Tubular
## 5  820053   05N02184     5.0 Very fine Dendritic tubular
## 6 1842112   10N04396     0.8    Medium           Tubular
## 
## $phstructure
##     phiid labsampnum structgrade structsize        structtype structid structpartsto
## 1  820051   05N02182    Moderate       Fine          Granular       NA            NA
## 2  820052   05N02183    Moderate     Medium Subangular blocky       NA            NA
## 3  820053   05N02184    Moderate     Medium Subangular blocky       NA            NA
## 4 1842112   10N04396    Moderate       Fine          Granular       NA            NA
## 5 1842114   10N04397    Moderate     Medium Subangular blocky       NA            NA
## 6 1842116   10N04398    Moderate     Medium Subangular blocky       NA            NA
# extract pedons into SoilProfileCollection
pedons <- s$SPC

# simplify color data
# note that both horizon ID and color "percentage" must be specified
s.colors <- simplifyColorData(s$morph$phcolor, id.var = 'labsampnum', wt='colorpct')

# merge color data into SPC
h <- horizons(pedons)
h <- join(h, s.colors, by='labsampnum', type='left', match='first')
horizons(pedons) <- h

# check
par(mar=c(0,0,0,0))
plot(pedons, color='moist_soil_color', print.id=FALSE, name='hzn_desgn')

# simplify fragment data
s.frags <- simplfyFragmentData(s$morph$phfrags, id.var='labsampnum')

# merge fragment data into SPC
h <- horizons(pedons)
h <- join(h, s.frags, by='labsampnum', type='left', match='first')
horizons(pedons) <- h


# check
par(mar=c(0,0,3,0))
plot(pedons, color='total_frags_pct', print.id=FALSE, name='hzn_desgn')
addVolumeFraction(pedons, 'total_frags_pct', pch=1, cex.min=0.25, cex.max = 0.5)

Fragment Volume >= 100%

There are cases where errors in the source data (NASIS Pedon Horizon Fragments table) suggest coarse fragment volume > 100%. As of soilDB 2.0.2 and aqp 1.15.4, the simplfyFragmentData() function will report such errors and the addVolumeFraction() function will fail gracefully by truncating fragment volume at 100%.

Give it a try with an example (thanks to Brian Gardner for finding it).

# get lab and morphologic data
s <- fetchKSSL(series='kettenbach', returnMorphologicData = TRUE)

# extract pedons into SoilProfileCollection
pedons <- s$SPC

# simplify fragment data
s.frags <- simplfyFragmentData(s$morph$phfrags, id.var='labsampnum')
## Warning: fragment volume >= 100%
## labsampnum:
## 82P02200
## 82P02204
# merge fragment data into SPC
h <- horizons(pedons)
h <- join(h, s.frags, by='labsampnum', type='left', match='first')
horizons(pedons) <- h

# flag horizons with >= 100% total fragment volume
pedons$fragment_errors <- factor(pedons$total_frags_pct >= 100)

# label pedons with user pedon ID
# label horizons with total fragment %
par(mar=c(0,0,3,5))
plot(pedons, color='fragment_errors', name='total_frags_pct', label='pedon_id', cex.name=0.85)
addVolumeFraction(pedons, 'total_frags_pct', pch=1, cex.min=0.25, cex.max = 0.5)
## Warning: 160 >= 100, likely a data error
## Warning: 180 >= 100, likely a data error

Check the original fragment data, based on the warnings issued by simplfyFragmentData(). Looks like multiple versions of the same data were accidentally entered. Errors like this can be corrected in NASIS, but may require assistance from the regional staff depending on KSSL pedon ownership. fetchKSSL() is based on a snapshot of NASIS and LIMS data, therefore corrections may not be available until after the next release cycle.

# identify and print affected rows
idx <- which(s$morph$phfrags$labsampnum %in% c('82P02200', '82P02204'))
s$morph$phfrags[idx, ]
##      phiid labsampnum fragvol         fragkind fragsize_l fragsize_r fragsize_h fragshp  fraground
## 8   570064   82P02200      10             <NA>         76         NA        250    <NA>       <NA>
## 9   570064   82P02200      70             <NA>          2         NA         76    <NA>       <NA>
## 10 4384673   82P02200      10 Basalt fragments         75         NA        250 Nonflat Subangular
## 11 4384673   82P02200      70 Basalt fragments          2         NA         75 Nonflat Subangular
## 24  571809   82P02204      40             <NA>          2         NA         76    <NA>       <NA>
## 25  571809   82P02204      50             <NA>         76         NA        250    <NA>       <NA>
## 26 4384681   82P02204      40 Basalt fragments          2         NA         75 Nonflat Subangular
## 27 4384681   82P02204      50 Basalt fragments         75         NA        250 Nonflat Subangular
##     fraghard
## 8       <NA>
## 9       <NA>
## 10 Indurated
## 11 Indurated
## 24      <NA>
## 25      <NA>
## 26 Indurated
## 27 Indurated

Fetch Data Associated with a Set of IDs

Fetching KSSL data for a set of IDs requires some additional "helper" functions. The following example fetches data according to a vector of "pedlabsampnum" IDs, extracts the horizon color data, filters missing data, and combines the results into a single table.

First, paste the following code into your script or console to initialize our helper functions.

# get data associated with a single ID
# no records -> NULL returned
getPedons <- function(x) {
  suppressMessages(fetchKSSL(pedlabsampnum=x, returnMorphologicData = TRUE))
}

# extract a morphologic table from each record set
# NULL data are filtered out
# converted to data.frame
extractMorphTable <- function(x, table='phcolor') {
  m <- lapply(x, function(i) i[['morph']][[table]])
  # index pedons with data
  idx <- which(sapply(m, length) > 0) 
  if(length(idx) < 1)
    return(NULL)
  # filter non-NULL and convert to DF
  m <- ldply(m[idx])
  return(m)
}

Next, use the getPedons() helper function to fetch lab and morphologic data according to a vector of pedlabsampnum IDs. Note that the results are a list of lists. Select morphologic data are extracted and combined with the extractMorphTable() function.

# pedons indexed by "pedlabsampnum" ID
pls <- c("04N0610", "04N0611", "04N0612", "04N0613")
# iterate over IDs and get data
res <- lapply(pls, getPedons)

# check: OK
str(res, 2)
## List of 4
##  $ :List of 2
##   ..$ SPC  :Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##   ..$ morph:List of 4
##  $ :List of 2
##   ..$ SPC  :Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##   ..$ morph:List of 4
##  $ :List of 2
##   ..$ SPC  :Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##   ..$ morph:List of 4
##  $ :List of 2
##   ..$ SPC  :Formal class 'SoilProfileCollection' [package "aqp"] with 7 slots
##   ..$ morph:List of 4
# extract phcolor data from list of records
# result is a data.frame with all non-NULL rows
phcolor <- extractMorphTable(res, table='phcolor')

# check: OK
head(phcolor)
##     phiid labsampnum colorpct colorhue colorvalue colorchroma colormoistst
## 1 1463128   04N03382      100     10YR          5           4        Moist
## 2 1463138   04N03383      100      5YR          4           6        Moist
## 3 1463139   04N03384      100      5YR          4           6        Moist
## 4 1463140   04N03385       33    2.5YR          4           8        Moist
## 5 1463140   04N03385       33    7.5YR          5           8        Moist
## 6 1463140   04N03385       34      10R          4           8        Moist

Finally, extract the SoilProfileCollection objects from res and combine into a single object. We can now join our combined color data with the combined pedon data.

# extract pedons from list
pedons <- lapply(res, function(i) i$SPC)
# combine into a single SPC object
pedons <- do.call('rbind', pedons)

# check: looks good
plot(pedons, color='estimated_ph_h2o', name='hzn_desgn')

# simplify combined color data
colors <- simplifyColorData(phcolor, id.var = 'labsampnum', wt='colorpct')

# join combined color data to SPC
h <- horizons(pedons)
h <- join(h, colors, by='labsampnum', type='left', match='first')
horizons(pedons) <- h

# check: looks good
plot(pedons, color='moist_soil_color', name='hzn_desgn')

Example Application

Data can be queried by "taxonname"" (typically a soil series), using a geographic bounding-box in WGS84 referenced coordinates, or by MLRA. Taxonname matching is case-insensitive and based on the most current taxonname in NASIS. See ?fetchKSSL() or this GitHub repository for details.

In this example, we are downloading KSSL data for the musick, chaix, and holland soils; commonly found on granitic rocks of the Sierra Nevada Mountain region of California.

# load libraries
library(soilDB)
library(plyr)
library(lattice)
library(maps)

# define plotting style
tps <- list(superpose.line=list(col=c('RoyalBlue', 'DarkRed', 'DarkGreen'), lwd=2))

# fetch KSSL data by fuzzy-matching of series name
musick <- fetchKSSL('musick')
holland <- fetchKSSL('holland')
chaix <- fetchKSSL('chaix')

# generate a basemap of northern California, with county boundaries
map('county', 'California', xlim=c(-123.25, -118.75), ylim=c(36.5, 41))
# add long/lat axes
map.axes()
# add locations of musick
points(y ~ x, data=site(musick), pch=21, bg='RoyalBlue')
# add locations of holland
points(y ~ x, data=site(holland), pch=21, bg='DarkRed')
# add locations of chaix
points(y ~ x, data=site(chaix), pch=21, bg='DarkGreen')
# add a simple legend
legend('topright', pch=21, pt.bg=c('RoyalBlue', 'DarkRed', 'DarkGreen'), legend=c('musick', 'holland', 'chaix'), bty='n')

Profile Plots

Generate "thematic" profile sketches by coloring horizons according to clay content.

par(mar=c(0,1,5,1))
plot(musick, name='hzn_desgn', color='clay')

plot(holland, name='hzn_desgn', color='clay')

plot(chaix, name='hzn_desgn', color='clay')

Combining Data

Since the results from fetchKSSL() function always returns data in the same format, it is possible to stack the results of several KSSL queries using rbind().

# check "correlated as" names from each object
table(musick$taxonname)
## 
## Musick MUSICK 
##     10      3
table(holland$taxonname)
## 
## Holland HOLLAND 
##      15       6
table(chaix$taxonname)
## 
## Chaix CHAIX 
##     8     4
# since there are multiple permuations of each, normalize soil series name by replacement
musick$taxonname <- 'musick'
holland$taxonname <- 'holl