Joins OSM road segments that share the same name

join_osm_roads.py

from psycopg2 import connect
from shapely.wkb import loads, dumps
from shapely.ops import linemerge
from osgeo import ogr
from osgeo.osr import SpatialReference
import os

db = connect('dbname=osm')
cur = db.cursor()
cur.execute('SELECT way, name, ST_SRID(way) FROM planet_osm_roads WHERE "highway" = \'secondary\'')

streetLines = {}
maxLen = 0
srid = None

for way, name, _srid in cur:
    if name is None:
        name = 'n/a'
    if name not in streetLines:
        streetLines[name] = []
        maxLen = max(maxLen, len(name))
    streetLines[name].append(loads(way.decode('hex')))

    if srid is None:
        srid = _srid
    elif srid != _srid:
        raise Exception('mulitple spatial references found in table')

for name in streetLines:
    geometry = linemerge(streetLines[name])
    geoms = hasattr(geometry, 'geoms') and geometry.geoms or [geometry]
    streetLines[name] = geoms

if os.path.exists('secondary-streets.shp'):
    os.remove('secondary-streets.shp')

# store as shapefile
driver = ogr.GetDriverByName('ESRI Shapefile')
source = driver.CreateDataSource('secondary-streets.shp')
srs = SpatialReference()
srs.ImportFromEPSG(srid)
layer = source.CreateLayer('default', srs, ogr.wkbLineString)

# add field definition for street names
nameField = ogr.FieldDefn('NAME', ogr.OFTString)
nameField.SetWidth(maxLen)
layer.CreateField(nameField)

feat = ogr.Feature(layer.GetLayerDefn())

for name in streetLines:
    for line in streetLines[name]:
        feat = ogr.Feature(layer.GetLayerDefn())
        feat.SetField(nameField.name, name)
        geom = ogr.CreateGeometryFromWkb(dumps(line))
        feat.SetGeometry(geom)
        layer.CreateFeature(feat)