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.

Ground area of a 0.5 degree buffer at three airports against a 50 km circle. area of a 0.5° buffer on the ground, km²SIN, 1.4°N9,649LHR, 51.5°N6,062ANC, 61.2°N4,701true 50 km circle7,84305000 km²10k km²Projected 50 km buffers: 7,842-7,845 km² at all three spatialworkflow.io
The area a degree buffer really covers.

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.

Three map panels drawn so a kilometre is the same in both directions, with a true 50 kilometre circle and a 0.5 degree buffer around Singapore, Heathrow and Anchorage: the buffer is a near-circle at Singapore but a narrow vertical ellipse at London and Anchorage, missing airfields to the east and west. Three map panels drawn so a kilometre is the same in both directions, with a true 50 kilometre circle and a 0.5 degree buffer around Singapore, Heathrow and Anchorage: the buffer is a near-circle at Singapore but a narrow vertical ellipse at London and Anchorage, missing airfields to the east and west.
Buffers seen on the ground.

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.

Chart of the ground area of a 0.5 degree buffer against latitude, falling from about 9,700 square kilometres at the equator through the 7,854 square kilometre line of a 50 kilometre circle near 36 degrees to under 2,000 at 80 degrees. Chart of the ground area of a 0.5 degree buffer against latitude, falling from about 9,700 square kilometres at the equator through the 7,854 square kilometre line of a 50 kilometre circle near 36 degrees to under 2,000 at 80 degrees.
Area against latitude.

Buffering is planar

Buffering points by a distance in kilometres. buffer in metres, not degreespointsEPSG:4326to_crsestimate_utm_crs()buffer50_000 mto_crs(4326)if neededtrue circlelooks stretched onthe mapWorldwide points: geodesic circles from Geod.fwd instead spatialworkflow.io
Buffer in metres, not degrees.

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

Checks before buffering coordinates. before bufferingwhat CRS?EPSG:4326 means degreesheed the warninggeographic CRS: likely incorrectdistance in metresprojected CRSglobal datageodesic circlescheck one areaπ·r² on the ellipsoidLHR: a 0.5° buffer missed 20 of 78 airfields within 50 km spatialworkflow.io
Before buffering.
  • 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.

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.