Unwrapped JTS with Python

While Justin has been hard at work wrapping JTS and GeoTools in clear and simple Python syntax, I’ve been tinkering with some JTS classes that aren’t wrapped in GeoScript Python. Well, not wrapped yet. One of the great things about using Jython (or any scripting language written for the Java platform) is that you can still explore JTS and GeoTools even if some functionality isn’t implemented in GeoScript.

Jared produced a convex hull and minimum bounding circle with GeoScript Groovy in a recent post. I tried to replicate this workflow but quickly realized certain classes and functions weren’t available in GeoScript Python. Did it matter? Not very much in this case. Here’s the convex hull and minimum bounding circle for a shapefile of Washington weather stations produced in a previous post. Pay attention to what’s being imported.

from geoscript.layer import Shapefile
from geoscript.viewer import plot
from com.vividsolutions.jts.geom import GeometryCollection
from com.vividsolutions.jts.geom import GeometryFactory
from com.vividsolutions.jts.geom import PrecisionModel
from com.vividsolutions.jts.algorithm import MinimumBoundingCircle

#Get the point shapefile
staLayer = Shapefile('/home/gregcorradini/wa_stations.shp')

#Create a list of feature geometries
pntGeoms = [f.geom for f in staLayer.features()]

#Create a GeometryCollection from our list of points
geomColl = GeometryCollection(pntGeoms,GeometryFactory(PreciscionModel(),4326))

#Get the geometry collection's convex hull
geomConvexH = geomColl.convexHull()

#Get the geometry collection's minimum bounding circle
mCircle = (MinimumBoundingCircle(geomColl)).getCircle()

# Plot the geometries, BANG
plot(pntGeoms + [geomConvexH] + [mCircle])

I imported GeoScript modules and JTS Java classes in the example above. Jython can interoperate with Java libraries and classes because under the hood it’s Java. But how did I know what parameters a GeometryCollection accepted? I took a peek at the JTS Java doc and followed the bouncing argument ball. The weaving of colorful geometries happens after that thanks to the great plotting functionality from Justin.

There are tricks and trapdoors to this unwrapped exploration - the syntax is not as clean and intuitive as working with GeoScript Python. But a lot of magic still happens. I try to use functions that return geometry objects so I can plot them with the geoscript.viewer.plot function. This way I can see each step visually. Another advantage is that JTS and GeoTools have very well organized Java docs (here and there). It’s hard to fail at finding your way with road maps this good.

For example, mkgeomatics was deciding how to visually represent the output of a least-cost analysis he performed on a Bellingham, WA dataset. Part of the cartographic workflow used Voronoi diagrams that were built in ArcGIS. Both of us wondered if we could leverage GeoScript Python and JTS to solve the Voronoi diagram work. After a little sandbox scripting I came up with the following test code. It outlines how to derive a Voronoi diagram from a Delaunay triangulation (A more succinct way to build Voronoi diagrams would use the com.vividsolutions.jts.triangulate.VoronoiDiagramBuilder class. But Dr JTS produced an informative post on the relationship between both that I wanted to explore).

from geoscript.layer import Shapefile
from geoscript.viewer import plot
from com.vividsolutions.jts.geom import GeometryCollection
from com.vividsolutions.jts.geom import GeometryFactory
from com.vividsolutions.jts.geom import PrecisionModel
from com.vividsolutions.jts.algorithm import MinimumBoundingCircle
from com.vividsolutions.jts.triangulate import DelaunayTriangulationBuilder

#Get a point shapefile
staLayer = Shapefile('/home/gregcorradini/wa_stations.shp')

#Create a list of feature geometries
pntGeoms = [f.geom for f in staLayer.features()]

#Create a GeometryCollection from our list of points
geomColl = GeometryCollection(pntGeoms,GeometryFactory(PreciscionModel(),4326))

#Create a DelaunayTriangulationBuilder object and pass our geometry collection
dtb = DelaunayTriangulationBuilder()
dtb.setSites(geomColl)

#Get the Delaunay triangle faces
triangles = dtb.getTriangles(GeometryFactory(PrecisionModel(),4326))

#Make a pretty triangle picture
plot(pntGeoms + [triangles])

#Get the object's corresponding voronoi polygons in a java array
voronoiArray = dtb.getSubdivision().getVoronoiCellPolygons(GeometryFactory(PrecisionModel(),4326))

#Make a pretty picture
plot(pntGeoms + [i for i in voronoiArray])

That last Voronoi picture is hard to see in detail. So let's save off the Voronoi polygons to a shapefile and give each one the attributes of the original point that falls within it. The function 'changeGeometryType' is from a previous post.

#Create a voronoi geom.Poly schema based on old layer schema
voronoiSchema = changeGeometryType('voronoiWAStations',geom.Polygon,staLayer)

#Create voronoi layer
vLayer = ws.create(schema=voronoiSchema)

#For each voronoi polygon test which point is within it and pass attributes
for poly in voronoiArray:
 for f in staLayer.features():
 
  #Is my point within the polygon?
  if f.geom.within(poly):
   
   #Copy all point feature attributes into a dictionary
   dictFeat = dict(f.iteritems())
   
   #Overwrite the_geom with voronoi polygon
   dictFeat['the_geom'] = poly
   
   #Create a new feature and add it to voronoi layer
   vLayer.add(voronoiSchema.feature(dictFeat))

The zoomed-in shot above shows the original points in pink, Delaunay triangles in gray and Voronoi diagram in green.

Buffering Features with GeoScript Python

In a previous post, Jared buffered points and added them to a polygon shapefile. Here's the workflow in Python GeoScript.

from geoscript import proj
from geoscript import geom
from geoscript.layer import Shapefile
from geoscript.feature import schema

def changeGeometryType(strName,geomType,shpParent):

   #Get all shapefile schema attributes except 'the_geom'
   parentS = [(str(i.name),i.typ) for i in shpParent.schema.fields if not issubclass(i.typ,geom.Geometry)]

   #Get shapefile Projection
   pPrj = shpParent.schema.get('the_geom').proj

   #Reinsert geomType into schema
   parentS.insert(0,('the_geom',geomType,pPrj))
             
   #Create centroid buffer schema
   centroidSchema = schema.Schema(strName,parentS)

   return centroidSchema


if __name__=='__main__':

   #Get centroid shapefile
   shp = Shapefile('/home/gregcorradini/GeoTools/geoscript_wrk/centroid_buffer/data/centroids.shp')

   #Create a geom.Polygon schema from centroid schema
   centroidSchema = changeGeometryType('centroid_buffer',geom.Polygon,shp)

   #Get shapefile workspace
   ws = shp.workspace

   #Create a buffered centroid layer based on centroid schema
   buffCentLayer = ws.create(schema=centroidSchema)

   #For each centroid feature, create a buffered centroid, add it to buffered layer
   for f in shp.features():
      
       #Copy all feature attributes into a dictionary
       dictFeat = dict(f.iteritems())
      
       #Buffer 'the_geom'  3 degrees
       dictFeat['the_geom'] = f.geom.buffer(3)

       print dictFeat
      
       #Create a buffered centroid feature with new dictionary attributes
       buffFeature = buffCentLayer.schema.feature(dictFeat,f.id)

       #Add new feature to buffered layer
       buffCentLayer.add(buffFeature)

Calculating Centroids with GeoScript Python

Jared and Tim have shown us how to create centroids from a polygon shapefile in two different GeoScript implementations. Now let's see what the code looks like in Python GeoScript.

from geoscript import proj
from geoscript import geom
from geoscript.layer import Shapefile
from geoscript.feature import schema

def changeGeometryType(strName,geomType,shpParent):

 #Put shapefile schema attributes except 'the geom' in a list of tuples
 parentS = [(str(i.name),i.typ) for i in shpParent.schema.fields if not issubclass(i.typ,geom.Geometry)]

 #Get shapefile Projection
 pPrj = shpParent.schema.field('the_geom').proj

 #reinsert geomType and proj into schema
 parentS.insert(0,('the_geom',geomType,pPrj))
  
 #create centroid schema with the list
 centroidSchema = schema.Schema(strName,parentS)

 return centroidSchema


if __name__=='__main__':

 #Get states shapefile
 shp = Shapefile('/home/gregcorradini/GeoTools/geoscript_wrk/centroid_buffer/data/US_shapefiles.shp')

 #Substitute geom.Point into schema
 centroidSchema = changeGeometryType('centroids',geom.Point,shp)

 #Get workspace
 ws = shp.workspace

 #Create a centroid layer based on centroidSchema
 centroidLayer = ws.create(schema=centroidSchema)

 #Get all state features
 stateFeatures = [f for f in shp.features()]

 #For each state polygon get its centroid and add it to our centroid layer
 for f in stateFeatures:
 
    #create dictionary to hold each feature's attributes
    newAtts = {}
    
    for key in f.attributes.keys():
  if issubclass(f.attributes[key].__class__,geom.Geometry):
   newAtts[str(key)] = f.geom.centroid
  else:
   newAtts[str(key)] = f.attributes[key]

    print newAtts

    #Create a feature based on dictionary
    centFeature = centroidLayer.schema.feature(newAtts,f.id)

    #Add it to centroid layer
    centroidLayer.add(centFeature)

Calculating Centroids with GeoScript Groovy

GeoScript uses the GeoTools and Java Topology Suite libraries to provide easy-to-use scripting APIs with implementations in Python, JavaScript, Scala and Groovy. The GeoScript web site provides simple code snippets in all four languages to help you get started. In this post, I address a real world example: calculating centroids from one GIS layer and saving them to another layer using Groovy. This example highlights the GeoScript Layer, Feature, and Geometry modules.

// import GeoScript modules
import geoscript.layer.*
import geoscript.feature.*
import geoscript.geom.*

// This is the Shapefile we want to extract centroids from
Shapefile shp = new Shapefile('states.shp')

// Create a new Schema (based on the above Shapefile, but with Point Geometry)
Schema schema = shp.schema.changeGeometryType('Point','states_centroids')

// Create the new centroid Layer
Layer centroidLayer = shp.workspace.create(schema)

// Use the Cursor to loop through each Feature in the Shapefile
Cursor cursor = shp.cursor
while(cursor.hasNext()) {

   // Get the next Feature
   Feature f = cursor.next()
 
   // Create a Map for the new attributes
   Map attributes = [:]

   // For each attribute in the shapefile
   f.attributes.each{k,v ->

      // If its Geometry, find the centroid
      if (v instanceof Geometry) {
          attributes[k] = v.centroid
      }
      // Else, they stay the same
      else {
         attributes[k] = v
      }
   }

   // Create a new Feature with the new attributes
   Feature feature = schema.feature(attributes, f.id)

   // Add it to the centroid Layer
   centroidLayer.add(feature)
}

// Always remember to close the cursor
cursor.close()

Here is what you get when you run the script against a shapefile of the 48 contiguous states (sorry Alaska and Hawaii):