Package {tectonicr}

Title:	Analyzing the Orientation of Maximum Horizontal Stress
Version:	0.4.9
Description:	Models the direction of the maximum horizontal stress using relative plate motion parameters. Statistical algorithms to evaluate the modeling results compared with the observed data. Provides plots to visualize the results. Methods described in Stephan et al. (2023) <doi:10.1038/s41598-023-42433-2> and Wdowinski (1998) <doi:10.1016/S0079-1946(98)00091-3>.
License:	GPL (≥ 3)
URL:	https://tobiste.github.io/tectonicr/
BugReports:	https://github.com/tobiste/tectonicr/issues
Depends:	R (≥ 4.1.0)
Imports:	boot, circular (≥ 0.5.0), dplyr, ggplot2, lifecycle, methods, RColorBrewer, sf, smoothr (≥ 1.0.1), spatstat.explore (≥ 3.2.7), spatstat.geom (≥ 3.2.9), spatstat.univar (≥ 2.0.3), spatstat.utils (≥ 3.0.4), terra, tidyr, viridis, zoo (≥ 1.8.12)
Suggests:	ggforce, grid, knitr, rmarkdown, roxygen2, testthat (≥ 3.0.0), tidyterra
VignetteBuilder:	knitr
Config/testthat/edition:	3
Encoding:	UTF-8
Language:	en-US
LazyData:	true
RoxygenNote:	7.3.3
NeedsCompilation:	no
Packaged:	2026-08-28 19:34:50 UTC; tstephan
Author:	Tobias Stephan [aut, cre]
Maintainer:	Tobias Stephan <tobias.stephan1@yahoo.com>
Repository:	CRAN
Date/Publication:	2026-08-29 12:50:02 UTC

library(tectonicr)

Description

Modeling the Direction of the Maximum Horizontal Stress using Relative Plate Motion

Details

Further details and theoretical background are provided by Wdowinski (1998) and Stephan et al. (2023).

Note

A list of documented functions may be viewed by typing help(package="tectonicr").

Author(s)

Tobias Stephan

References

Wdowinski (1998) "A theory of intraplate tectonics". JGR: Solid Earth, 103(3), 5037<U+2013>5059.

Stephan, T., Enkelmann, E., and Kroner, U. "Analyzing the horizontal orientation of the crustal stress adjacent to plate boundaries". Sci Rep 13. 15590 (2023).

See Also

Useful links:

• https://tobiste.github.io/tectonicr/

• Report bugs at https://github.com/tobiste/tectonicr/issues

Azimuth Conversion From PoR to Geographical Coordinate Reference System

Description

Conversion of PoR azimuths into geographical azimuths

Usage

PoR2Geo_azimuth(x, PoR, axial = TRUE)

Arguments

Value

numeric vector of transformed azimuths (in degrees)

References

Stephan, T., Enkelmann, E., and Kroner, U. "Analyzing the horizontal orientation of the crustal stress adjacent to plate boundaries". Sci Rep 13. 15590 (2023). doi:10.1038/s41598-023-42433-2.

See Also

PoR_shmax()

Examples

data("nuvel1")
# North America relative to Pacific plate:
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
data("san_andreas")
san_andreas$azi.PoR <- PoR_shmax(san_andreas, PoR)

# convert back to geo CRS
PoR2Geo_azimuth(san_andreas, PoR) |> head()

Azimuth Conversion from Geographical to PoR Coordinate Reference System

Description

Transforms azimuths and models the direction of maximum horizontal stress \sigma_\text{Hmax} in the Euler pole (Pole of Rotation) coordinate reference system. When type of plate boundary is given, it also gives the deviation from the theoretically predicted azimuth of \sigma_\text{Hmax}, the circular distance, and the normalized \chi^2 statistics.

Usage

PoR_azimuth(x, PoR, axial = TRUE)

PoR_shmax(x, PoR, type = c("none", "in", "out", "right", "left"), axial = TRUE)

Arguments

Details

The theoretical azimuth of \sigma_\text{Hmax} in the pole of rotation reference system is 0 (or 180), 45, 90, 135° if the stress is sourced by an outward, sinistral, inward, or dextral moving plate boundary, respectively. directions of \sigma_\text{Hmax} with respect to the four plate boundary types.

Value

PoR_azimuth returns numeric vector of the transformed azimuth in degrees. PoR_shmax returns either a numeric vector of the azimuths in the transformed coordinate system (in degrees), or a "data.frame" with

azi.PoR

the transformed azimuths (in degrees),

prd

the predicted azimuths (in degrees),

dev

the deviation between the transformed and the predicted azimuth in degrees (positive for counterclockwise deviation of observed azimuth wrt. predicted azimuth),

nchisq

the Norm \chi^2 test statistic, and

cdist

the angular distance between the transformed and the predicted azimuth.

References

Stephan, T., Enkelmann, E., and Kroner, U. "Analyzing the horizontal orientation of the crustal stress adjacent to plate boundaries". Sci Rep 13. 15590 (2023). doi:10.1038/s41598-023-42433-2.

See Also

model_shmax() to compute the theoretical direction of \sigma_\text{Hmax} in the geographical reference system. deviation_shmax() to compute the deviation of the modeled direction from the observed direction of \sigma_\text{Hmax}. norm_chisq() to calculate the normalized \chi^2 statistics. circular_distance() to calculate the angular distance.

Examples

data("nuvel1")
# North America relative to Pacific plate:
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")

data("san_andreas")
res <- PoR_shmax(san_andreas, PoR, type = "right")
head(res)

Coordinates of the Pole of Rotation Reference System

Description

Retrieve the PoR equivalent coordinates of an object

Usage

PoR_coordinates(x, PoR)

Arguments

Value

PoR_coordinates() returns data.frame with the PoR coordinates (lat.PoR, lon.PoR).

Examples

data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to Pacific plate
data("san_andreas")

# coordinates from sf object
san_andreas.por_sf <- PoR_coordinates(san_andreas, por)
head(san_andreas.por_sf)

# coordinates from data.frame
san_andreas.por_df <- PoR_coordinates(sf::st_drop_geometry(san_andreas), por)
head(san_andreas.por_df)

PoR coordinate reference system

Description

Create the reference system transformed in Euler pole coordinate

Usage

PoR_crs(x)

Arguments

Details

The PoR coordinate reference system is oblique transformation of the geographical coordinate system with the Euler pole coordinates being the the translation factors.

Value

Object of class crs

See Also

sf::st_crs()

Examples

data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to Pacific plate
PoR_crs(por)

Distance to Pole of Rotation

Description

Retrieve the (angular) great-circle distance between a point and the PoR (Euler pole).

Usage

PoR_distance(x, PoR, FUN = orthodrome)

Arguments

Value

numeric vector. Great-circle distance in degree

Examples

data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to Pacific plate
data("san_andreas")

# distance form sf object
PoR_distance(san_andreas, por) |> head()

# distance form data.frame
PoR_distance(sf::st_drop_geometry(san_andreas), por) |> head()
PoR_distance(sf::st_drop_geometry(san_andreas), por, FUN = orthodrome) |> head()
PoR_distance(sf::st_drop_geometry(san_andreas), por, FUN = vincenty) |> head()

Map of data in Pole of Rotation reference frame

Description

Transforms the spatial data and azimuths into the PoR reference frame and shows them in a map

Usage

PoR_map(
  x,
  PoR,
  pb = NULL,
  type = c("none", "in", "out", "right", "left"),
  show.deviation = FALSE,
  ...
)

Arguments

Value

plot

See Also

PoR_shmax(), axes(), tectonicr.colors()

Examples

data("nuvel1")
na_pa <- subset(nuvel1, nuvel1$plate.rot == "na")

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")

data("san_andreas")
PoR_map(san_andreas, PoR = na_pa, pb = plate_boundary, type = "right", show.deviation = TRUE)

Spatial Interpolation of SHmax in PoR Coordinate Reference System

Description

The data is transformed into the PoR system before the interpolation. The interpolation grid is returned in geographical coordinates and azimuths.

Usage

PoR_stress2grid(
  x,
  PoR,
  grid = NULL,
  PoR_grid = TRUE,
  lon_range = NULL,
  lat_range = NULL,
  gridsize = 2.5,
  remove_PoR = FALSE,
  ...
)

PoR_stress2grid_stats(
  x,
  PoR,
  grid = NULL,
  PoR_grid = TRUE,
  lon_range = NULL,
  lat_range = NULL,
  gridsize = 2.5,
  remove_PoR = FALSE,
  ...
)

Arguments

Details

Stress field and wavelength analysis in PoR system and back-transformed

Value

sf object containing

lon,lat

longitude and latitude in geographical CRS (in degrees)

lon.PoR,lat.PoR

longitude and latitude in PoR CRS (in degrees). Only if remove_PoR=TRUE

azi

geographical mean \sigma_\text{Hmax} in degree

azi.PoR

PoR mean \sigma_\text{Hmax} in degree. Only if remove_PoR=TRUE

sd

Standard deviation of \sigma_\text{Hmax} in degrees

R

Search radius in km

mdr

Mean distance of datapoints per search radius

N

Number of data points in search radius

See Also

stress2grid(), compact_grid()

Examples

data("san_andreas")
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
PoR_stress2grid(san_andreas, PoR) |> head()

## Not run: 
PoR_stress2grid_stats(san_andreas, PoR, mode = TRUE) |> head()

## End(Not run)

Centrically aligned geom_spoke marker

Description

"position" subclass "center_spoke" to center ggplot::geom_spoke() marker at its origin

Usage

PositionCenterSpoke

Format

An object of class PositionCenterSpoke (inherits from Position, ggproto, gg) of length 2.

Source

https://stackoverflow.com/questions/55474143/how-to-center-geom-spoke-around-their-origin/

Absolute Plate Velocity

Description

Calculates the absolute angular velocity of plate motion

Usage

abs_vel(w, alpha, r = earth_radius())

Arguments

Value

numeric (unit of velocity: km/Myr)

See Also

earth_radius()

Examples

abs_vel(0.21, 0)
abs_vel(0.21, 45)
abs_vel(0.21, 90)

Degrees to Radians

Description

Helper functions to transform between angles in degrees and radians.

Usage

rad2deg(rad)

deg2rad(deg)

Arguments

Value

numeric. angle in degrees or radians.

Examples

deg2rad(seq(-90, 90, 15))
rad2deg(seq(-pi / 2, pi / 2, length = 13))

Angle Between Two Vectors

Description

Calculates the angle between two vectors

Usage

angle_vectors(x, y)

Arguments

Value

numeric. angle in degrees

Examples

u <- c(1, -2, 3)
v <- c(-2, 1, 1)
angle_vectors(u, v) # 96.26395

Angular Randomisation Test of Homogeneity

Description

Performs Angular Randomisation Test for homogeneity on two samples of circular data after Ruxton et al. (2023). P-values are estimated using permutation.

Usage

ar_test(x, y, n_perm = 1000L, axial = FALSE, alpha = NULL)

Arguments

Details

Null Hypothesis (H_0): The two circular samples share an identical underlying probability distribution.

Alternative Hypothesis (H_{1}): The two samples come from different distributions.

Interpretation

• Small p-value (p < \alpha, e.g., <0.05): Reject the null hypothesis. This indicates strong evidence that the two samples come from different circular distributions (differing in central tendency/mean direction or shape).

• Large p-value (p \ge \alpha): Fail to reject the null hypothesis; there is insufficient evidence to claim the two circular samples differ.

Value

list containing the test statistic, the p-value, the significance value alpha and a logical decision whether to reject the null hypothesis or not.

Note

Concentration Differences: The test can suffer from markedly lower statistical power if the underlying unimodal distributions differ by concentration (dispersion/spread) rather than location—especially with small, uneven sample sizes where the smaller sample comes from the more concentrated distribution.

Axial/Multimodal Data: ART performs poorly and loses power when applied to axially symmetric or symmetrically multimodal distributions.

References

Ruxton, G.D., Malkemper, E.P. & Landler, L. Evaluating the power of a recent method for comparing two circular distributions: an alternative to the Watson U2 test. Sci Rep 13, 10007 (2023). https://doi.org/10.1038/s41598-023-36960-1

See Also

Other Tests: kuiper_test(), norm_chisq(), rayleigh-test, watson_test(), watson_two_sample, watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

set.seed(20250411)
x1 <- c(35, 45, 50, 55, 60, 70, 85, 95, 105, 120)
x2 <- c(75, 80, 90, 100, 110, 130, 135, 140, 150, 160, 165)
ar_test(x1, x2)

# San Andreas Fault Data:
data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
ar_test(sa.por$azi.PoR, rvm(100, 135, 10), axial = TRUE, alpha = 0.05)

Plot axes

Description

Show direction axes in a map

Usage

axes(
  x,
  y,
  angle,
  radius = 0.5,
  arrow.code = 1,
  arrow.length = 0,
  add = FALSE,
  ...
)

Arguments

Value

No return value, called for side effects

Examples

data("san_andreas")
axes(san_andreas$lon, san_andreas$lat, san_andreas$azi, add = FALSE)

Circular Mean Difference

Description

The circular mean difference is based on the sample circular distance

Usage

circular_mean_difference(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_mean_difference_alt(x, w = NULL, axial = TRUE, na.rm = TRUE)

Arguments

Value

numeric

References

Mardia, K.V., and Jupp, P.E (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA. doi:10.1002/9780470316979

See Also

sample_circular_distance()

Examples

data("san_andreas")
circular_mean_difference(san_andreas$azi)
circular_mean_difference(san_andreas$azi, weighting(san_andreas$unc))

circular_mean_difference_alt(san_andreas$azi)
circular_mean_difference_alt(san_andreas$azi, weighting(san_andreas$unc))

Summary Statistics of Circular Data

Description

Calculate the (weighted median) and standard deviation of orientation data.

Usage

circular_mean(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_var(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_sd(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_median(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_quantiles(x, w = NULL, axial = TRUE, na.rm = TRUE)

circular_IQR(x, w = NULL, axial = TRUE, na.rm = TRUE)

Arguments

Value

numeric vector

Note

Weighting may be the reciprocal of the data uncertainties.

Weightings have no effect on quasi-quantiles if length(x) %% 4 == 0.

References

Mardia, K.V. (1972). Statistics of Directional Data: Probability and Mathematical Statistics. London: Academic Press.

Mardia, K.V., and Jupp, P.E (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA. doi:10.1002/9780470316979

N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

Ziegler, M. O.; Heidbach O. (2019). Manual of the Matlab Script Stress2Grid v1.1. WSM Technical Report 19-02, GFZ German Research Centre for Geosciences. doi:10.2312/wsm.2019.002

Heidbach, O., Tingay, M., Barth, A., Reinecker, J., Kurfess, D., & Mueller, B. (2010). Global crustal stress pattern based on the World Stress Map database release 2008. Tectonophysics 482, 3<U+2013>15, doi:10.1016/j.tecto.2009.07.023

Examples

set.seed(1)
x <- rvm(10, 0, 100) %% 180
unc <- stats::runif(length(x), 0, 10)
w <- weighting(unc)
circular_mean(x, w)
circular_var(x, w)
circular_sd(x, w)
circular_median(x, w)
circular_quantiles(x, w)
circular_IQR(x, w)

data("san_andreas")
w2 <- weighting(san_andreas$unc)
circular_mean(san_andreas$azi)
circular_mean(san_andreas$azi, w2)
circular_median(san_andreas$azi)
circular_median(san_andreas$azi, w2)
circular_quantiles(san_andreas$azi)
circular_quantiles(san_andreas$azi, w2)
circular_var(san_andreas$azi)
circular_var(san_andreas$azi, w2)

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
circular_mean(sa.por$azi.PoR, w2)
circular_median(sa.por$azi.PoR, w2)
circular_var(sa.por$azi.PoR, w2)
circular_quantiles(sa.por$azi.PoR, w2)

Circular Kernel Density Estimation

Description

Kernel density estimates for circular data from a given kernel (von Mises, wrapped Cauchy, and wrapped Normal distribution and bandwidth

Usage

circular_density(
  x,
  z = NULL,
  bw = NULL,
  weights = NULL,
  na.rm = TRUE,
  from = 0,
  to = 360,
  n = 512L,
  axial = TRUE,
  kappa = NULL,
  rho = NULL,
  sd = NULL,
  c = NULL,
  kernel = c("vonmises", "wrappedcauchy", "wrappednormal", "wrappedlevy"),
  adjust = 1,
  subdensity = FALSE
)

Arguments

Value

Object of class "density"

See Also

stats::density(), dvm(), dwcauchy(), dwnorm(), and plot_density()

Examples

w <-  weighting(san_andreas$unc)

# von Mises kernel density
circular_density(san_andreas$azi, kappa = 100)
circular_density(san_andreas$azi, weights = w, kappa = 100)

# wrapped Cauchy kernel density
circular_density(san_andreas$azi, rho = 0.9, kernel = "wrappedcauchy")

# wrapped Normal kernel density
circular_density(san_andreas$azi, sd = 5, kernel = "wrappednormal")

Bootstrapped Estimates for Circular Dispersion

Description

Calculates bootstrapped estimates of the circular dispersion, its standard error and its confidence interval.

Usage

circular_dispersion_boot(
  x,
  y = NULL,
  w = NULL,
  w.y = NULL,
  R = 1000,
  conf.level = 0.95,
  ...
)

Arguments

Value

list containing:

MLE

the maximum likelihood estimate of the circular dispersion

sde

standard error of MLE

CI

lower and upper limit of the confidence interval of MLE

See Also

circular_dispersion()

Examples

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
circular_dispersion(sa.por$azi.PoR, y = 135, w = weighting(san_andreas$unc))
circular_dispersion_boot(sa.por$azi.PoR, y = 135, w = weighting(san_andreas$unc), R = 1000)

Circular Mode

Description

MLE angle (maximum density) using a circular distribution kernel with specified concentration

Usage

circular_mode(x, ...)

Arguments

Value

numeric

References

N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

Examples

set.seed(20250411)
x <- rvm(10, 0, 100)

# Mode of von Mises kernel density (the default)
circular_mode(x)

# Mode of wrapped Cauchy kernel density
circular_mode(x, kernel = "wrappedcauchy")

Circular plot

Description

Circular plot

Usage

circular_plot(
  main = NULL,
  labels = TRUE,
  at = seq(0, 360 - 45, 45),
  cborder = TRUE,
  ...
)

Arguments

Value

none

Note

Polar diagram where angles increase clockwise.

Quantile-Quantile Linearised Plot for Circular Distributions

Description

Uniformly distributed orientations should yield a straight line through the origin. Systematic departures from linearity will indicate preferred orientation.

Usage

circular_qqplot(
  x,
  axial = TRUE,
  xlab = paste("i/(n+1)"),
  ylab = NULL,
  main = "Circular Quantile-Quantile Plot",
  add_line = TRUE,
  col = "#B63679FF",
  ...
)

Arguments

Value

plot

References

Borradaile, G. J. (2003). Statistics of earth science data: their distribution in time, space, and orientation (Vol. 351, p. 329). Berlin: Springer.

See Also

cunif()

Other circ-qqplot: vm_qqplot()

Examples

set.seed(20250411)

# von Mises distribution
x_vm <- rvm(100, mean = 0, kappa = 2)
circular_qqplot(x_vm, pch = 20)

x_wcauchy <- rwcauchy(100, mean = 0, rho = 0.5)
circular_qqplot(x_wcauchy, pch = 20)

# circular uniform data
x_cunif <- rcunif(100)
circular_qqplot(x_cunif, pch = 20)

Circular Range

Description

Length of the smallest arc which contains all the observations. The circular range is based on the sample circular distance.

Usage

circular_range(x, axial = TRUE, na.rm = TRUE)

Arguments

Value

numeric. angle in degrees

References

Mardia, K.V., and Jupp, P.E (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA. doi:10.1002/9780470316979

See Also

sample_circular_distance()

Examples

roulette <- c(43, 45, 52, 61, 75, 88, 88, 279, 357)
circular_range(roulette, axial = FALSE)

data("san_andreas")
circular_range(san_andreas$azi)

Standard Error of Mean Direction of Circular Data

Description

Measure of the chance variation expected from sample to sample in estimates of the mean direction (after Mardia 1972). It is a parametric estimate of the the circular standard error of the mean direction by the particular form of the standard error for the von Mises distribution. The approximated standard error of the mean direction is computed by the mean resultant length and the MLE concentration parameter \kappa.

Usage

circular_sd_error(x, w = NULL, axial = TRUE, na.rm = TRUE)

Arguments

Value

numeric

References

• Batschelet, E. (1971). Recent statistical methods for orientation data. "Animal Orientation, Symposium 1970 on Wallops Island". Amer. Inst. Biol. Sciences, Washington.

• Mardia, K.V. (1972). Statistics of Directional Data: Probability and Mathematical Statistics. London: Academic Press.

• N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

• Davis (1986) Statistics and data analysis in geology. 2nd ed., John Wiley & Sons.

See Also

mean_resultant_length(), circular_mean()

Examples

# Example data from Davis (1986), pp. 316
finland_stria <- c(
  23, 27, 53, 58, 64, 83, 85, 88, 93, 99, 100, 105, 113,
  113, 114, 117, 121, 123, 125, 126, 126, 126, 127, 127, 128, 128, 129, 132,
  132, 132, 134, 135, 137, 144, 145, 145, 146, 153, 155, 155, 155, 157, 163,
  165, 171, 172, 179, 181, 186, 190, 212
)
circular_sd_error(finland_stria, axial = FALSE)

data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
circular_sd_error(sa.por$azi.PoR, w = weighting(san_andreas$unc))

Circular Summary Statistics

Description

Circular mean, standard deviation, variance, quasi-quantiles, mode, 95% confidence angle, standardized skewness and kurtosis

Usage

circular_summary(
  x,
  w = NULL,
  axial = TRUE,
  mode = FALSE,
  kappa = NULL,
  fisher.CI = FALSE,
  conf.level = 0.95,
  na.rm = FALSE
)

Arguments

Value

named vector

See Also

circular_mean(), circular_sd(), circular_var(), circular_quantiles(), confidence_angle(), second_central_moment(), circular_mode()

Examples

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
circular_summary(sa.por$azi.PoR)
circular_summary(sa.por$azi.PoR, w = weighting(san_andreas$unc))

Compact Smoothed Stress Field

Description

Filter smoothed stress field containing a range of search radii or kernel half widths to find shortest wavelength (R) with the least circular sd. or dispersion (or any statistic) for each coordinate, respectively.

Usage

compact_grid(x, type = c("stress", "dispersion"))

compact_grid2(x, ..., FUN = min)

Arguments

Value

sf object

See Also

stress2grid(), PoR_stress2grid(), kernel_dispersion(), stress2grid_stats(), dplyr::dplyr_tidy_select()

Examples

data("san_andreas")
res <- stress2grid(san_andreas)
compact_grid(res) |> head()

## Not run: 
res2 <- stress2grid_stats(san_andreas)
compact_grid2(res2, var, FUN = min)

## End(Not run)

Confidence Interval around the Mean Direction of Circular Data after Batschelet (1971)

Description

Probabilistic limit on the location of the true or population mean direction, assuming that the estimation errors are normally distributed.

Usage

confidence_angle(x, conf.level = 0.95, w = NULL, axial = TRUE, na.rm = TRUE)

confidence_interval(x, conf.level = 0.95, w = NULL, axial = TRUE, na.rm = TRUE)

Arguments

Details

The confidence angle gives the interval, i.e. plus and minus the confidence angle, around the mean direction of a particular sample, that contains the true mean direction under a given level of confidence.

Value

Angle in degrees

References

• Batschelet, E. (1971). Recent statistical methods for orientation data. "Animal Orientation, Symposium 1970 on Wallops Island". Amer. Inst. Biol. Sciences, Washington.

• Mardia, K.V. (1972). Statistics of Directional Data: Probability and Mathematical Statistics. London: Academic Press. (p. 146)

• Davis (1986) Statistics and data analysis in geology. 2nd ed., John Wiley & Sons.

• Jammalamadaka, S. Rao and Sengupta, A. (2001). Topics in Circular Statistics, Sections 3.3.3 and 3.4.1, World Scientific Press, Singapore.

See Also

mean_resultant_length(), circular_sd_error()

Examples

# Example data from Davis (1986), pp. 316
finland_stria <- c(
  23, 27, 53, 58, 64, 83, 85, 88, 93, 99, 100, 105, 113,
  113, 114, 117, 121, 123, 125, 126, 126, 126, 127, 127, 128, 128, 129, 132,
  132, 132, 134, 135, 137, 144, 145, 145, 146, 153, 155, 155, 155, 157, 163,
  165, 171, 172, 179, 181, 186, 190, 212
)
confidence_angle(finland_stria, axial = FALSE)
confidence_interval(finland_stria, axial = FALSE)

data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
confidence_angle(sa.por$azi.PoR, w = weighting(san_andreas$unc))
confidence_interval(sa.por$azi.PoR, w = weighting(san_andreas$unc))

Confidence Interval around the Mean Direction of Circular Data after Fisher (1993)

Description

For large samples (n >=25) i performs are parametric estimate based on sample_circular_dispersion(). For smaller size samples, it returns a bootstrap estimate.

Usage

confidence_interval_fisher(
  x,
  conf.level = 0.95,
  w = NULL,
  axial = TRUE,
  na.rm = TRUE,
  boot = FALSE,
  R = 1000L,
  quiet = FALSE
)

Arguments

Value

list

References

N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

Examples

# Example data from Davis (1986), pp. 316
finland_stria <- c(
  23, 27, 53, 58, 64, 83, 85, 88, 93, 99, 100, 105, 113,
  113, 114, 117, 121, 123, 125, 126, 126, 126, 127, 127, 128, 128, 129, 132,
  132, 132, 134, 135, 137, 144, 145, 145, 146, 153, 155, 155, 155, 157, 163,
  165, 171, 172, 179, 181, 186, 190, 212
)
confidence_interval_fisher(finland_stria, axial = FALSE)
confidence_interval_fisher(finland_stria, axial = FALSE, boot = TRUE)

data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
confidence_interval_fisher(sa.por$azi.PoR, w = weighting(san_andreas$unc))
confidence_interval_fisher(sa.por$azi.PoR, w = weighting(san_andreas$unc), boot = TRUE)

Coordinate Correction

Description

Corrects the longitudes or latitudes to value between -180&deg and 180&deg or -90&deg and 90°

Usage

longitude_modulo(x)

latitude_modulo(x)

Arguments

Value

numeric

Examples

longitude_modulo(-361 + 5 * 360) # -1
latitude_modulo(-91 + 5 * 180) # 89

Coordinate Transformations

Description

Converts vector between Cartesian and geographical coordinate systems

Usage

cartesian_to_geographical(n)

geographical_to_cartesian(p)

geographical_to_spherical(p)

Arguments

Value

Functions return a (2- or 3-dimensional) vector representing a point in the requested coordinate system.

See Also

cartesian_to_spherical() and spherical_to_cartesian() for conversions to spherical coordinates

Examples

n <- c(1, -2, 3)
cartesian_to_geographical(n)
p <- c(50, 10)
geographical_to_cartesian(p)

Coordinate Transformations

Description

Converts vector between Cartesian and spherical coordinate systems

Usage

cartesian_to_spherical(n)

spherical_to_cartesian(p)

spherical_to_geographical(p)

Arguments

Value

Functions return a (2- or 3-dimensional) vector representing a point in the requested coordinate system.

See Also

cartesian_to_geographical() and geographical_to_cartesian() for conversions to geographical coordinates

Examples

n <- c(1, -2, 3)
cartesian_to_spherical(n) # 36.699, -63.435
p <- c(50, 10)
spherical_to_cartesian(p) # 0.75, 0.13, 0.64

Global model of current plate motions

Description

Compilation of global models for current plate motions, including NUVEL1 (DeMets et al. 1990), NNR-NUVEL1A (DeMets et al., 1990), NNR-MORVEL56 (Argus et al., 2011), REVEL (Sella et al., 2002), GSRM2.1 (Kreemer et al., 2014), HS2-NUVEL1 (Gripp and Gordon, 1990), HS3-NUVEL1A (Gripp and Gordon, 2002), P073 (Chase 1978), AM1-2 (Minster and Jordan, 1978), ITRF2020-PPM (Altamimi et al. 2023), and PB2002 (Bird, 2003)

Usage

data('cpm_models')

Format

list containing objects of class data.frame

plate.name

The rotating plate

plate.rot

The abbreviation of the plate's name

lat,lon

Coordinates of the Pole of Rotation

angle

The amount of rotation (angle in 1 Myr)

plate.fix

The anchored plate, i.e. plate.rot moves relative to plate.fix

model

Model for current global plate motion

References

Altamimi, Z., Métivier, L., Rebischung, P., Collilieux, X., Chanard, K., Barnéoud, J., 2023. ITRF2020 Plate Motion Model. Geophys. Res. Lett. 50, 1–7. doi:10.1029/2023GL106373

Argus, D.F., Gordon, R.G., 1991. No-net-rotation model of current plate velocities incorporating plate motion model NUVEL-1. Geophys. Res. Lett. 18, 2039–2042. doi: 10.1029/91GL01532

Argus, D. F., Gordon, R. G., & DeMets, C. (2011). Geologically current motion of 56 plates relative to the no-net-rotation reference frame. Geochemistry, Geophysics, Geosystems, 12(11). 10.1029/2011GC003751.

Chase, C.G. (1978). Plate kinematics: The Americas, East Africa, and the rest of the world. Earth Planet. Sci. Lett. 37, 355–368. doi: doi:10.1016/0012-821X(78)90051-1

Bird, P. (2003), An updated digital model of plate boundaries, Geochem. Geophys. Geosyst., 4, 1027, doi: 10.1029/2001GC000252, 3.

DeMets, C., Gordon, R. G., Argus, D. F., & Stein, S. (1990). Current plate motions. Geophysical Journal International, 101(2), 425-478. doi:10.1111/j.1365-246X.1990.tb06579.x.

Gripp, A. E., & Gordon, R. G. (2002). Young tracks of hotspots and current plate velocities. Geophysical Journal International, 150(2), 321<U+2013>361. doi:10.1046/j.1365-246X.2002.01627.x.

Kreemer, C., Blewitt, G., & Klein, E. C. (2014). A geodetic plate motion and Global Strain Rate Model. Geochemistry, Geophysics, Geosystems, 15(10), 3849<U+2013>3889. doi: 10.1002/2014GC005407.

Minster, J. and Jorda, T. (1978). Present-day plate motions. Journal of Geophysical Research, 83, doi:10.1029/jb083ib11p05331.

Sella, G. F., Dixon, T. H., & Mao, A. (2002). REVEL: A model for Recent plate velocities from space geodesy. Journal of Geophysical Research: Solid Earth, 107(B4). doi: 10.1029/2000jb000033.

Examples

data("cpm_models")
head(cpm_models[[1]])

The Circular Uniform Distribution

Description

Density, probability distribution function, quantiles, and random generation for the circular uniform distribution.

Usage

rcunif(n, axial = FALSE)

dcunif(theta, axial = FALSE, log = FALSE)

pcunif(theta, axial = FALSE, lower.tail = TRUE, log.p = FALSE)

qcunif(p, axial = FALSE, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

dcunif gives the density, pcunif gives the probability of circular uniform distribution function, rcunif generates random deviates (in degrees), and qcunif provides quantiles (in degrees).

See Also

wnorm, wcauchy, and vonmises

Examples

set.seed(1)
x <- rcunif(5)

dcunif(x)
dcunif(x, axial = TRUE)

pcunif(x)
qcunif(c(.25, .5, .75))

Transforms coordinates and azimuths into PoR coordinates system

Description

Convenience function to add PoR coordinates and PoR azimuths to data

Usage

data2PoR(x, PoR)

Arguments

Value

sf object in PoR CRS with additional columns lon.PoR, lat.PoR, and azi.PoR

Examples

por <- subset(nuvel1, nuvel1$plate.rot == "na")
data2PoR(san_andreas, por)

Normalize Angle Between Two Directions

Description

Normalizes the angle between two directions to the acute angle in between, i.e. angles between 0 and 90°

Usage

deviation_norm(x, y = NULL)

Arguments

Value

numeric vector, acute angles between two directions, i.e. values between 0 and 90°

Author(s)

Tobias Stephan

Examples

deviation_norm(175, 5)
deviation_norm(c(175, 95, 0), c(5, 85, NA))
deviation_norm(c(-5, 85, 95, 175, 185, 265, 275, 355, 365))

Deviation of Observed and Predicted Directions of Maximum Horizontal Stress

Description

Calculate the angular difference between the observed and modeled direction of maximum horizontal stresses (\sigma_\text{Hmax}) along great circles, small circles, and loxodromes of the relative plate motion's Euler pole

Usage

deviation_shmax(prd, obs)

Arguments

Details

Deviation is positive for counterclockwise deviation of observed azimuth wrt. predicted azimuth.

Value

An object of class data.frame

dev.gc

Deviation of observed stress from modeled \sigma_\text{Hmax} following great circles

dev.sc

Small circles

dev.ld.cw

Clockwise loxodromes

dev.ld.ccw

Counter-clockwise loxodromes

Author(s)

Tobias Stephan

References

Stephan, T., Enkelmann, E., and Kroner, U. "Analyzing the horizontal orientation of the crustal stress adjacent to plate boundaries". Sci Rep 13. 15590 (2023). doi:10.1038/s41598-023-42433-2.

See Also

model_shmax() to calculate the theoretical direction of \sigma_\text{Hmax}.

Examples

data("nuvel1")
# North America relative to Pacific plate:
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")

# the point where we want to model the SHmax direction:
point <- data.frame(lat = 45, lon = 20)

prd <- model_shmax(point, PoR)
deviation_shmax(prd, obs = 90)

Circular Distance and Dispersion

Description

Circular distance between two angles and circular dispersion of angles about a specified angle.

Usage

circular_distance(x, y, axial = TRUE, na.rm = TRUE)

circular_dispersion(
  x,
  y = NULL,
  w = NULL,
  w.y = NULL,
  axial = TRUE,
  na.rm = TRUE
)

circular_sd2(x, y, w = NULL, axial = TRUE, na.rm = TRUE)

Arguments

Details

Circular dispersion is a measure for the spread of data like the variance. Dispersion measures the spread about a given angles, whereas the variance measures the spread about the mean (Mardia and Jupp, 1999). When y = NULL the dispersion is identical to the variance.

Circular standard deviation in circular_sd2() is the transformed dispersion instead of the variance as for circular_sd().

Value

circular_distance returns a numeric vector of positive numbers, circular_dispersion and circular_sd2() return a positive number.

Note

If y is NULL, than the circular variance is returned.

References

Mardia, K.V. (1972). Statistics of Directional Data: Probability and Mathematical Statistics. London: Academic Press.

Mardia, K.V., and Jupp, P.E (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA. doi:10.1002/9780470316979

See Also

circular_mean(), circular_var().

Examples

a <- c(0, 2, 359, 6, 354)
circular_distance(a, 10) # distance to single value

b <- a + 90
circular_distance(a, b) # distance to multiple values

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
circular_dispersion(sa.por$azi.PoR, y = 135)
circular_dispersion(sa.por$azi.PoR, y = 135, w = weighting(san_andreas$unc))
circular_sd2(sa.por$azi.PoR, y = 135, w = weighting(san_andreas$unc))

Distance between points

Description

Returns the great circle distance between a location and all grid point in km

Usage

dist_greatcircle(
  lat1,
  lon1,
  lat2,
  lon2,
  r = earth_radius(),
  method = c("haversine", "orthodrome", "vincenty", "euclidean")
)

Arguments

Value

numeric vector with length equal to length(lat1) or length(lat2)

See Also

orthodrome(), haversine(), vincenty()

Examples

# Haversine: (4149.157, 2296.583) km
dist_greatcircle(lat1 = 20, lon1 = 12, lat2 = c(50, 30), lon2 = c(40, 32))

# Orthodrome: (4149.157, 2296.583) km
dist_greatcircle(
  lat1 = 20, lon1 = 12, lat2 = c(50, 30), lon2 = c(40, 32),
  method = "orthodrome"
)

# Vincenty: (4149.157, 2296.583) km
dist_greatcircle(
  lat1 = 20, lon1 = 12, lat2 = c(50, 30), lon2 = c(40, 32),
  method = "vincenty"
)

# Euclidean (4076.220, 2284.169) km
dist_greatcircle(
  lat1 = 20, lon1 = 12, lat2 = c(50, 30), lon2 = c(40, 32),
  method = "euclidean"
)

Distance Binned Summary Statistics

Description

Circular summary statistics over intervals of distances.

Usage

distance_binned_stats(
  azi,
  distance,
  n.breaks = 10,
  width.breaks = NULL,
  unc = NULL,
  prd = NULL,
  prd.error = NULL,
  kappa = 2,
  R = 1000,
  conf.level = 0.95,
  ...
)

Arguments

Value

tibble containing the n values for aziin each bin, min/median/max distance of the bin, and the summary statistics for azi. If prd is specified, the normal Chi-squared statistic, dispersion and its standard error are returned as well.

See Also

circular_summary(), circular_dispersion(), and circular_dispersion_boot()

Examples

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")
data("san_andreas")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
san_andreas$distance <- distance_from_pb(
  x = san_andreas,
  PoR = PoR,
  pb = plate_boundary,
  tangential = TRUE
)
dat <- san_andreas |> cbind(PoR_shmax(san_andreas, PoR, "right"))

distance_binned_stats(dat$azi.PoR,
  distance = dat$distance, width.breaks = 1,
  unc = dat$unc, prd = 135
) |> head()

Distance from plate boundary

Description

Absolute distance of data points from the nearest plate boundary

Usage

distance_from_pb(x, PoR, pb, tangential = FALSE, km = FALSE, ...)

Arguments

Details

The distance to the plate boundary is the longitudinal or latitudinal difference between the data point and the plate boundary (along the closest latitude or longitude) for inward/outward or tangential plate boundaries, respectively.

Value

Numeric vector of the great circle distances in units defined by km.

Note

Stresses emanate from the plate boundary along great circles, small circles or loxodromes associated with the pole of rotation. Hence the emanation distance is not necessarily the shortest distance to the plate boundary, which is measured along a great circle unrelated to the pole of rotation. The differences are particularly notable when the plate boundary is kinked or for convergent and divergent plate boundaries.

References

Wdowinski, S. (1998). A theory of intraplate tectonics. Journal of Geophysical Research: Solid Earth, 103(3), 5037<U+2013>5059. http://dx.doi.org/10.1029/97JB03390

Examples

data("nuvel1")
na_pa <- subset(nuvel1, nuvel1$plate.rot == "na")

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")

data("san_andreas")
res <- distance_from_pb(
  x = san_andreas, PoR = na_pa, pb = plate_boundary, tangential = TRUE
)
head(res)

res.km <- distance_from_pb(
  x = san_andreas, PoR = na_pa, pb = plate_boundary, tangential = TRUE, km = TRUE
)
range(res.km)

Normalize angular distance on a sphere distance

Description

Helper function to express angular distance on the sphere in the range of 0 to 180 °

Usage

distance_mod(x)

Arguments

Value

numeric vector

Plate Stress Dummy Grid

Description

Helper functions to create a dummy grid for small circles, great circles, and loxodromes of an Euler pole

Usage

smallcircle_dummy(n)

greatcircle_dummy(n)

loxodrome_dummy(n, angle, cw)

Arguments

Value

data.frame

Earth's radius in km

Description

IERS mean radius of Earth in km (based on WGS 84)

Usage

earth_radius()

Value

numeric value

Equivalent rotation

Description

Transforms a sequence of rotations into a new reference system

Usage

equivalent_rotation(x, fixed, rot)

Arguments

Value

sequence of plate rotations in new reference system. Same object class as x

See Also

relative_rotation()

Examples

data(nuvel1) # load the NUVEL1 rotation parameters

# all nuvel1 rotation equivalent to fixed Africa:
equivalent_rotation(nuvel1, fixed = "af")
# relative plate motion between Eurasia and India:
equivalent_rotation(nuvel1, "eu", "in") # lat = 24.58, lon = 18.07, angle = 0.528

Concentration parameter of von Mises distribution

Description

Estimates the concentration parameter of a von Mises distribution, given a set of angular measurements.

Usage

est.kappa.MLE(x, w = NULL, bias = FALSE)

est.kappa(x, w = NULL, p = 2)

Arguments

Details

est.kappa.MLE() is the maximum likelihood estimate for MLE for \kappa.

est.kappa() uses an approximation based on the empirical equation:

\kappa = \frac{\bar{R}(p-\bar{R}^2)}{1-\bar{R}^2}

where \bar{R} is the mean resultant length and p is the dimensionality of the data (2 for circular data).

Value

numeric. Concentration of a von Mises distribution

See Also

vonmises()

Examples

set.seed(123)
x <- rvm(100, 90, 10)
w <- weighting(runif(100, 0, 10))

est.kappa(x, w)

est.kappa.MLE(x, w)

Euler pole object

Description

Creates an object of the orientation of the Euler pole axis

Usage

euler_pole(x, y, z = NA, geo = TRUE, angle = NA)

Arguments

Value

An object of class "euler.pole" containing the Euler pole axis in both geographical and Cartesian coordinates and the angle of rotation in radians.

Examples

euler_pole(90, 0, angle = 45)
euler_pole(0, 0, 1, geo = FALSE)

Azimuth visualization

Description

geom_azimuth() visualizes axial-directional vector fields using a geom to produce a new graphical layer, which allows aesthetic options. This layer can be overlaid on a map to improve visualisation of mapped data. The geom draws line segments (spokes) centered at (x, y) with a given orientation (angle in degrees) and length (radius). By default the spoke is centered using PositionCenterSpoke, so that the given coordinates mark the middle of the line. The azimuths are given as angles in degrees increasing clockwise from North.

Usage

geom_azimuth(
  mapping = NULL,
  data = NULL,
  stat = "azimuth",
  center = TRUE,
  radius = NULL,
  na.rm = FALSE,
  show.legend = NA,
  inherit.aes = TRUE,
  ...
)

Arguments

Value

A ggplot2 layer that adds axis-like spokes.

Aesthetics

geom_azimuth() understands the following aesthetics (required aesthetics in bold):

• x

• y

• angle (in degrees, transformed internally)

• radius

• colour

• alpha

• linewidth

• linetype

See Also

ggplot2::geom_spoke(), geom_azimuthpoint()

Examples

set.seed(20250411)
df <- data.frame(
  x = runif(5), y = runif(5),
  angle_deg = rvm(5, mean = 90, kappa = 10),
  radius = runif(5, 0.1, 2)
)

if (require("ggplot2")) {
  ggplot(df, aes(x, y)) +
    geom_azimuth(aes(angle = angle_deg), radius = .1, linewidth = 1.2, colour = "blue")
  if (require("grid")) {
    ggplot(df, aes(x, y, radius = radius)) +
      geom_azimuth(aes(angle = angle_deg), center = FALSE, colour = "red", arrow = grid::arrow())
  }
}

Azimuth + point visualization

Description

geom_azimuthpoint() draws line segments (spokes) like geom_azimuth(), but also places a point (marker) at the spoke's center ⁠(x, y)⁠.

Aesthetic rules:

• linewidth, linetype affect the spoke only

• shape affects the point only

• colour, alpha affect both spoke and point

• size sets the size of the point only

Usage

geom_azimuthpoint(
  mapping = NULL,
  data = NULL,
  stat = "identity",
  center = TRUE,
  na.rm = FALSE,
  show.legend = NA,
  inherit.aes = TRUE,
  size = 2,
  ...
)

Arguments

Value

A list of ggplot2 layers (spokes + points).

Aesthetics

geom_azimuthpoint() understands the following aesthetics (required aesthetics in bold):

• x

• y

• angle (in degrees, transformed internally; spoke only)

• radius (spoke only)

• colour (shared)

• alpha (shared)

• linewidth (spoke only)

• linetype (spoke only)

• shape (point only)

• size (point only, or via argument)

• fill (point only, for shapes that accept fill)

See Also

geom_azimuth(), ggplot2::geom_spoke(), ggplot2::geom_point()

Examples

set.seed(20250411)
df <- data.frame(
  x = runif(5), y = runif(5),
  angle_deg = rvm(5, mean = 90, kappa = 10),
  radius = runif(5, 0.1, 1),
  group = rep(1:2, length.out = 5)
)

if (require("ggplot2")) {
  ggplot(data = df, aes(x, y, angle = angle_deg, radius = radius)) +
    geom_azimuthpoint(aes(colour = factor(group), shape = factor(group)),
      linewidth = 1.1, linetype = "dashed",
      size = 3, alpha = 0.8
    )
}

Azimuth Between two Points

Description

Calculate initial bearing (or forward azimuth/direction) to go from point a to point b following great circle arc on a sphere.

Usage

get_azimuth(lat_a, lon_a, lat_b, lon_b)

Arguments

Details

get_azimuth() is based on the spherical law of tangents. This formula is for the initial bearing (sometimes referred to as forward azimuth) which if followed in a straight line along a great circle arc will lead from the start point a to the end point b.

\theta = \arctan2 (\sin \Delta\lambda \cos\psi_2, \cos\psi_1 \sin\psi_1-\sin\psi_1 \cos\psi_2 \cos\Delta\lambda)

where \psi_1, \lambda_1 is the start point, \psi_2, \lambda_2 the end point (\Delta\lambda is the difference in longitude).

Value

numeric. Azimuth in degrees

References

http://www.movable-type.co.uk/scripts/latlong.html

Examples

berlin <- c(52.517, 13.4) # Berlin
tokyo <- c(35.7, 139.767) # Tokyo
get_azimuth(berlin[1], berlin[2], tokyo[1], tokyo[2]) # 41.57361

Helper function to Distance from plate boundary

Description

Helper function to Distance from plate boundary

Usage

get_distance(lon, lat, pb.coords, tangential, km)

Arguments

See Also

distance_from_pb()

Helper function to get Distance from plate boundary

Description

Helper function to get Distance from plate boundary

Usage

get_projected_pb_strike(lon, lat, pb.coords, pb.bearing, tangential)

Arguments

See Also

projected_pb_strike()

Helper function to Equivalent rotation

Description

Helper function to Equivalent rotation

Usage

get_relrot(plate.rot, lat, lon, angle, fixed, fixed.ep)

Arguments

See Also

equivalent_rotation()

Homing directions of Pidgeons

Description

Directions of 10 pidgeons, example data from Mardia and Jupp (1999), pp. 93

Format

numeric vector of length 10

References

Mardia and Jupp (1999). Directional Statistics. John Wiley and Sons.

Examples

data("homing")
homing

World Stress Map Database (WSM)

Description

Download WSM2025 or WSM2016 database from the GFZ sever and applies optional filters. If destdir is specified, the data can be reloaded in a later R session using load_WSM() using the same arguments.

Usage

download_WSM(
  destdir = tempdir(),
  load = TRUE,
  version = c("2025", "2016"),
  ...
)

load_WSM(
  file,
  quality = c("A", "B", "C", "D", "E", "X"),
  lat_range = c(-90, 90),
  lon_range = c(-180, 180),
  depth_range = c(-Inf, Inf),
  type = c("BO", "BOC", "BOT", "BS", "DIF", "FMA", "FMF", "FMS", "GFI", "GFM", "GFS",
    "GVA", "HF", "HFG", "HFM", "HFH", "HFP", "HFS", "OC", "PC", "SWB", "SWL", "SWS"),
  regime = c("N", "NS", "T", "TS", "S", NA)
)

download_WSM2016(destdir = tempdir(), load = TRUE, ...)

load_WSM2016(
  file,
  quality = c("A", "B", "C", "D", "E"),
  lat_range = c(-90, 90),
  lon_range = c(-180, 180),
  depth_range = c(-Inf, Inf),
  type = c("BO", "BOC", "BOT", "BS", "DIF", "FMA", "FMF", "FMS", "GFI", "GFM", "GFS",
    "GVA", "HF", "HFG", "HFM", "HFP", "OC", "PC", "SWB", "SWL", "SWS"),
  regime = c("N", "NS", "T", "TS", "S", NA)
)

Arguments

Value

sf object of and the parsed numeric uncertainty (unc) based on the reported standard deviation and the quality rank. If load=FALSE, the path to the downloaded file is returned.

Note

Because of R-compatibility and easy readability reasons, the downloaded dataset is a modified version of the original, WSM server version: All column names have been changed from uppercase (in the original dataset) to lowercase characters. Unknown azimuth values are represented by NA values instead of 999 in the original. Unknown regimes are represented by NA instead of "U" in the original.

Source

https://datapub.gfz.de/download/10.5880.WSM.2025.001-Scbwez/WSM_Database_2025.csv

https://datapub.gfz-potsdam.de/download/10.5880.WSM.2016.001/wsm2016.csv

References

Heidbach, O., M. Rajabi, X. Cui, K. Fuchs, B. M<U+00FC>ller, J. Reinecker, K. Reiter, M. Tingay, F. Wenzel, F. Xie, M. O. Ziegler, M.-L. Zoback, and M. D. Zoback (2018): The World Stress Map database release 2016: Crustal stress pattern across scales. Tectonophysics, 744, 484-498, doi:10.1016/j.tecto.2018.07.007.

Heidbach, Oliver; Rajabi, Mojtaba; Di Giacomo, Domenico; Harris, James; Lammers, Steffi; Morawietz, Sophia; Pierdominici, Simona; Reiter, Karsten; von Specht, Sebastian; Storchak, Dmitry; Ziegler, Moritz O. (2025): World Stress Map Database Release 2025. GFZ Data Services. doi:10.5880/WSM.2025.001

Examples

## Not run: 
download_WSM(
  quality = c("A", "B", "C"), lat_range = c(51, 72),
  lon_range = c(-180, -130), depth_range = c(0, 10), type = "FMS"
)

## End(Not run)

Check if object is euler.pole

Description

Check if object is euler.pole

Usage

is.euler(x)

Arguments

Value

logical

Adaptive Kernel Dispersion

Description

Stress field and wavelength analysis using circular dispersion (or other statistical estimators for dispersion)

Usage

kernel_dispersion(
  x,
  stat = c("dispersion", "nchisq", "rayleigh"),
  grid = NULL,
  lon_range = NULL,
  lat_range = NULL,
  gridsize = 2.5,
  min_data = 3L,
  max_data = Inf,
  min_dist_threshold = 200,
  dist_threshold = 0.1,
  stat_threshold = Inf,
  R_range = seq(100, 2000, 100),
  ...
)

dispersion_grid(...)

Arguments

Value

sf object containing

lon,lat

longitude and latitude in degree

stat

output of function defined in stat

R

The rearch radius in km.

mdr

Mean distance of datapoints per search radius

N

Number of data points in search radius

Note

dispersion_grid() was renamed to kernel_dispersion() to create a more consistent API.

See Also

circular_dispersion(), norm_chisq(), weighted_rayleigh()

Examples

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
san_andreas_por <- data2PoR(san_andreas, PoR)
san_andreas_por$prd <- 135
kernel_dispersion(san_andreas_por) |> head()

Kuiper Test of Circular Uniformity

Description

A statistical test used to determine whether a set of angular or circular data points (such as times of day, compass directions, or degrees) are spread out evenly around a circle or if they cluster in some way.

Usage

kuiper_test(x, alpha = 0, axial = TRUE, quiet = FALSE)

Arguments

Details

The Null Hypothesis (H_0): The data are distributed completely uniformly (randomly and evenly) around the circle.

The Alternative Hypothesis (H_1): The data are not uniform and show a preference, clustering, or pattern somewhere on the circle.

The Test Statistic (V or D^{+} + D^{-}): It measures the greatest positive and negative differences between your data's empirical cumulative distribution and a theoretical uniform distribution.

Interpreting the Results

• High Test Statistic / Low p-value (p < \alpha, typically 0.05): You reject the null hypothesis. This means your data are not uniform; they have a significant preferred direction, grouping, or non-random pattern on the circle.

• Low Test Statistic / High p-value (p \ge 0.05): You fail to reject the null hypothesis. There is no strong evidence to say the data are different from a flat, uniform distribution. The points appear random across the circle.

Value

list containing the test statistic statistic and the significance level p.value.

Note

Kuiper's test statistic is a rotation-invariant Kolmogorov-type test statistic. The critical values of a modified Kuiper's test statistic are used according to the tabulation given in Stephens (1970).

See Also

Other Tests: ar_test(), norm_chisq(), rayleigh-test, watson_test(), watson_two_sample, watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

# Example data from Mardia and Jupp (1999), pp. 93
kuiper_test(homing, alpha = .05)

# San Andreas Fault Data:
data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
kuiper_test(sa.por$azi.PoR, alpha = .05)

Extract azimuths of line segments

Description

Extract azimuths of line segments

Usage

line_azimuth(x, warn = TRUE)

lines_azimuths(x)

Arguments

Details

It is recommended to perform line_azimuth() on single line objects, i.e. type "LINESTRING", instead of "MULTILINESTRING". This is because the azimuth of the last point of a line will be calculated to the first point of the next line otherwise. This will cause a warning message (if warn = TRUE). For "MULTILINESTRING" objects, use lines_azimuths().

Value

sf object of type "POINT" with the columns and entries of the first row of x

Examples

data("plates")

# one line:
subset(plates, pair == "af-eu") |>
  smoothr::densify() |>
  line_azimuth() |>
  head()

# multiple lines:
lines_azimuths(plates[1:5, ]) |> head()

Mean Cosine and Sine

Description

Mean Cosine and Sine

Usage

mean_SC(x, w = NULL, na.rm = TRUE)

Arguments

Value

named two element vector

Examples

## Not run: 
set.seed(1)
x <- rvm(100, 0, 5)
mean_SC(x)

## End(Not run)

Mean Resultant Length

Description

Measure of spread around the circle. It should be noted that: If R=0, then the data is completely spread around the circle. If R=1, the data is completely concentrated on one point.

Usage

mean_resultant_length(x, w = NULL, na.rm = TRUE)

Arguments

Value

numeric.

References

Mardia, K.V. (1972). Statistics of Directional Data: Probability and Mathematical Statistics. London: Academic Press.

Examples

# Example data from Davis (1986), pp. 316
finland_stria <- c(
  23, 27, 53, 58, 64, 83, 85, 88, 93, 99, 100, 105, 113,
  113, 114, 117, 121, 123, 125, 126, 126, 126, 127, 127, 128, 128, 129, 132,
  132, 132, 134, 135, 137, 144, 145, 145, 146, 153, 155, 155, 155, 157, 163,
  165, 171, 172, 179, 181, 186, 190, 212
)
mean_resultant_length(finland_stria, w = NULL, na.rm = FALSE) # 0.800

Theoretical Direction of Maximum Horizontal Stress in the geographical reference system.

Description

Models the direction of maximum horizontal stress \sigma_\text{Hmax} along great circles, small circles, and loxodromes at a given point or points according to the relative plate motion in the geographical coordinate reference system.

Usage

model_shmax(df, euler)

Arguments

Details

\sigma_\text{Hmax} following great circles is the (initial) bearing between the given point and the pole of relative plate motion. \sigma_\text{Hmax} along small circles, clockwise, and counter-clockwise loxodromes is 90°, +45°, and +135° (or -45°) to this great circle bearing, respectively.

Value

data.frame

gc

Azimuth of modeled \sigma_{Hmax} following great circles

sc

Small circles

ld.cw

Clockwise loxodromes

ld.ccw

Counter-clockwise loxodromes

Author(s)

Tobias Stephan

References

Stephan, T., Enkelmann, E., and Kroner, U. "Analyzing the horizontal orientation of the crustal stress adjacent to plate boundaries". Sci Rep 13. 15590 (2023). doi:10.1038/s41598-023-42433-2.

See Also

deviation_shmax() to compute the deviation of the modeled direction from the observed direction of \sigma_\text{Hmax}. PoR_shmax() to calculate the azimuth of \sigma_\text{Hmax} in the pole of rotation reference system.

Examples

data("nuvel1")
# North America relative to Pacific plate:
euler <- subset(nuvel1, nuvel1$plate.rot == "na")

# the point where we mant to model the SHmax direction:
point <- data.frame(lat = 45, lon = 20)

model_shmax(point, euler)

Normalized Chi-Squared Test for Circular Data

Description

A quantitative comparison between the predicted and observed directions of \sigma_\text{Hmax} is obtained by the calculation of the average azimuth and by a normalized \chi^2 test (Wdowsinki, 1998)

Usage

norm_chisq(obs, prd, unc)

Arguments

Details

The normalized \chi^2 test is

\text{Norm} \chi^2_i = \frac{ \sum^M_{i = 1} \left( \frac{\alpha_i - \alpha_{{predict}}}{\sigma_i} \right) ^2} {\sum^M_{i = 1} \left( \frac{90}{\sigma_i} \right) ^2 }

The value of the chi-squared test statistic is a number between 0 and 1 indicating the quality of the predicted \sigma_\text{Hmax} directions. Low values (\le 0.15) indicate good agreement, high values (> 0.7) indicate a systematic misfit between predicted and observed \sigma_\text{Hmax} directions.

Value

Numeric vector

References

Wdowinski, S., 1998, A theory of intraplate tectonics. Journal of Geophysical Research: Solid Earth, 103, 5037-5059, doi: 10.1029/97JB03390.

See Also

Other Tests: ar_test(), kuiper_test(), rayleigh-test, watson_test(), watson_two_sample, watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to
# Pacific plate
data(san_andreas)
point <- data.frame(lat = 45, lon = 20)
prd <- model_shmax(point, PoR)
norm_chisq(obs = c(50, 40, 42), prd = prd$sc, unc = c(10, NA, 5))

data(san_andreas)
prd2 <- PoR_shmax(san_andreas, PoR, type = "right")
norm_chisq(obs = prd2$azi.PoR, 135, unc = san_andreas$unc)

NUVEL-1 Global model of current plate motions

Description

NNR-NUVEL-1 global model of current plate motions by DeMets et al. 1990

Usage

data('nuvel1')

Format

An object of class data.frame

plate.name

The rotating plate

plate.rot

The abbreviation of the plate's name

lat,lon

Coordinates of the Pole of Rotation

angle

The amount of rotation (angle in 1 Myr)

plate.fix

The anchored plate, i.e. plate.rot moves relative to plate.fix

source

Reference to underlying study

References

DeMets, C., Gordon, R. G., Argus, D. F., & Stein, S. (1990). Current plate motions. Geophysical Journal International, 101(2), 425-478. doi:10.1111/j.1365-246X.1990.tb06579.x.

Examples

data("nuvel1")
head("nuvel1")

Plate Boundaries on the Earth

Description

Global set of present plate boundaries on the Earth based on NUVEL-1 model by DeMets et al. 1990

Usage

data('nuvel1_plates')

Format

An object of class sf

References

DeMets, C., Gordon, R. G., Argus, D. F., & Stein, S. (1990). Current plate motions. Geophysical Journal International, 101(2), 425-478. doi:10.1111/j.1365-246X.1990.tb06579.x.

Examples

data("nuvel1_plates")
head("nuvel1_plates")

Decomposition of Orientation Tensor in 2D

Description

Spectral decomposition of the 2D orientation tensor into two Eigenvectors and corresponding Eigenvalues provides provides a measure of location and a corresponding measure of dispersion, respectively.

Usage

ot_eigen2d(x, w = NULL, scale = FALSE)

principal_direction(x, w = NULL)

axial_strength(x, w = NULL)

axial_dispersion(x, w = NULL)

Arguments

Details

The Eigenvalues (\lambda_1 > \lambda_2) can be interpreted as the fractions of the variance explained by the orientation of the associated Eigenvectors. The two perpendicular Eigenvectors (a_1, a_2) are the "principal directions" with respect to the highest and the lowest concentration of orientation data.

The strength of the orientation is the largest eigenvalue \lambda_1 normalized by the sum of the eigenvalues (scale=TRUE). Then \lambda_2 = 1-\lambda_1 is a measure of dispersion of 2D orientation data with respect to a_1.

Value

ot_eigen2d returns a list of the Eigenvalues and the axial angles corresponding to the Eigenvectors. principal_direction(), axial_strength() and axial_dispersion() are convenience functions to return the orientation of the largest eigenvalue, the orientation strength, the axial dispersion respectively.

Note

Eigenvalues and Eigenvectors of the orientation tensor (inertia tensor) are also called "principle moments of inertia" and "principle axes of inertia", respectively.

See Also

ortensor2d()

Examples

test <- rvm(100, mean = 0, k = 10) / 2
ot_eigen2d(test)

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
sa_eig <- ot_eigen2d(sa.por$azi.PoR, w = weighting(san_andreas$unc), scale = TRUE)
print(sa_eig)

rose(sa.por$azi.PoR, muci = FALSE)
rose_line(sa_eig$vectors,
  col = c("red", "green"),
  radius = sa_eig$values, lwd = 2
)
graphics::legend("topright",
  legend = round(sa_eig$values, 2),
  col = c("red", "green"), lty = 1
)

principal_direction(sa.por$azi.PoR)

axial_strength(sa.por$azi.PoR)

axial_dispersion(sa.por$azi.PoR)

Orientation Tensor

Description

2D orientation tensor characterizes distribution of axial angles using the Eigenvalue method (Watson 1966, Scheidegger 1965).

Usage

ortensor2d(x, w = NULL, norm = FALSE)

Arguments

Details

The moment of inertia can be minimized by calculating the Cartesian coordinates of the orientation data, and calculating their covariance matrix. This yields

I = x \cdot x^\intercal

where x is the Cartesian vector of the orientations. Orientation tensor T and the inertia tensor I are related by

I = E - T

where E denotes the unit matrix, so that

T = \frac{1}{n} \sum_{i=i}^{n} x_i \cdot x_i^\intercal

Value

2x2 matrix

References

Watson, G. S. (1966). The Statistics of Orientation Data. The Journal of Geology, 74(5), 786–797.

Scheidegger, A. E. (1964). The tectonic stress and tectonic motion direction in Europe and Western Asia as calculated from earthquake fault plane solutions. Bulletin of the Seismological Society of America, 54(5A), 1519–1528. doi:10.1785/BSSA05405A1519

Bachmann, F., Hielscher, R., Jupp, P. E., Pantleon, W., Schaeben, H., & Wegert, E. (2010). Inferential statistics of electron backscatter diffraction data from within individual crystalline grains. Journal of Applied Crystallography, 43(6), 1338–1355. https://doi.org/10.1107/S002188981003027X

See Also

ot_eigen2d()

Examples

test <- rvm(100, mean = 0, k = 10)
ortensor2d(test)

Numerical values to World Stress Map Quality Ranking

Description

Assigns numeric values of the precision (sd.) of each measurement to the categorical quality ranking of the World Stress Map (A, B, C, D, E, X).

Usage

parse_wsm_quality(x)

quantise_wsm_quality(x)

Arguments

Value

"numeric". the standard deviation of stress azimuth

References

Heidbach, O., Barth, A., M<U+00FC>ller, B., Reinecker, J., Stephansson, O., Tingay, M., Zang, A. (2016). WSM quality ranking scheme, database description and analysis guidelines for stress indicator. World Stress Map Technical Report 16-01, GFZ German Research Centre for Geosciences. doi:10.2312/wsm.2016.001

Examples

parse_wsm_quality(c("A", "B", "C", "D", NA, "E", "X"))
data("san_andreas")
head(parse_wsm_quality(san_andreas$quality))

Global model of current plate motions

Description

PB2002 global model of current plate motions by Bird 2003

Usage

data('pb2002')

Format

An object of class data.frame

plate.name

The rotating plate

plate.rot

The abbreviation of the plate's name

lat,lon

Coordinates of the Pole of Rotation

angle

The amount of rotation (angle in 1 Myr)

plate.fix

The anchored plate, i.e. plate.rot moves relative to plate.fix

source

Reference to underlying study

References

Bird, P. (2003), An updated digital model of plate boundaries, Geochem. Geophys. Geosyst., 4, 1027, doi: 10.1029/2001GC000252, 3.

Examples

data("pb2002")
head("pb2002")

Plate Boundaries on the Earth

Description

Global set of present plate boundaries on the Earth based on PB2002 model by Bird (2003). Contains the plate boundary displacement types such as inward, outward, or tangentially displacement.

Usage

data('plates')

Format

An object of class sf

References

Bird, P. (2003), An updated digital model of plate boundaries, Geochem. Geophys. Geosyst., 4, 1027, doi: 10.1029/2001GC000252, 3.

Examples

data("plates")
head("plates")

Circular Kernel Density Plot

Description

Plots multiples of a von Mises, wrapped Cauchy, and wrapped Normal density distribution in a circular plot

Usage

plot_density(
  x,
  bw = NULL,
  kernel = c("vonmises", "wrappedcauchy", "wrappednormal"),
  weights = NULL,
  axial = TRUE,
  n = 512L,
  norm.density = TRUE,
  kappa = NULL,
  rho = NULL,
  sd = NULL,
  c = NULL,
  fill = FALSE,
  scale = 0,
  shrink = 1,
  add = TRUE,
  main = NULL,
  labels = TRUE,
  at = seq(0, 360 - 45, 45),
  cborder = TRUE,
  grid = FALSE,
  ...
)

Arguments

Value

plot or calculated densities

See Also

circular_density()

Other circular-plot: plot_points(), rose(), rose_geom, rose_stats()

Examples

# Filled von Mises kernel density curve inside the plot
plot_density(san_andreas$azi,
  kappa = 100,
  fill = TRUE, col = "#51127C80", border = "#51127CFF",
  grid = TRUE,
  add = FALSE
)

# Superimpose a wrapped Cauchy kernel distribution curve
plot_density(san_andreas$azi,
  rho = 0.9, kernel = "wrappedcauchy",
  fill = FALSE, col = "#FB8861FF",
  add = TRUE
)

# Superimpose a wrapped Normal kernel distribution curve
plot_density(san_andreas$azi,
  sd = 2, kernel = "wrappednormal",
  fill = FALSE, col = "#E65164FF",
  add = TRUE
)

# Superimpose a von Mises kernel density curve on a rose diagram:
rose(san_andreas$azi, grid = TRUE)
plot_density(san_andreas$azi,
  bw = 100, col = "#51127CFF",
  add = TRUE, lwd = 3
)

# Corona plot (density curve outside of a rose diagram plot):
w <- weighting(san_andreas$unc)
rose(san_andreas$azi, weights = w, dots = TRUE, stack = TRUE, dot_cex = 0.5, dot_pch = 21)
plot_density(san_andreas$azi, weights = w,
  bw = 100,
  scale = 1.1, shrink = 3, xpd = NA,
  col = "#51127CFF"
)

Add Points to a Circular Plot

Description

Add points to a plot of circular data points on the current graphics device.

Usage

plot_points(
  x,
  axial = TRUE,
  stack = FALSE,
  binwidth = 1,
  cex = 1,
  sep = 0.025,
  jitter_factor = 0,
  ...,
  scale = 1.1,
  add = TRUE,
  main = NULL,
  labels = TRUE,
  at = seq(0, 360 - 45, 45),
  cborder = TRUE
)

Arguments

Value

A list with information on the plot

See Also

Other circular-plot: plot_density(), rose(), rose_geom, rose_stats()

Examples

x <- rvm(100, mean = 90, k = 5)

# plot poinit without jitter
plot_points(x, add = FALSE)

# with some jitter
plot_points(x, jitter_factor = .2, add = FALSE)

# stacked dots:
plot_points(x, stack = TRUE, binwidth = 3, add = FALSE, xpd = TRUE)

Conversion between spherical PoR to geographical coordinate system

Description

Transformation from spherical PoR to geographical coordinate system and vice versa

Usage

geographical_to_PoR(x, PoR)

PoR_to_geographical(x, PoR)

Arguments

Value

object of same type of x with the transformed coordinates. If x was a data.frame, transformed coordinates are named lat.PoR and lon.PoR for PoR CRS, or lat and lon for geographical CRS).

Examples

data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to Pacific plate
data("san_andreas")
san_andreas.por <- geographical_to_PoR(san_andreas, por)
head(san_andreas.por)
head(PoR_to_geographical(san_andreas.por, por))

Conversion between spherical PoR to geographical coordinate system of data.frames

Description

Transformation from spherical PoR to geographical coordinate system and vice versa

Usage

geographical_to_PoR_df(x, PoR)

PoR_to_geographical_df(x, PoR)

Arguments

Value

"data.frame" with the transformed coordinates (lat.PoR and lon.PoR for PoR CRS, or lat and lon for geographical CRS).

Conversion between PoR to geographical coordinate system using quaternions

Description

Helper function for the transformation from PoR to geographical coordinate system or vice versa

Usage

geographical_to_PoR_quat(x, PoR)

PoR_to_geographical_quat(x, PoR)

Arguments

Value

two-element numeric vector

Conversion between PoR to geographical coordinates of sf data

Description

Transform spatial objects from PoR to geographical coordinate reference system and vice versa.

Usage

PoR_to_geographical_sf(x, PoR)

geographical_to_PoR_sf(x, PoR)

Arguments

Details

The PoR coordinate reference system is oblique transformation of the geographical coordinate system with the Euler pole coordinates being the translation factors.

Value

sf or SpatRast object of the data points in the transformed geographical or PoR coordinate system

Error of Model's Prediction

Description

The maximum error in the model's predicted azimuth given the Pole of rotations uncertainty and distance of the data point to the pole.

Usage

prd_err(dist_PoR, sigma_PoR = 1)

Arguments

Value

numeric vector. The maximum error for azimuths prediction (in degree)

References

Ramsay, J.A. Folding and fracturing of rocks. McGraw-Hill, New York, 1967.

See Also

PoR_shmax() and model_shmax() for the model's prediction, and orthodrome() for great circle distances.

Examples

prd_err(67, 1)

# San Andreas example:
data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to Pacific plate
data("san_andreas")
d <- PoR_distance(san_andreas, por)
prd_err(d)

Strike of the plate boundary projected on data point

Description

The fault's strike in the PoR CRS projected on the data point along the predicted stress trajectories.

Usage

projected_pb_strike(x, PoR, pb, tangential = FALSE, ...)

Arguments

Details

Useful to calculate the beta angle, i.e. the angle between \sigma_\text{Hmax} direction (in PoR CRS!) and the fault's strike (in PoR CRS). The beta angle is the same in geographical and PoR coordinates.

Value

Numeric vector of the strike direction of the plate boundary (in degree)

Note

The algorithm calculates the great circle bearing between line vertices. Since transform plate boundaries represent small circle lines in the PoR system, this great-circle azimuth is only a approximation of the true (small-circle) azimuth.

Examples

data("nuvel1")
na_pa <- subset(nuvel1, nuvel1$plate.rot == "na")

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")

data("san_andreas")
res <- projected_pb_strike(
  x = san_andreas, PoR = na_pa, pb = plate_boundary, tangential = TRUE
)
head(res)
head(san_andreas$azi - res) # beta angle

Plotting Stress Analysis Results

Description

Creates a set of plots including the azimuth as a function of the distance to the plate boundary, the Norm Chi-squared as a function of the distance to the plate boundary, the circular distance (and dispersion) a function of the distance to the plate boundary, a von-Mises QQ plot, and a rose diagram of the quality-weighted frequency distribution of the azimuths.

Usage

quick_plot(azi, distance, prd, unc = NULL, regime, width = 51)

Arguments

Details

Plot 1 shows the transformed azimuths as a function of the distance to the plate boundary. The red line indicates the rolling circular mean, stippled red lines indicate the 95% confidence interval about the mean.

Plot 2 shows the normalized \chi^2 statistics as a function of the distance to the plate boundary. The red line shows the rolling \chi^2 statistic.

Plot 3 shows the circular distance of the transformed azimuths to the predicted azimuth, as a function of the distance to the plate boundary. The red line shows the rolling circular dispersion about the prediction.

Plot 4 give the rose diagram of the transformed azimuths.

Value

four R base plots

See Also

PoR_shmax(), distance_from_pb(), circular_mean(), circular_dispersion(), confidence_interval_fisher(), norm_chisq(), weighted_rayleigh(), vm_qqplot()

Examples

data("nuvel1")
na_pa <- subset(nuvel1, nuvel1$plate.rot == "na")

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")

data("san_andreas")
res <- PoR_shmax(san_andreas, na_pa, "right")
d <- distance_from_pb(san_andreas, na_pa, plate_boundary, tangential = TRUE)
quick_plot(res$azi.PoR,
  distance = d, prd = res$prd, unc = san_andreas$unc,
  regime = san_andreas$regime
)

Conversion between PoR to geographical coordinate reference system of raster data

Description

Helper function to transform raster data set from PoR to geographical coordinates

Usage

geographical_to_PoR_raster(x, PoR)

PoR_to_geographical_raster(x, PoR)

Arguments

Value

terra "SpatRaster" object

Rayleigh Test of Circular Uniformity

Description

A test to determine whether a sample of circular or directional data is evenly spread out or clustered around a single specific direction. The test assesses the significance of the mean resultant length. rayleight_test_perm() uses permutation to estimate p-values.

Usage

rayleigh_test(x, mu = NULL, axial = TRUE, alpha = 0.05, quiet = FALSE)

rayleigh_test_perm(x, mu = NULL, axial = TRUE, n_perm = 1000L)

Arguments

Details

Hypotheses

Null Hypothesis H_0: The population is distributed uniformly (randomly) around the circle with no preferred direction.

Alternative Hypothesis H_1: The population is not uniform and has a unimodal (single-peaked) concentration in a preferred direction. When mu is specified), angles are non-uniformly distributed around the specified direction.

Mean Resultant Length (\bar{R} or R): A value between 0 and 1 that measures how concentrated the data points are.

• \bar{R} = 0: The data is completely spread out around the circle

• \bar{R} = 1: All data points point in the exact same direction.

p-value The probability of seeing data this clustered purely by chance under the assumption of uniformity.

Interpretation

• Small p-value (p < 0.05): Reject the null hypothesis. The length of the mean resultant differs significantly from zero, and the angles are not randomly distributed. You have strong evidence that the data points share a significant preferred or mean direction (unimodal clustering).

• Large p-value (p \ge 0.05): Fail to reject the null hypothesis. There is not enough evidence to claim a preferred direction, meaning the data looks random or uniform around the circle.

Value

a list with the components:

R or C

mean resultant length or the dispersion (if mu is specified). Small values of R (large values of C) will reject uniformity. Negative values of C indicate that vectors point in opposite directions (also lead to rejection).

statistic

test statistic

p.value

significance level of the test statistic

Note

Although the Rayleigh test is consistent against (non-uniform) von Mises alternatives, it is not consistent against alternatives with p = 0 (in particular, distributions with antipodal symmetry, i.e. axial data). Tests of non-uniformity which are consistent against all alternatives include Kuiper's test (kuiper_test()) and Watson's U^2 test (watson_test()).

Limitations

• The test assumes a unimodal alternative (one main peak).

• If your data has two opposite clusters (bimodal or axial data, like a bi-directional line trend), the Rayleigh test can yield a high/non-significant p-value because the opposing vectors cancel each other out.

References

Fisher, N. I. (1993) Statistical Analysis of Circular Data, Cambridge University Press.

See Also

mean_resultant_length(), circular_mean()

Other Tests: ar_test(), kuiper_test(), norm_chisq(), watson_test(), watson_two_sample, watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

# Example data from Mardia and Jupp (1999), pp. 93
rayleigh_test(homing, axial = FALSE) # Do not reject null hypothesis.
# R = 0.22; stat = 0.497, p = 0.62
rayleigh_test_perm(homing, axial = FALSE)

# Example data from Davis (1986), pp. 316
rayleigh_test(striae, axial = FALSE) # reject null hypothesis
rayleigh_test_perm(striae, axial = FALSE) # reject null hypothesis

rayleigh_test(striae, mu = 105, axial = FALSE) # reject null hypothesis
rayleigh_test_perm(striae, mu = 105, axial = FALSE) # reject null hypothesis

# Example data from Mardia and Jupp (1999), pp. 99
atomic_weight <- c(
  rep(0, 12), rep(3.6, 1), rep(36, 6), rep(72, 1),
  rep(108, 2), rep(169.2, 1), rep(324, 1)
)
rayleigh_test(atomic_weight, 0, axial = FALSE) # reject null hypothesis

# San Andreas Fault Data:
data(san_andreas)
rayleigh_test(san_andreas$azi) # reject null hypothesis
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
rayleigh_test(sa.por$azi.PoR, mu = 135) # reject null hypothesis
rayleigh_test_perm(sa.por$azi.PoR, mu = 135, n_perm = 1e3) # reject null hypothesis

Relative rotation between two rotations

Description

Calculates the relative rotation between two rotations, i.e. the difference from rotation 1 to rotation 2.

Usage

relative_rotation(r1, r2)

Arguments

Value

list. Euler axes (geographical coordinates) and Euler angles (in degrees)

References

Schaeben, H., Kroner, U. and Stephan, T. (2021). Euler Poles of Tectonic Plates. In B. S. Daza Sagar, Q. Cheng, J. McKinley and F. Agterberg (Eds.), Encyclopedia of Mathematical Geosciences. Encyclopedia of Earth Sciences Series (pp. 1–7). Springer Nature Switzerland AG 2021. doi: 10.1007/978-3-030-26050-7_435-1.

See Also

euler_pole() for class "euler.pole"

Examples

a <- euler_pole(90, 0, angle = 45)
b <- euler_pole(0, 0, 1, geo = FALSE, angle = -15)
relative_rotation(a, b) # axis: -90, -180; angle: 60
relative_rotation(b, a) # axis: 90, 0; angle: 60

Apply Rolling Functions using Circular Statistics

Description

A generic function for applying a function to rolling margins of an array.

Usage

roll_circstats(
  x,
  w = NULL,
  FUN,
  axial = TRUE,
  na.rm = TRUE,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

Arguments

Value

numeric vector with the results of the rolling function.

Note

If the rolling statistics are applied to values that are a function of distance it is recommended to sort the values first.

Examples

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")
data("san_andreas")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
distance <- distance_from_pb(
  x = san_andreas,
  PoR = PoR,
  pb = plate_boundary,
  tangential = TRUE
)
dat <- san_andreas[order(distance), ]
roll_circstats(dat$azi, w = 1 / dat$unc, circular_mean, width = 51) |> head()

Apply Rolling Functions using Circular Statistical Tests for Uniformity

Description

A generic function for applying a function to rolling margins of an array.

Usage

roll_normchisq(
  obs,
  prd,
  unc = NULL,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

roll_rayleigh(
  obs,
  prd,
  unc = NULL,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

roll_dispersion(
  x,
  y,
  w = NULL,
  w.y = NULL,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

roll_confidence(
  x,
  conf.level = 0.95,
  w = NULL,
  axial = TRUE,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

roll_dispersion_CI(
  x,
  y,
  w = NULL,
  w.y = NULL,
  R,
  conf.level = 0.95,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

roll_dispersion_sde(
  x,
  y,
  w = NULL,
  w.y = NULL,
  R,
  conf.level = 0.95,
  width = NULL,
  by.column = FALSE,
  partial = TRUE,
  fill = NA,
  ...
)

Arguments

Value

numeric vector with the test statistic of the rolling test. roll_dispersion_CI returns a 2-column matrix with the lower and the upper confidence limits

Note

If the rolling functions are applied to values that are a function of distance it is recommended to sort the values first.

Examples

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")
data("san_andreas")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
distance <- distance_from_pb(
  x = san_andreas,
  PoR = PoR,
  pb = plate_boundary,
  tangential = TRUE
)
dat <- san_andreas[order(distance), ]
dat.PoR <- PoR_shmax(san_andreas, PoR, "right")
roll_normchisq(dat.PoR$azi.PoR, 135, dat$unc) |> head()
roll_rayleigh(dat.PoR$azi.PoR, prd = 135, unc = dat$unc) |> head()
roll_dispersion(dat.PoR$azi.PoR, y = 135, w = 1 / dat$unc) |> head()
roll_confidence(dat.PoR$azi.PoR, w = 1 / dat$unc) |> head()

roll_dispersion_CI(dat.PoR$azi.PoR, y = 135, w = 1 / dat$unc, R = 10) |> head()

Apply Rolling Functions using Circular Statistics

Description

A generic function for applying a function to rolling margins of an array along an additional value.

Usage

distroll_circstats(
  x,
  distance,
  FUN,
  width = NULL,
  min_n = 2,
  align = c("right", "center", "left"),
  w = NULL,
  sort = TRUE,
  ...
)

distroll_confidence(
  x,
  distance,
  w = NULL,
  width = NULL,
  min_n = 2,
  align = c("right", "center", "left"),
  sort = TRUE,
  ...
)

distroll_dispersion(
  x,
  y,
  w = NULL,
  w.y = NULL,
  distance,
  width = NULL,
  min_n = 2,
  align = c("right", "center", "left"),
  sort = TRUE,
  ...
)

distroll_dispersion_sde(
  x,
  y,
  w = NULL,
  w.y = NULL,
  distance,
  width = NULL,
  min_n = 2,
  align = c("right", "center", "left"),
  sort = TRUE,
  ...
)

Arguments

Value

two-column vectors of (sorted) x and the rolled statistics along distance.

Note

distroll_circstats() and friends are complete, and for new code it is recommended switching to distance_binned_stats(), which is fasrter, easier to use, more featureful, and still under active development.

Examples

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")
data("san_andreas")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
san_andreas$distance <- distance_from_pb(
  x = san_andreas,
  PoR = PoR,
  pb = plate_boundary,
  tangential = TRUE
)
dat <- san_andreas |> cbind(PoR_shmax(san_andreas, PoR, "right"))

distroll_circstats(dat$azi.PoR,
  distance = dat$distance,
  w = 1 / dat$unc, FUN = circular_mean
) |> head()
distroll_confidence(dat$azi.PoR, distance = dat$distance, w = 1 / dat$unc) |> head()
distroll_dispersion(dat$azi.PoR,
  y = 135,
  distance = dat$distance, w = 1 / dat$unc
) |> head()
distroll_dispersion_sde(dat$azi.PoR,
  y = 135,
  distance = dat$distance, w = 1 / dat$unc, R = 100
) |> head()

# New functions
distance_binned_stats(
  dat$azi.PoR,
  distance = dat$distance, width.breaks = 1, unc = dat$unc, prd = 135
) |> head()

Rose Diagram

Description

Plots a rose diagram (rose of directions), the analogue of a histogram or density plot for angular data.

Usage

rose(
  x,
  weights = NULL,
  binwidth = NULL,
  bins = NULL,
  axial = TRUE,
  equal_area = TRUE,
  muci = TRUE,
  round_binwidth = 0,
  origin.text = "N",
  main = NULL,
  sub = NULL,
  at = seq(0, 360 - 45, 45),
  cborder = TRUE,
  labels = TRUE,
  col = "grey",
  dots = FALSE,
  dot_pch = 1,
  dot_cex = 1,
  dot_col = "slategrey",
  stack = FALSE,
  jitter_factor = 0,
  grid = FALSE,
  grid.lines = seq(0, 135, 45),
  grid.circles = seq(0.2, 1, 0.2),
  add = FALSE,
  ...
)

Arguments

Value

A window (class "owin") containing the plotted region or a list of the calculated frequencies.

Note

If bins and binwidth are NULL, an optimal bin width will be calculated using Scott (1979):

w_b = \frac{R}{n^{\frac{1}{3}}}

with n being the length of x, and the range R being either 180 or 360 degree for axial or directional data, respectively.

If "axial" == TRUE, the binwidth is adjusted to guarantee symmetrical fans.

See Also

Other circular-plot: plot_density(), plot_points(), rose_geom, rose_stats()

Examples

x <- rvm(100, mean = 90, k = 5)
rose(x, axial = FALSE, border = TRUE, grid = TRUE)

data("san_andreas")
rose(san_andreas$azi, main = "equal area")
rose(san_andreas$azi, equal_area = FALSE, main = "equal angle")

# weighted frequencies:
rose(san_andreas$azi, weights = 1 / san_andreas$unc, main = "weighted")

# add dots:
rose(san_andreas$azi, dots = TRUE, main = "dot plot", jitter = .2)

# stack dots:
rose(san_andreas$azi,
  dots = TRUE, stack = TRUE, dot_cex = 0.5, dot_pch = 21,
  main = "stacked dot plot"
)

Selecting optimal number of bins and width for rose diagrams

Description

Selecting optimal number of bins and width for rose diagrams

Usage

rose_bins(n, round = FALSE)

rose_binwidth(n, axial = TRUE, ...)

Arguments

Direction Lines and Fans in Circular Diagram

Description

Direction Lines and Fans in Circular Diagram

Usage

rose_line(x, radius = 1, axial = TRUE, add = TRUE, ...)

rose_fan(x, d, radius = 1, axial = TRUE, add = TRUE, ...)

Arguments

Value

No return value, called for side effects

See Also

Other circular-plot: plot_density(), plot_points(), rose(), rose_stats()

Examples

angles <- c(0, 10, 45)
radius <- c(.7, 1, .2)
lwd <- c(2, 1, .75)
col <- c(1, 2, 3)
rose_line(angles, radius = radius, axial = FALSE, add = FALSE, lwd = lwd, col = col)

Show Average Direction and Spread in Rose Diagram

Description

Adds the average direction (and its spread) to an existing rose diagram.

Usage

rose_stats(
  x,
  weights = NULL,
  axial = TRUE,
  avg = c("mean", "median", "sample_median"),
  spread = c("CI", "fisher", "sd", "IQR", "mdev"),
  avg.col = "#B63679FF",
  avg.lty = 2,
  avg.lwd = 1.5,
  spread.col = ggplot2::alpha("#B63679FF", 0.2),
  spread.border = FALSE,
  spread.lty = NULL,
  spread.lwd = NULL,
  add = TRUE,
  ...
)

Arguments

Value

plot or a two-element vector containing the calculated average and spread when assigned.

See Also

circular_mean(), circular_median(), circular_sample_median(), confidence_interval(), confidence_interval_fisher(), circular_sd(), circular_IQR(), circular_sample_median_deviation() for statistical parameters.

Other circular-plot: plot_density(), plot_points(), rose(), rose_geom

Examples

data("san_andreas")
rose(san_andreas$azi, weights = 1 / san_andreas$unc, muci = FALSE)
rose_stats(san_andreas$azi, weights = 1 / san_andreas$unc, avg = "sample_median", spread = "mdev")

Rotate Lines

Description

Rotates a set of straight lines around an angle

Usage

rotate_lines(theta, p, centre)

Arguments

Value

matrix

Sample circular dispersion

Description

Alternative versions of variance, dispersion a distance (Mardia and Jupp, 1999; pp. 19-20). These alternative dispersion has a minimum at the sample median.

Usage

sample_circular_variance(x, w = NULL, axial = TRUE)

sample_circular_distance(x, y, axial = TRUE, na.rm = TRUE)

sample_circular_dispersion(
  x,
  y = NULL,
  w = NULL,
  w.y = NULL,
  axial = TRUE,
  na.rm = TRUE
)

Arguments

References

N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

Mardia, K.V., and Jupp, P.E (1999). Directional Statistics, Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., Hoboken, NJ, USA. doi:10.1002/9780470316979

Examples

a <- c(0, 2, 359, 6, 354)
sample_circular_distance(a, 10) # distance to single value

b <- a + 90
sample_circular_distance(a, b) # distance to multiple values

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
sample_circular_variance(sa.por$azi.PoR)
sample_circular_dispersion(sa.por$azi.PoR, y = 135)
sample_circular_dispersion(sa.por$azi.PoR, y = 135, w = weighting(san_andreas$unc))

Sample Circular Median and Deviation

Description

Sample median direction for a vector of circular data

Usage

circular_sample_median(x, axial = TRUE, na.rm = TRUE)

circular_sample_median_deviation(x, axial = TRUE, na.rm = TRUE)

Arguments

Value

numeric

References

N.I. Fisher (1993) Statistical Analysis of Circular Data, Cambridge University Press.

Examples

set.seed(1)
x <- rvm(n = 100, mean = 0, kappa = 10)
circular_sample_median(x)
circular_sample_median_deviation(x)

data("san_andreas")
circular_sample_median(san_andreas$azi)
circular_sample_median_deviation(san_andreas$azi)

Second Central Momentum

Description

Measures the skewness (a measure of the asymmetry of the probability distribution) and the kurtosis (measure of the "tailedness" of the probability distribution). Standardized versions are the skewness and kurtosis normalized by the mean resultant length (Mardia 1972).

Usage

second_central_moment(x, w = NULL, axial = TRUE, na.rm = FALSE)

Arguments

Details

Negative values of skewness indicate skewed data in counterclockwise direction.

Large kurtosis values indicate tailed, values close to 0 indicate packed data.

Value

list containing

skewness

second central sine momentum, i.e. the skewness

std_skewness

standardized skewness

kurtosis

second central cosine momentum, i.e. the kurtosis

std_kurtosis

standardized kurtosis

Examples

data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
second_central_moment(sa.por$azi.PoR)
second_central_moment(sa.por$azi.PoR, w = weighting(san_andreas$unc))

Shortest distance between pairs of geometries

Description

The shortest Great Circle distance between pairs of geometries

Usage

shortest_distance_to_line(x, line, ellipsoidal = FALSE)

Arguments

Value

numeric. Shortest distance in meters

Examples

plate_boundary <- subset(plates, plates$pair == "na-pa")
shortest_distance_to_line(san_andreas, plate_boundary) |>
  head()

Quadrant-specific inverse of the tangent

Description

Returns the quadrant specific inverse of the tangent

Usage

atan2_spec(x, y)

atan2d_spec(x, y)

Arguments

Value

numeric.

References

Jammalamadaka, S. Rao, and Ambar Sengupta (2001). Topics in circular statistics. Vol. 5. world scientific.

Angle along great circle on spherical surface

Description

Smallest angle between two points on the surface of a sphere, measured along the surface of the sphere

Usage

orthodrome(lat1, lon1, lat2, lon2)

haversine(lat1, lon1, lat2, lon2)

vincenty(lat1, lon1, lat2, lon2)

Arguments

Details

"orthodrome"

based on the spherical law of cosines

"haversine"

uses haversine formula that is optimized for 64-bit floating-point numbers

"vincenty"

uses Vincenty formula for an ellipsoid with equal major and minor axes

Value

numeric. Angle in radians

References

Imboden, C. & Imboden, D. (1972). Formel fuer Orthodrome und Loxodrome bei der Berechnung von Richtung und Distanz zwischen Beringungs- und Wiederfundort. Die Vogelwarte 26, 336-346.

Sinnott, Roger W. (1984). Virtues of the Haversine. Sky and telescope 68(2), 158.

Vincenty, T. (1975). Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations. Survey Review, 23(176), 88<U+2013>93. doi:10.1179/sre.1975.23.176.88.

http://www.movable-type.co.uk/scripts/latlong.html

http://www.edwilliams.org/avform147.htm

Examples

berlin <- c(52.52, 13.41) |> deg2rad()
calgary <- c(51.04, -114.072) |> deg2rad()
orthodrome(berlin[1], berlin[2], calgary[1], calgary[2]) # 1.176406
haversine(berlin[1], berlin[2], calgary[1], calgary[2]) # 1.176406
vincenty(berlin[1], berlin[2], calgary[1], calgary[2]) # 1.176406

Spatial Interpolation of SHmax

Description

Stress field interpolation and wavelength analysis using a kernel (weighted) mean/median and standard deviation/IQR of stress data. Parameters can be adjusted to have inverse-distance-weighting (IDW) or nearest-neighbor interpolations (NN).

Usage

stress2grid(
  x,
  stat = c("mean", "median", "tensor"),
  grid = NULL,
  lon_range = NULL,
  lat_range = NULL,
  gridsize = 2,
  min_data = 3L,
  max_data = Inf,
  max_sd = Inf,
  threshold = deprecated(),
  min_dist_threshold = 200,
  arte_thres = deprecated(),
  method_weighting = FALSE,
  quality_weighting = TRUE,
  dist_weighting = c("inverse", "linear", "none"),
  idp = 1,
  qp = 1,
  mp = 1,
  dist_threshold = 0.1,
  R_range = seq(50, 1000, 50),
  axial = TRUE
)

stress2grid_stats(
  x,
  grid = NULL,
  lon_range = NULL,
  lat_range = NULL,
  gridsize = 2,
  min_data = 4L,
  max_data = Inf,
  threshold = deprecated(),
  min_dist_threshold = 200,
  arte_thres = deprecated(),
  method_weighting = FALSE,
  quality_weighting = TRUE,
  dist_weighting = c("inverse", "linear", "none"),
  idp = 1,
  qp = 1,
  mp = 1,
  dist_threshold = 0.1,
  R_range = seq(50, 1000, 50),
  mode = FALSE,
  kappa = 10,
  axial = TRUE,
  ...
)

Arguments

Details

stress2grid() is originally based on the MATLAB script "stress2grid" by Ziegler and Heidbach (2019): https://github.com/MorZieg/Stress2Grid. The tectonicr version has been significantly modified to provide better performance and more flexibility.

stress2grid_stats() is based on stress2grid() but calculates circular summary statistics (see circular_summary()).

Value

sf object containing

lon,lat

longitude and latitude in degrees

azi

Circular mean od median SHmax in degree

sd

Circular standard deviation or Quasi-IQR on the Circle of \sigma_\text{Hmax} in degrees

R

Search radius in km

mdr

Mean distance between grid point and datapoints per search radius

N

Number of data points in search radius

When stress2grid_stats(), azi and sd are replaced by the output of circular_summary().

Note

Although specialized for stress fields, this spatial interpolation algorithm can be applied to type of data. Please adjust the axial argument when applied to different, directional data sets.

References

Ziegler, M. and Heidbach, O. (2019). Matlab Script Stress2Grid v1.1. GFZ Data Services. doi:10.5880/wsm.2019.002

See Also

dist_greatcircle(), PoR_stress2grid(), compact_grid(), circular_mean(), circular_median(), circular_sd(), circular_summary()

Examples

data("san_andreas")

# Inverse Distance Weighting interpolation:
stress2grid(san_andreas, stat = "median") |> head()

stress2grid(san_andreas, stat = "tensor") |> head()

# Nearest Neighbor interpolation:
stress2grid(san_andreas, stat = "median", max_data = 5) |> head()

## Not run: 
stress2grid_stats(san_andreas, mode = TRUE) |> head()

## End(Not run)

Quick analysis of a stress data set

Description

Returns the converted azimuths, distances to the plate boundary, statistics of the model, and some plots.

Usage

stress_analysis(
  x,
  PoR,
  type = c("none", "in", "out", "right", "left"),
  pb,
  plot = TRUE,
  ...
)

Arguments

Value

list containing the following values:

results

data.frame showing the the coordinate and azimuth conversions (lat.PoR, lon.PoR, and azi.PoR), the predicted azimuths (prd), deviation angle from predicted (dev), circular distance (cdist), misfit to predicted stress direction (nchisq) and, if given, distance to tested plate boundary (distance)

stats

array with circular (weighted) mean, circular standard deviation, circular variance, circular median, skewness, kurtosis, the 95% confidence angle, circular dispersion, and the normalized Chi-squared test statistic

test

list containing the test results of the (weighted) Rayleigh test against the uniform distribution about the predicted orientation.

See Also

PoR_shmax(), distance_from_pb(), norm_chisq(), quick_plot(), circular_summary()

Examples

data("nuvel1")
na_pa <- subset(nuvel1, nuvel1$plate.rot == "na")

data("plates")
plate_boundary <- subset(plates, plates$pair == "na-pa")

data("san_andreas")
stress_analysis(san_andreas, na_pa, type = "right", plate_boundary, plot = TRUE)

Color palette for stress regime

Description

Color palette for stress regime

Usage

stress_colors()

Value

function

Examples

stress_colors()

Example crustal stress dataset

Description

Subsets of the World Stress Map (WSM) compilation of information on the crustal present-day stress field (Version 1.1. 2019).

Usage

data('san_andreas')

data('tibet')

data('iceland')

Format

A sf object / data.frame with 10 columns. Each row represents a different in-situ stress measurement:

id

Measurement identifier

lat

Latitude in degrees

lon

Longitude in degrees

azi

SHmax azimuth in degrees

unc

Measurement standard deviation (in degrees)

type

Type of measurement

depth

Depth in km

quality

WSM quality rank

regime

Stress regime

An object of class sf (inherits from data.frame) with 1126 rows and 10 columns.

An object of class sf (inherits from data.frame) with 1165 rows and 10 columns.

An object of class sf (inherits from data.frame) with 490 rows and 10 columns.

Details

⁠'san_andreas"⁠

contains 407 stress data adjacent to the San Andreas Fault to be tested against a tangentially displaced plate boundary.

"tibet"

contains 947 stress data from the Himalaya and Tibetan plateau to be tested against an inward-moving displaced plate boundary.

⁠'iceland⁠

contains 201 stress data from Iceland to be tested against a outward-moving displaced plate boundary.

Source

https://www.world-stress-map.org/

References

Heidbach, O., Barth, A., Müller, B., Reinecker, J., Stephansson, O., Tingay, M., & Zang, A. (2016). WSM quality ranking scheme, database description and analysis guidelines for stress indicator. WSM Technical Report; 16-01. GFZ German Research Centre for Geosciences. doi:10.2312/WSM.2016.001

See Also

download_WSM() for description of columns and stress regime acronyms

Examples

data("san_andreas")
head(san_andreas)

data("tibet")
head(tibet)

data("iceland")
head(iceland)

Theoretical Plate Tectonic Stress Paths

Description

Construct \sigma_{Hmax} lines that are following small circles, great circles, or loxodromes of an Euler pole for the relative plate motion.

Usage

eulerpole_paths(x, type = c("sc", "gc", "ld"), n = 10L, angle = 45, cw)

eulerpole_smallcircles(x, n = 10L)

eulerpole_greatcircles(x, n = 10L)

eulerpole_loxodromes(x, n = 10L, angle = 45, cw)

Arguments

Details

Maximum horizontal stress can be aligned to three types of curves related to relative plate motion:

Small circles

Lines that have a constant distance to the Euler pole. If x contains angle, output additionally gives absolute velocity on small circle (°/Myr -> km/Myr).

Great circles

Paths of the shortest distance between the Euler pole and its antipodal position.

Loxodromes

Lines of constant bearing, i.e. curves cutting small circles at a constant angle.

Value

sf object

Author(s)

Tobias Stephan

Examples

data("nuvel1")
por <- subset(nuvel1, nuvel1$plate.rot == "na") # North America relative to
# Pacific plate

eulerpole_smallcircles(por)
eulerpole_greatcircles(por)
eulerpole_loxodromes(x = por, angle = 45, n = 10, cw = FALSE)
eulerpole_loxodromes(x = por, angle = 30, cw = TRUE)
eulerpole_smallcircles(data.frame(lat = 30, lon = 10))

Glacial striae from Finland

Description

Directions of 51 glacial striations, example data from Davis (1986), pp. 316

Format

numeric vector of length 51

References

Davis, J.C. (1986) Statistics and Data Analysis in Geology. 2nd Edition, John Wiley & Sons, New York.

Examples

data("striae")
striae

SHmax direction resulting from multiple plate boundaries

Description

Calculates a \sigma_\text{Hmax} direction at given coordinates, sourced by multiple plate boundaries. This first-order approximation is the circular mean of the superimposed theoretical directions, weighted by the rotation rates of the underlying PoRs.

Usage

superimposed_shmax(df, PoRs, types, absolute = TRUE, PoR_weighting = NULL)

Arguments

Value

two column vector. azi is the resultant azimuth in degrees / geographical CRS), R is the resultant length.

See Also

model_shmax()

superimposed_shmax_PB() for considering distances to plate boundaries

Examples

data(san_andreas)
data(nuvel1)
pors <- subset(nuvel1, plate.rot %in% c("eu", "na"))
res <- superimposed_shmax(san_andreas, pors, types = c("in", "right"), PoR_weighting = c(2, 1))
head(res)

SHmax direction resulting from multiple plate boundaries considering distance to plate boundaries

Description

Calculates a \sigma_{Hmax} direction at given coordinates, sourced by multiple plate boundaries. This first-order approximation is the circular mean of the superimposed theoretical directions, weighted by the rotation rates of the underlying PoRs, the inverse distance to the plate boundaries, and the type of plate boundary.

Usage

superimposed_shmax_PB(
  x,
  pbs,
  model,
  rotation_weighting = TRUE,
  type_weights = c(divergent = 1, convergent = 3, transform_L = 2, transform_R = 2),
  idp = 1
)

Arguments

Value

two-column matrix. azi is the resultant azimuth (in degrees), R is the resultant length.

See Also

superimposed_shmax()

Examples

na_grid <- sf::st_make_grid(san_andreas, what = "centers", cellsize = 1)
na_plate <- subset(plates, plateA == "na" | plateB == "na")
cpm <- cpm_models[["NNR-MORVEL56"]]

# make divergent to ridge-push:
na_plate <- transform(na_plate, type = ifelse(na_plate$pair == "eu-na", "convergent", type))

res <- superimposed_shmax_PB(na_grid, na_plate, model = cpm, idp = 2)
head(res)

Colors for input variables

Description

assigns colors to continuous or categorical values for plotting

Usage

tectonicr.colors(
  x,
  n = 10,
  pal = NULL,
  categorical = FALSE,
  na.value = "grey",
  ...
)

Arguments

Value

named color vector

Examples

val1 <- c("N", "S", "T", "T", NA)
tectonicr.colors(val1, categorical = TRUE)
tectonicr.colors(val1, pal = stress_colors(), categorical = TRUE)

val2 <- runif(10)
tectonicr.colors(val2, n = 5)

Trigonometric Functions in Degrees

Description

Trigonometric functions expecting input in degrees.

Usage

sind(x)

cosd(x)

tand(x)

asind(x)

acosd(x)

atand(x)

atan2d(x1, x2)

cot(x)

cotd(x)

Arguments

Value

scalar or vector of numeric values.

Vector cross product

Description

Vector or cross product

Usage

vcross(x, y)

Arguments

Value

numeric vector of length 3

Examples

vcross(c(1, 2, 3), c(4, 5, 6))

von Mises Quantile-Quantile Plot

Description

Produces a Q-Q plot of the data against a specified von Mises distribution to graphically assess the goodness of fit of the model.

Usage

vm_qqplot(
  x,
  w = NULL,
  axial = TRUE,
  mean = NULL,
  kappa = NULL,
  xlab = "von Mises quantile function",
  ylab = "Empirical quantile function",
  main = "von Mises Q-Q Plot",
  col = "#B63679FF",
  add_line = TRUE,
  ...
)

Arguments

Value

plot

See Also

vonmises()

Other circ-qqplot: circular_qqplot()

Examples

set.seed(20250411)

# von Mises distribution
x_vm <- rvm(100, mean = 0, kappa = 2)
vm_qqplot(x_vm, pch = 20)

x_wcauchy <- rwcauchy(100, mean = 0, rho = 0.5)
vm_qqplot(x_wcauchy, pch = 20)

# circular uniform data
x_cunif <- rcunif(100)
vm_qqplot(x_cunif, pch = 20)

The von Mises Distribution

Description

Density, probability distribution function, quantiles, and random generation for the circular normal distribution with mean \mu and kappa \kappa.

Usage

rvm(n, mean, kappa)

dvm(theta, mean, kappa, axial = FALSE, log = FALSE)

pvm(theta, mean, kappa, from = NULL, tol = 1e-20)

qvm(p, mean, kappa, from = NULL, tol = .Machine$double.eps^(0.6), ...)

Arguments

Value

dvm gives the density, pvm gives the probability of the von Mises distribution function, rvm generates random deviates (in degrees), and qvm provides quantiles (in degrees).

See Also

cunif, wnorm, and wcauchy

Examples

set.seed(1)
x <- rvm(5, mean = 90, kappa = 2)

dvm(x, mean = 90, kappa = 2)
dvm(x, mean = 90, kappa = 2, axial = TRUE)

pvm(x, mean = 90, kappa = 2)
qvm(c(.25, .5, .75), mean = 90, kappa = 2)

Watson's U^2 Test for Goodness-of-Fit against known Distribution

Description

A non-parametric statistical test used for circular data to determine whether a sample fits a specified theoretical distribution

Usage

watson_test(
  x,
  alpha = NULL,
  dist = c("uniform", "vonmises"),
  axial = TRUE,
  quiet = FALSE
)

Arguments

Details

Hypotheses

Null Hypothesis (H_0): The circular sample comes from a specified theoretical distribution (such as a uniform distribution or a specific von Mises distribution).

Alternative Hypothesis (H_1): The circular sample does not follow the specified theoretical distribution.

Interpretation

To interpret the output of Watson's U^2 test, compare your calculated U^2 test statistic to the critical value from Watson's goodness-of-fit/homogeneity tables at your chosen significance level (\alpha, commonly set to 0.05), or check the resulting p-value:

• If U^2_{\text{calculated}} > U^2_{\text{critical}} (or p < \alpha): Reject the null hypothesis (H_0). Conclude that the data significantly deviates from the theoretical distribution.

• If U^2_{\text{calculated}} \le U^2_{\text{critical}} (or p \ge \alpha): Fail to reject the null hypothesis (H_0). There is not enough evidence to claim the data deviates from the expected model.

Value

list containing the test statistic statistic, the significance level p.value, the critical value critical.value, whether to reject the null hypothesis, the significance level alpha, the tested distribution dist, and the number of data n

Note

Watson's test statistic is a rotation-invariant Cramer - von Mises test. non-parametric, rank-based alternative to one-sample

References

Mardia and Jupp (1999). Directional Statistics. John Wiley and Sons.

See Also

vonmises and cunif

Other Tests: ar_test(), kuiper_test(), norm_chisq(), rayleigh-test, watson_two_sample, watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

# Example data from Mardia and Jupp (1999), pp. 93
watson_test(homing, axial = FALSE, alpha = .05)

# San Andreas Fault Data:
data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
watson_test(sa.por$azi.PoR, alpha = .05)
watson_test(sa.por$azi.PoR, alpha = .05, dist = "vonmises")
watson_test(sa.por$azi.PoR, alpha = .05, dist = "vonmises")

Watson's Two-Sample Test of Homogeneity

Description

Performs Watson's test for homogeneity on two samples of circular data. watson_two_test_perm() uses permutation to estimate p-values.

Usage

watson_two_test(x, y, alpha = NULL, axial = TRUE, quiet = FALSE)

watson_two_test_perm(x, y, axial = TRUE, n_perm = 1000L, alpha = NULL)

Arguments

Details

A two-sample Watson's U^2 permutation test determines whether two independent groups of circular data (angles or directions) come from the same underlying distribution.

Hypotheses

Null Hypothesis (H_0): The two samples come from the same circular distribution (the two groups of angles share a common distribution around the circle).

Alternative Hypothesis (H_1): The two samples come from different circular distributions (the groups are oriented or dispersed differently around the circle).

Interpretation

The Test Statistic (U^2) measures the distance between the cumulative distribution functions of the two circular samples. A larger U^2 value means the two sets of angles look more different from each other.

The P-Value represents the probability of getting a U² value as large as (or larger than) your observed value purely by chance, assuming the null hypothesis is true. It is calculated by shuffling the group labels across your data many times to build a reference permutation distribution.

Making a Decision

• Low p-value (p \le \alpha, usually 0.05): Reject the null hypothesis. Conclude that the two groups have significantly different circular distributions.

• High p-value (p > \alpha): Fail to reject the null hypothesis. There is not enough evidence to say the two groups are distributed differently around the circle.

Value

list. Watson's two-sample test of homogeneity is performed, and the results are printed. If alpha is specified and non-zero, the test statistic is printed along with the critical value and decision. If alpha is omitted, the test statistic is printed and a range for the p-value of the test is given.

Note

Critical values for the test statistic are obtained using the asymptotic distribution of the test statistic. It is recommended to use the obtained critical values and ranges for p-values only for combined sample sizes in excess of 17. Tables are available for smaller sample sizes and can be found in Mardia (1972) for instance.

See Also

Other Tests: ar_test(), kuiper_test(), norm_chisq(), rayleigh-test, watson_test(), watson_wheeler_test_perm(), weighted-rayleigh-test

Examples

set.seed(20250411)
x1 <- c(35, 45, 50, 55, 60, 70, 85, 95, 105, 120)
x2 <- c(75, 80, 90, 100, 110, 130, 135, 140, 150, 160, 165)
watson_two_test(x1, x2, axial = FALSE)
watson_two_test_perm(x1, x2, axial = FALSE)

data1 <- rvm(n=20, mean = 0, kappa=3)
data2 <- rvm(n=20, mean = 90, kappa=2)
watson_two_test(data1, data2, axial = FALSE)
watson_two_test_perm(data1, data2, axial = FALSE)


# San Andreas Fault Data:
data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
watson_two_test(sa.por$azi, 135, alpha = 0.05)
watson_two_test_perm(sa.por$azi, rvm(100, 135, 10), alpha = 0.05)

Watson-Wheeler Test of Homogeneity of Means

Description

A a non-parametric statistical test used to determine whether two or more independent samples of circular data (angles, directions, or periodic times) come from the same underlying population distribution. The difference between the samples can be in either the mean or the variance.

Usage

watson_wheeler_test_perm(x, y, axial = TRUE, n_perm = 1000L, alpha = NULL)

Arguments

Details

Hypotheses

Null Hypothesis (H_0) The samples come from identical populations (meaning both the mean direction and the dispersion/variance are homogeneous across groups).

Alternative Hypothesis (H_1): At least one sample comes from a different population distribution, which can be due to a difference in the mean direction, a difference in variance/concentration, or both.

Interpretation

Test Statistic (W): This value follows an approximate \chi^2 distribution. Higher values of W indicate larger discrepancies between the angular distributions of your groups.

• If the p-value is less than your significance level (commonly \alpha = 0.05), you reject the null hypothesis. This means you have strong evidence that the groups differ significantly in their central direction or spread around the circle.

• If the p-value is greater than 0.05, you fail to reject the null hypothesis, meaning there is no statistically significant evidence of difference among the groups.

Value

list

Note

Important Considerations & Limitations:

• Sensitivity to both mean and variance: Because it detects differences in either mean or variance, a significant result does not automatically mean the mean angles are different; it could be driven entirely by differences in concentration (variance).

• Sample size requirement: The chi-squared approximation requires each group to have a minimum sample size (typically at least 10 elements per group) to remain valid.

See Also

Other Tests: ar_test(), kuiper_test(), norm_chisq(), rayleigh-test, watson_test(), watson_two_sample, weighted-rayleigh-test

Examples

set.seed(20250411)
x1 <- c(35, 45, 50, 55, 60, 70, 85, 95, 105, 120)
x2 <- c(75, 80, 90, 100, 110, 130, 135, 140, 150, 160, 165)
watson_wheeler_test_perm(x1, x2, axial = FALSE)

data1 <- rvm(n=20, mean = 0, kappa=3)
data2 <- rvm(n=20, mean = 90, kappa=2)
watson_wheeler_test_perm(data1, data2, axial = FALSE)

# San Andreas Fault Data:
data(san_andreas)
data("nuvel1")
PoR <- subset(nuvel1, nuvel1$plate.rot == "na")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
watson_wheeler_test_perm(sa.por$azi.PoR, rvm(100, 135, 10))

The Wrapped Cauchy Distribution

Description

Density, probability distribution function, quantiles, and random generation for the circular wrapped Cauchy distribution with mean \mu and rho \rho

Usage

rwcauchy(n, mean, rho)

dwcauchy(theta, mean, rho, axial = FALSE, log = FALSE)

pwcauchy(
  theta,
  mean,
  rho,
  axial = FALSE,
  from = NULL,
  lower.tail = TRUE,
  log.p = FALSE
)

qwcauchy(
  p,
  mean,
  rho,
  axial = FALSE,
  from = NULL,
  lower.tail = TRUE,
  log.p = FALSE
)

Arguments

Value

dwcauchy gives the density, pwcauchy gives the probability of the wrapped Cauchy distribution function, rwcauchy generates random deviates (in degrees), and qrwcauchy provides quantiles (in degrees).

See Also

cunif, wnorm, and vonmises

Examples

set.seed(1)
x <- rwcauchy(5, mean = 90, rho = exp(-1))

dwcauchy(x, mean = 90, rho = exp(-1))
dwcauchy(x, mean = 90, rho = exp(-1), axial = TRUE)

pwcauchy(x, mean = 90, rho = exp(-1))
qwcauchy(c(.25, .5, .75), mean = 90, rho = exp(-1))

Weighted Goodness-of-fit Test for Circular Data

Description

Weighted version of the Rayleigh test (or V0-test) for uniformity against a distribution with a priori expected von Mises concentration. weighted_rayleight_test_perm() uses permutation to estimate p-values.

Usage

weighted_rayleigh(x, mu = NULL, w = NULL, axial = TRUE, quiet = FALSE)

weighted_rayleigh_perm(x, mu = NULL, w = NULL, axial = TRUE, n_perm = 1000L)

Arguments

Details

The Null hypothesis is uniformity (randomness). The alternative is a distribution with a (specified) mean direction (mu). If statistic >= p.value, the null hypothesis of randomness is rejected and angles derive from a distribution with a (or the specified) mean direction.

Value

a list with the components:

R or C

mean resultant length or the dispersion (if mu is specified). Small values of R (large values of C) will reject uniformity. Negative values of C indicate that vectors point in opposite directions (also lead to rejection).

statistic

Test statistic

p.value

significance level of the test statistic

See Also

Other Tests: ar_test(), kuiper_test(), norm_chisq(), rayleigh-test, watson_test(), watson_two_sample, watson_wheeler_test_perm()

Examples

# Load data
data("cpm_models")
data(san_andreas)
PoR <- equivalent_rotation(cpm_models[["NNR-MORVEL56"]], "na", "pa")
sa.por <- PoR_shmax(san_andreas, PoR, "right")
data("iceland")
PoR.ice <- equivalent_rotation(cpm_models[["NNR-MORVEL56"]], "eu", "na")
ice.por <- PoR_shmax(iceland, PoR.ice, "out")
data("tibet")
PoR.tib <- equivalent_rotation(cpm_models[["NNR-MORVEL56"]], "eu", "in")
tibet.por <- PoR_shmax(tibet, PoR.tib, "in")

# GOF test:
weighted_rayleigh(tibet.por$azi.PoR, mu = 90, w = 1 / tibet$unc)
weighted_rayleigh(ice.por$azi.PoR, mu = 0, w = 1 / iceland$unc)
weighted_rayleigh(sa.por$azi.PoR, mu = 135, w = 1 / san_andreas$unc)

weighted_rayleigh_perm(tibet.por$azi.PoR, mu = 90, w = 1 / tibet$unc)
weighted_rayleigh_perm(ice.por$azi.PoR, mu = 0, w = 1 / iceland$unc)
weighted_rayleigh_perm(sa.por$azi.PoR, mu = 135, w = 1 / san_andreas$unc)

Weighting Factors

Description

Helper function to transform uncertainty angles into weighting factors

Usage

weighting(
  x,
  method = c("linear-inverse", "inverse", "cosine", "none"),
  max.err = 90,
  norm = TRUE
)

Arguments

Details

Linear inverse: w = 1 - x/\sigma, where \sigma is the maximum error expected for x (e.g. 90^\circ).

Inverse: w = 1/x

Cosine: w = \cos{x}

Value

numeric

Examples

x <- seq(0, 90, 1)

plot(x, weighting(x, "inverse"),
  col = 1, type = "l",
  xlab = "Uncertainty angle in degrees", ylab = "weight"
)
lines(x, weighting(x, "cosine"), col = 2)
lines(x, weighting(x, "linear-inverse"), col = 3)
legend("topright",
  col = 1:3, lty = 1,
  legend = c("inverse", "cosine", "linear-inverse")
)

Indices of n smallest values in array

Description

Indices of n smallest values in array

Usage

which.nsmallest(x, n)

The Wrapped Normal Distribution

Description

Density, probability distribution function, quantiles, and random generation for the circular normal distribution with mean \mu and standard deviation \sigma.

Usage

rwnorm(n, mean = 0, sd = 1)

dwnorm(theta, mean = 0, sd = 1, axial = FALSE, log = FALSE)

pwnorm(theta, mean = 0, sd = 1, axial = FALSE, from = NULL, ...)

qwnorm(
  p,
  mean = 0,
  sd = 1,
  axial = FALSE,
  from = NULL,
  tol = .Machine$double.eps^(0.6),
  ...
)

Arguments

Value

dwnorm gives the density, pwnorm gives the probability of the wrapped normal distribution function, rwnorm generates random deviates (in degrees), and qwnorm provides quantiles (in degrees).

See Also

cunif, wnorm, wcauchy, and vonmises

Examples

set.seed(1)
x <- rwnorm(5, mean = 90, sd = 5)

dwnorm(x, mean = 90, sd = 5, axial = FALSE)
dwnorm(x, mean = 90, sd = 5, axial = TRUE)

pwnorm(x, mean = 90, sd = 5)
qwnorm(c(.25, .5, .75), mean = 90, sd = 5)