python中的测地线缓冲

Question

鉴于陆地多边形为Shapely MultiPolygon,我想找到代表海岸线周围12海里缓冲区的(多)多边形.

使用Shapely buffer方法不起作用,因为它使用欧几里德计算.

有人能告诉我如何计算python中的测地缓冲区吗？

Answer 1

这不是一个形状问题,因为在其文档中明确地说明该库仅用于平面计算.但是,为了回答您的问题,您应该指定用于多面体的坐标系.假设您正在使用WGS84投影(lat,lon),这是我在另一个SO问题中找到的一个配方(固定形状多边形对象,当不连续地在地图投影后).你需要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投影中使用的距离.
• 小心不要包括两极的缓冲(参见提到的SO问题).
• 我们使用多边形的质心来使投影居中这一事实应该保证答案是足够好的,但是如果你有特定的精度要求你不应该使用这个解决方案,或者至少对典型的多边形进行误差表征您正在使用.