Join OSM road segments that share the same reference number

gistfile1.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-czech')
cur = db.cursor()
types = "'motorway'"
cur.execute('SELECT way, ref, highway, ST_SRID(way) FROM planet_osm_roads WHERE "highway" IN (%s)' % types)

streetLines = {}
streetTypes = {}
maxLen = 0
srid = None

for way, ref, highway, _srid in cur:
    if ref is None:
        ref = 'n/a'
    ref = ref.decode('ISO-8859-1')
    if ref not in streetLines:
        streetLines[ref] = []
        streetTypes[ref] = highway
        maxLen = max(maxLen, len(ref))
    streetLines[ref].append(loads(way.decode('hex')))

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

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

if os.path.exists('czech-highways.shp'):
    os.remove('czech-highways.shp')

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

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

# add field definition for street names
typeField = ogr.FieldDefn('highway', ogr.OFTString)
typeField.SetWidth(10)
layer.CreateField(typeField)

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

for ref in streetLines:
    for line in streetLines[ref]:
        feat = ogr.Feature(layer.GetLayerDefn())
        feat.SetField(nameField.name, ref.encode('utf-8'))
        feat.SetField(typeField.name, streetTypes[ref])
        geom = ogr.CreateGeometryFromWkb(dumps(line))
        feat.SetGeometry(geom)
        layer.CreateFeature(feat)