My last post dealt with a similar subject, but it didn't quite tell the whole story. No sooner had I posted it and attempted to implement it in a project, I ran into a problem. Granted, it was on a Java GUI in Java code (not a Python project) - but the concept is universal. When is a floating point number small enough that another mathematical operation on it will take it to zero or to the realm of infinity?
Python (I'm using CPython 2.5) is pretty forgiving in this sense.
>>> 1.0/1.0e-308
1e+308
>>> 1.0/1.0e-309
inf
There comes a point when even a representable number equals zero.
>>> 0.0 == 1.0e-10000
True
With some trigonometric functions from the math module:
>>> math.atan2(1.0, 0.0000000000000001)
1.5707963267948966
>>> 2 * _
3.1415926535897931
>>> _ == math.pi
True
>>> math.cos(1.0000000000000000001e-250)
1.0
The equation for the slope of a line in the Cartesian plane is (y2 - y1)/(x2 - x1) where y2, y1, x2, and x1 are the x and y coordinates of the enpoints of the line.
One way I've handled the infinite slope problem in the past is by setting a value as infinity:
# one million
>>> INFINITY = 1.0e+6
>>> slope = 1000000.0
>>> slope == INFINITY
True
This works sometimes, but it's not totally safe from error.
What I ended up doing this time was checking for zero or near zero values in the numerator and denominator of the slope equation.
NEARZERO = 1.0e-6
if math.abs(y2 - y1) < NEARZERO:
# set y2 == y1 or vice versa
# and proceed with the calculations
# for a zero slope
if math.abs(x2 - x1) < NEARZERO:
# set x2 == x1 or vice versa
# and proceed with the calculations
# for an infinite slope
This solved my problem in the sense that it gave reasonable results that served the problem domain well.
Lastly, you could use a similar method to check for duplicate points. If both numerator and denominator of the slope equation are near zero, the line segment is too small and the point should be done away with:
if (math.abs(y2 - y1) < NEARZERO and
math.abs(x2 - x1) < NEARZERO):
# duplicate point
# delete it
Disclaimer: I am not a floating point math expert. There may be better solutions out there. These are a couple practical ones that have worked for me.
Saturday, April 2, 2011
Sunday, March 27, 2011
Polygon Offset Revisited
I had previously written a post concerning polygon offset that was incomplete in its scope. Specifically, it dodged the problem of horizontal and vertical lines in polygons.
Since then I've had to confront this problem again. Shown below is the code I'm using for the more complete solution. It's a bit verbose, but it works.
import math
def calcoffsetpoint(pt1, pt2, offset):
"""
Get a point offset from the line
segment pt1-pt2 distance "offset".
Return two tuple of coordinates.
"""
theta = math.atan2(pt2[1] - pt1[1],
pt2[0] - pt1[0])
theta += math.pi/2.0
return (pt1[0] - math.cos(theta) * offset,
pt1[1] - math.sin(theta) * offset)
def getoffsetintercept(pt1, pt2, m, offset):
"""
From points pt1 and pt2 defining a line
in the Cartesian plane, the slope of the
line m, and an offset distance,
calculates the y intercept of
the new line offset from the original.
"""
x, y = calcoffsetpoint(pt1, pt2, offset)
return y - m * x
def getpt(pt1, pt2, pt3, offset):
"""
Gets intersection point of the two
lines defined by pt1, pt2, and pt3;
offset is the distance to offset
the point from the polygon.
Valid for lines with slopes other
than zero or infinity.
Returns a two tuple of coordinates.
"""
# get first offset intercept
m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
boffset = getoffsetintercept(pt1, pt2, m, offset)
# get second offset intercept
mprime = (pt3[1] - pt2[1])/(pt3[0] - pt2[0])
boffsetprime = getoffsetintercept(pt2, pt3, mprime, offset)
# get intersection of two offset lines
newx = (boffsetprime - boffset)/(m - mprime)
newy = m * newx + boffset
return newx, newy
def getslopeandintercept(pt1, pt2, offset):
"""
Gets the slope and the intercept of the
offset line.
Result returned as a two tuple.
"""
m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
b = getoffsetintercept(pt1, pt2, m, offset)
return m, b
def getoffsetcornerpoint(pt1, pt2, pt3, offset):
"""
Gets intersection point of the two
lines defined by pt1, pt2, and pt3;
offset is the distance to offset
the point from the polygon.
Returns a two tuple of coordinates.
"""
# starting out with horizontal line
if (pt2[1] - pt1[1]) == 0.0:
ycoord = pt1[1] - math.cos(math.atan2(0.0, pt2[0] - pt1[0])) * offset
# a vertical line follows
if (pt3[0] - pt2[0]) == 0.0:
xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
# a sloped line follows
else:
m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
# calculate for x with ycoord
xcoord = (ycoord - offsetintercept)/m
# starting out with a vertical line
if (pt2[0] - pt1[0]) == 0.0:
xcoord = pt1[0] + math.sin(math.atan2(pt2[1] - pt1[1], 0.0)) * offset
# a horizontal line follows
if (pt3[1] - pt2[1]) == 0.0:
ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
# a sloped line follows
else:
m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
# calculate for y with xcoord
ycoord = m * xcoord + offsetintercept
# starting out with sloped line
if (pt2[1] - pt1[1]) != 0.0 and (pt2[0] - pt1[0]) != 0.0:
# if second line is horizontal
if (pt3[1] - pt2[1]) == 0.0:
ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
# calculate for x with y coord
xcoord = (ycoord - offsetintercept)/m
# if second line is vertical
elif (pt3[0] - pt2[0]) == 0.0:
xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
# solve for y with x coordinate
ycoord = m * xcoord + offsetintercept
# if both lines are sloped
else:
xcoord, ycoord = getpt(pt1, pt2, pt3, offset)
return xcoord, ycoord
def offsetpolygon(polyx, offset):
"""
Offsets a clockwise list of coordinates
polyx distance offset to the inside of
the polygon.
Returns list of offset points.
"""
polyy = []
# need three points at a time
for counter in range(0, len(polyx) - 3):
# get first offset intercept
pt = getoffsetcornerpoint(polyx[counter],
polyx[counter + 1],
polyx[counter + 2],
offset)
# append new point to polyy
polyy.append(pt)
# last three points
pt = getoffsetcornerpoint(polyx[-3], polyx[-2], polyx[-1], offset)
polyy.append(pt)
pt = getoffsetcornerpoint(polyx[-2], polyx[-1], polyx[0], offset)
polyy.append(pt)
pt = getoffsetcornerpoint(polyx[-1], polyx[0], polyx[1], offset)
polyy.append(pt)
return polyy
I used this polygon to test it.
This is a fairly simple geometric problem, but it can be confusing. Hopefully this will save someone the trouble of working through a solution.
Notes:
1) If you're using Python 2.x, feed your coordinates in as floats. This will ensure the absence of integer division, which often yields zeros, which in turn yields ZeroDivision errors.
2) The polygon fed to the offsetpolygon function should not have a closing point. The closing point will cause ZeroDivision errors. Instead, append the point before drawing the polygon in gnuplot or some other plotting tool.
Since then I've had to confront this problem again. Shown below is the code I'm using for the more complete solution. It's a bit verbose, but it works.
import math
def calcoffsetpoint(pt1, pt2, offset):
"""
Get a point offset from the line
segment pt1-pt2 distance "offset".
Return two tuple of coordinates.
"""
theta = math.atan2(pt2[1] - pt1[1],
pt2[0] - pt1[0])
theta += math.pi/2.0
return (pt1[0] - math.cos(theta) * offset,
pt1[1] - math.sin(theta) * offset)
def getoffsetintercept(pt1, pt2, m, offset):
"""
From points pt1 and pt2 defining a line
in the Cartesian plane, the slope of the
line m, and an offset distance,
calculates the y intercept of
the new line offset from the original.
"""
x, y = calcoffsetpoint(pt1, pt2, offset)
return y - m * x
def getpt(pt1, pt2, pt3, offset):
"""
Gets intersection point of the two
lines defined by pt1, pt2, and pt3;
offset is the distance to offset
the point from the polygon.
Valid for lines with slopes other
than zero or infinity.
Returns a two tuple of coordinates.
"""
# get first offset intercept
m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
boffset = getoffsetintercept(pt1, pt2, m, offset)
# get second offset intercept
mprime = (pt3[1] - pt2[1])/(pt3[0] - pt2[0])
boffsetprime = getoffsetintercept(pt2, pt3, mprime, offset)
# get intersection of two offset lines
newx = (boffsetprime - boffset)/(m - mprime)
newy = m * newx + boffset
return newx, newy
def getslopeandintercept(pt1, pt2, offset):
"""
Gets the slope and the intercept of the
offset line.
Result returned as a two tuple.
"""
m = (pt2[1] - pt1[1])/(pt2[0] - pt1[0])
b = getoffsetintercept(pt1, pt2, m, offset)
return m, b
def getoffsetcornerpoint(pt1, pt2, pt3, offset):
"""
Gets intersection point of the two
lines defined by pt1, pt2, and pt3;
offset is the distance to offset
the point from the polygon.
Returns a two tuple of coordinates.
"""
# starting out with horizontal line
if (pt2[1] - pt1[1]) == 0.0:
ycoord = pt1[1] - math.cos(math.atan2(0.0, pt2[0] - pt1[0])) * offset
# a vertical line follows
if (pt3[0] - pt2[0]) == 0.0:
xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
# a sloped line follows
else:
m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
# calculate for x with ycoord
xcoord = (ycoord - offsetintercept)/m
# starting out with a vertical line
if (pt2[0] - pt1[0]) == 0.0:
xcoord = pt1[0] + math.sin(math.atan2(pt2[1] - pt1[1], 0.0)) * offset
# a horizontal line follows
if (pt3[1] - pt2[1]) == 0.0:
ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
# a sloped line follows
else:
m, offsetintercept = getslopeandintercept(pt2, pt3, offset)
# calculate for y with xcoord
ycoord = m * xcoord + offsetintercept
# starting out with sloped line
if (pt2[1] - pt1[1]) != 0.0 and (pt2[0] - pt1[0]) != 0.0:
# if second line is horizontal
if (pt3[1] - pt2[1]) == 0.0:
ycoord = pt2[1] - math.cos(math.atan2(0.0, pt3[0] - pt2[0])) * offset
m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
# calculate for x with y coord
xcoord = (ycoord - offsetintercept)/m
# if second line is vertical
elif (pt3[0] - pt2[0]) == 0.0:
xcoord = pt2[0] + math.sin(math.atan2(pt3[1] - pt2[1], 0.0)) * offset
m, offsetintercept = getslopeandintercept(pt1, pt2, offset)
# solve for y with x coordinate
ycoord = m * xcoord + offsetintercept
# if both lines are sloped
else:
xcoord, ycoord = getpt(pt1, pt2, pt3, offset)
return xcoord, ycoord
def offsetpolygon(polyx, offset):
"""
Offsets a clockwise list of coordinates
polyx distance offset to the inside of
the polygon.
Returns list of offset points.
"""
polyy = []
# need three points at a time
for counter in range(0, len(polyx) - 3):
# get first offset intercept
pt = getoffsetcornerpoint(polyx[counter],
polyx[counter + 1],
polyx[counter + 2],
offset)
# append new point to polyy
polyy.append(pt)
# last three points
pt = getoffsetcornerpoint(polyx[-3], polyx[-2], polyx[-1], offset)
polyy.append(pt)
pt = getoffsetcornerpoint(polyx[-2], polyx[-1], polyx[0], offset)
polyy.append(pt)
pt = getoffsetcornerpoint(polyx[-1], polyx[0], polyx[1], offset)
polyy.append(pt)
return polyy
I used this polygon to test it.
This is a fairly simple geometric problem, but it can be confusing. Hopefully this will save someone the trouble of working through a solution.
Notes:
1) If you're using Python 2.x, feed your coordinates in as floats. This will ensure the absence of integer division, which often yields zeros, which in turn yields ZeroDivision errors.
2) The polygon fed to the offsetpolygon function should not have a closing point. The closing point will cause ZeroDivision errors. Instead, append the point before drawing the polygon in gnuplot or some other plotting tool.
Thursday, March 24, 2011
jython + joda-time
joda-time is a date-time library written in Java. It has advantages over Java's built in time utilities in terms of power and ease of use. It has some functionality that Python's own datetime library does not. Here is a sampling from within jython:
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> # make sure the jar is in the CLASSPATH
>>> sys.path.append('/home/carl/Downloads/joda-time-1.6.2/joda-time-1.6.2.jar')
>>> from org.joda import time as joda
>>> rightnow = joda.DateTime()
>>> rightnow
2011-03-24T19:54:45.544-06:00
>>> rightnow.getZone()
America/Edmonton
>>> # I'm actually in New Mexico - close enough
>>> rightnow.dayOfWeek().getAsText()
u'Thursday'
>>> rightnow.monthOfYear().getAsText()
u'March'
>>> rightnow.getDayOfMonth()
24
>>> # days are one indexed along with months
>>> # a date from the past
>>> declofindependence = joda.DateTime(1776, 7, 4, 0, 0, 0, 0)
The Period object allows you to count down to dates in the future of count from dates in the past in terms of years, months, and days:
>>> americasage = joda.Period(declofindependence, rightnow)
>>> americasage.getYears()
234
>>> americasage.getMonths()
8
>>> americasage.getWeeks()
2
>>> americasage.getDays()
6
joda-time can handle dates far into the future as well as those in the far distant past:
>>> farfuture = joda.DateTime(10000, 1, 1, 1, 0, 0, 0)
>>> farfuture.getYear()
10000
>>> farfuture
10000-01-01T01:00:00.000-07:00
>>>
Like most things Java, joda-time can be a bit more verbose than its Python equivalent. Still, the ability to get dates beyond 9999 and some of the functionality may make it well worth the trouble.
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> # make sure the jar is in the CLASSPATH
>>> sys.path.append('/home/carl/Downloads/joda-time-1.6.2/joda-time-1.6.2.jar')
>>> from org.joda import time as joda
>>> rightnow = joda.DateTime()
>>> rightnow
2011-03-24T19:54:45.544-06:00
>>> rightnow.getZone()
America/Edmonton
>>> # I'm actually in New Mexico - close enough
>>> rightnow.dayOfWeek().getAsText()
u'Thursday'
>>> rightnow.monthOfYear().getAsText()
u'March'
>>> rightnow.getDayOfMonth()
24
>>> # days are one indexed along with months
>>> # a date from the past
>>> declofindependence = joda.DateTime(1776, 7, 4, 0, 0, 0, 0)
The Period object allows you to count down to dates in the future of count from dates in the past in terms of years, months, and days:
>>> americasage = joda.Period(declofindependence, rightnow)
>>> americasage.getYears()
234
>>> americasage.getMonths()
8
>>> americasage.getWeeks()
2
>>> americasage.getDays()
6
joda-time can handle dates far into the future as well as those in the far distant past:
>>> farfuture = joda.DateTime(10000, 1, 1, 1, 0, 0, 0)
>>> farfuture.getYear()
10000
>>> farfuture
10000-01-01T01:00:00.000-07:00
>>>
Like most things Java, joda-time can be a bit more verbose than its Python equivalent. Still, the ability to get dates beyond 9999 and some of the functionality may make it well worth the trouble.
Saturday, March 19, 2011
WKT (Well Known Text) + JTS Topology + jython
The last few posts have dealt with the JTS Topology library. This one will wrap up that topic with a method for converting JTS's Geometry objects to and from text, WKT.
WKT is a standard of the Open Geospatial Consortium (OGC). Objects expressed in WKT are not much different than their string representations in JTS. Below is a code example of working back and forth between JTS and WKT.
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06)
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')
>>> from com.vividsolutions.jts.geom import GeometryFactory
>>> from com.vividsolutions.jts.geom import Coordinate
>>> from com.vividsolutions.jts.geom import LinearRing
>>> from com.vividsolutions.jts.io import WKTWriter
>>> from com.vividsolutions.jts.io import WKTReader
>>> coord1 = Coordinate(0, 0)
>>> coord2 = Coordinate(0, 1)
>>> coord3 = Coordinate(1, 1)
>>> coord4 = Coordinate(1, 0)
>>> coord5 = Coordinate(0, 0)
>>> coords = [coord1, coord2, coord3, coord4, coord5]
>>> from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
>>> cas = CoordinateArraySequence(coords)
>>> gf = GeometryFactory()
>>> lr = LinearRing(cas, gf)
>>> asquare = gf.createPolygon(lr, None)
>>> f = open('wkttest.xml', 'w')
>>> wktwriter = WKTWriter()
>>> # this is just my own personal preference
>>> # it is inefficient for storage, but
>>> # sometimes easier on the eyes
>>> wktwriter.setMaxCoordinatesPerLine(1)
>>> text = wktwriter.writeFormatted(asquare)
>>> f.write(text)
>>> f.close()
>>> print text
POLYGON ((0 0,
0 1,
1 1,
1 0,
0 0))
>>> wktreader = WKTReader()
>>> f = open('wkttest.xml', 'r')
>>> text = f.read()
>>> asquare2 = wktreader.read(text)
>>> asquare2
POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))
It's fairly easy to store JTS Geometry objects in human readable (WKT) format and reuse them.
There is a good bit more to WKT than what's shown here. This was just a simple example for getting started.
WKB (Well Known Binary) is also available in JTS - this is used for efficient storage of Geometry objects, mainly in databases.
JTS also has a WKTFileReader class which can read in multiple geometries from a file.
WKT is a standard of the Open Geospatial Consortium (OGC). Objects expressed in WKT are not much different than their string representations in JTS. Below is a code example of working back and forth between JTS and WKT.
Jython 2.5.2 (Release_2_5_2:7206, Mar 2 2011, 23:12:06)
[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal
Type "help", "copyright", "credits" or "license" for more information.
>>> import sys
>>> sys.path.append('/home/carl/Downloads/jts/lib/jts-1.11.jar')
>>> from com.vividsolutions.jts.geom import GeometryFactory
>>> from com.vividsolutions.jts.geom import Coordinate
>>> from com.vividsolutions.jts.geom import LinearRing
>>> from com.vividsolutions.jts.io import WKTWriter
>>> from com.vividsolutions.jts.io import WKTReader
>>> coord1 = Coordinate(0, 0)
>>> coord2 = Coordinate(0, 1)
>>> coord3 = Coordinate(1, 1)
>>> coord4 = Coordinate(1, 0)
>>> coord5 = Coordinate(0, 0)
>>> coords = [coord1, coord2, coord3, coord4, coord5]
>>> from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
>>> cas = CoordinateArraySequence(coords)
>>> gf = GeometryFactory()
>>> lr = LinearRing(cas, gf)
>>> asquare = gf.createPolygon(lr, None)
>>> f = open('wkttest.xml', 'w')
>>> wktwriter = WKTWriter()
>>> # this is just my own personal preference
>>> # it is inefficient for storage, but
>>> # sometimes easier on the eyes
>>> wktwriter.setMaxCoordinatesPerLine(1)
>>> text = wktwriter.writeFormatted(asquare)
>>> f.write(text)
>>> f.close()
>>> print text
POLYGON ((0 0,
0 1,
1 1,
1 0,
0 0))
>>> wktreader = WKTReader()
>>> f = open('wkttest.xml', 'r')
>>> text = f.read()
>>> asquare2 = wktreader.read(text)
>>> asquare2
POLYGON ((0 0, 0 1, 1 1, 1 0, 0 0))
It's fairly easy to store JTS Geometry objects in human readable (WKT) format and reuse them.
There is a good bit more to WKT than what's shown here. This was just a simple example for getting started.
WKB (Well Known Binary) is also available in JTS - this is used for efficient storage of Geometry objects, mainly in databases.
JTS also has a WKTFileReader class which can read in multiple geometries from a file.
Friday, March 18, 2011
Polygon Buffering with JTS Topology and jython
Previously, I had done a post regarding polygon offset. It turns out that JTS Topology has polygon buffering capability. I wanted to try this out and compare it to my results from the polygon offset post.
As with most things with a library like JTS, the software does most of the work for you:
# path to jar into 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
# LinearRing is for the creation of the Polygon
from com.vividsolutions.jts.geom import LinearRing
# CoordinateSequenceArray is for the
# creation of the LinearRing
from com.vividsolutions.jts.geom.impl import CoordinateArraySequence
# CAP_ROUND is for the Polygon.buffer operation
from com.vividsolutions.jts.operation.buffer.BufferOp import CAP_ROUND
# coordinates
# PT 1
MONASTERY = [(1.1, 0.75),
# PT 2
(1.2, 1.95),
# PT 3
(1.9, 1.96),
.
.
.
# PT 21
(1.1, 0.75)]
gf = GeometryFactory()
coords = [Coordinate(*coord) for coord in MONASTERY]
cas = CoordinateArraySequence(coords)
lr = LinearRing(cas, gf)
polyx = gf.createPolygon(lr, None)
# offset of 0.15 inward
# rounded edges
# 10 points per quarter circle of rounded edges
polyinset = polyx.buffer(-0.15, 10, CAP_ROUND)
The result is shown below:
The rounded edges give an entirely different look to the shape. My version is shown below:
I'm glad I did the egg the way I did it (it looks more realistic). Nonetheless, the polygon buffering capability of JTS is a useful tool, particularly for calculating offsets and distances for scientific or geographic purposes.
Notes: Shapely is a Python package (CPython) with much of the same functionality as JTS.
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.
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.
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.
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.
Subscribe to:
Posts (Atom)