[GIS] How buffering works in GeoPandas

buffergeopandasparametersshapely

I try to add a buffer on a set of lines via the buffer attribute of geopandas.

But by testing on a simple case (see code below) I realize that the size of the buffer is not kept according to the points latitude.

The buffer build at point A is well 2° but at point B is only 0.35° (checked via QGIS).

Could someone please explain to me how this buffer works under geopandas?

And if you have a solution to generate a buffer that is robust to the change of latitude I'm interested!

import geopandas as gpd

from shapely.geometry import LineString, Point

# 2 points creation
A = Point(0.0, 0.0)
B = Point(0.0, 80.0)

# Linestring creation
lines_geom = LineString([A, B])
crs = {'init': 'epsg:4326'} # WGS84 CRS
lines = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[lines_geom])       

# Buffer generation
lines_buffer_1deg = lines.buffer(1.0, resolution=1)
print(lines_buffer_1deg)

lines_buffer_1deg.to_file("test.shp")

As proposed in the comment I also tried with a projected CRS, EPSG:3395 world mercator with the same result (see code below and caption from QGIS).

import geopandas as gpd

from shapely.geometry import mapping, Polygon, LineString, Point

# 2 points creation
A = Point(0.0,0.0)
B = Point(0.0,9000000.0)

# Linestring creation
lines_geom = LineString([A,B])
crs = {'init': 'epsg:3395'} 
lines = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[lines_geom])       

# Buffer generation
lines_buffer_100km = lines.buffer(100000.0,resolution=1)
print(lines_buffer_100km)

lines_buffer_100km.to_file("test.shp")

EPSG:3395 around 0° latitude:
EPSG:3395 at around 0° latitude

EPSG:3395 around 62° latitude:
EPSG:3395 around 62° latitude

Best Answer

You are generating your buffer in an unprojected (geographic) CRS EPSG:4326, not "world mercator" as commented in your script.

Your buffers are generated correctly (if it can be said that using degrees as a measurement is "correct"...), with 2 degree width at all latitudes. However you are viewing them in EPSG:3857 a projected coordinate system known as "Web Mercator" or "Pseudo-Mercator" so they appear squashed at higher latitudes.

Basically, you are buffering in one coordinate system, if you display and measure in that coordinate system, you will get the expected buffer width. If you display/measure in a different coordinate system, you'll get a different answer.

Below are a couple of screenshots of the same shapefile (EPSG:4326 & 3395) output by your code, but using different CRS in QGIS to visualise/measure.

EPSG:4326

enter image description here

EPSG:3857

enter image description here

EPSG:3395

enter image description here

Note that using EPSG:3857 to do any sort of measurement is a bad idea. The Web Mercator projection is designed for one thing only, chopping up the globe into 256x256 tiles for web mapping. From an Esri article:

The modified Mercator projection used by Google, Bing, and ArcGIS Online is not designed to minimize distortion at all. Instead, it was engineered for convenience in working with cached map tiles. This projection fit the entire globe (well, most of the latitudes anyway) into a square area that could be covered by 256 x 256 pixel tiles. The projection sacrifices some accuracy because it is based on a perfect sphere (the earth is better approximated by a spheroid), but the biggest problem is the heavy vertical and horizontal stretching at extreme latitudes. This is evident from the enormous dimensions of Greenland and Antarctica relative to land masses closer to the Equator:

enter image description here