# Average a group of snow profiles 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. ```r 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", n_cores = NULL, ... ) 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 sarp.snowprofile::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]`. - `n_cores`: number of nodes to create for a cluster using the parallel package to speed up calculations (default = NULL) - `...`: 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 sarp.snowprofile::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'` ## Returns 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. ## 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 monotonic order. That is, of course non-physical, 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 ## Functions * `averageSP()`: convenient wrapper function * `dbaSP()`: DTW barycenter averaging of snow profiles (low level worker function) ## Examples ```r ## 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 } ``` ## 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 ## Author(s) fherla

averageSP function