I'm writing a routine for creating the convex hull envelop around shapefiles of islands. I could read the data, extract polygons vertices into a geometry and then compute the convex hull.
Now I need to save it to a shapefile of polygons.
Is there a way to transform the geometry (named convexHull
in my code) returned by the convexHull method into a polygon geometry; or am I forced to parse all vertices of convexHull
and create a ring of a polygon?
def doConvexHull(infile, outfile):
inH = ogr.Open(infile, 0)
if inH is None:
usage("Could not open file {0}. Exit.".format(infile))
layer = inH.GetLayer()
# get all polygons
thisGeometry = ogr.Geometry(ogr.wkbPoint)
for index in xrange(layer.GetFeatureCount()):
feature = layer.GetFeature(index)
geometry = feature.GetGeometryRef()
ring = geometry.GetGeometryRef(0)
points = ring.GetPointCount()
for p in xrange(points):
lon, lat, z = ring.GetPoint(p)
thisGeometry.AddPoint(lon, lat)
convexHull = thisGeometry.ConvexHull()
drv = ogr.GetDriverByName( "ESRI Shapefile" )
ds = drv.CreateDataSource( outfile )
if ds is None:
usage("Could not create file {0}".format( outfile) )
lyrname = "convexHull_${0}".format( layer.GetName() )
lyr = ds.CreateLayer( lyrname, layer.GetSpatialRef(), ogr.wkbPolygon )
thisFeature = ogr.Feature( layer.GetLayerDefn() )
thisFeature.SetGeometry( convexHull )
lyr.CreateFeature( thisFeature ): error from geometry
No comments:
Post a Comment