Dissolve small & residential roads into urban areas.

dissolver.py

from sys import stderr
from osgeo import ogr
from time import time
from math import sqrt

def not_big(feature):
    highway = feature.GetFieldAsString(feature.GetFieldIndex('type'))
    
    return highway in ('residential', 'unclassified', 'road', 'minor')

start = time()

ff = 100 # "football field"

fn = '/Users/migurski/Documents/MapBox/project/osm-carto/shps/areas_line_z10/areas_line_z10.shp'
fn = '/tmp/manila-z12.shp'
fn = '/tmp/uppsala-z12.shp'
fn = '/tmp/port-au-prince.osm-roads/port-au-prince.osm-roads.shp'

ds = ogr.Open(fn)
layer = ds.GetLayer()

ds2 = ds.GetDriver().CreateDataSource('/tmp/dissolved.shp')
layer2 = ds2.CreateLayer('dissolved', layer.GetSpatialRef(), ogr.wkbMultiPolygon)

features = [layer.GetFeature(fid) for fid in range(layer.GetFeatureCount())]
features = filter(not_big, features)
geometries = [feature.GetGeometryRef() for feature in features]
geometries.sort(key=lambda geom: geom.Centroid().GetY())

print len(geometries), 'geometries in %.1f seconds' % (time() - start)

start = time()

buffers = [geom.Buffer(5*ff, 3) for geom in geometries]

print len(buffers), 'buffers in %.1f seconds' % (time() - start)

start = time()

def cascaded_union(polys):
    '''
    '''
    if len(polys) == 2:
        return polys[0].Union(polys[1])

    if len(polys) == 1:
        return polys[0]
    
    if len(polys) == 0:
        return None
    
    half = len(polys) / 2
    poly1 = cascaded_union(polys[:half])
    poly2 = cascaded_union(polys[half:])
    
    return poly1.Union(poly2)

dissolved = cascaded_union(buffers)

print 'dissolved in %.1f seconds' % (time() - start)

start = time()

shrunken = dissolved.Buffer(-3*ff, 3).Buffer(-4*ff, 3)
reinflated = shrunken.Buffer(4*ff, 3)

print 'reinflated in %.1f seconds' % (time() - start)

result = ogr.Geometry(ogr.wkbMultiPolygon)

for i in range(reinflated.GetGeometryCount()):
    geom = reinflated.GetGeometryRef(i)
    
    if geom.GetGeometryCount() == 1 and geom.Area() > (10*ff)**2:
        print int(sqrt(geom.Area()))
        result.AddGeometry(geom)
    
    elif geom.GetGeometryCount() > 1:
        print 'here it is...', 
        part = ogr.Geometry(ogr.wkbPolygon)
        part.AddGeometry(geom.GetGeometryRef(0))
        
        for j in range(1, geom.GetGeometryCount()):
            hole = geom.GetGeometryRef(j)
            
            if hole.Area() > (10*ff)**2:
                part.AddGeometry(hole)
        
        result.AddGeometry(part)

feature = ogr.Feature(ogr.FeatureDefn())
feature.SetGeometry(result)
layer2.CreateFeature(feature)