This vignette demonstrates some additional spatially interpolated statistics of a stress field.
data("san_andreas")
data("cpm_models")
por <- cpm_models |>
subset(model == "NNR-MORVEL56") |>
equivalent_rotation("na", "pa")
plate_boundary <- subset(plates, plates$pair == "na-pa")circular_summary() yields several statistics estimates
for a given vector of angles, such as mean, median, standard deviation,
quasi-quantiles, mode, 95% confidence angle, as wells as the moments
(2nd = variance, 3rd = skewness, 4th = kurtosis):
circular_summary(san_andreas$azi, w = 1 / san_andreas$unc)
#> n mean sd var 25% quasi-median
#> 1126.0000000 9.9840096 19.2389339 0.2018830 15.0000000 35.6250119
#> 75% median mode 95%CI skewness kurtosis
#> 160.0000000 9.0000000 9.1406250 2.1700917 -0.3002548 2.4245403
#> R
#> 0.7981170Spatial analysis
Spatial analysis and interpolation of stress data using
stress2grid_stats() or PoR_stress2grid_stats()
(analysis in the PoR coordinate system) uses a moving window with a user
defined cell-size (im km) and calculates the summary statistics within
each cell:
spatial_stats_R <- PoR_stress2grid_stats(san_andreas, PoR = por, gridsize = 1, R_range = 100)
head(spatial_stats_R)
#> Simple feature collection with 6 features and 25 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: -108.3822 ymin: 36.60408 xmax: -106.3458 ymax: 38.89699
#> Geodetic CRS: WGS 84
#> # A tibble: 6 × 26
#> lon.PoR lat.PoR mean.PoR sd var `25%.PoR` `quasi-median.PoR` `75%.PoR`
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 76.8 -62.2 NA NA NA NA NA NA
#> 2 77.8 -62.2 NA NA NA NA NA NA
#> 3 79.8 -62.2 NA NA NA NA NA NA
#> 4 80.8 -62.2 NA NA NA NA NA NA
#> 5 81.8 -62.2 NA NA NA NA NA NA
#> 6 82.8 -62.2 NA NA NA NA NA NA
#> # ℹ 18 more variables: median.PoR <dbl>, mode.PoR <dbl>, `95%CI` <dbl>,
#> # skewness <dbl>, kurtosis <dbl>, meanR <dbl>, R <dbl>, N <int>, mdr <dbl>,
#> # geometry <POINT [°]>, lat <dbl>, lon <dbl>, mean <dbl>, `25%` <dbl>,
#> # `quasi-median` <dbl>, `75%` <dbl>, median <dbl>, mode <dbl>One can also specify a range of cell-sizes for a wavelength analysis:
spatial_stats <- PoR_stress2grid_stats(san_andreas, PoR = por, gridsize = 1, R_range = seq(50, 350, 100), kappa = 10)The mean azimuth for each grid cell:
trajectories <- eulerpole_loxodromes(x = por, n = 40, cw = FALSE)
ggplot(spatial_stats) +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(data = san_andreas, aes(lon, lat, angle = deg2rad(90 - azi)), radius = .5, color = "grey30", position = "center_spoke") +
geom_spoke(aes(lon, lat, angle = deg2rad(90 - mean), alpha = sd, color = mdr), radius = 1, position = "center_spoke", lwd = 1) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat)) +
scale_alpha(name = "Standard deviation", range = c(1, .25)) +
scale_color_viridis_c(
"Wavelength\n(R-normalized mean distance)",
limits = c(0, 1),
breaks = seq(0, 1, .25)
) +
facet_wrap(~R)
To filter the range of search windows to only keep the shortest wavelength (R) with the least variance for each grid cell, use compact_grid2().
spatial_stats_comp <- spatial_stats |>
compact_grid2(var)Interpolated median stress field color-coded by the skewness within each search window:
ggplot(spatial_stats_comp) +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(data = san_andreas, aes(lon, lat, angle = deg2rad(90 - azi)), radius = .5, color = "grey30", position = "center_spoke") +
geom_spoke(aes(lon, lat, angle = deg2rad(90 - median), alpha = `95%CI`, color = skewness), radius = 1, position = "center_spoke", lwd = 1) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat)) +
scale_alpha(name = "95% CI", range = c(1, .25)) +
scale_color_viridis_c(
"Skewness"
)
Interpolated mode of the stress field color-coded by the absolute kurtosis within each search window:
ggplot(spatial_stats_comp) +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(data = san_andreas, aes(lon, lat, angle = deg2rad(90 - azi)), radius = .5, color = "grey30", position = "center_spoke") +
geom_spoke(aes(lon, lat, angle = deg2rad(90 - mode), alpha = `95%CI`, color = abs(kurtosis)), radius = 1, position = "center_spoke", lwd = 1) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat)) +
scale_alpha(name = "95% CI", range = c(1, .25)) +
scale_color_viridis_c(
"|Kurtosis|"
)
Heat maps for the spatial statistics
PoR_stress2grid_stats() and
stress2grid_stats() allow to create heatmaps showing the
spatial patterns of any desired statistical estimate (from
circular_summary()). Some examples:
Spatial central moments
Spatial variance
ggplot(spatial_stats_comp) +
ggforce::geom_voronoi_tile(
aes(lon, lat, fill = var),
max.radius = .7, normalize = FALSE
) +
scale_fill_viridis_c(limits = c(0, 1)) +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(
aes(lon, lat, angle = deg2rad(90 - mean)),
radius = .5, position = "center_spoke", lwd = .2, colour = "white"
) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat))
Skewness:
Skewness is a measure for the asymmetry of the probability distribution. It can be either counterclockwise or clockwise skewed, hence values can range between negative and positive numbers, respectively. This can be best visualized in a diverging color-sequence:
ggplot(spatial_stats_comp) +
ggforce::geom_voronoi_tile(
aes(lon, lat, fill = skewness),
max.radius = .7, normalize = FALSE
) +
scale_fill_gradient2("|Skewness|", low = "#001260", mid = "#EBE5E0", high = "#590007") +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(
aes(lon, lat, angle = deg2rad(90 - median), alpha = var),
radius = .5, position = "center_spoke", lwd = .2, colour = "white"
) +
scale_alpha("variance", range = c(1, 0)) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat))
Kurtosis
Kurtosis is a measure of the “tailedness” of the probability distribution. Here, colors are in a square-root scale:
ggplot(spatial_stats_comp) +
ggforce::geom_voronoi_tile(
aes(lon, lat, fill = abs(kurtosis)),
max.radius = .7, normalize = FALSE
) +
scale_fill_viridis_c("|Kurtosis|", transform = "sqrt") +
geom_sf(data = plate_boundary, color = "red") +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(
aes(lon, lat, angle = deg2rad(90 - mode), alpha = var),
radius = .5, position = "center_spoke", lwd = .2, colour = "white"
) +
scale_alpha("variance", range = c(1, 0)) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat))
Kernel dispersion
Another way to analyse spatial misfits is the kernel dispersion,
i.e. the local dispersion within a user-defined window (kernel). The
kernel´s half width can be a single number (km) or a range of widths.
The latter requires to compact the grid result (x) to find
the smallest kernel size containing the the least dispersion
(compact_grid(x, 'dispersion')).
It is recommended to calculate the kernel dispersion on PoR transformed data to avoid angle distortions due to projections.
san_andreas_por <- san_andreas
san_andreas_por$azi <- PoR_shmax(san_andreas, por, "right")$azi.PoR # transform to PoR azimuth
san_andreas_por$prd <- 135 # test direction
san_andreas_kdisp <- kernel_dispersion(san_andreas_por, gridsize = 1, R_range = seq(50, 350, 100))
san_andreas_kdisp <- compact_grid(san_andreas_kdisp, "dispersion")
ggplot(san_andreas_kdisp) +
ggforce::geom_voronoi_tile(
aes(lon, lat, fill = stat),
max.radius = .7, normalize = FALSE
) +
scale_fill_viridis_c("Dispersion", limits = c(0, 1)) +
geom_sf(data = trajectories, lty = 2) +
geom_spoke(
data = san_andreas,
aes(lon, lat, angle = deg2rad(90 - azi), alpha = unc),
radius = .5, position = "center_spoke", lwd = .2, colour = "white"
) +
scale_alpha("Standard deviation", range = c(1, .25)) +
coord_sf(xlim = range(san_andreas$lon), ylim = range(san_andreas$lat))