Nothing
## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
out.width = "100%"
)
## ----setup--------------------------------------------------------------------
library(sarp.snowprofile.alignment)
## ----avgSPdisp, eval=FALSE, echo=TRUE-----------------------------------------
# ## compute the average profile (here with default settings)
# avgSP <- averageSP(SPgroup2)
# ## plot profile set and resulting average profile
# plot(SPgroup2, SortMethod = 'hs', xticklabels = 'originalIndices')
# plot(avgSP$avg)
## ----avgSPcalc, eval=TRUE, echo=FALSE, fig.asp=0.7, dpi = 300-----------------
avgSP <- averageSP(SPgroup2, n = 1, progressbar = FALSE)
layout(matrix(c(1, 2), ncol = 2), widths = c(2.5, 1))
par(cex = 0.3)
plot(SPgroup2, SortMethod = 'hs', xticklabels = 'originalIndices', ylim = c(0, 150), main = "individual profiles of example set 'SPgroup2'")
par(mar = c(5.1, 2.1, 2.1, 1.1))
plot(avgSP$avg, main = "average profile", ymax = 150)
## ---- echo=FALSE--------------------------------------------------------------
knitr::include_graphics("figures/averageSPdistributions.png")
## ---- echo=TRUE, eval=TRUE, fig.asp=0.5, dpi = 300----------------------------
## identify layer of interest: we need its row index
deepestSH_index <- min(findPWL(avgSP$avg, pwl_gtype = "SH")) # this is the deepest SH of the profile, in our specific case exactly what we want
## alternatively, you can (additionally) query for its date,
## or you can identify it manually by investigating the profile layers data frame:
# View(avgSP$avg$layers)
## backtrack this one layer
backtrackedlayers <- backtrackLayers(avgSP$avg, layer = deepestSH_index, profileSet = avgSP$set)
## analyze result:
str(backtrackedlayers)
## the backtrackedLayers object is a list of data frames, one data frame for each averaged layer (in our case 1);
## list elements are named by the height (cm) of the averaged layer
## ---- fig.asp=0.3, dpi = 300--------------------------------------------------
## compute whatever distribution you're interested: e.g., depth histogram
par(cex = 0.3)
hist(backtrackedlayers[[1]]$depth,
main = "Depth distribution of deepest SH layer", xlab = "Depth (cm)")
## ---- echo=TRUE, eval=TRUE----------------------------------------------------
## backtrack all layers by not providing a row index
backtrackedlayers <- backtrackLayers(avgSP$avg, profileSet = avgSP$set)
## -----------------------------------------------------------------------------
## identify proportion of underlying layers with good/fair/poor stability
## (here with RTA, but feel free to choose any other index or thresholds!)
transitional <- sapply(backtrackedlayers, function(bti) {
sum(bti$rta >= 0.6)/ length(bti$rta)
}) # proportion of layers fair stability
poor <- sapply(backtrackedlayers, function(bti) {
sum(bti$rta >= 0.8)/ length(bti$rta)
}) # proportion of layers poor stability
## ---- echo=TRUE, eval=TRUE, fig.asp=0.7, dpi = 300----------------------------
## visualize profile
layout(matrix(c(2, 1), 1, 2, byrow = TRUE), c(1.2, 1.8))
par(mar = c(9.1, 0, 2.1, 2.1), bg = "transparent", cex = 0.3)
plot(avgSP$avg, axes = FALSE, xlab = "")
axis(1, at = seq(5), labels = c("F", "4F", "1F", "P", "K"))
mtext("Hardness", side = 1, line = 5.5, cex = 0.3)
##
## stacked barplot for stability distributions
par(mar = c(9.1, 5.1, 2.1, 0))
x0 <- 0
ygrid <- as.numeric(names(backtrackedlayers)) # the names of the list define the height grid
cols <- c("gray90", "gray70", "gray20")
stackedhist <- rbind(data.frame(xleft = 0, xright = 1-transitional, ytop = ygrid,
ybottom = c(0, ygrid[1:(length(ygrid)-1)]), col = cols[1]),
data.frame(xleft = 1-transitional, xright = (1-transitional) + transitional,
ytop = ygrid, ybottom = c(0, ygrid[1:(length(ygrid)-1)]), col = cols[2]),
data.frame(xleft = 1-poor, xright = 1, ytop = ygrid,
ybottom = c(0, ygrid[1:(length(ygrid)-1)]), col = cols[3]))
plot(ygrid, xlim = c(0, 1), ylim = c(0, max(ygrid)), frame.plot = FALSE, type = 'n',
axes = FALSE, xlab = "", ylab = "Height (cm)")
mtext("Proportion of \nindividual profiles", side = 1, line = 5.5, cex = 0.3)
rect(xleft = x0+stackedhist$xleft, ybottom = stackedhist$ybottom,
xright = x0+stackedhist$xright, ytop = stackedhist$ytop,
col = stackedhist$col, border = NA)
xaxisticks <- c(0, 0.2, 0.5, 0.8, 1)
axis(1, at = xaxisticks, labels = rev(xaxisticks))
axis(2, at = pretty(ygrid))
abline(v = 0.2, lty = "dotted", col = "gray")
abline(v = 0.5, lty = "dotted", col = "gray")
abline(v = 0.8, lty = "dotted", col = "gray")
## ----timeseries---------------------------------------------------------------
## labeling of weak layers; again you can choose your own rules and thresholds!
SPspacetime <- snowprofileSet(lapply(SPspacetime, function(sp) {
labelPWL(sp, pwl_gtype = c("SH", "DH", "FC", "FCxr"), threshold_RTA = 0.8)
})) # label weak layers in each profile of the profile set 'SPspacetime'
## average along several days
avgSP <- averageSPalongSeason(SPspacetime)
## ---- dpi=300, fig.asp=0.7----------------------------------------------------
## explore the average timeseries object:
names(avgSP)
avgSP$call
avgSP$meta
## visualize the time series
par(cex = 0.3)
plot(avgSP$avgs, main = "Timeseries of average profile with median HS and median new snow amounts highligted", xlab = "Daily progression")
lines(avgSP$meta$date, avgSP$meta$hs_median)
lines(avgSP$meta$date, avgSP$meta$hs - avgSP$meta$thicknessPPDF_median, lty = "dashed")
## ---- fig.asp=0.7, dpi = 300--------------------------------------------------
## brief helper function for median vertical locations of specific layers (i.e., height or depth)
medianVLOC <- function(avgObj, pwldate, vloc = "Depth", pwlgt = c("SH", "DH"), date_range_earl = -5, draw = TRUE) {
mvl <- unname(do.call("c", lapply(avgObj$avgs, function(avg) {
median(avg$layers[, paste0("medianPredominant", vloc)][findPWL(avg, pwl_gtype = pwlgt, pwl_date = pwldate,
date_range_earlier = as.difftime(date_range_earl, units = "days"))])
})))
if (draw) {
if (vloc == "Depth") lines(avgObj$meta$date, avgObj$meta$hs_median - mvl, lty = "dotted", lwd = 2)
else if (vloc == "Height") lines(avgObj$meta$date, mvl, lty = "dotted", lwd = 2)
}
return(mvl)
}
## plot time series...
par(cex = 0.3)
plot(avgSP$avgs, main = "Time series with median depth of middle Nov 22 DH layer highlighted", xlab = "Daily progression")
## ... and apply above function to the Nov 22 weak layer
medianDepth_NOV22 <- medianVLOC(avgSP, "2018-11-23", vloc = "Depth")
## ---- fig.asp=0.7, dpi=300----------------------------------------------------
## rename the variable ppu_all to 'percentage' for subsequent plotting
avgSP$avgs <- snowprofileSet(lapply(avgSP$avgs, function(avg) {
avg$layers$percentage <- avg$layers$ppu_all
avg
}))
## overplot the grain type time series with the stability distribution
par(cex = 0.3)
plot(avgSP$avgs[avgSP$meta$date>="2018-09-20"], colAlpha = 0.5, main = "time series of average profile overplotted with stability distribution", xlab = "Daily progression")
plot(avgSP$avgs[avgSP$meta$date>="2018-09-20"], ColParam = "percentage", add = TRUE)
## ----echo=FALSE---------------------------------------------------------------
knitr::include_graphics("figures/averageSPseason_stability.png")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.