Thursday, March 17, 2011

Clipping Voronoi Polygons with JTS Topology and jython

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.

No comments:

Post a Comment