I am interested in finding all the vertices of all polygons which are not on the inside the most outer polygon that can be formed in a shapefile. Here is how I get all the vertices:

def __init__(self,shapeFilePath):
    driver = ogr.GetDriverByName('ESRI Shapefile')
    self.LandUses = []
    self.Codes = [] #Store Categories
    self.NPolygons = 1 # number of polygons in layer
    self.dataSource = driver.Open(shapeFilePath, 0)
    self.layer = self.dataSource.GetLayer()
    self.NPolygons = self.layer.GetFeatureCount()
    self.Boundaries = [""]*self.NPolygons   # store vertices of each polygon      
    for polygon in range(self.NPolygons):
        area = self.layer.GetFeature(polygon)
        Name = area.GetField(0)
        Kategorie = area.GetField(1)
        geometry = area.GetGeometryRef()
        boundary_raw = str(geometry.GetBoundary())
        #remove tail and head
        boundary = boundary_raw[12:-1]
        boundary = boundary.split(',')
        #print boundary
        #lenboundary = len(boundary)
        #convert each coordinate from string to float
        for x, point in enumerate(boundary):
            boundary[x] = point.split()
        boundingvertices=np.asarray(boundary, dtype=np.float64)
        self.Boundaries[polygon] = boundingvertices

But now I am stuck with a lot of vertices, and I don't know how to get the "mantle", here is a picture displaying what I need:

A Sample layer showing multiple polygons

Before I go hacking my own solution, I'd like to know if the GDAL Gods have been kind and made something like this already.

UPDATE: Found something call convex hull... not sure, but testing...
UPDATE #2: Of course one could just join all the polygons ...

  • 3
    There is a typical GIS process called "dissolve" that is likely what you are after (i.e. joins all of the polygons). Finding the convex (or alpha) hull I don't think is appropriate. The convex hull won't work for concave polygons, and finding the alpha hull to the desired resolution seems like over-kill compared to just dissolving the polygon. – Andy W Jun 12 '12 at 19:11
  • @AndyW, thanks for the hint! I found a pretty good solution. See my answer below. – Oz123 Jun 12 '12 at 22:36

I found a decent solution IMHO using Shapely. I am not limited in the amount of dependencies for that package I am writing (Hence, GDAL, is not my only tool). It also seems that Shapely can be used for many other things in my code.

A quick way to "dissolve" the polygons would be to use Shapely's

from shapely.ops import cascaded_union

enter image description here

See the full example of dissolving polygons.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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