liminfoStructural Geology Analysis with Stereonets
A comprehensive guide to projecting and analyzing planar structures (fault planes, bedding planes) and linear structures (fold axes, lineations) on stereonets. Covers Schmidt equal-area projection, great circle/pole plotting, mplstereonet Python code, and Kamb contour diagram analysis.
Problem
Required Tools
A matplotlib-based stereonet library. Generates great circles, poles, and contour diagrams in Python. Install with pip install mplstereonet.
A free stereonet program developed by Prof. Rick Allmendinger. GUI-based application supporting data input, great circle/pole plotting, and statistical analysis. Available for Windows/Mac/Linux.
Field instruments for measuring strike and dip of planar structures. Brunton compass is the standard; smartphone apps (e.g., FieldMove Clino) are also available.
Python visualization library. The foundation of mplstereonet, used for stereonet diagram customization and publication-quality figure output.
Solution Steps
Understanding Stereonet Fundamentals (Schmidt vs Wulff)
A stereonet is a method for projecting 3D orientation data onto a 2D circular diagram. It is one of the most important tools in structural geology. [Two Projection Methods] - Equal-Area (Schmidt Net): Preserves area, ideal for density analysis. Most commonly used in structural geology - Equal-Angle (Wulff Net): Preserves angles, primarily used in crystallography [Lower Hemisphere Projection] Structural geology always uses lower hemisphere projection. It represents looking down on the lower hemisphere below the horizontal plane. A plane is shown as a great circle, and the normal to a plane is shown as a pole (point). [Key Terminology] - Great Circle: The arc representing a plane projected onto the stereonet. Visually shows strike/dip orientation - Pole: A point representing the direction perpendicular to a plane. More useful than great circles for large datasets - Small Circle: A circle representing a cone projected onto the stereonet. Used for confidence intervals and rotation analysis
# Equal-Area vs Equal-Angle comparison
# Structural geology always uses Equal-Area (Schmidt) projection
import mplstereonet
import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5),
subplot_kw=dict(projection='stereonet'))
# Compare the same data with both projection methods
strikes = [120, 130, 125, 140, 135, 128, 132]
dips = [45, 50, 42, 55, 48, 44, 52]
for s, d in zip(strikes, dips):
ax1.plane(s, d, 'b-', linewidth=0.5)
ax1.pole(s, d, 'ro', markersize=4)
ax2.plane(s, d, 'b-', linewidth=0.5)
ax2.pole(s, d, 'ro', markersize=4)
ax1.set_title('Equal-Area (Schmidt)', fontsize=11)
ax2.set_title('Equal-Angle (Wulff)', fontsize=11)
ax1.grid(True)
ax2.grid(True)
plt.tight_layout()
plt.savefig('projection_comparison.png', dpi=150)
plt.show()Plotting Strike/Dip Data -- Great Circles and Poles
This step involves representing field-measured planar structure data as great circles and poles on a stereonet. [Strike/Dip Notation Systems] - Right-Hand Rule (RHR): With your right hand along the strike direction, your thumb points toward the dip direction. Example: N30E/45SE -> 030/45 - Dip-Direction convention: Written as dip-direction/dip-angle. Example: 120/45 (dip direction 120 degrees, dip angle 45 degrees) - Caution: Mixing these two conventions causes a 90-degree error! [Great Circles vs Poles] - Great Circles: Visually show the attitude of individual planes. Useful when data count is small (fewer than 10) - Poles: Points representing the perpendicular direction to a plane. Better for viewing large datasets (tens to hundreds) at a glance - Typical workflow: pole plot -> contour diagram
import mplstereonet
import matplotlib.pyplot as plt
import numpy as np
# Field measurement data (Strike/Dip, Right-Hand Rule convention)
# Example: bedding plane measurements
strikes = [30, 35, 28, 40, 32, 45, 38, 25, 42, 33,
130, 128, 135, 140, 125, 132, 138, 127, 136, 131]
dips = [45, 50, 42, 55, 48, 52, 47, 40, 53, 46,
60, 58, 65, 62, 55, 63, 59, 57, 64, 61]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5),
subplot_kw=dict(projection='stereonet'))
# --- Left: Great Circle plot ---
for s, d in zip(strikes, dips):
ax1.plane(s, d, 'b-', linewidth=0.8, alpha=0.6)
ax1.set_title('Great Circles', fontsize=12, fontweight='bold')
ax1.grid(True)
# --- Right: Pole plot ---
for s, d in zip(strikes, dips):
ax1.pole(s, d, 'ro', markersize=5, alpha=0.7)
ax2.pole(s, d, 'ro', markersize=6)
ax2.set_title('Poles to Planes', fontsize=12, fontweight='bold')
ax2.grid(True)
plt.suptitle('Bedding Plane Data (n=20)', fontsize=13)
plt.tight_layout()
plt.savefig('strike_dip_plot.png', dpi=150)
plt.show()
# To convert Dip-Direction to Strike:
# strike = dip_direction - 90 (RHR convention)
dip_direction = 120 # dip direction
dip_angle = 45 # dip angle
strike_rhr = dip_direction - 90 # -> 030
print(f"Dip-Direction {dip_direction}/{dip_angle} -> Strike/Dip {strike_rhr:03d}/{dip_angle}")Plotting Linear Structures (Trend/Plunge) -- Fold Axes, Lineations
Linear structures (lineations) are represented as points on a stereonet. These include fold axes, mineral lineations, intersection lineations, and similar features. [Trend/Plunge Notation] - Trend: The direction of the lineation projected onto the horizontal plane (0-360 degrees) - Plunge: The angle the lineation makes with the horizontal plane (0-90 degrees) - Example: 30/245 -> plunge 30 degrees, trend 245 degrees [Fold Axis Determination Methods] - Pi Diagram: The pole of the great circle (Pi girdle) formed by fold limb poles equals the fold axis - Beta Diagram: The intersection of adjacent fold limb planes equals the fold axis - mplstereonet provides the fit_girdle() function for automatic calculation
import mplstereonet
import matplotlib.pyplot as plt
import numpy as np
# Fold structure measurement data from both limbs
# Limb 1 (NW-dipping)
strikes_1 = [310, 315, 308, 320, 312, 318, 305, 322, 311, 316]
dips_1 = [35, 38, 32, 40, 36, 42, 30, 44, 34, 39]
# Limb 2 (SE-dipping)
strikes_2 = [130, 128, 135, 125, 132, 138, 127, 140, 133, 129]
dips_2 = [50, 48, 55, 45, 52, 58, 47, 60, 53, 49]
fig, ax = plt.subplots(figsize=(7, 7),
subplot_kw=dict(projection='stereonet'))
# Plot poles from both limbs
ax.pole(strikes_1, dips_1, 'bo', markersize=6, label='Limb 1 (NW)')
ax.pole(strikes_2, dips_2, 'rs', markersize=6, label='Limb 2 (SE)')
# Fit a Pi girdle to all poles (determines fold axis)
all_strikes = strikes_1 + strikes_2
all_dips = dips_1 + dips_2
# fit_girdle: Best-fit great circle through poles -> its pole = fold axis
fit_strike, fit_dip = mplstereonet.fit_girdle(all_strikes, all_dips)
ax.plane(fit_strike, fit_dip, 'g-', linewidth=2, label='Pi girdle')
# Fold axis = pole of the Pi girdle
fold_lon, fold_lat = mplstereonet.pole(fit_strike, fit_dip)
ax.line(fold_lat, fold_lon, 'g^', markersize=12, label='Fold axis')
# Calculate Trend/Plunge of fold axis
plunge, bearing = mplstereonet.geographic2plunge_bearing(fold_lon, fold_lat)
ax.set_title(f'Pi Diagram - Fold Axis: {plunge[0]:.0f}/{bearing[0]:.0f}',
fontsize=12, fontweight='bold')
ax.grid(True)
ax.legend(loc='upper left', fontsize=9)
plt.tight_layout()
plt.savefig('fold_axis_analysis.png', dpi=150)
plt.show()
print(f"Fold Axis -> Plunge: {plunge[0]:.1f}, Bearing: {bearing[0]:.1f}")Building Publication-Quality Stereonets with mplstereonet
A complete workflow for generating publication-quality stereonets using the mplstereonet library. [Installation] pip install mplstereonet matplotlib numpy [Key Functions] - ax.plane(strike, dip): Plot great circles - ax.pole(strike, dip): Plot poles - ax.line(plunge, bearing): Plot linear structures - ax.density_contourf(strikes, dips): Filled contour density diagram - ax.density_contour(strikes, dips): Line contour density diagram - mplstereonet.fit_girdle(): Best-fit great circle - mplstereonet.fit_pole(): Best-fit pole (mean orientation) [Working with CSV Data] Field survey data organized in CSV format can be bulk-processed using pandas.
import mplstereonet
import matplotlib.pyplot as plt
import numpy as np
# ===== Reading data from CSV files =====
# bedding_data.csv format:
# station,strike,dip,type
# ST001,035,45,bedding
# ST002,128,60,fault
import csv
def read_structural_data(filename):
"""Read structural data from CSV file and return as dictionary"""
data = {}
with open(filename, 'r') as f:
reader = csv.DictReader(f)
for row in reader:
dtype = row['type']
if dtype not in data:
data[dtype] = {'strikes': [], 'dips': []}
data[dtype]['strikes'].append(float(row['strike']))
data[dtype]['dips'].append(float(row['dip']))
return data
# ===== Comprehensive stereonet generation =====
# Example data (direct input instead of CSV)
np.random.seed(42)
bedding_s = np.random.normal(40, 12, 50) % 360
bedding_d = np.clip(np.random.normal(45, 8, 50), 5, 85)
fault_s = np.random.normal(160, 15, 20) % 360
fault_d = np.clip(np.random.normal(70, 10, 20), 10, 89)
fig, axes = plt.subplots(1, 3, figsize=(16, 5),
subplot_kw=dict(projection='stereonet'))
# (1) Great circles + poles
ax1 = axes[0]
ax1.plane(bedding_s, bedding_d, 'b-', linewidth=0.5, alpha=0.4)
ax1.pole(bedding_s, bedding_d, 'bo', markersize=3, alpha=0.6)
ax1.plane(fault_s, fault_d, 'r-', linewidth=0.5, alpha=0.4)
ax1.pole(fault_s, fault_d, 'rs', markersize=3, alpha=0.6)
ax1.set_title('Great Circles + Poles', fontsize=11)
ax1.grid(True)
# (2) Contour diagram - bedding planes
ax2 = axes[1]
cax = ax2.density_contourf(bedding_s, bedding_d,
measurement='poles', cmap='Reds', levels=8)
ax2.pole(bedding_s, bedding_d, 'ko', markersize=2, alpha=0.5)
fig.colorbar(cax, ax=ax2, shrink=0.6, label='Density (%)')
ax2.set_title(f'Bedding Poles (n={len(bedding_s)})', fontsize=11)
ax2.grid(True)
# (3) Contour diagram - fault planes
ax3 = axes[2]
cax2 = ax3.density_contourf(fault_s, fault_d,
measurement='poles', cmap='Blues', levels=8)
ax3.pole(fault_s, fault_d, 'ko', markersize=2, alpha=0.5)
fig.colorbar(cax2, ax=ax3, shrink=0.6, label='Density (%)')
ax3.set_title(f'Fault Poles (n={len(fault_s)})', fontsize=11)
ax3.grid(True)
plt.suptitle('Structural Data Analysis', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig('comprehensive_stereonet.png', dpi=200, bbox_inches='tight')
plt.show()Contour Diagram Analysis -- Kamb Contours and Preferred Orientation
Contour diagrams visualize the density distribution of poles as contour lines, enabling quantitative identification of preferred orientations. This is the core analytical method. [Kamb Contour Method] - Kamb (1959): Statistically significant density calculation. The counting area automatically adjusts based on sample size - Specified with the method='kamb' parameter in mplstereonet - Contour values are displayed in standard deviation (sigma) units [Interpreting Results] - Unimodal distribution: Single maximum -> one preferred planar orientation - Bimodal distribution: Two maxima -> two sets of planar structures (e.g., both limbs of a fold) - Girdle distribution: Distributed along a great circle -> fold structure present. Pole of the girdle = fold axis - Random distribution: No preferred orientation [Fold Axis Determination Procedure] 1. Plot all bedding plane poles 2. Generate contour diagram (Kamb method) 3. If maxima are distributed along a great circle (girdle), a fold exists 4. Fit the Pi girdle using fit_girdle() 5. The pole of the Pi girdle = fold axis (Trend/Plunge)
import mplstereonet
import matplotlib.pyplot as plt
import numpy as np
# Synthetic data representing a fold structure
np.random.seed(123)
# Limb 1: NW-dipping bedding
s1 = np.random.normal(315, 10, 40) % 360
d1 = np.clip(np.random.normal(35, 6, 40), 5, 85)
# Limb 2: SE-dipping bedding
s2 = np.random.normal(135, 10, 35) % 360
d2 = np.clip(np.random.normal(50, 7, 35), 5, 85)
all_strikes = np.concatenate([s1, s2])
all_dips = np.concatenate([d1, d2])
fig, axes = plt.subplots(1, 3, figsize=(16, 5),
subplot_kw=dict(projection='stereonet'))
# (1) Poles only
ax1 = axes[0]
ax1.pole(s1, d1, 'bo', markersize=4, alpha=0.6, label='Limb 1')
ax1.pole(s2, d2, 'rs', markersize=4, alpha=0.6, label='Limb 2')
ax1.set_title(f'Poles (n={len(all_strikes)})', fontsize=11)
ax1.grid(True)
ax1.legend(fontsize=8)
# (2) Kamb contour diagram
ax2 = axes[1]
cax = ax2.density_contourf(all_strikes, all_dips,
measurement='poles', method='kamb',
cmap='YlOrRd', levels=10)
ax2.density_contour(all_strikes, all_dips,
measurement='poles', method='kamb',
colors='black', linewidths=0.5, levels=10)
ax2.pole(all_strikes, all_dips, 'k.', markersize=2, alpha=0.3)
fig.colorbar(cax, ax=ax2, shrink=0.6, label='Kamb contour (sigma)')
ax2.set_title('Kamb Contour Diagram', fontsize=11)
ax2.grid(True)
# (3) Pi girdle + fold axis determination
ax3 = axes[2]
ax3.density_contourf(all_strikes, all_dips,
measurement='poles', method='kamb',
cmap='YlOrRd', levels=10, alpha=0.5)
# Fit Pi girdle
fit_s, fit_d = mplstereonet.fit_girdle(all_strikes, all_dips)
ax3.plane(fit_s, fit_d, 'g-', linewidth=2.5, label='Pi girdle')
# Fold axis (pole of the Pi girdle)
fold_lon, fold_lat = mplstereonet.pole(fit_s, fit_d)
ax3.line(fold_lat, fold_lon, 'g^', markersize=14,
markeredgecolor='black', label='Fold axis')
plunge, bearing = mplstereonet.geographic2plunge_bearing(fold_lon, fold_lat)
ax3.set_title(f'Fold Axis: {plunge[0]:.0f}/{bearing[0]:.0f}', fontsize=11)
ax3.grid(True)
ax3.legend(fontsize=8)
plt.suptitle('Fold Analysis: Pi Diagram with Kamb Contours',
fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('kamb_fold_analysis.png', dpi=200, bbox_inches='tight')
plt.show()
# ===== Statistical summary output =====
print("=" * 50)
print("Structural Analysis Summary")
print("=" * 50)
print(f"Total measurements: {len(all_strikes)}")
print(f" Limb 1: {len(s1)} measurements")
print(f" Limb 2: {len(s2)} measurements")
print(f"Pi girdle: {fit_s:.1f}/{fit_d:.1f}")
print(f"Fold axis: Plunge {plunge[0]:.1f}, Bearing {bearing[0]:.1f}")
# Calculate mean orientation
mean_s, mean_d = mplstereonet.fit_pole(all_strikes, all_dips)
print(f"Mean pole: {mean_s:.1f}/{mean_d:.1f}")Common Mistakes
Confusing Right-Hand Rule and Dip-Direction conventions, causing a 90-degree error
Strike/Dip (RHR) and Dip-Direction/Dip notations differ by 90 degrees. In RHR, Strike 030/Dip 45 is equivalent to Dip-Direction 120/45. Always verify which notation was used before entering data, and use the formula strike = dip_direction - 90 for conversion. mplstereonet uses the RHR notation by default.
Mixing equal-area (Schmidt) and equal-angle (Wulff) projections, distorting density analysis
Contour diagrams in structural geology must always be performed on equal-area (Schmidt net) projections. Equal-angle (Wulff net) projections do not preserve area, causing density distortion. mplstereonet uses equal-area projection by default, so do not change this setting. In Allmendinger Stereonet, verify that the Equal-Area mode is selected.
Plotting on upper hemisphere projection, making results incomparable with published literature
Structural geology universally uses lower hemisphere projection as the standard. Upper hemisphere projection is only used in some mineralogy/crystallography contexts. Always verify that your stereonet software is set to "Lower Hemisphere." mplstereonet uses lower hemisphere projection by default.
Over-relying on contour diagrams with insufficient data, leading to incorrect interpretations
Contour diagrams require a minimum of 25-30 measurement points to be statistically meaningful. With fewer than 10 data points, use only great circle and pole plots, and treat contour diagrams as supplementary. The Kamb method auto-corrects for sample size, but artifacts can still appear with small datasets.
Computing only the mean pole instead of fitting a girdle when determining fold axes
fit_pole() calculates the mean direction of a single cluster, while fit_girdle() fits a great circle (girdle) distribution. To determine a fold axis, you must use fit_girdle() to find the Pi girdle and then compute its pole. Using fit_pole() on all data returns an intermediate value between the two limbs, which is completely different from the actual fold axis.