contour_angelier_data.py
¶
Reproduce Figure 5 from Vollmer, 1995 to illustrate different density contouring methods.
import matplotlib.pyplot as plt
import mplstereonet
import parse_angelier_data
def plot(ax, strike, dip, rake, **kwargs):
ax.rake(strike, dip, rake, 'ko', markersize=2)
ax.density_contour(strike, dip, rake, measurement='rakes', linewidths=1,
cmap='jet', **kwargs)
# Load data from Angelier, 1979
strike, dip, rake = parse_angelier_data.load()
# Setup a subplot grid
fig, axes = mplstereonet.subplots(nrows=3, ncols=4)
# Hide azimuth tick labels
for ax in axes.flat:
ax.set_azimuth_ticks([])
contours = [range(2, 18, 2), range(1, 21, 2), range(1, 22, 2)]
# "Standard" Kamb contouring with different confidence levels.
for sigma, ax, contour in zip([3, 2, 1], axes[:, 0], contours):
# We're reducing the gridsize to more closely match a traditional
# hand-contouring grid, similar to Kamb's original work and Vollmer's
# Figure 5. `gridsize=10` produces a 10x10 grid of density estimates.
plot(ax, strike, dip, rake, method='kamb', sigma=sigma,
levels=contour, gridsize=10)
# Kamb contouring with inverse-linear smoothing (after Vollmer, 1995)
for sigma, ax, contour in zip([3, 2, 1], axes[:, 1], contours):
plot(ax, strike, dip, rake, method='linear_kamb', sigma=sigma,
levels=contour)
template = r'$E={}\sigma$ Contours: ${}\sigma,{}\sigma,\ldots$'
ax.set_xlabel(template.format(sigma, *contour[:2]))
# Kamb contouring with exponential smoothing (after Vollmer, 1995)
for sigma, ax, contour in zip([3, 2, 1], axes[:, 2], contours):
plot(ax, strike, dip, rake, method='exponential_kamb', sigma=sigma,
levels=contour)
# Title the different methods
methods = ['Kamb', 'Linear\nSmoothing', 'Exponential\nSmoothing']
for ax, title in zip(axes[0, :], methods):
ax.set_title(title)
# Hide top-right axis... (Need to implement Diggle & Fisher's method)
axes[0, -1].set_visible(False)
# Schmidt contouring (a.k.a. 1%)
plot(axes[1, -1], strike, dip, rake, method='schmidt', gridsize=25,
levels=range(3, 20, 3))
axes[1, -1].set_title('Schmidt')
axes[1, -1].set_xlabel(r'Contours: $3\%,6\%,\ldots$')
# Raw data.
axes[-1, -1].set_azimuth_ticks([])
axes[-1, -1].rake(strike, dip, rake, 'ko', markersize=2)
axes[-1, -1].set_xlabel('N={}'.format(len(strike)))
plt.show()