如何从 OpenStreetMap 获取关系的几何? [英] How do I get the geometry of a relation from OpenStreetMap?

问题描述

我有一个 numpy.ndarray 地理坐标，我想看看其中哪些位于阿拉斯加境内.为此，我想从 OpenStreetMap 获取阿拉斯加州的多边形，然后使用一些形状库(可能是 Shapely)来查询哪些点位于里面.但是，我被困在第 1 步:我无法获得多边形的几何形状.我安装了 OSMPythonTools(但如果有更好的工具，我很乐意切换)，我可以像这样查询阿拉斯加

I have a numpy.ndarray of geocoordinates and I want to see which of them lie inside Alaska. For that, I want to get the multipolygon of the state of Alaska from OpenStreetMap and then use some shape library (probably Shapely) to query which of the points lie inside. However, I am stuck at step 1: I cannot get the geometry of the multipolygon. I have the OSMPythonTools installed (but if there is a better tool for the job, I'm happy to switch) and I can query them for Alaska like this

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() 获取这个对象的几何形状，但它只返回

I then want to get the geometry of this object using alaska.geometry(), but that returns only

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

引发此异常是因为alaska.__members() 中构成阿拉斯加外边界的方式不包含几何图形，然后 API 假定已遇到关系并引发令人困惑的异常.我假设我需要运行一个中间步骤，从 OSM 查询所有这些成员并加载它们的几何图形，我该怎么做?

This exception is raised because the ways constituting the outer boundary of Alaska in alaska.__members() do not contain geometry, and then the API assumes that a relation has been encountered and raises a confusing exception. I assume that I need to run an intermediate step that queries all these members from OSM and loads their geometry, how do I do that?

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

Alternatively, I know that the Overpass API can return geometries, so I assume something like

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

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

might work, but this specific query is empty and using the Overpass API for a single Relation object whose ID I know feels very wrong, is it not?

推荐答案

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

I found a GIS question on SX describing how to convert an overpass query result into a multipolygon – well, actually just a list of polygons, but I do know how to convert those to a multipolygon.

使用 按元素 ID 进行 Overpass 查询 我实际上可以得到一个单个对象，所以 Overpass 不是这个任务的糟糕 API.

Using the Overpass query by element ID I can actually get just a single object, so Overpass is not a bad API for this task.

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

That linked question uses overpy instead of OSMPythonTools, but OSMPythonTools insists on a boundary box or area to limit the search, and it also applies some magic to build a query from its parameters instead of just taking the provided query, so switching libraries may be the right thing to do.

结果代码对于应该是一个简单查询的内容来说非常长，并且将我的 ndarray 中的每个坐标对转换为 shapely.geometry.Point 听起来效率很低，但至少这是有效的.

The resulting code is then surprisingly long for what should be a simple query, and converting every coordinate pair in my ndarray into a shapely.geometry.Point sounds inefficient, but at least this works.

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))

这篇关于如何从 OpenStreetMap 获取关系的几何?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！