Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question

3 Answers 3

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.

share|improve this answer
1  
Many thanks for you comment (I saw it after posting an answer to my own question). I will use your snippet in my code. Thanks again! –  Bruno von Paris Feb 1 '13 at 11:32
    
Minor note: convex hull returns the minimum geometry dimension that covers the input geometry. –  om_henners Feb 1 '13 at 11:35
    
A correct and fair point, but that will be a polygon for all practical real-world purposes :) –  MappaGnosis Feb 1 '13 at 11:43

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()
share|improve this answer
    
Thanks for the detailed comments and code snippets! –  Bruno von Paris Feb 1 '13 at 13:27

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()
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.