Generating Soil Component Relation Diagrams

Details pending. A recent discussion on the topic.

Simple Example

# load libraries
library(soilDB)
library(sharpshootR)
library(plyr)
library(igraph)
library(RColorBrewer)

# load sample data associated with the amador soil series
data(amador)

# map unit keys, component names, component percentages
head(amador)
##    mukey   compname comppct_r
## 1 461845     amador        45
## 2 461845  gillender        40
## 3 461845     pardee         3
## 4 461845     peters         3
## 5 461845 ranchoseco         3
## 6 461845      vleck         3
# convert into adjacency matrix, based on no. times component co-occur
m.1 <- component.adj.matrix(amador, method='occurrence')
print(m.1)
##            amador auburn bellota exchequer gillender hornitos pardee pentz peters ranchoseco vleck
## amador          0      2       1         4         2        4      2     7      4          2     2
## auburn          0      0       0         0         0        0      0     2      2          0     0
## bellota         0      0       0         0         0        0      0     1      0          0     0
## exchequer       0      0       0         0         0        4      0     4      0          0     0
## gillender       0      0       0         0         0        0      2     0      2          2     2
## hornitos        0      0       0         0         0        0      0     4      0          0     0
## pardee          0      0       0         0         0        0      0     0      2          2     2
## pentz           0      0       0         0         0        0      0     0      2          0     0
## peters          0      0       0         0         0        0      0     0      0          2     2
## ranchoseco      0      0       0         0         0        0      0     0      0          0     2
## vleck           0      0       0         0         0        0      0     0      0          0     0
## attr(,"method")
## [1] "occurrence"
# convert into adjacency matrix, weighted by component percent
m.2 <- component.adj.matrix(amador)
print(round(m.2, 2))
##            amador auburn bellota exchequer gillender hornitos pardee pentz peters ranchoseco vleck
## amador          0      0    0.01      0.03       0.2     0.03   0.01  0.04   0.02       0.01  0.01
## auburn          0      0    0.00      0.00       0.0     0.00   0.00  0.09   0.21       0.00  0.00
## bellota         0      0    0.00      0.00       0.0     0.00   0.00  0.22   0.00       0.00  0.00
## exchequer       0      0    0.00      0.00       0.0     1.00   0.00  0.61   0.00       0.00  0.00
## gillender       0      0    0.00      0.00       0.0     0.00   0.08  0.00   0.07       0.08  0.08
## hornitos        0      0    0.00      0.00       0.0     0.00   0.00  0.61   0.00       0.00  0.00
## pardee          0      0    0.00      0.00       0.0     0.00   0.00  0.00   0.79       1.00  1.00
## pentz           0      0    0.00      0.00       0.0     0.00   0.00  0.00   0.07       0.00  0.00
## peters          0      0    0.00      0.00       0.0     0.00   0.00  0.00   0.00       0.79  0.79
## ranchoseco      0      0    0.00      0.00       0.0     0.00   0.00  0.00   0.00       0.00  1.00
## vleck           0      0    0.00      0.00       0.0     0.00   0.00  0.00   0.00       0.00  0.00
## attr(,"standardization")
## [1] "max"
## attr(,"metric")
## [1] "jaccard"
## attr(,"method")
## [1] "community.matrix"
# compare two methods
par(mfcol=c(1,2), mar=c(0,0,1,0))
plotSoilRelationGraph(m.1, s='amador')
title('Occurence')
plotSoilRelationGraph(m.2, s='amador')
title('Community Matrix')

plot of chunk libraries-and-sample

# compute adjacency matrix via community matrix method
m <- component.adj.matrix(amador)
# try a bunch of layouts
par(mfrow=c(2,3), mar=c(0,0,1,0))
plotSoilRelationGraph(m, s='amador', main='fruchterman.reingold (default)')
plotSoilRelationGraph(m, s='amador', g.layout=layout.grid, main='grid')
plotSoilRelationGraph(m, s='amador', g.layout=layout.circle, main='circle')
plotSoilRelationGraph(m, s='amador', g.layout=layout.mds, main='MDS')
plotSoilRelationGraph(m, s='amador', g.layout=layout.spring, main='spring')
plotSoilRelationGraph(m, s='amador', g.layout=layout.star, main='star')

plot of chunk layout-examples

# setup page
par(mfcol=c(1,2), mar=c(0,0,0,0))
# plot as network diagram, with Amador soil highlighted
plotSoilRelationGraph(m, s='amador')
# plot as dendrogram, with Amador soil highlighted
plotSoilRelationGraph(m, s='amador', plot.style='dendrogram')

plot of chunk dendrogram-representation

# tighten margins
par(mfcol=c(1,2), mar=c(0,0,1,0))
# plot as network diagram, with Amador soil highlighted
plotSoilRelationGraph(m, s='amador', main='Full Graph')
plotSoilRelationGraph(m, s='amador', spanning.tree='max', main='Maximum Spanning Tree')

plot of chunk mst

par(mar=c(1,1,1,1), mfcol=c(1,5))
for(i in seq(0, 1, length.out = 5)) {
  plotSoilRelationGraph(m, vertex.scaling.factor=3, del.edges=i)
  title(paste0('Edge weight < ', i, '-tile pruned'), cex.main=0.9, line=-1)
  box()
}

plot of chunk prune-edges

par(mar=c(1,1,1,1), mfcol=c(1,5))
for(i in seq(0, 1, length.out = 5)) {
  plotSoilRelationGraph(m, vertex.scaling.factor=3, del.edges=i, spanning.tree='max')
  title(main=paste0('Edge weight < ', i, '-tile pruned'), sub='Max Spanning Tree', cex.main=0.9, line=-1)
  box()
}

plot of chunk prune-edges-mst

NASIS via local database

# get legend/DMU/component data from NASIS
# in this case, CA630 data
d <- get_component_data_from_NASIS_db()

# normalize component names
d$compname <- tolower(d$compname)

# remove misc. areas components
d <- subset(d, compkind != 'miscellaneous area')

# remove some higher-taxa and non-soil components
d <- subset(d, ! compname %in% c('lithic haploxeralfs', 'rock outcrop', 'aquic haploxeralfs', 'ultic haploxeralfs', 'riparian'))

# select only those columns that we really need
d <- d[, c('dmudesc', 'compname', 'comppct_r', 'majcompflag')]

# how frequently is any given component a "major" component?
maj.comp.freq <- ddply(d, 'compname', summarize, 1 - prop.table(table(majcompflag))[1])
names(maj.comp.freq) <- c('compname', 'freq')

# extract map unit symbols from DMU description labels
d$musym <- gsub('CA630', '', d$dmudesc)

# check top couple of rows
head(d)
##     dmudesc           compname comppct_r majcompflag musym
## 3 CA6301011           orthents        50           1  1011
## 5 CA6301012           orthents        50           1  1012
## 6 CA6301030  aquic argixerolls        95           1  1030
## 7 CA6301030 typic haploxeralfs         5           0  1030
## 8 CA6301031  aquic argixerolls         5           0  1031
## 9 CA6301031 typic haploxeralfs        95           1  1031
# convert into adjacency matrix, weighted by component percentage
m <- component.adj.matrix(d, mu='musym', co='compname', wt='comppct_r')
# tighten margins
par(mar=c(0,0,2,0))
plotSoilRelationGraph(m, vertex.scaling.factor=1, main='CA630 Components')