Saturday, June 12, 2010

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)




No comments:

Post a Comment

Introducing GeoScript

GeoScript adds geo capabilities to dynamic scripting languages such as JavaScript, Python, Scala and Groovy.