如何从 OpenStreetMap 获取关系的几何形状？答案

【问题标题】：How do I get the geometry of a relation from OpenStreetMap?如何从 OpenStreetMap 获取关系的几何形状？
【发布时间】：2020-03-13 12:06:56
【问题描述】：

我有一个numpy.ndarray 的地理坐标，我想看看其中哪些位于阿拉斯加境内。为此，我想从 OpenStreetMap 获取阿拉斯加州的多面体，然后使用一些形状库（可能是 Shapely）来查询其中的点。但是，我被困在第 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 查询所有这些成员并加载它们的几何图形，我该怎么做？

另外，我知道 Overpass API 可以返回几何图形，所以我假设类似

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

可能有效，但是这个特定的查询是空的，并且使用 Overpass API 处理一个我知道其 ID 的 Relation 对象感觉非常错误，不是吗？

【问题讨论】：

标签： python openstreetmap

【解决方案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))

【讨论】：

【解决方案2】：

我认为 OSMPythonTools 在幕后施展魔法的说法是错误的。如果您使用overpassQueryBuilder，OSMPythonTools 会为您组装查询，但您也可以提交一个字符串查询：overpass.query(...)。因此 OSMPythonTools 应该是适合此目的的工具。我们可能会就此询问作者。

【讨论】：