如何从OpenStreetMap获取关系的几何图形？

Question

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

或者，我知道OverpassAPI可以返回几何体，所以我假设

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

可能可以工作，但是这个特定的查询是空的，并且使用OverpassAPI来处理单个关系对象，我知道它的ID感觉非常错误，不是吗

Answer 1

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

Answer 2

我发现了一个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))