averageSP {sarp.snowprofile.alignment}

Average a group of snow profiles

Description

The functions dbaSP and averageSP implement Dynamic Time Warping Barycenter Averaging of snow profiles. The convenient wrapper averageSP takes care of choosing several appropriate initial conditions and picking the optimal end result (by minimizing the mean squared error between the average profile and the profile set). To pay appropriate attention to (thin) weak layers, weak layers need to be labeled in the profiles. You can either do that manually before calling this routine to suit your personal needs, or you can provide specific properties (in classifyPWLs) so that weak layers be labeled according to these properties by sarp.snowprofile::labelPWL. For more details, refer to the reference paper.

Usage

averageSP(
  SPx,
  n = 5,
  sm = summary(SPx),
  progressbar = requireNamespace("progress", quietly = TRUE),
  progressbar_pretext = NULL,
  classifyPWLs = list(pwl_gtype = c("SH", "DH")),
  classifyCRs = list(pwl_gtype = c("MFcr", "IF", "IFsc", "IFrc")),
  proportionPWL = 0.5,
  breakAtSim = 0.9,
  breakAfter = 2,
  verbose = FALSE,
  tz = "auto",
  ...
)

dbaSP(
  SPx,
  Avg,
  sm = summary(SPx),
  resamplingRate = 0.5,
  proportionPWL = 0.3,
  maxiter = 10,
  breakAtSim = 0.99,
  breakAfter = 1,
  plotChanges = FALSE,
  verbose = TRUE,
  tz = "auto",
  ...
)

Arguments

SPx

SPx a snowprofileSet object. Note that the profile layers need to contain a column called ⁠$layerOfInterest⁠ which classifies weak layers. While averageSP will label weak layers automatically if not done by the user beforehand, dbaSP won't do that but fail instead!; consider thinking about how you want to label weak layers, see Description, classifyPWLs below, and the references. Also note, that if you wish to average the rescaled profile set, do so manually before calling this function (see examples).

n

the number of initial conditions that will be used to run dbaSP; see also chooseICavg.

sm

a summary of SPx metadata

progressbar

should a progressbar be displayed (the larger n, the more meaningful the progressbar)

progressbar_pretext

a character string to be prepended to the progressbar (mainly used by higher level cluster function)

classifyPWLs

an argument list for a function call to sarp.snowprofile::findPWL which returns relevant PWLs for identifying initial conditions. Importantly, these arguments will also be used to label weak layers in the profiles, if these labels do not yet exist in the layers objects as column ⁠$layerOfInterest⁠. Check out the documentation of findPWL to familiarize yourself with your manifold options!

classifyCRs

an argument list for a function call to sarp.snowprofile::findPWL which returns relevant crusts for identifying initial conditions.

proportionPWL

decimal number that specifies the proportion required to average an ensemble of grain types as weak layer type. A value of 0.3, for example, means that layers will get averaged to a PWL type if 30% of the layers are of PWL type. Meaningful range is between ⁠[0.1, 0.5]⁠. Values larger than 0.5 get set to 0.5.

breakAtSim

stop iterations when simSP between the last average profiles is beyond that value. Can range between ⁠[0, 1]⁠. Default values differ between dbaSP and averageSP.

breakAfter

integer specifying how many values of simSP need to be above breakAtSim to stop iterating. Default values differ between dbaSP and averageSP.

verbose

print similarities between old and new average in between iterations?

tz

timezone of profiles; necessary for assigning the correct timezone to the average profile's ddate/bdate. Either 'auto' or a timezone known to ⁠[as.POSIXct]⁠.

...

alignment configurations which are passed on to dbaSP and then further to dtwSP. Note, that you can't provide rescale2refHS, which is always set to FALSE. If you wish to rescale the profiles, read the description of the SPx parameter and the examples.

Avg

the initial average snow profile: either a snowprofile object or an index to an initial average profile in SPx

resamplingRate

Resampling rate for a regular depth grid among the profiles

maxiter

maximum number of iterations

plotChanges

specify whether and how you want to plot the dba process: either FALSE, 'TRUE=='iterations'⁠, or ⁠'averages+last''

Details

Technical note: Since the layer characteristics of the average profile represent the median characteristics of the individual profiles, it can happen that ddates of the averaged layers are not in a monotonical order. That is, of course unphysical, but we specifically decided not to override these values to highlight these slight inconsistencies to users, so that they can decide how to deal with them. As a consequence, the function sarp.snowprofile::deriveDatetag does not work for these average profiles with ddate inconsistencies, but throws an error. The suggested workaround for this issue is to apply that function to all individual profiles before computing the average profile. This ensures that bdates or datetags are also included in the average profile.

For developers: Including new variables into the averaging/dba routines can be done easily by following commit #9f9e6f9

Value

A list of class avgSP that contains the fields

• ⁠$avg⁠: the resulting average profile

• ⁠$set⁠: the corresponding resampled profiles of the group

• ⁠$call⁠: (only with averageSP) the function call

• ⁠$prelabeledPWLs⁠: (only with averageSP) boolean scalar whether PWLs (or any other layers of interest) were prelabeled before this routine (TRUE) or labeled by this routine with the defaults specified in classifyPWLs (FALSE)

The profile layers of the average profile refer to the median properties of the predominant layers. For example, if you labeled all SH/DH layers as your 'layersOfInterest', and you find a SH or DH layer in the average profile, then it means that the predominant grain type is SH/DH (i.e., more profiles than specified in proportionPWL have that layer) and layer properties like hardness, p_unstable, etc refer to the median properties of these SH/DH layers. If you find a RG layer in your average profile, it means that most profiles have that RG layer and the layer properties refer to the median properties of all these RG layers. There are two exceptions to this rule, one for height/depth, and one for layer properties with the ending ⁠_all⁠, such as ppu_all:

• height and depth provide the vertical grid of the average profile, and for algorithmic reasons, this grid is not always equal to the actual median height or depth of the predominant layers. To account for that, two layer columns exist called medianPredominantHeight and medianPredominantDepth.

• Properties ending with ⁠_all⁠: For example, while ppu refers to the proportion of profiles, whose predominant layers are unstable (i.e., p_unstable >= 0.77), ppu_all refers to the the proportion of profiles, whose layers are unstable while taking into account all individual layers matched to this average layer (i.e., despite grain type, etc).

• Other layer properties specific to the average profile: distribution ranges between ⁠[0, 1]⁠ and specifies the proportion of profiles that contain the predominant layer described in the other properties.

Functions

• averageSP(): convenient wrapper function

• dbaSP(): DTW barycenter averaging of snow profiles (low level worker function)

Author(s)

fherla

References

Herla, F., Haegeli, P., and Mair, P. (2022). A data exploration tool for averaging and accessing large data sets of snow stratigraphy profiles useful for avalanche forecasting, The Cryosphere, 16(8), 3149–3162, https://doi.org/10.5194/tc-16-3149-2022

See Also

averageSPalongSeason

Examples

## EXAMPLES OF averageSP
this_example_runs_about_10s <- TRUE
if (!this_example_runs_about_10s) {  # exclude from cran checks

## compute the average profile of the demo object 'SPgroup'
## * by labeling SH/DH layers as weak layers,
##   - choosing 3 initial conditions with an above average number of weak layers
##   - in as many depth ranges as possible
## * and neglecting crusts for initial conditions

  avgList <- averageSP(SPgroup, n = 3,
                       classifyPWLs = list(pwl_gtype = c("SH", "DH")),
                       classifyCRs = NULL)

  opar <- par(mfrow = c(1, 2))
  plot(avgList$avg, ymax = max(summary(avgList$set)$hs))
  plot(avgList$set, SortMethod = "unsorted", xticklabels = "originalIndices")
  par(opar)


## compute the average profile of the demo object 'SPgroup'
## * by labeling SH/DH/FC/FCxr layers with an RTA threshold of 0.65 as weak layers,
## * otherwise as above

  SPx <- computeRTA(SPgroup)
  avgList <- averageSP(SPx, n = 3,
                       classifyPWLs = list(pwl_gtype = c("SH", "DH", "FC", "FCxr"),
                                           threshold_RTA = 0.65),
                       classifyCRs = NULL)

  opar <- par(mfrow = c(1, 2))
  plot(avgList$avg, ymax = max(summary(avgList$set)$hs))
  plot(avgList$set, SortMethod = "unsorted", xticklabels = "originalIndices")
  par(opar)

## compute the average profile of the other demo object 'SPgroup2', which
## contains more stability indices, such as SK38 or p_unstable
## * by labeling SH/DH/FC/FCxr layers that either
##   - have an SK38 below 0.95, *or*
##   - have a p_unstable above 0.77

  SPx <- snowprofileSet(SPgroup2)
  avgList <- averageSP(SPx,
                       classifyPWLs = list(pwl_gtype = c("SH", "DH", "FC", "FCxr"),
                                           threshold_SK38 = 0.95, threshold_PU = 0.77))

  opar <- par(mfrow = c(1, 2))
  plot(avgList$avg, ymax = max(summary(avgList$set)$hs))
  plot(avgList$set, SortMethod = "unsorted", xticklabels = "originalIndices")
  par(opar)

}



## EXAMPLES OF dbaSP
## either rescale profiles beforehand...
if (FALSE) {  # don't run in package check to save time
  SPx <- reScaleSampleSPx(SPgroup)$set          # rescale profiles
  SPx <- snowprofileSet(lapply(SPx, labelPWL))  # label PWLs
  DBA <- dbaSP(SPx, 5, plotChanges = TRUE)      # average profiles
}

## or use unscaled snow heights:
if (FALSE) {  # don't run in package check to save time
  SPx <- snowprofileSet(lapply(SPgroup, labelPWL))  # label PWLs
  DBA <- dbaSP(SPx, 5, plotChanges = TRUE)          # average profiles
}