How to Map Earthquake Focal Mechanisms in Python
On this page
Problem statement
An epicentre says where an earthquake happened; a focal mechanism says how the fault moved: strike-slip, thrust or normal, and in which direction. Seismologists draw it as a "beach ball", a stereographic picture of the radiation pattern, and maps of beach balls show the stress regime of a region at a glance. The USGS publishes a moment tensor for most earthquakes above magnitude 5, but turning the six tensor components into beach balls, faulting classes and stress directions takes some geometry that usually hides inside a dedicated library.
This page downloads the USGS moment tensors for the first month of the February 2023 Kahramanmaraş, Türkiye, earthquake sequence, computes each event's moment magnitude from its tensor, classifies the faulting style from the principal axes, draws beach balls from scratch with NumPy, and measures the direction of maximum compression.
Quick answer
Classify each mechanism by the plunges of its pressure (P), tension (T) and null (N) axes:
def style(p, t, n):
if t >= 50: return "thrust"
if p >= 52: return "normal"
if n >= 45: return "strike-slip"
return "oblique"
M["style"] = [style(*x) for x in M[["p_plunge", "t_plunge", "n_plunge"]].values]
print(f"{len(M)} moment tensors, M{M.mag.min()} to M{M.mag.max()}")
print(M["style"].value_counts().to_dict())58 moment tensors, M4.5 to M7.8
{'strike-slip': 35, 'normal': 14, 'oblique': 6, 'thrust': 3}
35 of the 58 mechanisms are strike-slip, as expected on the East Anatolian fault system. But 14 are normal-faulting events and 3 thrusts, and 6 fall between classes: a large sequence activates faults of several kinds around the main rupture.
From tensor to magnitude, beach balls, compression direction
1. From tensor to magnitude
The scalar seismic moment is the norm of the tensor, M0 = √(Σ Mij² / 2), in newton metres, and the moment magnitude is Mw = ⅔ (log10 M0 − 9.1). Check it against the catalogue magnitude:
comp = M[["mrr", "mtt", "mpp"]].pow(2).sum(axis=1) + 2 * M[["mrt", "mrp", "mtp"]].pow(2).sum(axis=1)
M["m0"] = np.sqrt(comp / 2)
M["mw_tensor"] = 2 / 3 * (np.log10(M.m0) - 9.1)
print(f"largest M0: {M.m0.max():.3e} N m (Mw {M.mw_tensor.max():.2f})")
print(f"Mw from the tensor minus catalogue magnitude: median {(M.mw_tensor - M.mag).median():+.2f}, largest {(M.mw_tensor - M.mag).abs().max():.2f}")
print(f"share of the total moment released by the two largest events: {M.m0.nlargest(2).sum() / M.m0.sum():.1%}")largest M0: 5.389e+20 N m (Mw 7.75)
Mw from the tensor minus catalogue magnitude: median -0.01, largest 0.07
share of the total moment released by the two largest events: 99.2%
The M7.8's tensor has a scalar moment of 5.389 × 10²⁰ N m, Mw 7.75. Across all 58 events the magnitude computed from the tensor matches the catalogue magnitude to within 0.07, with a median difference of -0.01. And the two largest events, the M7.8 and the M7.5, released 99.2% of all the moment of the month: the other 56 earthquakes together amount to less than one percent.
2. Draw a beach ball
A beach ball is the lower focal hemisphere seen from above in an equal-area projection, shaded where first motions are compressional. For each point on the hemisphere, the sign of nᵀMn, with n the ray direction in ComCat's Up-South-East coordinates, says which:
def beachball(row, n=201):
T = np.array([[row.mrr, row.mrt, row.mrp], [row.mrt, row.mtt, row.mtp], [row.mrp, row.mtp, row.mpp]])
x, y = np.meshgrid(np.linspace(-1, 1, n), np.linspace(-1, 1, n))
r = np.hypot(x, y); inside = r <= 1
inc = 2 * np.arcsin(np.clip(r / np.sqrt(2), 0, 1)) # equal-area: r = sqrt(2) sin(i/2)
az = np.arctan2(x, y) # azimuth clockwise from north
ray = np.stack([-np.cos(inc), -np.sin(inc) * np.cos(az), np.sin(inc) * np.sin(az)]) # up, south, east
amp = np.einsum("i...,ij,j...->...", ray, T, ray)
return np.where(inside, amp > 0, np.nan)
bb = beachball(M.loc[M.mag.idxmax()])
print(f"M7.8 beach ball: compressional share of the hemisphere {np.nanmean(bb):.2f}")
print(f"its nodal planes: strike/dip/rake {M.loc[M.mag.idxmax(), ['np1_strike', 'np1_dip', 'np1_rake']].astype(float).round(0).tolist()} and "
f"{M.loc[M.mag.idxmax(), ['np2_strike', 'np2_dip', 'np2_rake']].astype(float).round(0).tolist()}")M7.8 beach ball: compressional share of the hemisphere 0.46
its nodal planes: strike/dip/rake [318.0, 89.0, -179.0] and [228.0, 89.0, -1.0]
The M7.8's beach ball is the classic strike-slip cross: four quadrants, 46% of the hemisphere compressional, with two nearly vertical nodal planes. Drawn from the full tensor, which is 81% double couple, the nodal lines bend near the centre; a best-double-couple drawing, as on USGS event pages, has straight ones. One strikes 228° with a rake of -1°, left-lateral slip on a plane running north-east to south-west, parallel to the East Anatolian fault; the other, 318° with a rake of -179°, would be right-lateral slip on a plane at right angles to it.
3. The direction of compression
For strike-slip and thrust events the P axis is close to horizontal, and its azimuth is the direction of greatest compression. Average it as an axis, doubling the angles so 10° and 190° count as the same direction:
ss = M[M["style"] == "strike-slip"]
a = np.radians(2 * ss.p_az)
mean_p = (np.degrees(np.arctan2(np.sin(a).mean(), np.cos(a).mean())) / 2) % 180
R = np.hypot(np.sin(a).mean(), np.cos(a).mean())
print(f"strike-slip events: {len(ss)}; mean P-axis azimuth {mean_p:.0f}° (or {mean_p + 180:.0f}°); resultant length {R:.2f}")
print(f"P-axis plunges: median {ss.p_plunge.median():.0f}°, all within {ss.p_plunge.max():.0f}° of horizontal")strike-slip events: 35; mean P-axis azimuth 19° (or 199°); resultant length 0.71
P-axis plunges: median 15°, all within 44° of horizontal
For the 35 strike-slip events, the axial mean P-axis azimuth is 19°, compression from north-north-east to south-south-west, with a resultant length of 0.71 on a scale from 0 (no preferred direction) to 1 (all identical). Their P axes plunge a median of 15°, close to horizontal.
How much of each mechanism is a double couple
The data used on this page
The USGS moment tensors for every earthquake of magnitude 4.5 or more between 35° and 40°N and 35° and 40°E from 6 February to 6 March 2023 that has one, from the FDSN event service with producttype=moment-tensor, and each event's detail record for the tensor itself. The preferred moment-tensor product comes first in each record.
import json
from pathlib import Path
import numpy as np
import pandas as pd
import requests
Path("data/mt").mkdir(parents=True, exist_ok=True)
Q = "https://earthquake.usgs.gov/fdsnws/event/1/query"
def cached(path, **params):
p = Path(path)
if not p.exists():
r = requests.get(Q, params=params, timeout=300)
r.raise_for_status()
p.write_text(r.text)
return json.loads(p.read_text())
lst = cached("data/turkey_2023_mt_list.json", format="geojson", starttime="2023-02-06", endtime="2023-03-07", minmagnitude=4.5,
minlatitude=35, maxlatitude=40, minlongitude=35, maxlongitude=40, producttype="moment-tensor")
rows = []
for f in lst["features"]:
d = cached(f"data/mt/{f['id']}.json", format="geojson", eventid=f["id"])
q = d["properties"]["products"]["moment-tensor"][0]["properties"]
g = lambda k: float(q[k])
rows.append({"id": f["id"], "mag": d["properties"]["mag"], "lon": d["geometry"]["coordinates"][0], "lat": d["geometry"]["coordinates"][1],
"depth": d["geometry"]["coordinates"][2], "source": q.get("derived-magnitude-type"),
**{c: g(f"tensor-{c}") for c in ("mrr", "mtt", "mpp", "mrt", "mrp", "mtp")},
**{f"np{i}_{k}": g(f"nodal-plane-{i}-{k}") for i in (1, 2) for k in ("strike", "dip", "rake")},
"p_az": g("p-axis-azimuth"), "p_plunge": g("p-axis-plunge"), "t_plunge": g("t-axis-plunge"), "n_plunge": g("n-axis-plunge"),
"dc": g("percent-double-couple") if "percent-double-couple" in q else np.nan})
M = pd.DataFrame(rows)
print(f"{len(M)} events with a moment tensor; types {M.source.value_counts().to_dict()}")58 events with a moment tensor; types {'Mwr': 39, 'Mww': 19}
58 events in the month have a USGS moment tensor: 39 regional moment tensors (Mwr) for the smaller events and 19 W-phase solutions (Mww) for the larger.
Example 1: how much of each mechanism is a double couple
A slip on a single planar fault gives a pure double couple. The percent double couple measures how close each tensor is:
print(f"double couple: median {M.dc.median():.0%}; events below 70%: {(M.dc < 0.7).sum()} of {M.dc.notna().sum()}")
print(M.groupby("style").dc.median().round(2).to_string())double couple: median 78%; events below 70%: 22 of 58
style
normal 0.81
oblique 0.82
strike-slip 0.77
thrust 0.78
The median event is 78% double couple, and 22 of the 58 are below 70%. The share is similar for every faulting style, 77% to 82%, so in this sequence the double-couple share does not set one kind of faulting apart from another.
What a moment tensor describes
A moment tensor describes the forces at the source as a 3×3 symmetric matrix. For slip on a fault it is a double couple: two pairs of opposite forces, with no net torque, whose radiation pattern divides the focal sphere into four quadrants, two compressional and two dilatational, separated by two perpendicular nodal planes. One of those is the fault and the other the auxiliary plane; the tensor alone cannot say which, but aftershock alignments and mapped faults usually can. The P and T axes bisect the quadrants: the direction of greatest shortening and of greatest extension at the source.
The faulting style follows from which axis is vertical. A vertical T axis means the crust is shortened horizontally and thickened, a thrust; a vertical P axis means horizontal extension, a normal fault; a vertical null axis means both P and T are horizontal, strike-slip. The plunge thresholds on this page are a simplified form of the classes used for stress maps; borderline events are best called oblique. The 2023 sequence ruptured the East Anatolian fault zone, a left-lateral strike-slip boundary between the Anatolian and Arabian plates, and its mechanisms say so: nearly vertical nodal planes, horizontal P axes and one plane parallel to the fault.
Coordinates, fault planes and symbol size
- Coordinates: ComCat tensors are in Up-South-East (r, θ, φ) coordinates, the Harvard CMT convention; other sources use North-East-Down, and mixing them flips the beach ball.
- Which plane is the fault: choose with geology or aftershocks; the moment tensor treats both nodal planes equally.
- Beach ball size on maps: scale by magnitude, not moment, or the M7.8 hides everything else; and draw circles in screen units, not map units, so they stay round in a projected map.
- Several solutions: an event can carry W-phase, body-wave and regional tensors; the first product in the record is the preferred one.
Internal links
- Fixing USGS earthquake times that show up as 1970 — the same sequence, from the same service
- How to draw an earthquake depth cross-section in Python — hypocentres in depth
- How to map earthquakes against plate boundaries in Python — the plate boundary context
FAQ
What is a focal mechanism?
A description of the fault orientation and slip direction of an earthquake, derived from the radiation pattern of its seismic waves and usually drawn as a beach ball.
How do I get moment tensors from the USGS in Python?
Query the event service with producttype=moment-tensor, then read each event's detail GeoJSON and take the tensor-mrr to tensor-mtp properties of its first moment-tensor product.
How do I draw a beach ball without ObsPy?
Evaluate nᵀMn over a grid of ray directions on the lower hemisphere in an equal-area projection and shade the positive, compressional, values, as in step 2.
What does percent double couple mean?
How much of the moment tensor is explained by slip on a single planar fault; low values indicate complex ruptures, volume changes or poorly resolved solutions.