Python中的测地线缓冲

Question

假设已经有了一个由Shapely创建的MultiPolygon地理区域，现在需要找到代表海岸线周围12海里缓冲区的（多边形）区域。

使用Shapely的buffer方法无法实现，因为它使用欧几里得计算。

请问有谁能告诉我如何在Python中计算大圆距离缓冲区？

- ARF

---

我建议在 gis.SE 上移动 / 提问这个问题。 - Mike T

@MikeT 我同意将问题转移到gis.SE。请继续。谢谢。 - ARF

由于SO管理员负载过重，自己复制/粘贴/删除此问题可能更快。 - Mike T

@VikashPandey，请停止在与问题无关的地方垃圾邮件随机的GIS相关标签。 - sco1

Answer 1

这不是一个合适的问题，因为shapely在其文档中明确指出该库仅用于平面计算。然而，为了回答您的问题，您应该指定您用于多边形的坐标系统。假设您正在使用WGS84投影（纬度，经度），这是我在另一个SO问题中找到的解决方法（fix-up-shapely-polygon-object-when-discontinuous-after-map-projection）。您需要使用pyproj库。

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函数的功能如下：首先，在多边形重心处建立一个等角方位投影。然后，将多边形的坐标系更改为该投影的坐标系。之后，在该位置构建缓冲区，然后将缓冲多边形的坐标系更改为WGS84坐标系。

需要特别注意以下几点：

• 您需要找出如何将所需距离转换为aeqd投影中使用的距离。
• 小心不要将包括极点在内的缓冲区。
• 使用多边形重心作为投影中心点应该足够准确，但如果您有特定的精度要求，则不应使用此解决方案，或者至少对您正在使用的典型多边形进行误差表征。

- eguaio