sarp.snowprofile.alignment (version 2.0.2)

simSP: Similarity measure between snow profile pairs

Description

This function calculates a similarity measure for two snow profiles that have been aligned onto the same height grid (either through DTW or resampling). If one profile contains more layers than the other one, the layers with a non-matched height represent missing layers and will be treated accordingly. The similarity measure is compatible with top-down alignments and is symmetric with respect to its inputs, i.e. simSP(P1, P2) == simSP(P2, P1). Several different approaches of computing the measure have been implemented by now, see Details below.

Usage

simSP(
  ref,
  qw,
  gtype_distMat_simSP = sim2dist(grainSimilarity_evaluate(triag = FALSE)),
  simType = "HerlaEtAl2021",
  nonMatchedSim = 0,
  nonMatchedThickness = 10,
  verbose = FALSE,
  returnDF = FALSE,
  apply_scalingFactor = FALSE,
  simWeights = c(gtype = 1/3, hardness = 1/3, stability = 1/3),
  ...
)

Value

Either a scalar similarity between [0, 1] with 1 referring to the two profiles being identical, or (if returnDF is TRUE) a list with the elements $sim and $simDF.

Arguments

ref

snowprofile object 1

qw

snowprofile object 2 (matched layers need to be on the same height grid of ref)

gtype_distMat_simSP

a distance matrix that stores distance information of grain types (Be careful to convert similarities, as in grainSimilarity_evaluate, into dissimilarities with sim2dist.)

simType

the similarity measure can be computed in several different ways (of sophistication). See Details section. Possible choices

• simple

• HerlaEtAl2021 (= simple2)

• tsa_WLdetection & rta_WLdetection

• layerwise & rta_scaling

• remotesensing

• wsum_scaled

nonMatchedSim

sets the similarity value of non-matched layers [0, 1]. "indifference" = 0.5, penalty < 0.5. Note that dtwSP sets the same value and overrides the default value in this function!

nonMatchedThickness

If NA, every unique non-matched layer (i.e., contiguous resampled layers with identical properties) contributes to the overall similarity by 1 x nonMatchedSim. In that case, 5cm of non-matched new snow has the same effect on the overall similarity as 50cm of non-matched new snow. To make the effect of non-matched layers dependent on the layer thickness, provide a positive number to nonMatchedThickness. For nonMatchedThickness = 10, every 10cm of a unique non-matched layer contribute to the overall similarity by 1 x nonMatchedSim. So, 50cm of non-matched new snow would contribute 5 times stronger than 5cm of non-matched new snow. Note that dtwSP sets the same value and overrides the default value in this function!

verbose

print similarities of different grain classes to console? default FALSE

returnDF

additionally return the similarities of the grain classes as data.frame (analogously to verbose); the return object then has the fields $sim and $simDF

apply_scalingFactor

Only applicable to types layerwise and wsum_scaled: TRUE or FALSE, see Details.

simWeights

a numeric vector with exact names that specifies the weights for the weighted averaging in wsum_scaled

...

not used, but necessary to absorb unused inputs from dtwSP

Details

Several approaches of computing the similarity measure (simple, HerlaEtAl2021, tsa_WLdetection, rta_WLdetection) represent different flavours of the approach detailed in Herla et al (2021). In essence, they are a simple approach to incorporate avalanche hazard relevant characteristics into the score by computing the score as arithmetic mean of 4 different grain type classes:

• weak layers (wl): SH and DH

• new snow (pp): PP and DF

• crusts (cr): MFcr and IF

• bulk: the rest (i.e., predominantly RG, FC, FCxr --- MF falls also in here, will maybe be adjusted in future.)

Additionally, for classes wl and cr, vertical windows are computed to weigh layers more heavily that have no other wl or cr grain types in their neighborhood.

Type simple deviates from simple2 (= HerlaEtAl2021) by computing the aforementioned vertical windows based on heuristic depth ranges (i.e., Surface--30cm depth--80cm depth--150cm depth--Ground). It is otherwise identical to the simple2 type, which computes as many numbers of equidistant vertical windows as number of wl or cr are present in the profile.

Type tsa_WLdetection employs a similar approach as simple, but it identifies weak layers (wl) based on the Threshold Sum Approach (>= 5 TSA, lemons, German 'Nieten'). Therefore, the original profiles need to contain grain size information, which allows you to pre-compute the lemons for all layers (additionally to the otherwise necessary gain type and hardness information). It is thus more targeted to simulated profiles or detailed manual profiles of very high quality. While the former two types neglect hardness information of wl and cr classes, this type does not. Type rta_WLdetection works analogous, but uses RTA instead of TSA and a threshold of >= 0.8.

Unlike the former types, layerwise applies no weighting at all if used as per default. That means that the similarity of each individual layer contributes equally to the overall similarity measure. It is, however, very flexible in that any custom scaling factor can be applied to each layer. The resulting similarity score is then computed by

• sim = sim_gtype x sim_hardness (i.e., an array of similarities, one for each layer)

• simSP = sum(sim * scalingFactor) / sum(scalingFactor),

where the denominator ensures that the resulting score will be within [0, 1]. If you want to explore your own scaling approach, both input snow profiles need to contain a column called $layers$scalingFactor that store the desired factor. Type rta_scaling is a special case of layerwise, where the scaling is determined by the relative lemons of each layer (RTA, see Monti et al 2013). Type remotesensing makes use of the layerwise algorithm, but triggers an alternative similarity computation beforehand. Similarity is first computed from density and Optical Grain Size (ogs), and then the layerwise similarity is called upon to compute the global sim score.

The newest approach wsum_scaled differs from all approaches before on a foundational level. While all other approaches compute the similarity of two layers by multiplying their similarities in various layer properties (e.g., gtype, hardness), this approach computes a weighted sum of the similarities of three layer properties: gtype, hardness, layer stability. Differently than previous approaches, the layer stability is not only used for scaling purposes but also for the similarity calculation itself. By scaling the similarity with stability, unstable layers get more weight in the resulting score. By additionally including the similarity of layer stability in the similarity calculation, profiles with similar stability patterns get a higher score. By using a weighted sum to combine the three layer properties, the approach is identical to how the underlying alignment of the profiles is computed. The resulting similarity score is computed by

• sim = w1 x sim_gtype + w2 x sim_hardness + w3 x sim_stability (i.e., an array of similarities, one for each layer)

• simSP = sum(sim * stability) / sum(stability),

where layer stability defaults to p_unstable, or to scalingFactor (if apply_scalingFactor is TRUE).

NOTE that for all types that include stability indices (TSA, RTA, p_unstable, scalingFactor), these measures need to be computed prior to aligning the profiles (and therefore need to be present in the profiles provided to this function!)

References

Herla, F., Horton, S., Mair, P., & Haegeli, P. (2021). Snow profile alignment and similarity assessment for aggregating, clustering, and evaluating of snowpack model output for avalanche forecasting. Geoscientific Model Development, 14(1), 239–258. https://doi.org/10.5194/gmd-14-239-2021

Monti, F., & Schweizer, J. (2013). A relative difference approach to detect potential weak layers within a snow profile. Proceedings of the 2013 International Snow Science Workshop, Grenoble, France, 339–343. Retrieved from https://arc.lib.montana.edu/snow-science/item.php?id=1861

Examples

## first align two profiles, then assess the similarity of the aligned profiles
alignment <- dtwSP(SPpairs$A_modeled, SPpairs$A_manual)
SIM <- simSP(alignment$queryWarped, alignment$reference, verbose = TRUE)

## similarity of identical profiles
SIM <- simSP(alignment$queryWarped, alignment$queryWarped, verbose = TRUE)

## non-matched layers become apparent here:
alignment <- plotSPalignment(SPpairs$C_day1, SPpairs$C_day2, keep.alignment = TRUE,
                             rescale2refHS = FALSE, checkGlobalAlignment = FALSE)
simSP(alignment$queryWarped, alignment$reference, nonMatchedSim = 0.5)
## smaller similarity score due to 'penalty' of non-matched layers:
simSP(alignment$queryWarped, alignment$reference, nonMatchedSim = 0)
## even smaller similarity score due to higher impact of non-matched layer thickness:
simSP(alignment$queryWarped, alignment$reference, nonMatchedSim = 0, nonMatchedThickness = 1)

## detect WL based on lemons (instead of grain type alone):
P1 <- computeTSA(SPpairs$D_generalAlignment1)
P2 <- computeTSA(SPpairs$D_generalAlignment2)
alignment <- dtwSP(P1, P2, simType = "tsa_wldetection")
# sim based on WL-detection with TSA:
simSP(alignment$queryWarped, alignment$reference, type = "tsa_wldetection", verbose = TRUE)
# sim solely based on grain type, neglecting TSA information
simSP(alignment$queryWarped, alignment$reference, type = "simple", verbose = TRUE)

## RTA scaling type
P1 <- computeRTA(P1)
P2 <- computeRTA(P2)
alignment <- dtwSP(P1, P2, simType = "rta_scaling")
# sim based on scaling with RTA
simSP(alignment$queryWarped, alignment$reference, type = "rta_scaling")
# sim based on WL-detection with RTA
simSP(alignment$queryWarped, alignment$reference, type = "rta_wldetection")
# sim based on WL-detection with TSA
simSP(alignment$queryWarped, alignment$reference, type = "tsa_wldetection")

## layerwise similarity (i) unscaled...
simSP(alignment$queryWarped, alignment$reference, type = "layerwise", verbose = TRUE)

##... or (ii) with custom scaling factor (example only illustrative)
alignment$queryWarped$layers$scalingFactor <- 0.1
alignment$queryWarped$layers$scalingFactor[findPWL(alignment$queryWarped)] <- 1
alignment$reference$layers$scalingFactor <- 0.1
alignment$reference$layers$scalingFactor[findPWL(alignment$reference)] <- 1
simSP(alignment$queryWarped, alignment$reference, type = "layerwise",
      apply_scalingFactor = TRUE, verbose = TRUE)