Tuesday, March 15, 2011

jython, JTS Topology, and Voronoi Diagrams

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.

2 comments:

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

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

    ReplyDelete