如何从 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

我在 SX 上发现了一个 GIS 问题，描述了如何将立交桥查询结果转换为多边形\xe2\x80\x93 ，实际上只是一个多边形列表，但我确实知道如何将它们转换为多边形。

\n\n

通过元素 ID使用Overpass 查询实际上只能获取单个对象，因此 Overpass 对于此任务来说并不是一个糟糕的 API。

\n\n

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

\n\n

对于一个简单的查询来说，结果代码出人意料地长，并将 my 中的每个坐标ndarray对转换为shapely.geometry.Point听起来效率很低，但至少这是有效的。

\n\n

import overpy\nimport shapely.geometry as geometry\nfrom shapely.ops import linemerge, unary_union, polygonize\n\nquery = """[out:json][timeout:25];\nrel(1116270);\nout body;\n>;\nout skel qt; """\napi = overpy.Overpass()\nresult = api.query(query)\n\nlss = [] #convert ways to linstrings\n\nfor ii_w,way in enumerate(result.ways):\n    ls_coords = []\n\n    for node in way.nodes:\n        ls_coords.append((node.lon,node.lat)) # create a list of node coordinates\n\n    lss.append(geometry.LineString(ls_coords)) # create a LineString from coords\n\n\nmerged = linemerge([*lss]) # merge LineStrings\nborders = unary_union(merged) # linestrings to a MultiLineString\npolygons = list(polygonize(borders))\nalaska = geometry.MultiPolygon(polygons)\n\nassert alaska.contains(geometry.Point(-147.7798220, 64.8564400))\n

\n