Generating Soil Component Relation Diagrams
D.E. Beaudette
2017-06-14
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 |
|---|---|---|
| 461845 | amador | 45 |
| 461845 | gillender | 40 |
| 461845 | pardee | 3 |
| 461845 | peters | 3 |
| 461845 | ranchoseco | 3 |
| 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')
# compare vertex scaling methods
par(mfcol=c(1,2), mar=c(0,0,1,0))
plotSoilRelationGraph(m.2, s='amador')
title('Degree')
plotSoilRelationGraph(m.2, s='amador', vertex.scaling.method = 'distance')
title('Distance from "Amador"')
# 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_on_grid, main='grid')
plotSoilRelationGraph(m, s='amador', g.layout=layout_in_circle, main='circle')
plotSoilRelationGraph(m, s='amador', g.layout=layout_with_mds, main='MDS')
plotSoilRelationGraph(m, s='amador', g.layout=layout_with_gem, main='GEM')
plotSoilRelationGraph(m, s='amador', g.layout=layout_as_star, main='star')
# 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')
# 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')
# 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='min', main='Minumum Spanning Tree')
par(mar=c(1,1,1,1), mfcol=c(1,5))
for(i in seq(0, 1, length.out = 5)) {
plotSoilRelationGraph(m, s='amador', vertex.scaling.factor=3, del.edges=i)
title(paste0('Edge weight < ', i, '-tile pruned'), cex.main=0.9, line=-1)
box()
}
par(mar=c(1,1,1,1), mfcol=c(1,5))
for(i in seq(0, 1, length.out = 5)) {
plotSoilRelationGraph(m, s='amador', 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()
}
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')
plotSoilRelationGraph(m, spanning.tree='max', vertex.scaling.factor=1, edge.scaling.factor=8, s='amador', main='CA630 Components\n(max spanning tree)')
plotSoilRelationGraph(m, spanning.tree=0.75, vertex.scaling.factor=1, edge.scaling.factor=8, s='amador', main='CA630 Components\n(max spanning tree + edges with weights > 75th percentile)')
plotSoilRelationGraph(m, vertex.scaling.factor=1, edge.scaling.factor=8, s='loafercreek', del.edges=0.25, main='CA630 Components\n(edge weights < 25th percentile pruned)')
plotSoilRelationGraph(m, spanning.tree='max', vertex.scaling.factor=1, edge.scaling.factor=8, s='amador', main='CA630 Components\n(max spanning tree and maj. component flag label adjustment)', vertex.label.cex=sqrt(maj.comp.freq$freq+0.5))
SSURGO via SDA
# select by survey area(s)
q <- "SELECT
component.mukey, comppct_r, lower(compname) as compname
FROM legend
INNER JOIN mapunit ON mapunit.lkey = legend.lkey
LEFT OUTER JOIN component ON component.mukey = mapunit.mukey
WHERE legend.areasymbol IN ('CA654')
AND compkind IN ('Series', 'Taxadjunct')"
# run query, process results, and return as data.frame object
res <- SDA_query(q)
# convert into adj. matrix using component percentage weights
m <- component.adj.matrix(res)# tighten margins
par(mar=c(0,0,2,0))
plotSoilRelationGraph(m, vertex.scaling.factor=1, edge.scaling.factor=8, main='CA654 Components\nFull Graph')
plotSoilRelationGraph(m, vertex.scaling.factor=1, edge.scaling.factor=8, spanning.tree='max', main='CA654 Components\nMaximum Spanning Tree')
plotSoilRelationGraph(m, vertex.scaling.factor=1, edge.scaling.factor=8, spanning.tree=0.75, main='CA654 Components\nMaximum Spanning Tree + edges with weight > 75th percentile')
This document is based on aqp version 1.11, soilDB version 1.8.5, and sharpshootR version 1.3.