Where do earthquakes and volcanoes happen — and why there?¶
EPS 88 · PyEarth. Open your own copy on DataHub: click here.
Two files today. One holds every earthquake of magnitude 5.5 and above located from 2000 to the end of 2025; the other, every volcano the Smithsonian Institution lists as active in about the last twelve thousand years. Neither file is a map, and yet each one draws one.
You will draw those maps from longitude and latitude, lay the plate boundaries over them, and colour every earthquake by how deep it was. Then you will turn the Andes on their side and watch a sinking plate appear. Last, you will count earthquakes by size and compare them with eruptions.
Every place you write something opens with a pencil icon and the words Your turn, and is
followed by a cell that already has the shape of the answer written in, with ... wherever your
code goes. Replace every ..., run every cell so your answers and figures are
saved in the file, then download this notebook itself — the .ipynb file — and upload it
to Gradescope. Not a PDF: the marking reads your notebook, and a PDF cannot be read.
In JupyterLab: File ▸ Download, or right-click the file in the left-hand panel and choose
Download.
Two habits from the first minute. A cell runs when you press Shift+Enter, and the notebook remembers everything it has already run — so when something breaks and you cannot see why, Kernel → Restart Kernel and Run All Cells throws the memory away and rebuilds it from the top. That is never the wrong thing to do.
What you'll be able to do¶
The science. Say where earthquakes and volcanoes happen and which plate boundaries they sit on, show from the side that the deep earthquakes lie on a sinking plate, and read counts that span a factor of a thousand.
The code. Pick the rows of a table with a condition, and join several conditions with & ·
draw a map as plt.scatter of longitude against latitude, with layers on top · colour every dot by
a column with c= and add a key with plt.colorbar · draw a side view, with depth down the page ·
put two panels side by side with plt.subplot, and count an axis in factors of ten with
plt.yscale("log").
Eight places where you write something: five in class, three at home. Each is headed
Your turn. The cell under it has every line already started, with ... where your code goes
and a comment saying what goes there. Where the cell after it is a self-check, run it: it
tells you whether something is off, and you do not need to read it.
- Where do earthquakes and volcanoes happen?
- What draws those lines?
- Where are the deep earthquakes, and why there?
- How many small ones for each big one?
Setup¶
Run the next cell once; you are not expected to follow it. It reads the USGS earthquake catalogue and the Smithsonian's volcano and eruption tables from the live archives, falling back to the copy stored with the course, and three files that ship with the course: the coastlines, the plate boundaries, and a model of the sinking plate under South America.
Every count printed in this notebook came from the copy stored with the course. The archives are revised all the time, so a slightly different number is the archive having moved, not you having made a mistake.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
# house style, set once, so every plot cell below holds only what matters
plt.rcParams.update({"figure.dpi": 110,
"axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True})
CACHE = "https://raw.githubusercontent.com/AI4EPS/EPS88_PyEarth/main/data"
def load(url, cache_name):
"""Read the live source; fall back to the copy stored with the course."""
# Ask the live archive first. If it is down, or you are offline, read the copy stored with
# the course instead, so the notebook still runs.
try:
return pd.read_csv(url)
except Exception as e:
print("live source unreachable, using the cached copy:", type(e).__name__)
return pd.read_csv(CACHE + "/" + cache_name)
USGS = ("https://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&orderby=time-asc"
"&starttime=2000-01-01&endtime=2026-01-01&minmagnitude=5.5")
GVP = ("https://webservices.volcano.si.edu/geoserver/GVP-VOTW/ows?service=WFS&version=1.0.0"
"&request=GetFeature&typeName=GVP-VOTW:Smithsonian_VOTW_Holocene_")
quakes = load(USGS, "week04_2000-01-01_2026-01-01_M5.5.csv")
quakes = quakes[quakes["type"] == "earthquake"] # drops a 2017 nuclear test and a 2022 eruption near Tonga
volcanoes = load(GVP + "Volcanoes&outputFormat=csv", "week04_gvp_volcanoes.csv")
eruptions = load(GVP + "Eruptions&outputFormat=csv", "week04_gvp_eruptions.csv")
# the first words of each volcano's setting: "Subduction zone", "Rift zone" or "Intraplate"
volcanoes["setting"] = volcanoes["Tectonic_Setting"].str.split(" / ").str[0]
# these three ship with the course, so they are read straight from it
coast = pd.read_csv(CACHE + "/coastlines.csv")
boundaries = pd.read_csv(CACHE + "/plate_boundaries.csv")
slab = pd.read_csv(CACHE + "/slab2_south_america.csv")
print("earthquakes:", len(quakes), " volcanoes:", len(volcanoes), " eruptions:", len(eruptions))
live source unreachable, using the cached copy: URLError live source unreachable, using the cached copy: URLError
live source unreachable, using the cached copy: URLError
earthquakes: 12847 volcanoes: 1214 eruptions: 11089
1. Where do earthquakes and volcanoes happen?¶
Longitude runs from -180 to 180 around the world and latitude from -90 to 90 up it, so a scatter
plot of one against the other is already a map. plt.xlim and plt.ylim fix the edges at the
whole world, and plt.gca().set_aspect("equal") draws a degree across the same length as a degree
up. The coastline goes on with one plt.plot, as in the first week.
plt.figure(figsize=(9, 4.5))
plt.scatter(quakes["longitude"], quakes["latitude"], s=1, color="crimson")
plt.plot(coast["lon"], coast["lat"], color="0.6", lw=0.6)
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title(str(len(quakes)) + " earthquakes of magnitude 5.5 and above, 2000-2025")
plt.show()
The dots are not spread over the planet; they are lines. One rings the Pacific: Alaska, Japan, Indonesia, New Zealand, the whole west coast of the Americas. Another runs down the middle of the Atlantic, far from any coast, drawn by nothing but earthquakes.
This map stretches the far north and south: it draws a degree of longitude the same width at every latitude, while on the ground the degrees close in towards the poles. For finding where things are, that is a fair trade.
2. What draws those lines?¶
The lines are the edges of the tectonic plates. plate_boundaries.csv holds those edges as
longitudes and latitudes, and its kind column says how the two plates meet:
- ridge — they pull apart, and new ocean floor forms between them;
- transform — they slide past each other;
- trench — one plate bends down and sinks under the other. The file also files the collision fronts of the Himalaya and the Zagros under this word, where no ocean floor sinks.
Each kind is the rows of the table whose kind is that word. boundaries["kind"] == "ridge" asks
every row the same question and hands back a mask of True and False, and boundaries[mask] keeps
the rows where it is True — last week's move, on a table.
ridges = boundaries[boundaries["kind"] == "ridge"]
transforms = boundaries[boundaries["kind"] == "transform"]
trenches = boundaries[boundaries["kind"] == "trench"]
plt.figure(figsize=(9, 4.5))
plt.scatter(quakes["longitude"], quakes["latitude"], s=1, color="0.75")
plt.plot(coast["lon"], coast["lat"], color="0.6", lw=0.5)
plt.plot(ridges["lon"], ridges["lat"], color="tab:blue", lw=1, label="ridge")
plt.plot(transforms["lon"], transforms["lat"], color="tab:green", lw=1, label="transform")
plt.plot(trenches["lon"], trenches["lat"], color="black", lw=1, label="trench")
plt.legend(loc="lower left")
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title(str(len(quakes)) + " earthquakes and the plate boundaries")
plt.show()
The boundaries land on the dots: the blue ridge down the Atlantic is the line the earthquakes drew on their own, and black trenches trace most of the rim of the Pacific.
Do the volcanoes sit on the same lines, and on which kind? The Smithsonian classifies the setting
of every volcano, and the setup cell kept the first words of it in a column called setting:
"Subduction zone" where a plate sinks at a trench, "Rift zone" where plates pull apart, and
"Intraplate" for a volcano far from any boundary.
✏️ Your turn 1¶
Which kind of plate boundary holds most of the volcanoes?
Make three tables from volcanoes, one for each setting, the same way ridges was made from
boundaries above. Then draw each one as triangles in its own colour, on top of the ridges and
trenches. The lines at the end print what percent of all the volcanoes each setting holds.
Use these names, because the self-check looks for them: subduction, rift and intraplate.
subduction = volcanoes[volcanoes["setting"] == "Subduction zone"]
rift = volcanoes[volcanoes["setting"] == "Rift zone"]
intraplate = volcanoes[volcanoes["setting"] == "Intraplate"]
plt.figure(figsize=(9, 4.5))
plt.plot(coast["lon"], coast["lat"], color="0.7", lw=0.5)
plt.plot(ridges["lon"], ridges["lat"], color="0.4", lw=0.8, ls="--", label="ridge")
plt.plot(trenches["lon"], trenches["lat"], color="black", lw=0.8, label="trench")
plt.scatter(subduction["Longitude"], subduction["Latitude"], marker="^", color="crimson", label="subduction zone")
plt.scatter(rift["Longitude"], rift["Latitude"], marker="^", color="dodgerblue", label="rift zone")
plt.scatter(intraplate["Longitude"], intraplate["Latitude"], marker="^", color="orange", label="intraplate")
plt.legend(loc="lower left")
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title(str(len(volcanoes)) + " volcanoes, by setting")
plt.show()
print("subduction zone:", round(100 * len(subduction) / len(volcanoes)), "percent")
print("rift zone: ", round(100 * len(rift) / len(volcanoes)), "percent")
print("intraplate: ", round(100 * len(intraplate) / len(volcanoes)), "percent")
subduction zone: 69 percent rift zone: 18 percent intraplate: 12 percent
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (subduction["setting"] == "Subduction zone").all() and (rift["setting"] == "Rift zone").all() and (intraplate["setting"] == "Intraplate").all(), "each table should keep only the rows whose setting is that text: compare volcanoes[\"setting\"] with it"
assert 800 < len(subduction) < 880 and 200 < len(rift) < 240 and 130 < len(intraplate) < 170, "each table should keep every volcano with that setting: pick the rows from volcanoes, the whole table"
print(f"✓ Your turn 1 — {len(subduction)} volcanoes at subduction zones, {len(rift)} at rifts, {len(intraplate)} intraplate")
✓ Your turn 1 — 838 volcanoes at subduction zones, 221 at rifts, 149 intraplate
69 percent of the catalogued volcanoes sit at subduction zones, and there that is a mechanism rather than a coincidence. The sinking plate carries water down, locked inside its minerals. Once the plate is roughly a hundred kilometres deep, the water is driven out into the hot mantle above it, and water lowers the temperature at which rock melts, so mantle that would have stayed solid melts. The melt rises and builds a line of volcanoes set back from the trench.
The blue rift volcanoes (18 percent) sit where plates pull apart: on the ridge through Iceland and along the rift down East Africa, and the orange ones (12 percent) sit far from any line, like Hawaii in the middle of the Pacific. The catalogue lists volcanoes people can find and name, mostly on land, while most of the magma Earth erupts comes out along the mid-ocean ridges, under kilometres of water, where eruptions are rarely seen.
3. Where are the deep earthquakes, and why there?¶
Every earthquake in the catalogue also has a depth, in kilometres. Seismologists sort them into
three classes: shallow, down to 70 km; intermediate, 70 to 300 km; and
deep, below 300 km. Deep earthquakes are strange: a few hundred kilometres down, rock is hot
and under enormous pressure, and it should not break the way rock near the surface does. So where
they happen is a clue.
# ── Checkpoint ── run this if you restarted the kernel or fell behind ──
trenches = boundaries[boundaries["kind"] == "trench"]
✏️ Your turn 2¶
How many earthquakes are deep, and where are they?
A map can show depth if each dot's colour carries it. plt.scatter takes c= for the column to
colour by and cmap= for the colour scheme; vmin and vmax set which depths the two ends of the
colour scale mean. Draw every earthquake coloured by its depth, from 0 km to 600 km, add a key to
the colours, and lay the trenches on top as black lines.
plt.figure(figsize=(9, 4.5))
plt.scatter(quakes["longitude"], quakes["latitude"], s=2,
c=quakes["depth"], cmap="plasma_r", vmin=0, vmax=600)
plt.colorbar(label="depth (km)")
plt.plot(trenches["lon"], trenches["lat"], color="black", lw=0.8)
plt.plot(coast["lon"], coast["lat"], color="0.6", lw=0.6)
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title(str(len(quakes)) + " earthquakes, coloured by depth")
plt.show()
deep = quakes[quakes["depth"] > 300]
print(len(deep), "of", len(quakes), "earthquakes are deeper than 300 km:", round(100 * len(deep) / len(quakes)), "percent")
785 of 12847 earthquakes are deeper than 300 km: 6 percent
Only 6 percent of the earthquakes are deep, and they are not scattered: nearly every cluster of dark dots lies beside a black trench line, and along the ridge down the Atlantic every dot is pale. But the dark dots are not on the trench lines. They sit off to one side of them.
A map cannot show which way that is going, so turn one stretch of trench on its side. First, pick the stretch: the earthquakes under southern Peru, Chile, Bolivia and Argentina.
✏️ Your turn 3¶
Pick the earthquakes in one box.
A box needs four conditions, one for each edge: latitude greater than -30 and less than -15,
longitude greater than -80 and less than -60. Each condition is a mask; join them with &, which
keeps a row only when all four are True. Put each condition in its own brackets: Python does &
before > and <, so without them it would join the wrong pieces and stop with an error. The map under the cell shows the box you picked.
Use these names, because the self-check looks for them: andes.
andes = quakes[(quakes["latitude"] > -30) & (quakes["latitude"] < -15)
& (quakes["longitude"] > -80) & (quakes["longitude"] < -60)]
print(len(andes), "earthquakes in the box")
plt.figure(figsize=(9, 4.5))
plt.scatter(quakes["longitude"], quakes["latitude"], s=1, color="0.8")
plt.scatter(andes["longitude"], andes["latitude"], s=6, c=andes["depth"], cmap="plasma_r", vmin=0, vmax=600)
plt.colorbar(label="depth (km)")
plt.plot(coast["lon"], coast["lat"], color="0.6", lw=0.6)
plt.plot(trenches["lon"], trenches["lat"], color="black", lw=1)
plt.plot([-80, -60, -60, -80, -80], [-30, -30, -15, -15, -30], color="black", ls="--")
plt.xlim(-95, -30)
plt.ylim(-60, 15)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.title("the box: " + str(len(andes)) + " earthquakes, coloured by depth")
plt.show()
437 earthquakes in the box
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (andes["latitude"] > -30).all() and (andes["latitude"] < -15).all() and (andes["longitude"] > -80).all() and (andes["longitude"] < -60).all(), "every row of andes should sit inside the box: check each of the four conditions"
assert len(andes) == ((quakes["latitude"] > -30) & (quakes["latitude"] < -15) & (quakes["longitude"] > -80) & (quakes["longitude"] < -60)).sum(), "andes should keep every earthquake inside the box: join all four conditions with &"
print(f"✓ Your turn 3 — {len(andes)} earthquakes in the box")
✓ Your turn 3 — 437 earthquakes in the box
✏️ Your turn 4¶
Which side of the coast are the deep earthquakes on, and how far from the shallow ones?
Draw the box from the side: longitude across, depth down the page. plt.ylim(700, 0) puts 0 km at
the top, as in the Earth. Split the box into shallow earthquakes (no deeper than 70 km)
and deep ones (deeper than 300 km), draw them in blue and red over the grey of all of them,
and mark the middle of each colour with a dashed line at its median longitude. Read the answer off
your figure; the lines at the end print the two medians and the distance between them.
Use these names, because the self-check looks for them: shallow and deep_andes.
shallow = andes[andes["depth"] <= 70]
deep_andes = andes[andes["depth"] > 300]
plt.figure(figsize=(9, 4.5))
plt.scatter(andes["longitude"], andes["depth"], s=5, color="0.8", label="all")
plt.scatter(shallow["longitude"], shallow["depth"], s=8, color="tab:blue", label="shallow, 0-70 km")
plt.scatter(deep_andes["longitude"], deep_andes["depth"], s=12, color="tab:red", label="deep, below 300 km")
plt.axvline(shallow["longitude"].median(), color="tab:blue", ls="--")
plt.axvline(deep_andes["longitude"].median(), color="tab:red", ls="--")
plt.axvline(-70.2, color="black", lw=1.5, label="the coast at 22.5°S")
plt.ylim(700, 0)
plt.xlim(-80, -60)
plt.xlabel("longitude (degrees east)")
plt.ylabel("depth (km)")
plt.title("The Andes from the side: " + str(len(andes)) + " earthquakes")
plt.legend(loc="lower left")
plt.show()
print("shallow median longitude:", round(shallow["longitude"].median(), 1))
print("deep median longitude: ", round(deep_andes["longitude"].median(), 1))
print("deep minus shallow:", round(deep_andes["longitude"].median() - shallow["longitude"].median(), 1), "degrees")
shallow median longitude: -71.0 deep median longitude: -63.3 deep minus shallow: 7.7 degrees
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert shallow["depth"].max() <= 70 and deep_andes["depth"].min() > 300, "shallow keeps the rows no deeper than 70 km, and deep_andes the rows deeper than 300 km"
assert len(shallow) + len(deep_andes) < len(andes) and len(deep_andes) > 20, "take shallow and deep_andes from andes, the box, not from the whole catalogue"
print(f"✓ Your turn 4 — {len(shallow)} shallow and {len(deep_andes)} deep earthquakes in the box")
✓ Your turn 4 — 214 shallow and 44 deep earthquakes in the box
The shallow earthquakes crowd the coast, and the deep ones sit about 8 degrees further east, under the continent. In between, the grey dots are not a cloud: they lie along a line that starts near the coast and slopes down to the east.
Seismologists have mapped that surface. Slab2, a model published by the US Geological Survey,
gives the depth of the top of the sinking plate under the world's major subduction zones, built from earthquake
locations like these together with other measurements. The setup cell loaded its South American
part as slab, with columns lon, lat and depth. One latitude of it, picked with a mask, is a
line across the box.
profile = slab[slab["lat"] == -22.5] # the top of the sinking plate, along 22.5 degrees south
plt.figure(figsize=(9, 4.5))
plt.scatter(andes["longitude"], andes["depth"], s=6, c=andes["depth"], cmap="plasma_r", vmin=0, vmax=600)
plt.plot(profile["lon"], profile["depth"], color="black", lw=2, label="top of the sinking plate (Slab2), at 22.5°S")
plt.axvline(-70.2, color="0.4", ls=":", label="the coast at 22.5°S")
plt.ylim(700, 0)
plt.xlim(-80, -60)
plt.xlabel("longitude (degrees east)")
plt.ylabel("depth (km)")
plt.title("The Andes from the side, with the plate they sit on: " + str(len(andes)) + " earthquakes")
plt.legend(loc="lower left")
plt.show()
The earthquakes below 70 km roughly follow the black line: under the coast the plate is some 37 km down, and it keeps sinking eastwards to more than 500 km. The deepest cluster lies a little further south, near 27°S, where Slab2 takes the plate to about 620 km. This is the Nazca plate, the floor of the eastern Pacific, going down beneath South America. A slab is a plate that has sunk into the mantle, usually at a trench; the deep earthquakes happen inside slabs. A slab stays colder than the mantle around it for millions of years, which is thought to be why it can still break at those depths; exactly how rock breaks 600 km down is still argued about.
4. How many small ones for each big one?¶
That was where; now how big. Volcanologists score an eruption on the Volcanic Explosivity
Index (VEI), a whole number from 0 up, and the eruption table carries it in ExplosivityIndexMax.
An index counts steps, not amounts: from VEI 2 up, each step is about ten times the ash and pumice
thrown out, and each step of magnitude about 32 times the energy. So Tambora's great eruption
(VEI 7) threw out about a hundred times the ash and pumice of Mount St Helens in 1980 (VEI 5). The
tenfold step is how Newhall and Self defined the VEI in 1982; the 32 comes from Gutenberg and
Richter's relation between magnitude and energy.
The cell below counts the eruptions in each class, and draws the counts twice, side by side.
plt.subplot(1, 2, k) splits a figure into one row of two panels and draws in panel k. The left
panel counts its axis the ordinary way; plt.yscale("log") makes the right one count in factors
of ten.
vei_counts = eruptions["ExplosivityIndexMax"].value_counts().sort_index()
plt.figure(figsize=(9, 4.5))
plt.subplot(1, 2, 1)
plt.bar(vei_counts.index, vei_counts.values, color="darkorange")
plt.xlabel("VEI (an index, no units)")
plt.ylabel("number of eruptions")
plt.title("ordinary axis")
plt.subplot(1, 2, 2)
plt.bar(vei_counts.index, vei_counts.values, color="darkorange")
plt.yscale("log")
plt.xlabel("VEI (an index, no units)")
plt.title("log axis, " + str(vei_counts.sum()) + " eruptions with a VEI")
plt.show()
print(vei_counts)
ExplosivityIndexMax 0.0 1019 1.0 1441 2.0 4030 3.0 1172 4.0 516 5.0 181 6.0 52 7.0 7 Name: count, dtype: int64
On the left the axis climbs to fit the 4,030 eruptions at VEI 2, and the big ones vanish: the 7 at VEI 7 are a hairline. On the right every class is readable. When the values span factors of a thousand, plot the exponents instead and a curve becomes a line. From VEI 2 to VEI 6 the bars step down in a roughly straight line, each class between 2.3 and 3.5 times smaller than the one before. The two smallest classes fall short of that line, most likely because a small eruption long ago, far from anyone, left nothing to record.
✏️ Your turn 5¶
Do earthquakes fall off the same way, and how many times rarer is each step of magnitude?
Make the same two panels for the earthquakes' magnitudes, with plt.hist in place of plt.bar:
edges gives one bar for every half step of magnitude. np.histogram, from last week, hands back
the same counts as numbers. A whole step of magnitude is two bars, so per_step is the first
bar's count divided by the third's.
Use these names, because the self-check looks for them: counts and per_step.
edges = [5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5] # half a step of magnitude per bar
counts, edges = np.histogram(quakes["mag"], bins=edges)
plt.figure(figsize=(9, 4.5))
plt.subplot(1, 2, 1)
plt.hist(quakes["mag"], bins=edges, color="crimson")
plt.xlabel("magnitude")
plt.ylabel("number of earthquakes")
plt.title("ordinary axis")
plt.subplot(1, 2, 2)
plt.hist(quakes["mag"], bins=edges, color="crimson")
plt.yscale("log")
plt.xlabel("magnitude")
plt.title("log axis, " + str(len(quakes)) + " earthquakes")
plt.show()
per_step = counts[0] / counts[2]
print("magnitude 5.5 to 6:", counts[0], " 6.5 to 7:", counts[2], " so each whole step is", round(per_step, 1), "times rarer")
magnitude 5.5 to 6: 8973 6.5 to 7: 799 so each whole step is 11.2 times rarer
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(counts) == 8 and counts[0] > 8000 and counts[-1] < 10, "counts should hold the eight bar heights of quakes[\"mag\"], cut at edges"
assert 9 < per_step < 13, "per_step is the first bar's count divided by the third's: counts[0] / counts[2]"
print(f"✓ Your turn 5 — each whole step of magnitude is {round(per_step, 1)} times rarer")
✓ Your turn 5 — each whole step of magnitude is 11.2 times rarer
Only 28 of 12,847 reach magnitude 8. On the log axis the earthquake bars fall in a roughly straight line too: below magnitude 8, each whole step holds about 7 to 11 times fewer earthquakes than the step below. Eruptions and earthquakes are rare in the same way: every step up the scale is several times rarer, by a roughly steady factor. Since each step multiplies the size, the count falls as a power of the size: a power law. On an ordinary axis it hides everything big; on a log axis it is a line, and in week 7 that line becomes a forecast of how often the biggest ones come.
The question, answered¶
Earthquakes and volcanoes happen along the edges of the tectonic plates: the earthquakes drew those lines before the boundary file did, and 69 percent of the catalogued volcanoes sit at subduction zones, where water carried down by a sinking plate lowers the melting point of the mantle above it. Turned on its side, the Andes show why the deep earthquakes sit inland of the trench: they happen inside the Nazca plate as it sinks beneath South America. And the big events are rare in a regular way — each step up in size several times rarer — which only a log axis shows.
Week 4 summary¶
The science¶
Where do earthquakes and volcanoes happen — and why there?
Earthquakes and volcanoes sit on plate boundaries; the deep earthquakes lie on slabs sinking beneath the trenches. Pick the plot for the question: a map for where, a side view for how deep, a log axis for how many. Big eruptions and big earthquakes are rare in the same way: each step up the scale is several times rarer, so their counts fall as a roughly straight line on a log axis.
- Log axes. When the values span factors of a thousand, plot the exponents instead and a curve becomes a line.
- Index scale. An index counts steps, not amounts: from VEI 2 up, each step is about ten times the ash and pumice thrown out, and each step of magnitude about 32 times the energy.
- Slab. A slab is a plate that has sunk into the mantle, usually at a trench; the deep earthquakes happen inside slabs.
The code¶
| Function | What it does |
|---|---|
| NumPy | |
mask & other_mask |
True only where both masks are True, the array version of and |
np.histogram(values, bins) |
the counts a histogram would draw, handed back as numbers instead of drawn |
| pandas | |
table[mask] |
keep only the rows where the mask is True |
column.count() / column.median() |
how many, and the middle value |
| Matplotlib | |
plt.figure(figsize=(w, h)) |
start a new figure, w inches wide and h tall — the first line of every plot cell |
plt.subplot(rows, cols, k) |
draw several panels in one figure; k picks which one |
plt.bar(positions, heights) |
one bar per category — the right chart when the bottom axis is a label |
plt.yscale("log") |
count the axis in factors of ten, so classes a thousandfold apart are all readable |
plt.scatter(x, y, c=values, cmap=...) |
colour every dot by a third column, with vmin/vmax fixing what the ends mean |
plt.colorbar(label=...) |
the key that says what the colours mean |
plt.hist(values, bins) |
draw how the numbers are spread; bins is either how many to use, or the edges to use |
plt.axvline(x) |
a vertical line at x, to mark a boundary on a plot |
Homework¶
Class drew the Andes from the side and found the deep earthquakes on a sinking plate. The first two parts draw two more places from the side; the third maps the volcanoes erupting now.
If you have restarted the kernel, run the setup cell at the top and then the checkpoint below.
# ── Checkpoint ── run this if you restarted the kernel or fell behind ──
trenches = boundaries[boundaries["kind"] == "trench"]
✏️ Your turn 6¶
Which way does the plate go down under Japan?
Do for Japan what Your turns 3 and 4 did for the Andes: pick the earthquakes in the box from latitude 30 to 45 and longitude 128 to 148, draw them from the side with the shallow ones blue and the deep ones red, and mark each median with a dashed line, with depth running down the page.
So: under Japan, which side of the shallow earthquakes are the deep ones on? Read it off your figure; the line at the end prints the distance between the two medians.
Use these names, because the self-check looks for them: japan, shallow_japan and
deep_japan.
japan = quakes[(quakes["latitude"] > 30) & (quakes["latitude"] < 45)
& (quakes["longitude"] > 128) & (quakes["longitude"] < 148)]
shallow_japan = japan[japan["depth"] <= 70]
deep_japan = japan[japan["depth"] > 300]
plt.figure(figsize=(9, 4.5))
plt.scatter(japan["longitude"], japan["depth"], s=5, color="0.8", label="all")
plt.scatter(shallow_japan["longitude"], shallow_japan["depth"], s=8, color="tab:blue", label="shallow, 0-70 km")
plt.scatter(deep_japan["longitude"], deep_japan["depth"], s=12, color="tab:red", label="deep, below 300 km")
plt.axvline(shallow_japan["longitude"].median(), color="tab:blue", ls="--")
plt.axvline(deep_japan["longitude"].median(), color="tab:red", ls="--")
plt.ylim(700, 0)
plt.xlabel("longitude (degrees east)")
plt.ylabel("depth (km)")
plt.title("Japan from the side: " + str(len(japan)) + " earthquakes")
plt.legend(loc="lower right")
plt.show()
print("deep minus shallow median longitude:", round(deep_japan["longitude"].median() - shallow_japan["longitude"].median(), 1), "degrees")
deep minus shallow median longitude: -5.5 degrees
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (japan["latitude"] > 30).all() and (japan["latitude"] < 45).all() and (japan["longitude"] > 128).all() and (japan["longitude"] < 148).all() and 780 < len(japan) < 900, "japan should keep every earthquake inside the box, and no others: four conditions joined with &"
assert shallow_japan["depth"].max() <= 70 and deep_japan["depth"].min() > 300 and len(deep_japan) < len(japan), "shallow_japan and deep_japan come from japan: no deeper than 70 km, and deeper than 300 km"
print(f"✓ Homework 1 — {len(shallow_japan)} shallow and {len(deep_japan)} deep earthquakes under Japan")
✓ Homework 1 — 758 shallow and 45 deep earthquakes under Japan
What the numbers mean. Under Japan the deep earthquakes sit about 5.5 degrees west of the shallow ones, where under the Andes they sat about 8 degrees east. The side view shows why the sign flips: the deep dots slope down towards the west, so the plate under Japan sinks westwards, where the Nazca plate sinks eastwards under South America. In both places the deep earthquakes lie on the side away from the ocean the plate came from. The median is a blunt tool, though: this box holds more than one stretch of trench, and a different box could move the gap.
✏️ Your turn 7¶
Do the earthquakes under the Himalaya go as deep as Japan's?
The Himalaya are where India runs into Asia: a collision of two continents, not ocean floor sinking
at a trench. Pick the earthquakes in the box from latitude 26 to 36 and longitude 78 to 97,
and draw them from the side next to Japan, on the same depth scale, with plt.subplot. The
Himalaya run from west to east, so their side view puts latitude across.
So: do the earthquakes under the Himalaya go as deep as the ones under Japan?
Use these names, because the self-check looks for them: himalaya.
himalaya = quakes[(quakes["latitude"] > 26) & (quakes["latitude"] < 36)
& (quakes["longitude"] > 78) & (quakes["longitude"] < 97)]
plt.figure(figsize=(9, 4.5))
plt.subplot(1, 2, 1)
plt.scatter(japan["longitude"], japan["depth"], s=5, color="tab:red")
plt.ylim(700, 0)
plt.xlabel("longitude (degrees east)")
plt.ylabel("depth (km)")
plt.title("Japan: " + str(len(japan)) + " earthquakes")
plt.subplot(1, 2, 2)
plt.scatter(himalaya["latitude"], himalaya["depth"], s=5, color="tab:brown")
plt.ylim(700, 0)
plt.xlabel("latitude (degrees north)")
plt.title("the Himalaya: " + str(len(himalaya)) + " earthquakes")
plt.show()
print("deepest under Japan:", round(japan["depth"].max()), "km under the Himalaya:", round(himalaya["depth"].max()), "km")
deepest under Japan: 586 km under the Himalaya: 69 km
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (himalaya["latitude"] > 26).all() and (himalaya["latitude"] < 36).all() and (himalaya["longitude"] > 78).all() and (himalaya["longitude"] < 97).all() and 60 < len(himalaya) < 110, "himalaya should keep every earthquake inside its box: four conditions joined with &"
print(f"✓ Homework 2 — the deepest earthquake under the Himalaya is {round(himalaya['depth'].max())} km down, against {round(japan['depth'].max())} km under Japan")
✓ Homework 2 — the deepest earthquake under the Himalaya is 69 km down, against 586 km under Japan
What the numbers mean. No. The deepest earthquake under the Himalaya is 69 km down, against 586 km under Japan, and none of the 83 is deeper than 70 km, while Japan's band runs down to nearly 600 km. That fits the Himalaya being a collision of continents: continental rock is light, and does not sink into the mantle the way ocean floor does, so there is no deep slab to break. The box leaves out the Hindu Kush and the Pamir, west of it, where the same collision does make earthquakes more than 200 km down, so this is one stretch, not every collision. What the earthquakes cannot show is how far north, and how deep, India's plate reaches beneath Tibet; seismologists image that from how seismic waves travel through the rock.
✏️ Your turn 8¶
Are the volcanoes erupting now on the same boundaries as all the volcanoes?
Last_Eruption_Year says when each volcano last erupted. Keep the ones that have erupted in
2000 or later as recent, and draw a world map: every volcano as a small grey triangle, the
recent ones on top as larger red triangles, and the trenches as black lines. The line at the end
prints what percent of the recent volcanoes are at subduction zones; compare it with the
69 percent of all volcanoes from Your turn 1.
Use these names, because the self-check looks for them: recent.
recent = volcanoes[volcanoes["Last_Eruption_Year"] >= 2000]
plt.figure(figsize=(9, 4.5))
plt.scatter(volcanoes["Longitude"], volcanoes["Latitude"], s=6, marker="^", color="0.7", label="all volcanoes")
plt.scatter(recent["Longitude"], recent["Latitude"], s=25, marker="^", color="crimson", label="since 2000")
plt.plot(trenches["lon"], trenches["lat"], color="black", lw=0.8, label="trench")
plt.plot(coast["lon"], coast["lat"], color="0.7", lw=0.5)
plt.xlim(-180, 180)
plt.ylim(-90, 90)
plt.gca().set_aspect("equal")
plt.xlabel("longitude (degrees east)")
plt.ylabel("latitude (degrees north)")
plt.legend(loc="lower left")
plt.title(str(len(recent)) + " of " + str(len(volcanoes)) + " volcanoes have erupted since 2000")
plt.show()
recent_subduction = recent[recent["setting"] == "Subduction zone"]
print(len(recent), "volcanoes have erupted since 2000;", len(recent_subduction), "of them,",
round(100 * len(recent_subduction) / len(recent)), "percent, are at subduction zones")
231 volcanoes have erupted since 2000; 189 of them, 82 percent, are at subduction zones
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert (recent["Last_Eruption_Year"] >= 2000).all() and 200 < len(recent) < 270, "recent should keep every volcano that last erupted in 2000 or later"
print(f"✓ Homework 3 — {len(recent)} volcanoes have erupted since 2000")
✓ Homework 3 — 231 volcanoes have erupted since 2000
What the numbers mean. 189 of the 231 volcanoes that have erupted since 2000, 82 percent, are at subduction zones, more than the 69 percent of the whole catalogue: the red triangles crowd the rim of the Pacific. The numbers leave open how much of that is the Earth and how much the record. An eruption on the sea floor, kilometres down, is rarely seen or dated, and the ridges, where most magma comes out, barely appear in either list. Watching the sea floor itself would settle it.