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.

## Sunday, March 27, 2011

## 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.

## Sunday, March 13, 2011

### jython + the JTS geometry library

My favorite library for dealing with two dimensional points, lines, and polygons is Polygon by Jörg Rädler. In a Java/jython environment, though, that's not available. This led me to try my hand at JTS.

Greg Corradini has a post that highlights JTS use with jython and geoscript. His overview was very helpful, but I wanted to get to know JTS more thoroughly.

The first thing I came up against when using JTS was the documentation. It's all javadocs. This isn't necessarily a bad thing, but, along with the library's java-centric design, it takes a little getting used to. If concepts like abstract class, iterface, extending, and implementing don't immediately ring a bell to the user, use of the javadocs will make them clear over time.

That shouldn't discourage us, though. Let's fire up the interpreter and try something simple.

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')

[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')

I usually store things I'm testing out in my Downloads folder (I'm using KDE on OpenBSD). This may not be the best practice. In any case, to import from jts I need access to the jar in my sys.path.

Making a Point object in JTS in a bit complicated. There are Coordinate objects and there are Point objects. The Coordinate is a more basic class with less functionality than the Point class. Point inherits from the abstract class Geometry. This leads it to have a number of methods, some meaningless, like convexHull. Shown below is the process required to make a Point object.

>>> from com.vividsolutions.jts.geom import Point

>>> from com.vividsolutions.jts.geom import GeometryFactory

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

>>> # get CoordinateArraySequence of size one

>>> cas = CoordinateArraySequence(1)

>>> factoryx = GeometryFactory()

>>> ptx = Point(cas, factoryx)

>>> from com.vividsolutions.jts.geom import GeometryFactory

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

>>> # get CoordinateArraySequence of size one

>>> cas = CoordinateArraySequence(1)

>>> factoryx = GeometryFactory()

>>> ptx = Point(cas, factoryx)

We have a Point object. A couple explanatory notes:

1) the documentation for the default GeometryFactory reads thus: "Constructs a GeometryFactory that generates Geometries having a floating PrecisionModel and a spatial-reference ID of 0." For what I'm doing this is fine.

2) The preferred constructor for the Point object takes a CoordinateArraySequence of length one instead of a single Coordinate object. This is analogous to a character in Python being a string of length one.

We don't know yet what our values are for our Point's coordinates. Let's have a look.

>>> ptx

POINT (0 0)

POINT (0 0)

OK, zero-zero. Let's assign something to those coordinates away from the origin.

>>> ptx.getCoordinates()[0].x = 44

>>> ptx.getCoordinates()[0].y = 22

>>> ptx

POINT (44 22)

>>> ptx.getCoordinates()[0].y = 22

>>> ptx

POINT (44 22)

Even though the Point has only one coordinate, we still need to specify that with the index 0.

Well, that was a lot of work for just one point. Let's try something more efficient and fun.

>>> cas2 = CoordinateArraySequence(6)

>>> somecoordinates = [(1, 5), (7, 14), (22, 44), (36, 12), (19, 1), (12, 4)]

>>> for numx in range(len(somecoordinates)):

... cas2.getCoordinate(numx).x = somecoordinates[numx][0]

... cas2.getCoordinate(numx).y = somecoordinates[numx][1]

...

>>> cas2

((1.0, 5.0, NaN), (7.0, 14.0, NaN), (22.0, 44.0, NaN), (36.0, 12.0, NaN), (19.0, 1.0, NaN), (12.0, 4.0, NaN))

>>> from com.vividsolutions.jts.geom import LineString

>>> ls = LineString(cas2, factoryx)

>>> ch = ls.convexHull()

>>> ch

POLYGON ((19 1, 1 5, 22 44, 36 12, 19 1))

>>>

LineString is one of the classes that can be instantiated from a CoordinateArraySequence. I went with that mainly for ease of use.

There is a great deal more to JTS. This post was mainly for getting started with the library.

>>> somecoordinates = [(1, 5), (7, 14), (22, 44), (36, 12), (19, 1), (12, 4)]

>>> for numx in range(len(somecoordinates)):

... cas2.getCoordinate(numx).x = somecoordinates[numx][0]

... cas2.getCoordinate(numx).y = somecoordinates[numx][1]

...

>>> cas2

((1.0, 5.0, NaN), (7.0, 14.0, NaN), (22.0, 44.0, NaN), (36.0, 12.0, NaN), (19.0, 1.0, NaN), (12.0, 4.0, NaN))

>>> from com.vividsolutions.jts.geom import LineString

>>> ls = LineString(cas2, factoryx)

>>> ch = ls.convexHull()

>>> ch

POLYGON ((19 1, 1 5, 22 44, 36 12, 19 1))

>>>

LineString is one of the classes that can be instantiated from a CoordinateArraySequence. I went with that mainly for ease of use.

There is a great deal more to JTS. This post was mainly for getting started with the library.

## Wednesday, March 2, 2011

### jython + Google Guava's base.CaseFormat

I was trying to research Google Collections for use in jython and stumbled upon Solid Craft's blog and his description of Guava base's CaseFormat utility. Where I work we're a Java shop, and use of camelCase infiltrates all code, even jython code. CaseFormat provides functionality for changing these variable names to a number of different formats.

$ /usr/local/jdk-1.7.0/bin/java -jar jython.jar

Jython 2.5.1 (Release_2_5_1:6813, Sep 26 2009, 13:47:54)

[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal

Type "help", "copyright", "credits" or "license" for more information.

>>> import sys

>>> # get downloaded jar into sys.path

>>> sys.path.append('/home/carl/Downloads/guava-r08/guava-r08.jar')

>>> from com.google.common.base import CaseFormat as cf

>>> dir(cf)

['LOWER_CAMEL', 'LOWER_HYPHEN', 'LOWER_UNDERSCORE', 'UPPER_CAMEL', . . .]

>>> varname = 'someVariable'

>>> cf.LOWER_CAMEL.to(cf.LOWER_UNDERSCORE, varname)

u'some_variable'

>>> cf.LOWER_CAMEL.to(cf.UPPER_UNDERSCORE, varname)

u'SOME_VARIABLE'

>>> cf.LOWER_CAMEL.to(cf.UPPER_CAMEL, varname)

u'SomeVariable'

>>>

For cycling through a file full of camel case or underscored variable names, even by hand in the interpreter, this could save some time and typing.

$ /usr/local/jdk-1.7.0/bin/java -jar jython.jar

Jython 2.5.1 (Release_2_5_1:6813, Sep 26 2009, 13:47:54)

[OpenJDK Client VM (Sun Microsystems Inc.)] on java1.7.0-internal

Type "help", "copyright", "credits" or "license" for more information.

>>> import sys

>>> # get downloaded jar into sys.path

>>> sys.path.append('/home/carl/Downloads/guava-r08/guava-r08.jar')

>>> from com.google.common.base import CaseFormat as cf

>>> dir(cf)

['LOWER_CAMEL', 'LOWER_HYPHEN', 'LOWER_UNDERSCORE', 'UPPER_CAMEL', . . .]

>>> varname = 'someVariable'

>>> cf.LOWER_CAMEL.to(cf.LOWER_UNDERSCORE, varname)

u'some_variable'

>>> cf.LOWER_CAMEL.to(cf.UPPER_UNDERSCORE, varname)

u'SOME_VARIABLE'

>>> cf.LOWER_CAMEL.to(cf.UPPER_CAMEL, varname)

u'SomeVariable'

>>>

For cycling through a file full of camel case or underscored variable names, even by hand in the interpreter, this could save some time and typing.

Subscribe to:
Posts (Atom)