Python OGR: convert a point geometry into polygon

Question

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

Answer 1

The OGR convexHull method returns a polygon so the answer to your question is 'no'. You do not need to parse the vertices or convert the geometry.

I notice that you are defining thisGeometry to be of type ogr.wkbPoint. wkbPoint geometry is for storing single points and not multiple points or polygons. You also don't need to parse every single point but can add entire geometries to your collection wholesale. For your use-case I suggest you try:

thisGeometry = ogr.Geometry(ogr.wkbGeometryCollection)
for index in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()
    thisGeometry.AddGeometry(geometry)

convexHull = thisGeometry.ConvexHull()

The result of this should be a normal convex hull polygon.

Answer 2

In this case what is happening is that you are creating a geometry of type Point. Point geometries can only contain one actual point.

thisGeometry.GetGeometryCount() #1

ConvexHull() actully creates a geometry of the minimum number of dimensions that contains all of the geometry of the parent object. So if there is one point, the convex hull will be a point, for two points it will be a line (assuming no coincident points), and for three or more it will be a polygon (assuming the points aren't in a line and aren't coincident).

To check this the easiest way is to use the ogr.GeometryTypeToName method:

ogr.GeometryTypeToName(convexHull.GetGeometryType())

The easiest way for you to get the convex hull for a group of islands would be to create a GeometryCollection and get the convex hull from that, like so:

geom_collection = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in layer: #layers (and datasets) have a list interface
    geometry = feature.GetGeometryRef()
    geom_collection.AddGeometry(geometry)

convex_hull = geom_collection.ConvexHull()

And this convex hull should be a polygon (as it's created from a collection of polygons) which means that you should be able to use it as input for your output dataset.

If you still wanted to get the convex hull of multiple points you would have to create a geometry of type MultiPoint and then add the points to that MultiPoint object individually.

thisGeometry = ogr.Geometry(ogr.wkbMultiPoint)
for feature in layer:
    geometry = feature.GetGeometryRef()
    for ring in geometry: #This is assuming that the geometry is a polygon
        for point in ring:
            thisGeometry.AddGeometry(point)

convexHull = thisGeometry.ConvexHull()

Answer 3

I realised that I did not use the proper type for the initial collection of points. I set it to wkbLineString, which allowed me to use the AddPoint method on the collection (in the initial code, I had only 1 point at the end in my collection.

I also made some corrections, by adding a field (that is mandatory) and destroying the data source to force writing to file.

Hereafter the corrected code:

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.wkbLineString)
    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" )
    if os.path.exists(outfile):
        drv.DeleteDataSource(outfile)
    ds = drv.CreateDataSource( outfile )
    if ds is None:
        usage("Could not create file {0}".format( outfile) )

    # fields
    fldDfn = ogr.FieldDefn('id', ogr.OFTInteger)
    fldDfn.SetWidth(4)

    lyrname = "convexHull_${0}".format( layer.GetName() )
    lyr = ds.CreateLayer( lyrname, layer.GetSpatialRef(), ogr.wkbPolygon )
    lyr.CreateField(fldDfn)

    thisFeature = ogr.Feature( lyr.GetLayerDefn() )
    thisFeature.SetGeometry( convexHull )
    thisFeature.SetField('id',1)
    lyr.CreateFeature( thisFeature )

    ds.Destroy()