Last time I barely scratched the surface of the JTS library. As it turns out, through ignorance I made a few things harder than they needed to be. The GeometryFactory class is quite versatile and has a number of methods for creating objects from simple Coordinate objects. In the example below, I've managed to create a clipped Voronoi diagram without too much code or difficulty.

import sys

# this line will be different depending on where

# the jts jar is stored and may be unnecessary if

# the jar is in your CLASSPATH

sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')

from com.vividsolutions.jts.geom import Coordinate

from com.vividsolutions.jts.geom import GeometryFactory

from com.vividsolutions.jts.triangulate import VoronoiDiagramBuilder

from com.vividsolutions.jts.geom import Envelope

# first some coordinates

coord1 = Coordinate(5, 7)

coord2 = Coordinate(3, 6)

coord3 = Coordinate(1, 4)

coord4 = Coordinate(2, 5)

coord5 = Coordinate(4, 2)

coord6 = Coordinate(7, 2)

coord7 = Coordinate(5, 5)

coord8 = Coordinate(2, 6)

coord5 = Coordinate(4, 2)

coords = [coord1, coord2,

coord3, coord4,

coord5, coord6,

coord7, coord8]

# the handy GeometryFactory

geomfactx = GeometryFactory()

# a MultiPoint Geometry object

mp = geomfactx.createMultiPoint(coords)

# the diagram builder class

vdb = VoronoiDiagramBuilder()

# it would be nice to clip the

# diagram with a rectangle;

# this is what the Envelope

# object is for

lowerleft = Coordinate(0, 0)

upperright = Coordinate(8, 8)

env = Envelope(lowerleft, upperright)

vdb.setClipEnvelope(env)

# load our sites (points) from the

# MultiPoint object

vdb.setSites(mp)

# this is where the builder object clips the

# the polygons with the Envelope rectangle

diagramx = vdb.getDiagram(geomfactx)

print 'Number of polygons = %d' % diagramx.getNumGeometries()

polygonx = diagramx.getGeometryN(0)

print 'One polygon:'

for coordsetx in polygonx.getCoordinates():

print coordsetx

Output:

Number of polygons = 8

One polygon:

(-5.0, -4.0, NaN)

(-5.0, 8.25, NaN)

(0.5, 5.5, NaN)

(2.7, 3.3, NaN)

(-2.166666666666666, -4.0, NaN)

(-5.0, -4.0, NaN)

There's more that can be done here (for instance, clipping the polygons with a non-rectangular shape). As was true of the last post, this is just a small fraction of what the JTS Topology library can do.

You might want to use the VoronoiDiagramBuilder getDiagram method - that's a bit simpler than getting the Subdivision

ReplyDelete@DrJTS Thanks so much for the tip. I've changed the code accordingly.

ReplyDelete