Как получить геометрию отношения из OpenStreetMap? - PullRequest
Новая
       10

Как получить геометрию отношения из OpenStreetMap?

0 голосов
/ 13 марта 2020

У меня есть numpy.ndarray геокоординат, и я хочу посмотреть, какие из них ie на Аляске. Для этого я хочу получить мультиполигон штата Аляска из OpenStreetMap, а затем использовать некоторую библиотеку форм (возможно, Shapely) для запроса, какая из точек l ie внутри. Тем не менее, я застрял на шаге 1: я не могу получить геометрию мультиполигона. У меня установлен OSMPythonTools (но если есть более подходящий инструмент для работы, я с удовольствием переключусь), и я могу запросить их на Аляске следующим образом: 

from OSMPythonTools.nominatim import Nominatim
from OSMPythonTools.api import Api

nominatim = Nominatim()
api = Api()

alaska_id = nominatim.query("Alaska, United States of America").areaId()

alaska = api.query('relation/{:}'.format(alaska_id - 3600000000))

Затем я хочу получить геометрия этого объекта с использованием alaska.geometry(), но которая возвращает только 

Exception: [OSMPythonTools.Element] Cannot build geometry: geometry information not included. (way/193430587)

Это исключение возникает из-за того, что пути, составляющие внешнюю границу Аляски в alaska.__members(), не содержат геометрию, а затем API предполагает, что отношение было встречено и вызывает запутанное исключение. Я предполагаю, что мне нужно выполнить промежуточный шаг, который запрашивает все эти элементы из OSM и загружает их геометрию, как мне это сделать?

В качестве альтернативы, я знаю, что API Overpass может возвращать геометрии, поэтому я предполагаю что-то например, 

query = overpassQueryBuilder(
    area=alaska_id,
    elementType=['relation'],
    selector='"id"="1116270"',
    includeGeometry=True)

может работать, но этот конкретный запрос c пустой и использует API-интерфейс Overpass для одного объекта Relation, ID которого, как я знаю, кажется очень неправильным, не так ли?

1 Ответ

0 голосов
/ 13 марта 2020

Я нашел вопрос ГИС на SX, описывающий, как преобразовать результат запроса обхода в мультиполигон - ну, на самом деле это просто список полигонов, но я знаю, как преобразовать их в мультиполигон.

Использование запроса Overpass по идентификатору элемента Я действительно могу получить только один объект, поэтому Overpass - неплохой API для этой задачи.

Этот связанный вопрос использует overpy вместо OSMPythonTools, но OSMPythonTools настаивает на граничной рамке или области, чтобы ограничить поиск, и он также применяет некоторые маги c, чтобы построить запрос из своих параметров, вместо того, чтобы просто брать предоставленный запрос, поэтому переключение библиотек это может быть правильным решением.

Результирующий код удивительно длинен для того, что должно быть простым запросом, и преобразование каждой пары координат в моем ndarray в shapely.geometry.Point звучит неэффективно, но по крайней мере это работает. 

import overpy
import shapely.geometry as geometry
from shapely.ops import linemerge, unary_union, polygonize

query = """[out:json][timeout:25];
rel(1116270);
out body;
>;
out skel qt; """
api = overpy.Overpass()
result = api.query(query)

lss = [] #convert ways to linstrings

for ii_w,way in enumerate(result.ways):
    ls_coords = []

    for node in way.nodes:
        ls_coords.append((node.lon,node.lat)) # create a list of node coordinates

    lss.append(geometry.LineString(ls_coords)) # create a LineString from coords


merged = linemerge([*lss]) # merge LineStrings
borders = unary_union(merged) # linestrings to a MultiLineString
polygons = list(polygonize(borders))
alaska = geometry.MultiPolygon(polygons)

assert alaska.contains(geometry.Point(-147.7798220, 64.8564400))
...