# setup environment
library(reshape)
library(lattice)
# load example data
data(sp1)
# test that mean of 1cm slotted property is equal to the
# hz-thickness weighted mean value of that property
sp1.sub <- subset(sp1, sub=id == 'P009')
hz.wt.mean <- with(sp1.sub,
sum((bottom - top) * prop) / sum(bottom - top)
)
a <- soil.slot(sp1.sub)
all.equal(mean(a$p.mean), hz.wt.mean)
#
# 1 standard usage, and plotting example
#
# slot at two different segment sizes
a <- soil.slot(sp1)
b <- soil.slot(sp1, seg_size=5)
# stack into long format
ab <- make.groups(a, b)
ab$which <- factor(ab$which, levels=c('a','b'),
labels=c('1-cm Interval', '5-cm Interval'))
# manually add mean +/- SD
ab$upper <- with(ab, p.mean+p.sd)
ab$lower <- with(ab, p.mean-p.sd)
# use mean +/- 1 SD
# custom plotting function for uncertainty viz.
xyplot(top ~ p.mean | which, data=ab, ylab='Depth',
xlab='mean bounded by +/- 1 SD',
lower=ab$lower, upper=ab$upper, ylim=c(250,-5), alpha=0.5,
panel=panel.depth_function,
prepanel=prepanel.depth_function,
layout=c(2,1), scales=list(x=list(alternating=1))
)
# use median and IQR
# custom plotting function for uncertainty viz.
xyplot(top ~ p.q50 | which, data=ab, ylab='Depth',
xlab='median bounded by 25th and 75th percentiles',
lower=ab$p.q25, upper=ab$p.q75, ylim=c(250,-5), alpha=0.5,
panel=panel.depth_function,
prepanel=prepanel.depth_function,
layout=c(2,1), scales=list(x=list(alternating=1))
)
#
# 1.1 try slotting categorical variables
#
# normalize horizon names:
sp1$name[grep('O', sp1$name)] <- 'O'
sp1$name[grep('A[0-9]', sp1$name)] <- 'A'
sp1$name[grep('AB', sp1$name, ignore.case=TRUE)] <- 'A'
sp1$name[grep('BA', sp1$name)] <- 'B'
sp1$name[grep('Bt', sp1$name)] <- 'B'
sp1$name[grep('Bw', sp1$name)] <- 'B'
sp1$name[grep('C', sp1$name)] <- 'C'
sp1$name[grep('R', sp1$name)] <- 'R'
# generate new data for testing soil.slot()
y <- with(sp1, data.frame(id=id, top=top, bottom=bottom, prop=name))
# convert name to a factor
y$prop <- factor(y$prop)
# fix factor levels
y$id <- factor(y$id)
# default slotting-- 1cm intervals
a <- soil.slot(y)
# reshape into long format for plotting
a.long <- melt(a, id.var=c('top','bottom'))
a.long$variable <- factor(a.long$variable, levels=c('O','A','B','C','R'))
# plot horizon type proportions
xyplot(top ~ value | variable, data=a.long, subset=value > 0,
ylim=c(150, -5), type=c('S','g'), horizontal=TRUE, layout=c(4,1), col=1 )
#
# 2. depth probability via contributing fraction
# note that this assumes that we are not missing data in 'prop'
# get around NA by making makeing a fake column filled with 1
# like this:
# sp1$prop <- 1
#
a <- soil.slot(sp1)
xyplot(top ~ contributing_fraction, data=a, ylim=c(250, -5), type='S', horizontal=TRUE, asp=4)
#
# 3.1 standard aggregation
#
a <- soil.slot(sp1)
# manually add mean +/- SD
a$upper <- with(a, p.mean+p.sd)
a$lower <- with(a, p.mean-p.sd)
# use custom plotting function for uncertainty viz.
xyplot(top ~ p.mean, data=a,
lower=a$lower, upper=a$upper, ylim=c(250,-5), alpha=0.5,
panel=panel.depth_function,
prepanel=prepanel.depth_function
)
#
# 3.2 standard aggregation, by grouping variable
#
require(plyr)
# try str(sp1) for hints
sp1$group <- with(sp1, factor(group))
a <- ddply(sp1, .(group), .fun=soil.slot)
# manually add mean +/- SD
a$upper <- with(a, p.mean+p.sd)
a$lower <- with(a, p.mean-p.sd)
# use custom plotting function for uncertainty viz.
xyplot(
top ~ p.mean, data=a, groups=group, subscripts=TRUE,
lower=a$lower, upper=a$upper, ylim=c(250,-5), alpha=0.5,
panel=panel.depth_function,
prepanel=prepanel.depth_function
)
#
# 3.3 use of weights
#
data(sp1)
# some fake weights
wts <- data.frame(id=levels(sp1$id), wt=c(3,1,2,1,2,1,4,1,2))
# merge wtih original data
g <- merge(sp1, wts, by='id')
# generate horizon mid points
g$mid <- with(g, (bottom + top) / 2)
# aggregate and add upper/lower intervals via SD
a <- soil.slot(g, use.wts=TRUE)
a$upper <- with(a, p.mean + p.sd)
a$lower <- with(a, p.mean - p.sd)
a$wt.upper <- with(a, p.wtmean + p.wtsd)
a$wt.lower <- with(a, p.wtmean - p.wtsd)
# check influence of weights
plot(mid ~ prop, data=g, ylim=c(240,0), cex=sqrt(g$wt), xlim=c(-5,60))
lines(top ~ p.mean, data=a, lwd=2)
lines(top ~ upper, data=a, lty=2)
lines(top ~ lower, data=a, lty=2)
lines(top ~ p.wtmean, data=a, col=2, lwd=2)
lines(top ~ wt.upper, data=a, col=2, lty=2)
lines(top ~ wt.lower, data=a, col=2, lty=2)
# annotate with explanation
legend('bottomright',
legend=c('wt = 1','wt = 2','wt = 3','wt = 4','un-weighted','weighted'),
pch=c(1,1,1,1,NA,NA), lwd=c(1,1,1,1,2,2), lty=c(NA,NA,NA,NA,1,1),
pt.cex=sqrt(c(1,2,3,4,1,1)), col=c(1,1,1,1,1,2), bty='n')
#
# try again with larger segment sizes
#
# aggregate and add upper/lower intervals via SD
a <- soil.slot(g, use.wts=TRUE, seg_size=10)
a$upper <- with(a, p.mean + p.sd)
a$lower <- with(a, p.mean - p.sd)
a$wt.upper <- with(a, p.wtmean + p.wtsd)
a$wt.lower <- with(a, p.wtmean - p.wtsd)
# convert to long format
library(reshape)
a.long <- melt(a, id.var=c('top','bottom'), measure.var=c('p.mean','p.wtmean','upper','lower','wt.upper','wt.lower'))
# red lines are weighted
# point symbols are sized proportional to their weights
xyplot(cbind(top, bottom) ~ value, groups=variable, data=a.long,
ylim=c(260,-10), ylab='Depth', xlab='Property',
par.settings=list(superpose.line=list(
col=c('black','red','black','black','red','red'),
lwd=c(2,2,1,1,1,1),
lty=c(1,1,2,2,2,2))),
panel=function(...) {
panel.points(g$prop, g$mid, cex=sqrt(g$wt), col=1)
panel.depth_function(...)
}
)Run the code above in your browser using DataLab