4

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

6
  • pushing the question down the stack, GeoPandas uses shapely for this operation, so the question really is what is shapely doing?
    – Llaves
    Commented Jun 25, 2020 at 22:35
  • 2
    You have "World mercator projection" written beside a command that specifies your EPSG as 4326. One of these is wrong (either the comment or the code). Shapely doesn't care about your units or projection at all, so my guess is that your issue lies in your projection. It's almost always a bad idea to buffer shapes that are in degree coordinates (i.e. 4326). Convert them to a projected CRS in meters first.
    – Jon
    Commented Jun 25, 2020 at 22:55
  • @Jon - I, too, was going to comment about using a projected CRS as well until I tried his example in QGIS in EPSG:4326. I was slightly surprised that it produced the result he was looking for - a buffer a uniform 2 degrees wide running from the equator to 80 degrees north. I can imagine use cases where this is the desired result.
    – Llaves
    Commented Jun 26, 2020 at 3:13
  • Yes you right @Jon about the CRS, it's a copying error and edit the question. I also tried with projected EPSG:3395 and have the same result (see the new question editing).
    – Robin
    Commented Jun 26, 2020 at 6:54
  • @LLaves As it should. Like I said, shapely doesn't care about your units or projection so if you ask for a buffer of "2", you'll get a buffer of "2." The reason it's a bad idea (usually) to buffer in degrees is that your buffer width/length will not be the same everywhere as the length of a degree of longitude shrinks towards the poles, or put another way, 2 degrees longitude is a much smaller length in Alaska than at the equator.
    – Jon
    Commented Jun 26, 2020 at 15:00

1 Answer 1

2

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

9
  • Thank you for the answer. You're right, it's a copying error in the commentary. I edit the question. But still strange to me, with the same CRS chosen in QGIS I can't have 2 degrees at the top polygon.
    – Robin
    Commented Jun 26, 2020 at 7:18
  • Yes I agree with you but if you look at the new editing above I'm working with the same CRS EPSG:3395 for building and viewing the shapefile. And stil have trouble with the EPSG;4326 in QGIS that you didn't have (same CRS used for building and viewing).
    – Robin
    Commented Jun 26, 2020 at 9:05
  • I've no idea what you're doing then, see my edit. I'm getting exactly as expected when buffering using your code and displaying in 3395...
    – user2856
    Commented Jun 26, 2020 at 9:16
  • Is it possible to have access your version number of geopandas and shapely? I tried with another SIG (global mapper) and have the same trouble. I also tried to build the buffer of 2 different point (lat/long 0.0 and lat 80.0 long 0.0), and have the same buffer.area information...
    – Robin
    Commented Jun 26, 2020 at 14:14
  • shapely == '1.6.4', geopandas == '0.7.0'
    – user2856
    Commented Jun 26, 2020 at 21:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.