data(sp1)
depths(sp1) <- id ~ top + bottom
# scale properties within each profile
# scaled = (x - mean(x)) / sd(x)
sp1$d <- profileApply(sp1, FUN=function(x) round(scale(x$prop), 2))
plot(sp1, name='d')
# compute depth-wise differencing by profile
# note that our function expects that the column 'prop' exists
f <- function(x) { c(x$prop[1], diff(x$prop)) }
sp1$d <- profileApply(sp1, FUN=f)
plot(sp1, name='d')
# compute depth-wise cumulative sum by profile
# note the use of an anonymous function
sp1$d <- profileApply(sp1, FUN=function(x) cumsum(x$prop))
plot(sp1, name='d')
# compute profile-means, and save to @site
# there must be some data in @site for this to work
site(sp1) <- ~ group
sp1$mean_prop <- profileApply(sp1, FUN=function(x) mean(x$prop, na.rm=TRUE))
# re-plot using ranks defined by computed summaries (in @site)
plot(sp1, plot.order=rank(sp1$mean_prop))
##
## helper functions: these must be modified to suit your own data
##
# compute the weighted-mean of some property within a given diagnostic horizon
# note that this requires conditional eval of data that may contain NA
# see ?slab for details on the syntax
# note that function expects certain columns within 'x'
f.diag.wt.prop <- function(x, d.hz, prop) {
# extract diagnostic horizon data
d <- diagnostic_hz(x)
# subset to the requested diagnostic hz
d <- d[d$diag_kind == d.hz, ]
# if missing return NA
if(nrow(d) == 0)
return(NA)
# extract depths and check for missing
sv <- c(d$featdept, d$featdepb)
if(any(is.na(sv)))
return(NA)
# create formula from named property
fm <- as.formula(paste('~', prop))
# return just the (weighted) mean, accessed from @horizons
s <- slab(x, fm, seg_vect=sv)$p.mean
return(s)
}
# conditional eval of thickness of some diagnostic feature or horizon
# will return a vector of length(x), you can save to @site
f.diag.thickness <- function(x, d.hz) {
# extract diagnostic horizon data
d <- diagnostic_hz(x)
# subset to the requested diagnostic hz
d <- d[d$diag_kind == d.hz, ]
# if missing return NA
if(nrow(d) == 0)
return(NA)
# compute thickness
thick <- d$featdepb - d$featdept
return(thick)
}
# conditional eval of property within particle size control section
# makes assumptions about the SPC that is passed-in
f.psc.prop <- function(x, prop) {
# these are accessed from @site
sv <- c(x$psctopdepth, x$pscbotdepth)
# test for missing PCS data
if(any(is.na(sv)))
return(NA)
# this should never happen... unless someone made a mistake
# check to make sure that the lower PSC boundary is shallower than the depth
if(sv[2] > max(x))
return(NA)
# create formula from named property
fm <- as.formula(paste('~', prop))
# return just the (weighted) mean, accessed from @horizons
s <- slab(x, fm, seg_vect=sv)$p.mean
return(s)
}
Run the code above in your browser using DataLab