435 lines
15 KiB
Python
435 lines
15 KiB
Python
|
#!/usr/bin/python
|
||
|
|
|
||
|
|
"""
|
||
|
|
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:
|
||
|
|
|
||
|
|
DO NOT CHANGE AFTER THIS LINE.
|
||
|
|
|
||
|
|
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 <crschmidt@crschmidt.net>"
|
||
|
|
__version__ = "$Id: polyshp2osm-WaterBody.py,v 1.3 2009/03/30 20:54:13 tbook Exp $"
|
||
|
|
|
||
|
|
gdal_install = """
|
||
|
|
Installing GDAL depends on your platform. Information is available at:
|
||
|
|
|
||
|
|
http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries
|
||
|
|
|
||
|
|
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
|
||
|
|
|
||
|
|
try:
|
||
|
|
try:
|
||
|
|
from osgeo import ogr
|
||
|
|
except ImportError:
|
||
|
|
import ogr
|
||
|
|
except ImportError:
|
||
|
|
__doc__ += gdal_install
|
||
|
|
if DONT_RUN:
|
||
|
|
print __doc__
|
||
|
|
sys.exit(2)
|
||
|
|
print "OGR Python Bindings not installed.\n%s" % gdal_install
|
||
|
|
sys.exit(1)
|
||
|
|
|
||
|
|
def close_file():
|
||
|
|
""" Internal. Close an open file."""
|
||
|
|
global open_file
|
||
|
|
if not open_file.closed:
|
||
|
|
open_file.write("</osm>")
|
||
|
|
open_file.close()
|
||
|
|
|
||
|
|
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:
|
||
|
|
close_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("&", "&").replace("'", """).replace("<", "<").replace(">", ">").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."
|
||
|
|
return
|
||
|
|
for count in range(ring.GetPointCount() - 1):
|
||
|
|
ids.append(id_counter)
|
||
|
|
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):
|
||
|
|
ringways.append(ids)
|
||
|
|
ids = []
|
||
|
|
ids.append(id_counter - 1)
|
||
|
|
ids.append(firstnode)
|
||
|
|
ringways.append(ids)
|
||
|
|
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):
|
||
|
|
ids.append(id_counter)
|
||
|
|
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>"
|
||
|
|
ringways.append(way_id)
|
||
|
|
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>"
|
||
|
|
ringways.append(way_id)
|
||
|
|
|
||
|
|
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
|
||
|
|
singleton."""
|
||
|
|
global id_counter, file_counter, counter, file_name, open_file, namespace
|
||
|
|
|
||
|
|
if no_source:
|
||
|
|
namespace=None
|
||
|
|
|
||
|
|
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.ResetReading()
|
||
|
|
l.SetSpatialFilterRect(extent[0] + slice_width * i, extent[2], extent[0] + (slice_width * (i + 1)), extent[3])
|
||
|
|
|
||
|
|
start_new_file()
|
||
|
|
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()
|
||
|
|
continue
|
||
|
|
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)
|
||
|
|
start_new_file()
|
||
|
|
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()
|
||
|
|
continue
|
||
|
|
# 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
|
||
|
|
outerways.append(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
|
||
|
|
tags={}
|
||
|
|
#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]
|
||
|
|
else:
|
||
|
|
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:
|
||
|
|
innerways.append(way)
|
||
|
|
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)
|
||
|
|
|
||
|
|
close_file()
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
if DONT_RUN:
|
||
|
|
print __doc__
|
||
|
|
sys.exit(2)
|
||
|
|
|
||
|
|
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",
|
||
|
|
dest="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.",
|
||
|
|
default="poly_output")
|
||
|
|
(options, args) = parse.parse_args()
|
||
|
|
|
||
|
|
if not len(args):
|
||
|
|
print "No shapefile name given!"
|
||
|
|
parse.print_help()
|
||
|
|
sys.exit(3)
|
||
|
|
|
||
|
|
kw = {}
|
||
|
|
for key in ('slice_count', 'obj_count', 'output_location', 'no_source'):
|
||
|
|
kw[key] = getattr(options, key)
|
||
|
|
|
||
|
|
try:
|
||
|
|
run(args[0], **kw)
|
||
|
|
except AppError, E:
|
||
|
|
print "An error occurred: \n%s" % E
|