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 am interested in finding all the vertices of all polygons which are not on the inside the most outer polygoin that can be formed in a shape file: 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)
        self.LandUses.append(Name)
        self.Codes.append(Kategorie)
        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.

Thanks a bunch! Oz

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

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

1 Answer

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
...
cascaded_union(polygons)

enter image description here

See the full example of dissolving polygons.

share|improve this answer
add comment

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.