在地球投影上相交两个形状良好的多边形

Question

我知道，shapely只使用笛卡尔坐标系。地球上有两个点是纬度和经度坐标。我需要创建一个半径为1km的缓冲区，并找到多边形，这个缓冲区相交的地方。在

但是建筑

buffer=点（54.4353,65.87343）。buffer（0.001）创建了一个简单的圆，但在地球上的投影中它变成了椭圆，但我需要两个半径为1km的实圆。在

我想，我需要把我的缓冲区转换成新的投影，然后和它相交，但现在不知道怎么做才正确。在

Answer 1

你要照你说的做。为此，您需要使用一个处理投影的库（pyproj是这里的选择）。在Geodesic buffering in python中有一个类似的问题

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

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

def point_buff_on_globe(lat, lon, radius):
    #First, you build the Azimuthal Equidistant Projection centered in the 
    # point given by WGS84 lat, lon coordinates
    aeqd = pyproj.Proj(proj='aeqd', ellps='WGS84', datum='WGS84',
                       lat_0=lat, lon_0=lon)
    #You then transform the coordinates of that point in that projection
    project_coords = pyproj.transform(wgs84_globe, aeqd,  lon, lat)
    # Build a shapely point with that coordinates and buffer it in the aeqd projection
    aeqd_buffer = Point(project_coords).buffer(radius) 
    # Transform back to WGS84 each coordinate of the aeqd buffer.
    # Notice the clever use of sh_transform with partial functor, this is 
    # something that I learned here in SO. A plain iteration in the coordinates
    # will do the job too.
    projected_pol = sh_transform(partial(pyproj.transform, aeqd, wgs84_globe),
                          aeqd_buffer)
    return projected_pol

函数point_buff_on_globe将给你一个纬度上的多边形，这是在以该点为中心的方位角等距投影中缓冲给定点的结果（您可以根据您的要求做的最好）。两个观察结果：

1. 我不记得radius参数的单位。我想是以米为单位的，所以如果你需要一个10公里的缓冲区，你需要通过它10e3。但是请检查一下！在
2. 注意在半径为宽或点之间相距较远时使用此选项。当点靠近投影的中心点时，投影效果很好。在