6. Spatial Analysis

This vignette demonstrates some additional spatially interpolated statistics of a stress field.

library(tectonicr)
library(ggplot2) # load ggplot library

data("san_andreas")

data("cpm_models")
por <- cpm_models[["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 (, i.e. 2nd moment = 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           CI     skewness     kurtosis            R 
#>  160.0000000    9.0000000    5.4361579   -0.3002548    2.4245403    0.7981170

Spatial 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)
subset(spatial_stats_R, !is.na(mean)) |> head()
#> Simple feature collection with 6 features and 23 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -110.8569 ymin: 34.99224 xmax: -107.9222 ymax: 39.97615
#> Geodetic CRS:  WGS 84
#>     lon.PoR   lat.PoR  mean.PoR       sd       var   25%.PoR quasi-median.PoR
#> 24 84.83339 -61.15394  62.75215 32.70755 0.4788668  45.75289         70.68846
#> 25 85.83339 -61.15394  84.28615 35.03544 0.5266049  70.68846         74.88230
#> 26 86.83339 -61.15394 115.54625 49.91392 0.7808168  54.91857         90.72089
#> 35 75.83339 -60.15394 148.00270 22.18433 0.2590570 148.14101        150.69395
#> 42 82.83339 -60.15394  72.66492 14.94546 0.1272296  67.77922         75.91797
#> 43 83.83339 -60.15394  81.60114 22.00234 0.2554180  47.91126         76.92189
#>      75%.PoR median.PoR        CI   skewness   kurtosis     meanR   R N
#> 24  72.28118   70.68846 180.00000 -1.5897927  1.1003251 0.5211332 100 5
#> 25 127.81138   80.70468 180.00000 -1.0181445 -1.6083568 0.4733951 100 4
#> 26 111.23202   90.72089 180.00000  0.6657998 -0.7271910 0.2191832 100 5
#> 35 157.91252  151.01426 132.03984  1.7075040  3.5695991 0.7409430 100 4
#> 42  95.51644   72.35056  73.76277  1.0204726  4.1055756 0.8727704 100 4
#> 43  78.01755   76.92189  90.71189 -0.1367145  0.1865552 0.7445820 100 5
#>          mdr                   geometry      lat       lon      mean       25%
#> 24 0.6719937 POINT (-110.1257 39.15542) 39.15542 -110.1257 120.36224 103.36298
#> 25 0.1758576 POINT (-110.4569 39.56431) 39.56431 -110.4569 142.56213 128.96443
#> 26 0.6370203 POINT (-110.7839 39.97615) 39.97615 -110.7839 174.48900 113.86132
#> 35 0.6239140 POINT (-107.9222 34.99224) 34.99224 -107.9222  19.09291  19.23121
#> 42 0.6932675 POINT (-110.5034 37.78685) 37.78685 -110.5034 128.29658 123.41089
#> 43 0.6433647 POINT (-110.8569 38.19922) 38.19922 -110.8569 137.88256 104.19268
#>    quasi-median        75%    median
#> 24    128.29855 129.891269 128.29855
#> 25    133.15827   6.087356 138.98065
#> 26    149.66364 170.174771 149.66364
#> 35     21.78415  29.002725  22.10446
#> 42    131.54964 151.148106 127.98222
#> 43    133.20331 134.298968 133.20331

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), mode = TRUE)

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 = .3, linewisth = .5, color = "grey30", position = "center_spoke") +
  geom_spoke(aes(lon, lat, angle = deg2rad(90 - mean), alpha = sd, color = mdr), radius = .75, 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 = .3, color = "grey30", position = "center_spoke") +
  geom_spoke(aes(lon, lat, angle = deg2rad(90 - median), alpha = CI, color = skewness), radius = .5, 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 = .3, color = "grey30", position = "center_spoke") +
  geom_spoke(aes(lon, lat, angle = deg2rad(90 - mode), alpha = CI, color = abs(kurtosis)), radius = .5, 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))

Tobias Stephan