Shapely 1.2a6 with pictures

One thing that Shapely has lacked is one or two dirt simple example programs to keep the API real and help explain its use. I did something about this over the past couple of nights: 1.2a6 includes two easy to understand, easy to run scripts. I hope users profit from them. Myself, I found that they demanded a new and improved API feature. I'll explain.

First, here's an example of using Shapely to construct patches by growing buffer regions out from a set of points and dissolving those regions together as they intersect, and plotting the results with Matplotlib. This is run-of-the-mill GIS stuff, yes, but done in style.

http://trac.gispython.org/lab/raw-attachment/wiki/Examples/dissolve.png

A plate of blue-speckled brains splattered on the floor, or is it just me?

The interesting part of the complete, amply-documented dissolve.py script is here:

import pylab
from shapely.ops import cascaded_union

patches = cascaded_union(spots)

pylab.figure(num=None, figsize=(4, 4), dpi=180)

for patch in patches.geoms:
    x, y = patch.exterior.xy
    pylab.fill(x, y, color='#cccccc', aa=True)
    pylab.plot(x, y, color='#666666', aa=True, lw=1.0)
    for hole in patch.interiors:
        x, y = hole.xy
        pylab.fill(x, y, color='#ffffff', aa=True)
        pylab.plot(x, y, color='#999999', aa=True, lw=1.0)

pylab.text(-25, 25,
    "Patches: %d, total area: %.2f" % (len(patches.geoms), patches.area))

pylab.savefig('dissolve.png')

The xy property is completely new in 1.2a6, inspired by how awkwardly I had to slice and dice coordinates when writing this example against 1.2a5. It provides two Python arrays that are immediately usable with Numpy or Matplotlib. Speaking of Matplotlib: I'd love to know how to fill a patch but not its holes (you'll notice that I'm faking the emptiness of the holes in this example).

What would would you have to go through to pyplot ArcGIS scripting results?

Shapely doesn't just make grey matter go splat, it can also toss brains in the air and pierce them with lasers:

http://trac.gispython.org/lab/raw-attachment/wiki/Examples/intersect.png

Or make a fair facsimile thereof. What's really going on in intersect.py is an analysis of a HTML5 geolocation (latitude, longitude, heading, and speed) trajectory's intersection with a cluster of patches. The intercepted patches are plotted in red and the intersecting segments of the trajectory itself are also plotted in red. Finally, scalar properties of different geometries are used in a text label. The example vector intercepts 2 of the 7 patches along 5 segments with a total length (to one decimal place) of 26.1:

import pylab
from shapely.geometry import LineString

# Represent the following geolocation parameters
#
# initial position: -25, -25
# heading: 45.0
# speed: 50*sqrt(2)
#
# as a line
vector = LineString(((-25.0, -25.0), (25.0, 25.0)))

# Find intercepted and missed patches. List the former so we can count them
intercepts = [patch for patch in patches.geoms if vector.intersects(patch)]
misses = (patch for patch in patches.geoms if not vector.intersects(patch))

pylab.figure(num=None, figsize=(4, 4), dpi=180)

for spot in misses:
    x, y = spot.exterior.xy
    pylab.fill(x, y, color='#cccccc', aa=True)
    pylab.plot(x, y, color='#999999', aa=True, lw=1.0)
    for hole in spot.interiors:
        x, y = hole.xy
        pylab.fill(x, y, color='#ffffff', aa=True)
        pylab.plot(x, y, color='#999999', aa=True, lw=1.0)

for spot in intercepts:
    x, y = spot.exterior.xy
    pylab.fill(x, y, color='red', alpha=0.25, aa=True)
    pylab.plot(x, y, color='red', alpha=0.5, aa=True, lw=1.0)
    for hole in spot.interiors:
        x, y = hole.xy
        pylab.fill(x, y, color='#ffffff', aa=True)
        pylab.plot(x, y, color='red', alpha=0.5, aa=True, lw=1.0)

pylab.arrow(-25, -25, 50, 50, color='#999999', aa=True,
    head_width=1.0, head_length=1.0)

intersection = vector.intersection(patches)
for segment in intersection.geoms:
    x, y = segment.xy
    pylab.plot(x, y, color='red', aa=True, lw=1.5)

pylab.text(-28, 25,
    "Patches: %d/%d (%d), total length: %.1f" \
     % (len(intercepts), len(patches.geoms),
        len(intersection.geoms), intersection.length))

pylab.savefig('intersect.png')

Install GEOS 3.2.0 (Windows users can get it from a PostGIS 1.5 installer, but will have to copy the DLLs to a location one can glean only from looking at shapely/geos.py. YMMV until we have Shapely 1.2 installers) then grab the new distribution with easy_install or pip (as well as Numpy and Matplotlib) and give them a try:

$ python /usr/local/bin/dissolve.py
$ python /usr/local/bin/intersect.py

I think this is pretty much the last 1.2 alpha.