Details pending. A recent discussion on the topic.
# 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')
# 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')
# 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')
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()
}
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()
}
# 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')