More IPython Notebook and Shapely

$ ipython notebook --pylab inline

from functools import reduce
from itertools import islice
from descartes import PolygonPatch
from shapely.geometry import Point
from shapely.ops import unary_union

BLUE = '#6699cc'

def halton(base):
    """Returns an iterator over an infinite Halton sequence."""
    def value(index):
        result = 0.0
        f = 1.0/base
        i = index
        while i > 0:
            result += f * (i % base)
            i = i/base
            f = f/base
        return result
    i = 1
    while i > 0:
        yield value(i)
        i += 1

def draw(axes, item):
    """Given matplotlib axes and a geometric item, adds the item as a patch
    and returns the axes so that reduce() can accumulate more patches."""
    axes.add_patch(
        PolygonPatch(item, fc=BLUE, ec=BLUE, alpha=0.5, zorder=2))
    return axes

# Zip together two 100-item sequences to make 100 pseudo-random points.
coords = zip(
    list(islice(halton(5), 20, 120)),
    list(islice(halton(7), 20, 120)) )

# Buffer the points to make overlapping patches.
patches = [Point(xy).buffer(0.06) for xy in coords]

# Note: with ipython's --pylab option, we've effectively imported
# all symbols from matplotlib's pylab module.
figsize(8, 8)

# Perform a left fold on the patches, applying the draw function (above)
# with the current axes as the accumulator. Aka "map rendering for hipsters."
reduce(draw, patches, gca())

# Dissolve continuously overlapping patches with the unary_union function.
u = unary_union(patches)
print u.geom_type
print u.area

# Output:
# MultiPolygon
# 0.87333863506

# A MultiPolygon is an iterator over its polygon parts, so we can perform
# a fold on it as well.
reduce(draw, u, gca())