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