1. Default algorithm
1.1 Averaging of profiles
The default averaging algorithm can be accessed conveniently through
the function averageSP and is implemented exactly as
described in the reference paper. You can also consult the documentation
(?averageSP) for more details.
## 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)
Note: To speed up the computation in this vignette, the above figure was
not computed with default settings in the background, but with a limited
number of initial conditions. The displayed result might not represent
the optimal average profile, so recompute it yourself with the default
settings displayed above, or experiment with other settings!
It might be worthwhile to reflect on which layers you consider weak
layers. With the function argument classifyPWLs you can
make the averaging algorithm consider your perspective on weak layers
and make them more likely to get included in the resulting profile. One
approach (the default) could be to consider all surface/depth hoar
layers as weak layers. Another approach could be to consider all
persistent grain types that have one particular stability index above or
below a certain threshold. Or you can also consider multiple different
stability indices. Make sure to consult the documentation for some more
details and tips.
1.2 Retrieval of distributions: Backtracking of layers
Since all layers of all individual profiles are matched against the layers of the average profile, we can backtrack each layer of the average profile to obtain all underlying layers that were matched against it. This allows us to compute various layer distributions or profile distributions as visualized in Figure 2 of the reference paper:
To compute an underlying distribution for the average example profile
from above, we have to call the function backtrackLayers.
You can either backtrack all layers of the average profile, or select
specific layers of your interest.
Backtracking individual layers
In the next example, let’s only backtrack the surface hoar layer that is buried just below the new snow in the average example profile:
## 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)
#> List of 1
#> $ 111:'data.frame': 5 obs. of 29 variables:
#> ..$ station : int [1:5] 77561 80858 75907 84723 82519
#> ..$ datetime : POSIXct[1:5], format: "2019-01-21 17:00:00" "2019-01-21 17:00:00" ...
#> ..$ elev : int [1:5] 1979 1859 1845 2037 1871
#> ..$ band : chr [1:5] "TL" "TL" "TL" "TL" ...
#> ..$ angle : num [1:5] 0 0 0 0 0
#> ..$ aspect : num [1:5] 0 0 0 0 0
#> ..$ hs : num [1:5] 111 114 142 131 118
#> ..$ station_id : chr [1:5] "VIR077561" "VIR080858" "VIR075907" "VIR084723" ...
#> ..$ date : Date[1:5], format: "2019-01-21" "2019-01-21" ...
#> ..$ depth : num [1:5] 6.5 8.5 9.5 9.5 8.5
#> ..$ height : num [1:5] 104 106 133 122 110
#> ..$ thickness : num [1:5] 0.5 0.5 0.5 0.5 0.5
#> ..$ ddate : POSIXct[1:5], format: "2019-01-17 15:30:00" "2019-01-16 09:30:00" ...
#> ..$ gtype : Factor w/ 9 levels "DF","DH","FC",..: 7 2 3 7 3
#> ..$ hardness : num [1:5] 1.36 1 1 1.12 1
#> ..$ gsize : num [1:5] 1.24 1.55 1.32 3.65 1.31
#> ..$ density : num [1:5] 171 113 100 140 95
#> ..$ sphericity : num [1:5] 0 0.09 0.27 0 0.25
#> ..$ v_strain_rate : num [1:5] -1 -0.2 -0.5 -0.4 -0.5
#> ..$ rta : num [1:5] 0.69 0.7 0.58 0.99 0.67
#> ..$ sk38 : num [1:5] 6 6 6 6 6
#> ..$ shear_strength : num [1:5] 0.31 0.22 0.17 0.32 0.16
#> ..$ slab_rho : num [1:5] 120 107 129 96 112
#> ..$ slab_rhogs : num [1:5] 241 211 280 199 241
#> ..$ crit_cut_length: num [1:5] 0.1 0.08 0.14 0.03 0.11
#> ..$ p_unstable : num [1:5] 0.64 0.6 0.78 0.63 0.7
#> ..$ bdate : POSIXct[1:5], format: "2019-01-18 10:30:00" "2019-01-18 05:15:00" ...
#> ..$ datetag : Date[1:5], format: "2019-01-18" "2019-01-18" ...
#> ..$ layerOfInterest: logi [1:5] TRUE TRUE FALSE TRUE FALSE
## 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## 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)")
Backtracking all layers simultaneously
This is most meaningful if you’re interested in the stability distributions of all layers.
## 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 stabilityCurrently, visualizing the stability distributions for the entire profile needs to be done by hand because no convenient wrapper function exists yet. In the following you get an idea of how that can be done:
## 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")
Admittedly, the distribution looks very pixelated, but that’s solely due
to the few number of underlying profiles (it’s 5 underlying profiles in
our case here). Grab a data set of yours and try it yourself with more
profiles!
