我如何从OpenStreetMap获取关系的几何？

Question

我有numpy.ndarray个地理坐标，我想看看其中哪个位于阿拉斯加内部。为此，我想从OpenStreetMap获取阿拉斯加州的多面体，然后使用一些形状库（可能是Shapely）来查询其中的哪些点。但是，我陷入了第一步：无法获得多边形的几何形状。我已经安装了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查询所有这些成员并加载其几何形状，我该怎么做？

或者，我知道Overpass API可以返回几何，所以我假设类似的东西

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

可能有效，但是此特定查询为空，并且对我知道的ID非常不正确的单个Relation对象使用Overpass API，不是吗？

Answer 1

我发现了GIS question on SX describing how to convert an overpass query result into a multipolygon –好的，实际上只是一个多边形列表，但是我确实知道如何将其转换为多多边形。

使用Overpass query by element ID实际上我只能得到一个对象，因此对于这个任务来说，Overpass并不是一个不错的API。

该链接的问题使用overpy而不是OSMPythonTools，但是OSMPythonTools坚持使用边界框或区域来限制搜索，并且它还运用了一些魔术来根据其参数构建查询，而不仅仅是采用提供的查询，因此切换库可能是正确的选择。

然后生成的代码对于应该是一个简单查询的长度出乎意料的长，并且将我的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))