[GIS] Python shapely intersection with buffer in meter

python 3shapely

Having a line and a polygon, I use shapely to determine the intersection between the two with a defined buffer. Both the line and the polygon are in WGS84 coordinates, while the buffer is in meter.

from shapely.geometry import LineString, Polygon
polygon = Polygon([
    (51.9981816, 4.3677729),
    (51.9983431, 4.3683897),
    (51.9980583, 4.3685846),
    (51.9977777, 4.3687773),
    (51.9976163, 4.3681605),
    (51.9978555, 4.3679964),
    (51.9981816, 4.3677729)
])
line = LineString([
    (51.9976836984107, 4.368050461186726),
    (51.99757198557364, 4.368489603711765),
    (51.9973971301288, 4.3684738261236475)
])

Finding the intersection is easy:

line.intersection(polygon)

Finding the intersection with a buffer in degrees (because coordinates of the objects are in degrees):

line.intersection(polygon.buffer(0.2))

But how to find the intersection with a buffer in meters? I understand that Shapely is not aware of any coordinate reference system; also no area metric CRS is known (OceanoGraphy is globally).

Best Answer

Shapely knows nothing about your units. Therefore, when you call polygon.buffer(0.2), shapely will buffer the coordinates in polygon by 0.2 units. In your case, these are not meters, but degrees.

If you want to buffer in meters, you first need to reproject your Polygon and Linestring into a CRS that uses meters. I would use geopandas for this, as it has a convenient .to_crs() method that will do the reprojection quickly. This also requires that you know the CRS of the coordinates you're providing (likely they are epsg:4326). You could also use pyproj, which might be easier if you're not dealing with layers with many polygons:

from pyproj import Proj, transform

inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3857')
x1,y1 = -11705274.6374,4826473.6922
x1p,y1p = transform(inProj,outProj,x1,y1)

I will leave the looping up to you (you can feed in lists rather than single points, I believe) for all of your coordinates.

Important: Choosing the proper CRS is not trivial, and depends on your application and precision needs. I have shown EPSG:3857, which is commonly used by services like Google/Bing maps, etc. However, all projections will have some distortions that vary (usually) with distance from the equator (or the projection's central latitude). One commonly-used option is to determine which UTM Zone your region of interest lies within, then choose the appropriate EPSG.

(EPSG is just a numeric code that maps to a projection.)

Related Question