Following up on the last post on Voronoi diagrams and the JTS topology library, I will show clipping of the Voronoi polygons with another polygon (a polygonal bounding region) in this post.

The plot below is what we have currently, a Voronoi diagram (with points) clipped along a rectangular boundary:

Shown below is what I'd like to accomplish - clipping the red polygons with the blue surrounding one:

Here is the code, continued from the last post:

# because JTS can handle holes and "donuts", it

# has a LinearRing structure for constructing

# polygons

from com.vividsolutions.jts.geom import LinearRing

# the LinearRing requires a CoordinateArraySequence

# (it will not accept a list)

from com.vividsolutions.jts.geom.impl import CoordinateArraySequence

boundingcoords = [(10, 0), (2, 0),

(-2, 4), (2, 12),

(10, 6), (10, 0)]

boundingcoords = [Coordinate(*coord) for coord in boundingcoords]

casboundingcoords = CoordinateArraySequence(boundingcoords)

lr = LinearRing(casboundingcoords, geomfactx)

# None here fills in the spot where any inner rings

# would be in the polygon

boundingpoly = geomfactx.createPolygon(lr, None)

newpolygons = []

numpolygons = diagramx.getNumGeometries()

for numx in range(numpolygons):

newpolygons.append(diagramx.getGeometryN(

numx).intersection(boundingpoly))

And the result:

Note: it is possible to get more than one polygon as the result of each intersection. Through visual inspection I determined there would only be one polygon per intersection with the bounding area and structured my list (newpolygons) accordingly.

## Thursday, March 17, 2011

Subscribe to:
Post Comments (Atom)

## No comments:

## Post a Comment