Python中的测地线缓冲

问题描述

将陆地多边形指定为形状 MultiPolygon，我想找到(Multi-代表例如的多边形海岸线周围有12海里的缓冲区.

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.

使用Shapely buffer方法不起作用，因为它使用了欧式计算.

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

有人可以告诉我如何在python中计算测地线缓冲区吗?

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

推荐答案

这不是一个形状问题，因为形状明确地在其文档中告知该库仅用于平面计算.但是，为了回答您的问题，您应该指定用于多面体的坐标系.假设您使用的是WGS84投影(纬度，经度)，这是我在另一个SO问题().您将需要pyproj库.

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函数执行以下操作.首先，建立一个以多边形质心为中心的方位角等距投影.然后，将多边形的坐标系更改为该投影.之后，在此处建立缓冲区，然后将缓冲的多边形的坐标系更改为WGS84坐标系.

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.

需要一些特殊的照顾:

• 您将需要了解如何将所需距离转换为等距投影中使用的距离.
• 请注意不要缓冲两极(请参见所提到的SO问题).
• 我们使用多边形的质心来使投影居中这一事实应保证答案足够好，但是如果您有特定的精度要求，则不应该使用此解决方案，或者至少要对误差进行表征您正在使用的典型多边形.

这篇关于Python中的测地线缓冲的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！