How to Find Cities Near Volcanoes in Python
On this page
Problem statement
Naples lies beside Vesuvius, Seattle within sight of Mount Rainier, and Mexico City between Popocatépetl and the volcanic fields of its own basin. A list of large cities within some distance of a volcano that has erupted in recent centuries is a common first step in hazard writing, and it is a nearest-neighbour search on the globe: two point sets in longitude and latitude, distances in kilometres. Done with plain Euclidean distance on degrees it is wrong by a factor that grows towards the poles. This page uses a ball tree with the haversine metric to find the places in Natural Earth's populated places within 30 and 100 km of a volcano that erupted since 1500.
Quick answer
Convert both sets to radians, build a BallTree with metric="haversine", and query each city's nearest volcano; distances come back in radians, so multiply by the Earth's radius:
c = eru[eru.Activity_Type == "Confirmed Eruption"]
v = vol[vol.Volcano_Number.isin(c[c.StartDateYear >= 1500].Volcano_Number)].reset_index(drop=True)
tree = BallTree(np.radians(np.c_[v.geometry.y, v.geometry.x]), metric="haversine")
dist, idx = tree.query(np.radians(np.c_[places.geometry.y, places.geometry.x]), k=1)
places["km"], places["volcano"] = dist[:, 0] * 6371.0, v.Volcano_Name.values[idx[:, 0]]
near = places[places.km <= 30].sort_values("pop_max", ascending=False)
print(near[["name", "adm0name", "pop_max", "volcano", "km"]].head(8).round(1).to_string(index=False))
print(f"volcanoes erupted since 1500: {len(v)}; places within 30 km: {len(near)}, within 100 km: {(places.km <= 100).sum()}")
name adm0name pop_max volcano km
Bandung Indonesia 2394000 Tangkuban Parahu 20.1
Naples Italy 2250000 Campi Flegrei 8.9
Haikou China 2046189 Hainan Volcanic Field 18.7
Quito Ecuador 1701000 Guagua Pichincha 11.7
San Salvador El Salvador 1433000 San Salvador 9.1
San José Costa Rica 1284000 Irazu 25.4
Guatemala City Guatemala 1024000 Pacaya 27.9
Bacolod Philippines 949354 Kanlaon 29.5
volcanoes erupted since 1500: 577; places within 30 km: 133, within 100 km: 505
577 volcanoes have had a confirmed eruption since 1500. 133 of Natural Earth's places lie within 30 km of one of them, and 505 within 100 km. The largest within 30 km is Bandung, 20.1 km from Tangkuban Parahu, then Naples, 8.9 km from the Campi Flegrei caldera it partly stands in, Haikou on Hainan, and Quito below Guagua Pichincha. San Salvador and Guatemala City complete a run of Central American capitals.
The distance, done right
1. Degrees are not kilometres
deg = np.hypot(places.geometry.x - v.geometry.x.values[idx[:, 0]], places.geometry.y - v.geometry.y.values[idx[:, 0]]) * 111.32
ratio = deg / places.km
print(f"Euclidean degrees x 111.32 against haversine, places within 100 km: median {ratio[places.km <= 100].median():.2f}, max {ratio[places.km <= 100].max():.2f}")
Euclidean degrees x 111.32 against haversine, places within 100 km: median 1.02, max 476.11
For most nearby pairs the shortcut of degrees times 111.32 km is within a few per cent, a median ratio of 1.02, because most of these places are in the tropics, where a degree of longitude is nearly as long as a degree of latitude. Towards 60° north it overstates distances by up to about twice, because a degree of longitude there is half as long. The maximum is 476 times too long: a place and a volcano either side of the 180° meridian, whose longitudes differ by almost 360 degrees, though they are close. Haversine distance has neither problem.
2. A geographic CRS warning
import warnings
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always")
places.geometry.head(1).distance(v.geometry.iloc[0])
print(w[0].message if w else "no warning")
Geometry is in a geographic CRS. Results from 'distance' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
GeoPandas warns when distance is called on longitude and latitude, because the result is in degrees. Reprojecting is the fix for a region. For the whole globe no single projection keeps distances right everywhere, which is why this page uses haversine distances on the sphere instead.
3. Check one pair
print(places[places.name == "Naples"][["name", "adm0name", "volcano", "km"]].round(1).to_string(index=False))
ves = vol[vol.Volcano_Name == "Vesuvius"].geometry.iloc[0]
nap = places[(places.name == "Naples") & (places.adm0name == "Italy")].geometry.iloc[0]
print(f"Naples to Vesuvius: {6371.0 * 2 * np.arcsin(np.sqrt(np.sin(np.radians(ves.y - nap.y) / 2) ** 2 + np.cos(np.radians(nap.y)) * np.cos(np.radians(ves.y)) * np.sin(np.radians(ves.x - nap.x) / 2) ** 2)):.1f} km")
name adm0name volcano km
Naples United States of America Chichon, El 1530.9
Naples Italy Campi Flegrei 8.9
Naples to Vesuvius: 15.6 km
Natural Earth has two Naples, and a lookup by name alone finds the one in Florida, 1,530.9 km from El Chichón in Mexico. Filtering by country gives the Italian city: its nearest volcano is Campi Flegrei at 8.9 km, while Vesuvius is 15.6 km away. Each place is matched to its single nearest volcano, so a city between two volcanoes is counted once.
People near volcanoes
The data used on this page
The Smithsonian Institution Global Volcanism Program's Volcanoes of the World database (version 5.3 is cited as doi:10.5479/si.GVP.VOTW5-2025.5.3): its Holocene volcano and eruption tables, read from the program's WFS service on 19 October 2026 and saved on this site, because the database is edited as eruptions are reported. Populated places come from Natural Earth's 1:10m populated places (public domain), whose pop_max is the population of the city's urban area, an estimate rather than a census figure.
from pathlib import Path
import zipfile
import geopandas as gpd
import numpy as np
import pandas as pd
import requests
Path("data").mkdir(exist_ok=True)
# Smithsonian Global Volcanism Program, Volcanoes of the World: Holocene volcanoes, Holocene eruptions and Pleistocene
# volcanoes as GeoJSON from its WFS service, with Natural Earth's populated places and 1:110m countries.
# Saved 19 October 2026, because the database is edited as eruptions are reported.
# Live: https://webservices.volcano.si.edu/geoserver/GVP-VOTW/ows?service=WFS&version=2.0.0&request=GetFeature
# &outputFormat=json&typeName=GVP-VOTW:Smithsonian_VOTW_Holocene_Volcanoes (and ..._Holocene_Eruptions)
ARCHIVE = "gvp-votw-2026-10-19.zip"
if not Path(f"data/{ARCHIVE}").exists():
Path(f"data/{ARCHIVE}").write_bytes(requests.get(f"https://data.spatialworkflow.io/{ARCHIVE}", timeout=600).content)
Z = zipfile.ZipFile(f"data/{ARCHIVE}")
vol = gpd.read_file(Z.open("Smithsonian_VOTW_Holocene_Volcanoes.json"))
eru = gpd.read_file(Z.open("Smithsonian_VOTW_Holocene_Eruptions.json"))
from sklearn.neighbors import BallTree
Path("data/ne_10m_populated_places_simple.zip").write_bytes(Z.read("ne_10m_populated_places_simple.zip"))
places = gpd.read_file("data/ne_10m_populated_places_simple.zip")
print(f"Volcanoes of the World: {len(vol):,} Holocene volcanoes, {len(eru):,} eruptions; Natural Earth places: {len(places):,}")
Volcanoes of the World: 1,214 Holocene volcanoes, 11,089 eruptions; Natural Earth places: 7,342
Example 1: by country
by = near.groupby("adm0name").agg(places=("name", "size"), people_m=("pop_max", lambda s: s.sum() / 1e6)).sort_values("people_m", ascending=False)
print(by.head(6).round(2).to_string())
places people_m
adm0name
Indonesia 19 7.63
Italy 3 3.17
Japan 11 2.74
El Salvador 10 2.37
Ecuador 7 2.29
China 1 2.05
Indonesia has 19 places within 30 km of a recently erupting volcano, with 7.63 million people between them. Italy's three are dominated by Naples, and El Salvador has ten in a small country. The populations are Natural Earth's metropolitan estimates, so the sums are rough, and they leave out everyone who doesn't live in a listed place.
Example 2: the volcanoes with the most people nearby
print((near.groupby("volcano").pop_max.sum() / 1e6).nlargest(6).round(2).to_string())
volcano
Tangkuban Parahu 2.39
Campi Flegrei 2.25
Hainan Volcanic Field 2.05
Guagua Pichincha 1.88
San Salvador 1.56
Irazu 1.48
Tangkuban Parahu, Campi Flegrei and the Hainan volcanic field each have more than two million people in places within 30 km, then Guagua Pichincha, San Salvador and Irazú. Distance is only a first screen: hazard depends on how a volcano erupts, the wind and the terrain. Lava rarely travels 30 km, while ash can fall hundreds of kilometres away.
Internal links
- Nearest-neighbour joins — the method in general
- Volcano types and settings — which volcanoes these are
- Seismic hazard and exposed population — the same exposure question for earthquakes
- Countries with erupting volcanoes — counting by country