How do I get the geometry of a relation from OpenStreetMap?

Question

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

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)

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?

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)

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?

Answer 1

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.

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.

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.

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