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