Component Interpretation Summary via SDA

D.E. Beaudette
2016-08-17
This document is based on soilDB version 1.8-1.

You probably need to install some packages before this will work. Have a look at this tutorial to get setup.

Lets get some interpretations that are relevant to MLRA 18 (Sierra Nevada Foothills) soils so that we can compare all components named Auburn and Dunstone. Just for fun lets throw in a very different soil from MLRA 17 (San Joaquin Valley); Hanford.

Note that SDA has a 100,000 row limit, so we have to be a little clever when writing the queries. In this case, letting the database filter out NULL fuzzy ratings and the reasons for ratings.

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

# beware, there are hard limits (10k rows) on what can be returned by SDA
q <- "SELECT component.cokey, compname, mrulename, interplr, interplrc
FROM component INNER JOIN cointerp ON component.cokey = cointerp.cokey
WHERE compname IN ('Auburn', 'Dunstone', 'Hanford')
AND seqnum = 0
AND mrulename IN ('ENG - Construction Materials; Topsoil', 
'ENG - Sewage Lagoons', 'ENG - Septic Tank Absorption Fields', 
'ENG - Unpaved Local Roads and Streets', 
'AGR - California Revised Storie Index (CA)', 
'AGR - Pesticide Loss Potential-Leaching')
AND interplr IS NOT NULL;"

# query and check
x <- SDA_query(q)
head(x)

OK, so how do the fuzzy numbers (interplr values) for these components compare? Box and whisker plots are the simplest, but sometimes density plots are helpful for viewing the complete distribution. Looks like Auburn and Dunstone are pretty similar, at least in terms of these interpretations.

# compute number of rules and soils
n.soils <- length(unique(x$compname))
n.rules <- length(unique(x$mrulename))

# compare population with box-whisker plot
bwplot(compname ~ interplr | mrulename, data=x, layout=c(1,n.rules))

densityplot(~ interplr | mrulename, groups=compname, data=x, layout=c(1,n.rules), auto.key=list(points=FALSE, lines=TRUE, columns=n.soils), scales=list(y=list(relation='free')))

Those figures were useful, sometimes it is nice to see the median values for each interpretation (rows) and component (columns).

s <- ddply(x, c('mrulename', 'compname'), summarize, 
           low=quantile(interplr, probs=0.05, na.rm=TRUE), 
           rv=quantile(interplr, probs=0.5, na.rm=TRUE), 
           high=quantile(interplr, probs=0.95, na.rm=TRUE))

knitr::kable(dcast(s, mrulename ~ compname, value.var = 'rv'), caption = "Median Fuzzy Ratings")

What about the categorical ratings? Note that the kable function from the knitr package makes nice HTML tables for us.

knitr::kable(xtabs(~ interplrc + compname , data=x, subset= mrulename == 'ENG - Construction Materials; Topsoil'), caption="ENG - Construction Materials; Topsoil")

knitr::kable(xtabs(~ interplrc + compname , data=x, subset= mrulename == 'AGR - California Revised Storie Index (CA)'), caption = "AGR - California Revised Storie Index (CA)'")

knitr::kable(xtabs(~ interplrc + compname , data=x, subset= mrulename == 'ENG - Septic Tank Absorption Fields'), caption = "ENG - Septic Tank Absorption Fields")

knitr::kable(xtabs(~ interplrc + compname , data=x, subset= mrulename == 'ENG - Sewage Lagoons'), caption = "ENG - Sewage Lagoons")

library(cluster)
library(ape)

# set list of soil series (component names)
soil.list <- c('Pardee', 'Yolo', 'Capay', 'Aiken', 'Amador', 'Pentz', 'Sobrante',
'Argonaut', 'Toomes', 'Jocal', 'Holland', 'Auburn', 'Dunstone', 
'Hanford', 'Redding', 'Columbia', 'San Joaquin', 'Fresno')

# set list of relevant interpretations
interp.list <- c('ENG - Construction Materials; Topsoil', 
'ENG - Sewage Lagoons', 'ENG - Septic Tank Absorption Fields', 
'ENG - Unpaved Local Roads and Streets', 
'AGR - California Revised Storie Index (CA)')

# compose query
q <- paste0("SELECT compname, mrulename, AVG(interplr) as interplr_mean
FROM component INNER JOIN cointerp ON component.cokey = cointerp.cokey
WHERE compname IN ", format_SQL_in_statement(soil.list), "
AND seqnum = 0
AND mrulename IN ", format_SQL_in_statement(interp.list), "
AND interplr IS NOT NULL
GROUP BY compname, mrulename;")

# send query
x <- SDA_query(q)

# reshape long -> wide
x.wide <- dcast(x, compname ~ mrulename, value.var = 'interplr_mean')
knitr::kable(x.wide, digits = 3, caption="Mean Fuzzy Ratings for Select Soil Series")

# create distance matrix
d <- daisy(x.wide[, -1])

# cluster via divisive hierachical method
h <- as.hclust(diana(d))

# transfer compname labels and convert to 'ape' class for plotting
h$labels <- x.wide$compname
h <- as.phylo(h)

# plot as dendrogram
par(mar=c(1,1,3,1))
plot(h)
title('Component Similarity via Select Interpretation Fuzzy Values')