2

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?

  • Extracting such amounts of data makes only sense if you run your own Overpass instance! For coastlines you should take a look at openstreetmapdata.com/data/coast – mmd Oct 24 '17 at 12:58
  • Thanks for the comment @mmd. I agree that's a lot of data if you do the whole world, but if you only care about a few countries, I think is a valid approach. That said, thanks for the reference to the coast data! – jorgeh Oct 24 '17 at 16:47
2

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.

1

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:

polygon visualisation

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.