Fixing Buffers Drawn in Degrees Instead of Kilometres
On this page
Problem statement
You buffered points by 0.5 to get "about 50 km" zones, or by 0.01 for "about a kilometre", and the results are wrong in ways that depend on where the points are: zones that are ellipses on the ground, too narrow east-west, counts of nearby features that change with latitude, and areas that make no sense. In a GeoDataFrame in EPSG:4326, buffer(0.5) means half a degree, and a degree of longitude shrinks from 111 km at the equator to 56 km at 60°N. GeoPandas warns about it, and the warning is easy to scroll past.
This page buffers real airports from OurAirports at three latitudes by half a degree, measures what the buffers really cover on the ground, counts the airfields each one catches against a true 50 km circle, and fixes the buffers by projecting to a local CRS or drawing geodesic circles.
Quick answer
Buffer in a projected CRS in metres, local to the points, and project back:
pts = gpd.GeoDataFrame(H, geometry=gpd.points_from_xy(H.longitude_deg, H.latitude_deg), crs=4326)
wrong = pts.buffer(0.5) # half a degree: warns, and is not 50 km
right = gpd.GeoSeries([gpd.GeoSeries([g], crs=4326).to_crs(gpd.GeoSeries([g], crs=4326).estimate_utm_crs()).buffer(50_000).to_crs(4326).iloc[0]
for g in pts.geometry], crs=4326)
for code, w, r in zip(H.index, wrong, right):
print(f"{code}: 0.5° buffer area {abs(geod.geometry_area_perimeter(w)[0]) / 1e6:7,.0f} km²; 50 km buffer {abs(geod.geometry_area_perimeter(r)[0]) / 1e6:7,.0f} km² (π·50² = {np.pi * 2500:,.0f})")
SIN: 0.5° buffer area 9,649 km²; 50 km buffer 7,845 km² (π·50² = 7,854)
LHR: 0.5° buffer area 6,062 km²; 50 km buffer 7,842 km² (π·50² = 7,854)
ANC: 0.5° buffer area 4,701 km²; 50 km buffer 7,843 km² (π·50² = 7,854)
Projected buffers of 50 km have the area of a 50 km circle, about 7,843 km², at all three airports. The 0.5° buffers do not: 9,649 km² at Singapore, where half a degree is about 55 km in every direction, 6,062 km² at London and only 4,701 km² at Anchorage, 40% less than the circle meant.
What a degree buffer reaches, airfields caught, geodesic circles
1. What a degree buffer reaches
Measure how far the edge of a 0.5° buffer is from the centre to the north and to the east, at each airport:
for code, r in H.iterrows():
n = geod.inv(r.longitude_deg, r.latitude_deg, r.longitude_deg, r.latitude_deg + 0.5)[2] / 1000
e = geod.inv(r.longitude_deg, r.latitude_deg, r.longitude_deg + 0.5, r.latitude_deg)[2] / 1000
print(f"{code} at {r.latitude_deg:5.1f}°: reaches {n:.1f} km north, {e:.1f} km east (ratio {e / n:.2f})")
SIN at 1.4°: reaches 55.3 km north, 55.6 km east (ratio 1.01)
LHR at 51.5°: reaches 55.6 km north, 34.7 km east (ratio 0.62)
ANC at 61.2°: reaches 55.7 km north, 26.9 km east (ratio 0.48)
Half a degree reaches about 55.5 km north everywhere, but east only 55.6 km at Singapore, 34.7 km at London and 26.9 km at Anchorage: the buffer is an ellipse on the ground, less than half as wide as it is tall in Alaska.
2. How many airfields each catches
Count OurAirports places inside the degree buffer and inside the true 50 km circle, measured with geodesic distance:
for (code, r), w in zip(H.iterrows(), wrong):
d = geod.inv(np.full(len(A), r.longitude_deg), np.full(len(A), r.latitude_deg), A.longitude_deg.values, A.latitude_deg.values)[2] / 1000
in_true = d <= 50
in_deg = A_pts.within(w).values
print(f"{code}: true 50 km {in_true.sum():>4}; 0.5° buffer {in_deg.sum():>4}; missed {(in_true & ~in_deg).sum():>3}; extra beyond 50 km {(in_deg & ~in_true).sum():>3}")
SIN: true 50 km 17; 0.5° buffer 17; missed 0; extra beyond 50 km 0
LHR: true 50 km 78; 0.5° buffer 63; missed 20; extra beyond 50 km 5
ANC: true 50 km 95; 0.5° buffer 84; missed 16; extra beyond 50 km 5
At Singapore the degree buffer happens to catch the same 17 airfields as a true 50 km circle. At London it misses 20 of the 78 within 50 km and includes 5 beyond it; at Anchorage it misses 16 of 95. The error depends on latitude and on where the neighbours happen to lie.
3. A geodesic circle without projecting
For points spread over the world, draw each circle from destination points at every bearing; no CRS choice is needed:
def geodesic_circle(lon, lat, km, n=180):
b = np.linspace(0, 360, n, endpoint=False)
lo, la, _ = geod.fwd(np.full(n, lon), np.full(n, lat), b, np.full(n, km * 1000))
return Polygon(zip(lo, la))
circles = gpd.GeoSeries([geodesic_circle(r.longitude_deg, r.latitude_deg, 50) for _, r in H.iterrows()], crs=4326)
for code, c, r in zip(H.index, circles, right):
print(f"{code}: geodesic circle {abs(geod.geometry_area_perimeter(c)[0]) / 1e6:,.0f} km²; area not shared with the UTM buffer {abs(geod.geometry_area_perimeter(c.symmetric_difference(r))[0]) / 1e6:.1f} km²")
SIN: geodesic circle 7,852 km²; area not shared with the UTM buffer 8.1 km²
LHR: geodesic circle 7,852 km²; area not shared with the UTM buffer 11.0 km²
ANC: geodesic circle 7,852 km²; area not shared with the UTM buffer 10.1 km²
A geodesic circle built from destination points covers 7,852 km², and it differs from the UTM buffer by only 8 to 11 km² of area, the slivers between a 180-sided polygon and a buffered curve and the small scale error of UTM. Either is a correct 50 km zone.
The warning you skipped
The data used on this page
Three large airports at different latitudes, Singapore Changi, London Heathrow and Anchorage, and every open place in OurAirports, a public-domain database, for the counts. Areas and distances are geodesic, from pyproj.Geod on WGS84.
import warnings
from pathlib import Path
import geopandas as gpd
import numpy as np
import pandas as pd
import requests
from pyproj import Geod
from shapely.geometry import Polygon
warnings.filterwarnings("ignore", message="Geometry is in a geographic CRS")
Path("data").mkdir(exist_ok=True)
p = Path("data/airports.csv")
if not p.exists():
r = requests.get("https://davidmegginson.github.io/ourairports-data/airports.csv", timeout=300)
r.raise_for_status()
p.write_bytes(r.content)
A = pd.read_csv(p, keep_default_na=False, na_values=[""])
A = A[A.type != "closed"].reset_index(drop=True)
A_pts = gpd.GeoSeries(gpd.points_from_xy(A.longitude_deg, A.latitude_deg), crs=4326)
H = A.dropna(subset=["iata_code"]).set_index("iata_code").loc[["SIN", "LHR", "ANC"]]
geod = Geod(ellps="WGS84")
print(H[["municipality", "latitude_deg", "longitude_deg"]].round(3).to_string())
municipality latitude_deg longitude_deg
iata_code
SIN Singapore 1.350 103.994
LHR London 51.471 -0.460
ANC Anchorage 61.179 -149.993
The three airports sit at 1.4°N, 51.5°N and 61.2°N.
Example 1: the warning you skipped
GeoPandas warns when buffering geographic coordinates. Capture it to see the message:
with warnings.catch_warnings(record=True) as w:
warnings.simplefilter("always")
gpd.GeoSeries(gpd.points_from_xy([0], [51.5]), crs=4326).buffer(0.5)
print(w[0].message if w else "no warning")
Geometry is in a geographic CRS. Results from 'buffer' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.
GeoPandas warns that results from buffer are likely incorrect in a geographic CRS. Python shows each warning only once per line of code by default, which is part of why it is easy to miss.
Buffering is planar
Buffering is a planar operation: shapely draws a circle of the given radius in whatever units the coordinates are in. In EPSG:4326 those units are degrees of latitude and longitude, which are not the same length on the ground. A degree of latitude is always about 111 km, but a degree of longitude is 111 km times the cosine of the latitude, so a 0.5° buffer reaches 55 km north and south everywhere, but only 55 km east at the equator, 35 km at London and 26 km at Anchorage. Its area on the ground shrinks by the same factor, and anything counted inside it changes with latitude for no reason in the data.
There are two fixes. For points in one region, project to a CRS in metres that suits the region, such as the local UTM zone, buffer in metres and project back if needed. For points spread across the globe, either project each point to its own zone or draw geodesic circles from destination points, which need no projection at all. Both give true circles on the ground, which look like ellipses stretched east-west on a latitude-longitude map. That stretch is the map, not the circle.
Approximate fixes, large radii and 180 degrees
- Approximate fix: dividing the radius by cos(latitude) for longitude only, a scaled ellipse, works for small radii away from the poles; geodesic circles are simpler and exact.
- Large radii: a projected buffer of several hundred kilometres accumulates the projection's distortion; use geodesic circles.
- Crossing 180°: geodesic circles near the antimeridian need splitting before becoming polygons.
- Unit sanity check: a buffer distance below 1 on projected data, or above 1 on geographic data, is almost always the wrong unit.
Internal links
- How to find a destination point from distance and bearing in Python — the geodesic circle in detail
- How to buffer in Python with GeoPandas — buffering in a projected CRS
- How to find all points within a radius in Python — counting without buffers
FAQ
Why is my buffer an ellipse?
On a latitude-longitude map a true circle looks like an ellipse; but a buffer drawn in degrees is a circle on the map and an ellipse on the ground, too narrow east-west.
How do I buffer points by kilometres in GeoPandas?
Project to a CRS in metres, for example with gdf.to_crs(gdf.estimate_utm_crs()), buffer by the distance in metres, and project back if needed.
How many degrees is 50 km?
About 0.45° of latitude everywhere, but 0.45° divided by cos(latitude) of longitude: about 0.72° at London.