data(ca630)
str(ca630)
# 2. promote to SoilProfileCollection
# combine site+horizon data into single DF
ca <- join(ca630$lab, ca630$site, type='inner')
# promote to SoilProfileCollection
depths(ca) <- pedon_key ~ hzn_top + hzn_bot
# extract site data
site(ca) <- ~ mlra + ssa + lon + lat + cntrl_depth_to_top + cntrl_depth_to_bot + sampled_taxon_name
# extract spatial data as SpatialPoints
coordinates(ca) <- ~ lon + lat
# assign CRS data
proj4string(ca) <- '+proj=latlong +datum=NAD83'
# check the result
ca
# 3. aggregate \%BS 7 for all profiles into 1 cm slices
a <- slab(ca, fm= ~ bs_7)
# 4. plot median & IQR by 1 cm slice
xyplot(
top ~ p.q50, data=a, lower=a$p.q25, upper=a$p.q75,
ylim=c(160,-5), alpha=0.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=list(superpose.line=list(col='black', lwd=2))
)
# 4. 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)
xyplot(
top ~ p.q50, groups=factor(mlra, levels=c('18','22')) , data=a, lower=a$p.q25, upper=a$p.q75,
ylim=c(160,-5), alpha=0.5, scales=list(y=list(tick.num=7, alternating=3), x=list(alternating=1)),
panel=panel.depth_function, prepanel=prepanel.depth_function,
ylab='Depth (cm)', xlab='Base Saturation at pH 8.2',
par.settings=list(superpose.line=list(col=c('black','blue'), lty=c(1,2), lwd=2)),
auto.key=list(columns=2, title='MLRA', points=FALSE, lines=TRUE)
)
# 5. safely compute hz-thickness weighted mean CEC (pH 7)
head(lab.agg.cec_7 <- ddply(ca630$lab, .(pedon_key),
.fun=summarise, CEC_7=wtd.mean(bs_7, weights=hzn_bot-hzn_top)))
# 6. extract a SPDF with horizon data along a slice at 25 cm
s.25 <- slice(ca, fm=25 ~ bs_7 + CEC7 + ex_acid)
spplot(s.25, zcol=c('bs_7','CEC7','ex_acid'))
# note that the ordering is preserved:
all.equal(s.25$pedon_key, profile_id(ca))
# 6.1 extract a data.frame with horizon data along a slices c(10,20,50)
s.multiple <- slice(ca, fm=c(10,20,50) ~ bs_7 + CEC7 + ex_acid)
# 7. Extract the 2nd horizon from all profiles as SPDF
ca.2 <- ca[, 2]
# 8. Extract the all horizons from profiles 1:10 as data.frame
ca.1.to.10 <- ca[1:10, ]
# basic plot method: profile plot
plot(ca.1.to.10, name='hzn_desgn')
Run the code above in your browser using DataLab