Skip to content

Instantly share code, notes, and snippets.

Last active June 20, 2017 08:16
Show Gist options
  • Save robintw/9366322 to your computer and use it in GitHub Desktop.
Save robintw/9366322 to your computer and use it in GitHub Desktop.
Code to convert OSM XML files with waypoints defining polygons into a shapefile with those polygons in it
import xml.etree.cElementTree as ET
from shapely.geometry import mapping, Polygon
import fiona
import os
input_filename = "villages.osm"
output_filename = "villages.shp"
# Parse the XML from the OSM file
tree = ET.ElementTree(file=input_filename)
# Get the root of the XML (the <osm> node)
r = tree.getroot()
# Get all of the individual points from the file into a dictionary,
# with the keys being the ids of the nodes
nodes = {}
# Put all nodes in a dictionary
for n in r.findall('node'):
nodes[n.attrib['id']] = (float(n.attrib['lon']), float(n.attrib['lat']))
# Create a dictionary to hold the polygons we create
polygons = {}
# For each 'way' in the file
for way in r.findall("way"):
coords = []
# Look at all of the children of the <way> node
for c in way.getchildren():
if c.tag == "nd":
# If it's a <nd> tag then it refers to a node, so get the lat-lon of that node
# using the dictionary we created earlier and append it to the co-ordinate list
ll = nodes[c.attrib['ref']]
if c.tag == "tag":
# If it's a <tag> tag (confusing, huh?) then it'll have some useful information in it for us
# Get out the information we want as an attribute (in this case village name - which could be
# in any of a number of fields, so we check)
if c.attrib['v'] not in ("residential", 'village'):
village_name = c.attrib['v']
# Take the list of co-ordinates and convert to a Shapely polygon
polygons[village_name] = Polygon(coords)
# Set up the schema for the shapefile
# In this case we have two columns: id and name
schema = {
'geometry': 'Polygon',
'properties': {'id': 'int', 'name': 'str'},
# Remove the shapefile if it exists already
i = 0
# Write the shapefile
with, 'w', 'ESRI Shapefile', schema) as c:
# For every polygon we stored earlier
for name, p in polygons.iteritems():
i += 1
print "Writing: %s" % name
# Write it to the shapefile
'geometry': mapping(p),
'properties': {'id': i, 'name': name},
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment