Structural Geology in R • structr

structr is a free and open-source package for R that provides tools for structural geology. The toolset includes

Analysis and visualization of orientation data of structural geology (including, stereographic projections, contouring, fabric plots, and statistics),
Statistical analysis: spherical mean and variance, confidence regions, hypothesis tests, cluster analysis of orientation data (sph_cluster(), and geodesic regression to find the best-fitting great circle or small circle through orientation data (regression_greatcircle() and regression_smallcircle()),
Reconstruction of fabric orientations in oriented drillcores by transforming the α, β, and γ angles (drillcore_transformation(),
Stress analysis: reconstruction of stress orientation and magnitudes from fault-slip data (stress inversion based on Michael, 1984: slip_inversion()), extracting the maximum horizontal stress of a 3D stress tensor (SH()), and visualization of magnitudes of stress in the Mohr circle (Mohr_plot()),
Calculation fault displacement components,
Strain analysis (R_f/ϕ), contouring on the unit hyperboloid, Fry plots and Hsu plots
Vorticity analysis using the Rigid Grain Net method (RGN_plot()), and
Direct import of your field data from StraboSpot projects (read_strabo_JSON()).

The {structr} package is all about structures in 3D. For analyzing orientations in 2D (statistics, rose diagrams, etc.), check out the tectonicr package!

Installation

You can install the development version of structr from GitHub with:

# install.packages("devtools")
devtools::install_github("tobiste/structr")

Documentation

The detailed documentation can be found at https://tobiste.github.io/structr/

Examples

These are some basic examples which shows you what you can do with {structr}. First we load the package

library(structr)

Stereographic and equal-area projection

Plot orientation data in equal-area, lower hemisphere projection:

data("example_planes")
data("example_lines")

par(xpd = NA)
stereoplot(title = "Lambert equal-area projection", sub = "Lower hemisphere")
points(example_lines, col = "#B63679", pch = 19, cex = .5)
points(example_planes, col = "#000004", pch = 1, cex = .5)

legend("topright", legend = c("Lines", "Planes"), col = c("#B63679", "#000004"), pch = c(19, 1), cex = 1)

Density

Density shown by contour lines and filled contours:

par(mfrow = c(1, 2))
contour(example_planes)
points(example_planes, col = "grey", cex = .5)
title(main = "Planes")

image(example_lines)
points(example_lines, col = "grey", cex = .5)
title(main = "Lines")

Spherical statistics

Calculation of arithmetic mean, geodesic mean, confidence cones and eigenvectors… and plotting them in the equal-area projection:

planes_mean <- sph_mean(example_planes)
planes_geomean <- geodesic_mean(example_planes)
planes_eig <- ot_eigen(example_planes)$vectors

par(mfrow = c(1, 2), xpd = NA)
stereoplot(title = "Planes", guides = FALSE)
points(example_planes, col = "lightgrey", pch = 1, cex = .5)
lines(planes_eig, col = c("#FB8861FF", "#FEC287FF", "#FCFDBFFF"), lty = 1:3)
points(planes_mean, col = "#B63679", pch = 19, cex = 1)
points(planes_geomean, col = "#E65164FF", pch = 19, cex = 1)
points(planes_eig, col = c("#FB8861FF", "#FEC287FF", "#FCFDBFFF"), pch = 19, cex = 1)
legend(
  0, -1.1,
  xjust = .5,
  legend = c("Arithmetic mean", "Geodesic mean", "Eigen 1", "Eigen 2", "Eigen 3"),
  col = c("#B63679", "#E65164FF", "#FB8861FF", "#FEC287FF", "#FCFDBFFF"),
  pch = 19, lty = c(NA, NA, 1, 2, 3),
  cex = .75
)

lines_mean <- sph_mean(example_lines)
lines_delta <- delta(example_lines)
lines_confangle <- confidence_ellipse(example_lines)

stereoplot(title = "Lines", guides = FALSE)
points(example_lines, col = "lightgrey", pch = 1, cex = .5)
points(lines_mean, col = "#B63679", pch = 19, cex = 1)
stereo_confidence(lines_confangle, col = "#E65164FF")
lines(lines_mean, ang = lines_delta, col = "#FB8861FF")
legend(
  0, -1.1,
  xjust = .5,
  legend = c("Arithmetic mean", "95% confidence cone", "63% data cone"),
  col = c("#B63679", "#E65164FF", "#FB8861FF"),
  pch = c(19, NA, NA), lty = c(NA, 1, 1), cex = .75
)

Orientation tensor and fabric plots

The shape parameters of the orientation tensor of the above examples planes and lines can be visualized in two ways:

par(mfrow = c(1, 2), xpd = NA)
vollmer_plot(example_planes, col = "#000004", pch = 16)
vollmer_plot(example_lines, col = "#B63679FF", pch = 16, add = TRUE)
title("Fabric plot of Vollmer (1990)")

hsu_fabric_plot(example_planes, col = "#000004", pch = 16)
hsu_fabric_plot(example_lines, col = "#B63679FF", pch = 16, add = TRUE)

legend(
  2.5, -.25,
  xjust = .5, horiz = TRUE, xpd = NA,
  legend = c("Planes", "Lines"), col = c("#000004", "#B63679FF"), pch = 16
)

Best-fit great and small-circles (geodesic regression)

Finds the best-fit great or small-circle for a given set of vectors by applying geodesic regression:

set.seed(20250411)
data("gray_example")
cleavage <- gray_example[1:8, ]
bedding <- gray_example[9:16, ]

cleavage_gc <- regression_greatcircle(cleavage)
bedding_gc <- regression_greatcircle(bedding)

cleavage_sc <- regression_smallcircle(cleavage)
bedding_sc <- regression_smallcircle(bedding)

par(mfrow = c(1, 2), xpd = NA)
stereoplot(title = "Best greatcircle", guides = FALSE)
lines(cleavage_gc$vec, col = "#000004FF")
lines(bedding_gc$vec, col = "#B63679")
points(cleavage, col = "#1D1147")
points(bedding, col = "#E65164", pch = 4)

legend(
  0, -1.1,
  xjust = .5,
  col = c("#000004FF", "#B63679"),
  lty = c(1, 1), legend = c("Cleavage greatcircle", "Bedding greatcircle"), bg = "white"
)

stereoplot(title = "Best smallcircle", guides = FALSE)
lines(cleavage_sc$vec, cleavage_sc$cone, col = "#000004FF")
lines(bedding_sc$vec, bedding_sc$cone, col = "#B63679")
points(cleavage, col = "#1D1147")
points(bedding, col = "#E65164", pch = 4)

legend(0, -1.1,
  xjust = .5,
  col = c("#000004FF", "#B63679"), lty = c(1, 1), legend = c("Cleavage smallcircle", "Bedding smallcircle"), bg = "white"
)

Fault plots

Graphical representation of fault-slip data using Angelier plot (slip vector on fault plane great circle) and Hoeppener plot (fault slip vector projected on pole to fault plane):

data("angelier1990")
faults <- angelier1990$TYM

par(mfrow = c(1, 2))
stereoplot(title = "Angelier plot")
angelier(faults)

stereoplot(title = "Hoeppener plot")
hoeppener(faults, points = FALSE)

Fault-slip inversion

Compute deviatoric stress tensor and calculate 95% confidence intervals using bootstrap samples:

set.seed(20250411)
faults_stress <- slip_inversion(faults, boot = 10)

Visualize the slip inversion results (orientation of principal stresses):

cols <- c("#000004FF", "#B63679FF", "#FEC287FF")
R_val <- round(faults_stress$R, 2)
R_CI <- round(faults_stress$R_conf, 2)

stereoplot(
  title = "Principal stress axes",
  sub = paste0("Relative stress magnitudes R = ", R_val, " | ", "95% CI: [", R_CI[1], ", ", R_CI[2], "]"),
  guides = FALSE
)
angelier(faults, col = "grey80")
stereo_confidence(faults_stress$principal_axes_conf$sigma1, col = cols[1])
stereo_confidence(faults_stress$principal_axes_conf$sigma2, col = cols[2])
stereo_confidence(faults_stress$principal_axes_conf$sigma3, col = cols[3])
text(faults_stress$principal_axes,
  label = rownames(faults_stress$principal_axes),
  col = cols, adj = -.25
)

Visualize the accuracy of the slip inversion by showing the deviation angle (β) between the theoretical slip and the actual slip vector:

beta <- faults_stress$fault_data$beta
beta_mean <- round(faults_stress$beta)
beta_CI <- round(faults_stress$beta_CI)

stereoplot(
  title = "Stress inversion accuracy",
  sub = bquote("Average deviation" ~ bar(beta) == .(beta_mean) * degree ~ "\U00B1" ~ .(beta_CI) * degree),
  guides = FALSE
)
angelier(faults, col = assign_col(beta))
legend_c(
  seq(min(beta), max(beta), 10),
  title = bquote("Deviation angle" ~ beta ~ "(" * degree * ")")
)

Azimuth of the maximum horizontal stress (in degrees) for the slip inversion result:

# Simply call
# faults_stress$SHmax
# faults_stress$SHmax_CI # confidence interval

SH(
  S1 = faults_stress$principal_axes[1, ],
  S2 = faults_stress$principal_axes[2, ],
  S3 = faults_stress$principal_axes[3, ],
  R = faults_stress$R
)
#> [1] 60.80844

Mohr circles

The Mohr circle for the slip inversion result:

Mohr_plot(
  sigma1 = faults_stress$principal_vals[1],
  sigma2 = faults_stress$principal_vals[2],
  sigma3 = faults_stress$principal_vals[3],
  unit = NULL, include.zero = FALSE
)
points(faults_stress$fault_data$sigma_n, abs(faults_stress$fault_data$sigma_s),
  col = assign_col(beta), pch = 16
)

Strain analysis

2D Strain

Aspect ratio of finite strain ellipses vs orientation of long-axis (Rf/ϕ)

data(ramsay)

par(mfrow = c(1, 2))
Rphi_plot(r = ramsay[, 1], phi = ramsay[, 2])
Rphi_polar_plot(ramsay[, 1], ramsay[, 2], proj = "eqd")

3D Strain

Finite strain ellipsoids plotted in Flinn diagram and Hsu diagram:

data("holst")
R_XY <- holst[, "R_XY"]
R_YZ <- holst[, "R_YZ"]

par(mfrow = c(1, 2))
flinn_plot(R_XY, R_YZ, log = TRUE, col = "#B63679", pch = 16)
hsu_plot(R_XY, R_YZ, col = "#B63679", pch = 16)

Vorticity analysis

Aspect ratio of finite strain ellipses of porphyroclasts vs orientation of long-axis with respect to foliation plotted in the Rigid Grain Net

data(shebandowan)
set.seed(20250411)

# Color code porphyroclasts by size of clast (area in log-scale):
RGN_plot(shebandowan$r, shebandowan$phi, col = assign_col(log(shebandowan$area)), pch = 16)

Author

Tobias Stephan (tstephan@lakeheadu.ca)

Feedback, issues, and contributions

I welcome feedback, suggestions, issues, and contributions! If you have found a bug, please file it here with minimal code to reproduce the issue.

License

MIT License