How to get a Multipolygon object from Overpass QL?

Question

I'm new to OSM/Overpass QL and there's a simple problem that I can't solve. I've been able to figure out the a query that will return the maritime boundary of a country. For example:

(rel[boundary="administrative"][name="Puerto Rico"][admin_level="4"];>;);way(r);(._;>;);out;

This returns the boundary as a set of ways. Using overpy I've been able to get the returned ways and display them using shapely/matplotlib, and I'm confident they have the data I need.

What I can't figure out is how can I reliably convert the returned set of ways to a MultiPolygon object, ideally a shapely.geometry.MyltiPolygon.

The above query is only an example that produces a fairly simple set of ways, but there are other queries that return much more complex results, such as

(rel[boundary="administrative"]["name:en"="French Polynesia"][admin_level="3"];>;);way(r);(._;>;);out;

I understand that Overpass doesn't handle (multi)polygons, but I assume there must be a simple way to achieve this using python/GeoPandas.shapely. I've tried some manual approaches, but they are error-prone and I'd be extremely surprised if there are no implemented librariessolutions to achieve this. Any ideas?

Answer 1

I was able to solve my problem using the function ways2poly function I wrote in this gist.

I also added code to run the Overpass queries (using overpy), in case anyone needs it.

Answer 2

I convert OSM ways to polygons using Shapely's linemerge, unary_union, and polygonize.

Like so:

import overpass
import shapely.geometry as geometry
from shapely.ops import linemerge, unary_union, polygonize

api = overpy.Overpass()

query = """ YOUR QUERY HERE """
response = api.query(query)

lss = [] #convert ways to linstrings

for ii_w,way in enumerate(response.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))

Running your query and plotting in matplotlib: