Пересечение двух фигурных многоугольников в проекции Земли.

насколько я знаю, shapely используют только декартову систему координат. У меня есть две точки на земле с координатами широты и долготы. Мне нужно создать буфер с радиусом 1 км вокруг этих двух точек и найти многоугольник, где эти буферы пересекаются.

Но конструкция

buffer = Point(54.4353,65.87343).buffer(0.001) создать простой круг, но в проекции на Землю он становится эллипсом, а мне нужно два реальных круга с радиусом 1 км.

Я думаю, мне нужно преобразовать мои буферы в новую проекцию, а затем пересечь ее, но не знаю, как правильно это сделать.


python geospatial shapely pyproj geos
person Михаил Агарков    schedule 16.06.2017    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Вы должны делать то, что вы говорите. Для этого вам нужно будет использовать библиотеку, которая обрабатывает проекции (здесь можно выбрать pyproj). Аналогичный вопрос есть в геодезическая буферизация в 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. Остерегайтесь использовать это с широким радиусом или точками, которые находятся далеко друг от друга. Проекции работают хорошо, когда точки находятся рядом с точкой, в которой вы центрируете проекцию.
person eguaio    schedule 04.07.2017