Skip to content

Instantly share code, notes, and snippets.

@jboynyc
Last active November 24, 2016 14:38
Show Gist options
  • Save jboynyc/4710bfdb61e0aeecf41a to your computer and use it in GitHub Desktop.
Save jboynyc/4710bfdb61e0aeecf41a to your computer and use it in GitHub Desktop.
Create a geometry from an OpenStreetMap element
import requests
from functools import partial
from collections import OrderedDict
from xml.dom.minidom import parseString as xml_parse
from shapely.geometry import MultiPolygon, Polygon, Point
def osm_getter(element_id, endpoint):
assert endpoint in ('node', 'way', 'relation')
OSM_API_URL = 'http://www.openstreetmap.org/api/0.6/{}/{}'
r = requests.get(OSM_API_URL.format(endpoint, element_id))
if r.status_code == 200:
return xml_parse(r.text)
elif r.status_code == 404:
raise ValueError('No such {}: {}.'.format(endpoint, element_id))
else:
raise ConnectionError('OSM API responded with {}.'.format(r.status_code))
get_node = partial(osm_getter, endpoint='node')
get_way = partial(osm_getter, endpoint='way')
get_relation = partial(osm_getter, endpoint='relation')
def node_to_point(node_id):
node = get_node(node_id).getElementsByTagName('node')[0]
return Point(map(float,
(node.attributes['lon'].value,
node.attributes['lat'].value)))
def way_to_polygon(way_id):
points = OrderedDict()
way_tree = get_way(way_id)
for nd in way_tree.getElementsByTagName('nd'):
node_id = nd.attributes['ref'].value
points[node_id] = node_to_point(node_id)
return Polygon(p.coords[0] for _, p in points.items())
def relation_to_multipolygon(relation_id):
polygons = OrderedDict()
relation_tree = get_relation(relation_id)
for member in relation_tree.getElementsByTagName('member'):
if member.attributes['type'].value == 'way':
way_id = member.attributes['ref'].value
polygons[way_id] = way_to_polygon(way_id)
return MultiPolygon([p for _, p in polygons.items()])
@jboynyc
Copy link
Author

jboynyc commented Aug 5, 2015

Is Strawberry Fields in Central Park?

>>> way_to_polygon(161202890).contains(node_to_point(3248816675))
True

Is Vienna in Australia?

>>> relation_to_multipolygon(80500).contains(node_to_point(17328659))
# sadly this doesn't work yet

NB: relation_to_multipolygon can take a long time to run, especially for relations with many members, because API calls are not made in parallel.

@jboynyc
Copy link
Author

jboynyc commented Nov 24, 2016

It's probably better to use this database.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment