CODEWITHSUNDEEP
pythongisshapelycartopy
ARF
ARF

Reputation: 7684

Geodesic buffering in python

Given land polygons as a Shapely MultiPolygon, I want to find the (Multi-)Polygon that represents the e.g. 12 nautical mile buffer around the coastlines.

Using the Shapely buffer method does not work since it uses euclidean calculations.

Can somebody tell me how to calculate geodesic buffers in python?

Upvotes: 5

Views: 3852

Answers (1)

eguaio
eguaio

Reputation: 3954

This is not a shapely problem, since shapely explicitly tells in its documentation that the library is for planar computation only. Nevertheless, in order to answer your question, you should specify the coordinate systems you are using for your multipolygons. Assuming you are using WGS84 projection (lat,lon), this is a recipe I found in another SO question (fix-up-shapely-polygon-object-when-discontinuous-after-map-projection). You will need pyproj library.

import pyproj
from shapely.geometry import MultiPolygon, Polygon
from shapely.ops import transform as sh_transform
from functools import partial

wgs84_globe = pyproj.Proj(proj='latlong', ellps='WGS84')

def pol_buff_on_globe(pol, radius):
    _lon, _lat = pol.centroid.coords[0]
    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84',
                       lat_0=_lat, lon_0=_lon)
    project_pol = sh_transform(partial(pyproj.transform, wgs84_globe, aeqd), pol)
    return sh_transform( partial(pyproj.transform, aeqd, wgs84_globe),
                          project_pol.buffer(radius))

def multipol_buff_on_globe(multipol, radius):
    return MultiPolygon([pol_buff_on_globe(g, radius) for g in multipol])

pol_buff_on_globe function does the following. First, build an azimuthal equidistant projection centered in the polygon centroid. Then, change the coordinate system of the polygon to that projection. After that, builds the buffer there, and then change the coordinate system of the buffered polygon to WGS84 coordinate system.

Some special care is needed:

  • You will need to find out how to translate the distance you want to the distance used in aeqd projection.
  • Be careful of not buffering including the poles (see the mentioned SO question).
  • The fact that we are using the centroid of the polygon to center the projection should guaranty the answer is good enough, but if you have specif precision requirements you should NOT USE this solution, or at least make a characterization of the error for the typical polygon you are using.

Upvotes: 7

Related Questions