Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Dissolve small & residential roads into urban areas.
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.