if (FALSE) {
library(tactile)
library(lattice)
library(Hmisc)
library(sp)
# check the data out:
data(ca630)
str(ca630)
# note that pedon_key is the link between the two tables
# make a copy of the horizon data
ca <- ca630$lab
# promote to a SoilProfileCollection class object
depths(ca) <- pedon_key ~ hzn_top + hzn_bot
# add site data, based on pedon_key
site(ca) <- ca630$site
# ID data missing coordinates: '|' is a logical OR
(missing.coords.idx <- which(is.na(ca$lat) | is.na(ca$lon)))
# remove missing coordinates by safely subsetting
if(length(missing.coords.idx) > 0)
ca <- ca[-missing.coords.idx, ]
# register spatial data
initSpatial(ca) <- ~ lon + lat
# assign a coordinate reference system
prj(ca) <- 'EPSG:4269'
# check the result
print(ca)
# aggregate %BS 7 for all profiles into 1 cm slices
a <- slab(ca, fm= ~ bs_7)
# plot median & IQR by 1 cm slice
xyplot(
top ~ p.q50,
data = a,
lower=a$p.q25,
upper=a$p.q75,
alpha=0.5,
ylim=c(160,-5),
scales = list(alternating = 1, y = list(tick.num = 7)),
panel = panel.depth_function,
prepanel = prepanel.depth_function,
ylab='Depth (cm)', xlab='Base Saturation at pH 7',
par.settings = tactile.theme(superpose.line = list(col = 'black', lwd = 2))
)
# aggregate %BS at pH 8.2 for all profiles by MLRA, along 1 cm slices
# note that mlra is stored in @site
a <- slab(ca, mlra ~ bs_8.2)
# keep only MLRA 18 and 22
a <- subset(a, subset=mlra %in% c('18', '22'))
# plot median & IQR by 1 cm slice, using different colors for each MLRA
xyplot(
top ~ p.q50,
groups = factor(mlra),
data = a,
lower=a$p.q25,
upper=a$p.q75,
alpha=0.25,
sync.colors = TRUE,
ylim=c(160,-5),
scales = list(alternating = 1, y = list(tick.num = 7)),
panel = panel.depth_function,
prepanel = prepanel.depth_function,
ylab='Depth (cm)', xlab='Base Saturation at pH 7',
par.settings = tactile.theme(superpose.line = list(lwd = 2)),
auto.key = list(lines = TRUE, points = FALSE, columns = 2)
)
# Extract the 2nd horizon from all profiles as SPDF
ca.2 <- ca[, 2]
# subset profiles 1 through 10
ca.1.to.10 <- ca[1:10, ]
# basic plot method: profile plot
par(mar = c(0, 0, 3, 1))
plotSPC(ca.1.to.10, name='hzn_desgn', color = 'CEC7')
}
Run the code above in your browser using DataLab