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.

World map with red triangles for 577 volcanoes erupting since 1500 and circles sized by population for 133 places within 30 km, largest in Indonesia, Italy, Central America and the Andes. World map with red triangles for 577 volcanoes erupting since 1500 and circles sized by population for 133 places within 30 km, largest in Indonesia, Italy, Central America and the Andes.
Places within 30 km of a volcano that has erupted since 1500.

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.

Scatter of the ratio of degree-based to haversine distance against latitude: close to 1 in the tropics, rising to about 2 near 60 degrees, with one outlier near 500 across the antimeridian. Scatter of the ratio of degree-based to haversine distance against latitude: close to 1 in the tropics, rising to about 2 near 60 degrees, with one outlier near 500 across the antimeridian.
The error of measuring in degrees, by latitude.

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.

Horizontal bars of people in places within 30 km by country: Indonesia 7.63 million in 19 places, Italy 3.17, Japan 2.74, El Salvador 2.37 and Ecuador 2.29. Horizontal bars of people in places within 30 km by country: Indonesia 7.63 million in 19 places, Italy 3.17, Japan 2.74, El Salvador 2.37 and Ecuador 2.29.
People living close to recently active volcanoes, by country.

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.