You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

434 lines
15 KiB

This script is designed to act as assistance in converting shapefiles
to OpenStreetMap data. This file is optimized and tested with MassGIS
shapefiles, converted to EPSG:4326 before being passed to the script.
You can perform this conversion with
ogr2ogr -t_srs EPSG:4326 new_file.shp old_file.shp
It is expected that you will modify the fixed_tags, tag_mapping, and
boring_tags attributes of this script before running. You should read,
or at least skim, the code up until it says:
to accomodate your own data.
Modified to support a maximum waylength and multiple outer ways in a
multipolygon relationship
adapted by devemux86
__author__ = "Christopher Schmidt <>"
__version__ = "$Id:,v 1.3 2009/03/30 20:54:13 tbook Exp $"
gdal_install = """
Installing GDAL depends on your platform. Information is available at:
For Debian-based systems:
apt-get install python-gdal
will usually suffice.
# These tags are attached to all exterior ways. You can put any key/value pairs
# in this dictionary.
fixed_tags = {
'natural': 'nosea',
'layer': '-5'
# Here are a number of functions: These functions define tag mappings. The API
# For these functions is that they are passed the attributes from a feature,
# and they return a list of two-tuples which match to key/value pairs.
# The following Ftypes are not imported
# 46100 - Submerged Stream
ignoreField = "tile_x"
ignoreValues = []
# The following Ftypes are imported, but don't have any OSM tags attached:
# 36100 - Playa - An area from which water evaporates
def ftype(data):
"""Type of body - From NHD Ftype"""
natural = {
'LakePond': 'water',
'SwampMarsh': 'wetland'
if 'ftype' in data:
if data['ftype'] in natural:
return [('natural', natural[data['ftype']])]
return None
def fcode(data):
"""For features where the ftype is not specific enough"""
landuse = {
'43600': 'reservoir', #Reservoir
'43601': 'reservoir', #Reservoir
'43602': 'reservoir', #Reservoir
'43603': 'reservoir', #Reservoir
'43604': 'reservoir', #Reservoir
'43605': 'reservoir', #Reservoir
'43606': 'reservoir', #Reservoir
'43607': 'reservoir', #Reservoir
'43609': 'reservoir', #Reservoir
'43610': 'reservoir', #Reservoir
'43611': 'reservoir', #Reservoir
'43612': 'reservoir', #Reservoir
'43614': 'reservoir', #Reservoir
'43615': 'reservoir', #Reservoir
'43616': 'reservoir', #Reservoir
'43617': 'reservoir', #Reservoir
'43618': 'reservoir', #Reservoir
'43619': 'reservoir', #Reservoir
'43620': 'reservoir' #Reservoir
manmade = {
'43613': 'reservoir_covered' #Covered Reservoir
leisure = {
'43608': 'swimming_pool' #Swimming Pool
if 'fcode' in data:
keys = []
if data['fcode'] in landuse:
keys.append( ('landuse', landuse[data['fcode']]) )
if data['fcode'] in manmade:
keys.append( ('man_made', manmade[data['fcode']]) )
if data['fcode'] in leisure:
keys.append( ('leisure', leisure[data['fcode']]) )
return keys
return None
# The most important part of the code: define a set of key/value pairs
# to iterate over to generate keys. This is a list of two-tuples: first
# is a 'key', which is only used if the second value is a string. In
# that case, it is a map of lowercased fieldnames to OSM tag names: so
# fee_owner maps to 'owner' in the OSM output.
# if the latter is callable (has a __call__; is a function), then that
# method is called, passing in a dict of feature attributes with
# lowercased key names. Those functions can then return a list of
# two-tuples to be used as tags, or nothin' to skip the tags.
tag_mapping = [
('ftype', ftype),
('fcode', fcode),
('gnis_name', 'name'),
('gnis_id', 'gnis:feature_id'),
('elevation', 'ele'),
# These tags are not exported, even with the source data; this should be
# used for tags which are usually calculated in a GIS. AREA and LEN are
# common.
boring_tags = [ 'error', 'tile_x', 'tile_y' ]
# Namespace is used to prefix existing data attributes. If 'None', or
# '--no-source' is set, then source attributes are not exported, only
# attributes in tag_mapping.
namespace = "NHD"
#namespace = None
# Uncomment the "DONT_RUN = False" line to get started.
#DONT_RUN = True
DONT_RUN = False
# Set the maximum length of a way (in nodes) before it is split into
# shorter ways
Max_Waylength = 1500000
# =========== DO NOT CHANGE AFTER THIS LINE. ===========================
# Below here is regular code, part of the file. This is not designed to
# be modified by users.
# ======================================================================
import sys
from osgeo import ogr
except ImportError:
import ogr
except ImportError:
__doc__ += gdal_install
print __doc__
print "OGR Python Bindings not installed.\n%s" % gdal_install
def close_file():
""" Internal. Close an open file."""
global open_file
if not open_file.closed:
def start_new_file():
""" Internal. Open a new file, closing existing file if neccesary."""
global open_file, file_counter
file_counter += 1
if open_file:
open_file = open("%s%s.osm" % (file_name, file_counter), "w")
print >>open_file, "<osm version='0.6'>"
def clean_attr(val):
"""Internal. Hacky way to make attribute XML safe."""
val = str(val)
val = val.replace("&", "&amp;").replace("'", "&quot;").replace("<", "&lt;").replace(">", "&gt;").strip()
return val
def add_ring_nodes(ring):
"""Internal. Write the outer ring nodes."""
global open_file, id_counter
ringways = []
ids = []
firstnode = id_counter
if range(ring.GetPointCount() - 1) == 0 or ring.GetPointCount() == 0:
print >>sys.stderr, "Degenerate ring."
for count in range(ring.GetPointCount() - 1):
print >>open_file, "<node timestamp='1969-12-31T23:59:59Z' changeset='-1' id='%s' version='1' lon='%s' lat='%s' />" % (id_counter, ring.GetX(count), ring.GetY(count))
id_counter += 1
if (count > 0) and ((count % (Max_Waylength - 1)) == 0):
ids = []
ids.append(id_counter - 1)
return ringways
def add_ring_way(ring):
"""Internal. write out the 'holes' in a polygon."""
global open_file, id_counter
ids = []
ringways = []
for count in range(ring.GetPointCount() - 1):
print >>open_file, "<node timestamp='1969-12-31T23:59:59Z' changeset='-1' version='1' id='%s' lon='%s' lat='%s' />" % (id_counter, ring.GetX(count), ring.GetY(count))
id_counter += 1
if len(ids) == 0:
return None
print >>open_file, "<way timestamp='1969-12-31T23:59:59Z' changeset='-1' id='%s' version='1'>" % id_counter
way_id = id_counter
id_counter += 1
count = 0
for i in ids:
print >>open_file, "<nd ref='%s' />" % i
count += 1
if count >= Max_Waylength - 1:
count = 0
print >>open_file, "</way>"
print >>open_file, "<way timestamp='1969-12-31T23:59:59Z' changeset='-1'id='%s' version='1'>" % id_counter
way_id = id_counter
id_counter += 1
print >>open_file, "<nd ref='%s' />" % i
print >>open_file, "<nd ref='%s' />" % ids[0]
print >>open_file, "</way>"
return ringways
open_file = None
file_name = None
id_counter = 22951459320
file_counter = 0
counter = 0
class AppError(Exception): pass
def run(filename, slice_count=1, obj_count=100000000, output_location=None, no_source=False):
"""Run the converter. Requires open_file, file_name, id_counter,
file_counter, counter to be defined in global space; not really a very good
global id_counter, file_counter, counter, file_name, open_file, namespace
if no_source:
if output_location:
file_name = output_location
ds = ogr.Open(filename)
if not ds:
raise AppError("OGR Could not open the file %s" % filename)
l = ds.GetLayer(0)
max_objs_per_file = obj_count
extent = l.GetExtent()
if extent[0] < -180 or extent[0] > 180 or extent[2] < -90 or extent[2] > 90:
raise AppError("Extent does not look like degrees; are you sure it is? \n(%s, %s, %s, %s)" % (extent[0], extent[2], extent[1], extent[3]))
slice_width = (extent[1] - extent[0]) / slice_count
seen = {}
print "Running %s slices with %s base filename against shapefile %s" % (
slice_count, file_name, filename)
for i in range(slice_count):
l.SetSpatialFilterRect(extent[0] + slice_width * i, extent[2], extent[0] + (slice_width * (i + 1)), extent[3])
f = l.GetNextFeature()
obj_counter = 0
last_obj_split = 0
while f:
start_id_counter = id_counter
if f.GetFID() in seen:
f = l.GetNextFeature()
seen[f.GetFID()] = True
# done = False
# while f.GetField(ignoreField) in ignoreValues:
# if l.GetNextFeature() == None:
# done = True
# break
# f = l.GetNextFeature()
# if done == True:
# break
outerways = []
innerways = []
geom = f.GetGeometryRef()
ring = geom.GetGeometryRef(0)
objcount = ring.GetPointCount()
for i in range(1, geom.GetGeometryCount()):
objcount += geom.GetGeometryRef(i).GetPointCount()
objcount += 1
if (obj_counter - last_obj_split + objcount) > max_objs_per_file:
print "Splitting file with %s objs" % (obj_counter - last_obj_split)
last_obj_split = obj_counter
if objcount > max_objs_per_file:
print "Warning: a feature contains %i objects which is more than the %i object limit. It will be placed in a file by itself." % (objcount, max_objs_per_file)
ringways = add_ring_nodes(ring)
if not ringways or len(ringways) == 0:
f = l.GetNextFeature()
# Write out the outer ways in the relation
for count in range( len(ringways) ):
ids = ringways[count]
if ids and len(ids) > 1:
print >>open_file, "<way timestamp='1969-12-31T23:59:59Z' changeset='-1' id='%s' version='1'>" % id_counter
id_counter += 1
for i in ids:
print >>open_file, "<nd ref='%s' />" % i
#Write out the fields for the way
field_count = f.GetFieldCount()
fields = {}
for field in range(field_count):
value = f.GetFieldAsString(field)
name = f.GetFieldDefnRef(field).GetName()
if namespace and name and value and name not in boring_tags:
print >>open_file, "<tag k='%s:%s' v='%s' />" % (namespace, name, clean_attr(value))
fields[name.lower()] = value
#Perform the specified field mappting
for tag_name, map_value in tag_mapping:
if hasattr(map_value, '__call__'):
tag_values = map_value(fields)
if tag_values:
for tag in tag_values:
tags[tag[0]] = tag[1]
if tag_name in fields:
tags[map_value] = fields[tag_name].title()
for key, value in tags.items():
if key and value:
print >>open_file, "<tag k='%s' v='%s' />" % (key, clean_attr(value))
#Write fixed tabs
for name, value in fixed_tags.items():
print >>open_file, "<tag k='%s' v='%s' />" % (name, clean_attr(value))
print >>open_file, "</way>"
if (geom.GetGeometryCount() > 1) or (len(ringways) > 1):
#add the inner ways
for i in range(1, geom.GetGeometryCount()):
ringways = add_ring_way(geom.GetGeometryRef(i))
for way in ringways:
print >>open_file, "<relation timestamp='1969-12-31T23:59:59Z' changeset='-1' id='%s' version='1' ><tag k='type' v='multipolygon' />" % id_counter
id_counter += 1
for way in outerways:
print >>open_file, '<member type="way" ref="%s" role="outer" />' % way
for way in innerways:
print >>open_file, '<member type="way" ref="%s" role="inner" />' % way
print >>open_file, "</relation>"
counter += 1
f = l.GetNextFeature()
obj_counter += (id_counter - start_id_counter)
if __name__ == "__main__":
print __doc__
from optparse import OptionParser
parse = OptionParser(usage="%prog [args] filename.shp", version=__version__)
parse.add_option("-s", "--slice-count", dest="slice_count",
help="Number of horizontal slices of data", default=1,
action="store", type="int")
parse.add_option("-o", "--obj-count",
help="Maximum number of objects in a single .osm file",
default=1000000, type="int")
parse.add_option("-n", "--no-source", dest="no_source",
help="Do not store source attributes as tags.",
action="store_true", default=False)
parse.add_option("-l", "--output-location",
dest="output_location", help="base filepath for output files.",
(options, args) = parse.parse_args()
if not len(args):
print "No shapefile name given!"
kw = {}
for key in ('slice_count', 'obj_count', 'output_location', 'no_source'):
kw[key] = getattr(options, key)
run(args[0], **kw)
except AppError, E:
print "An error occurred: \n%s" % E