Extracting vertices from polygon geometry

12.10.2010 18:01 ·  GIS  ·  gdal, python, tips

Imagine you have a vector layer containing polygons, both single- and multipart. Your task is to get vertices of each geometry, including all rings and parts. Here is how it can be done using GDAL and Python.

Let’s start with simplest case — singlepart polygons. In that case everything is straightforward — we need to find out the number of vertices, iterate over them and get coordinated. For demostration purposes the code will look like this:

def dump_geometry(geom):
    print "Exterior ring"
    ring = geom.GetGeometryRef(0)
    for i in range(ring.GetPointCount()):
        print i, ring.GetX(i), ring.GetY(i)

A polygon with “holes” is considered as a set of “rings”: one outer ring defining the boundary and several inner rings describing the “holes”. The outer ring always comes first. Thus, the algorithm will transform as follows: loop over all rings; if it is the first iteration - this is an outer ring, otherwise - inner ring; extract vertices from the ring and output their coordinates.

def dump_geometry(geom):
    nRings = geom.GetGeometryCount()
    for i in range(nRings):
        ring = geom.GetGeometryRef(i)
        if i == 0:
            print "Exterior ring"
        else:
            print "Interior ring", i

        for j in range(ring.GetPointCount()):
            print j, ring.GetX(j), ring.GetY(j)

Handling multipart polygons is more complicated: they can be either singlepart polygons or contain “holes”, “islands”, or even “holes” and “islands” at the same time. Let’s try to figure out how to deal with them. A polygon consisting of several “islands” can be considered as a set of independent polygons. In turn, each of those polygons can have one or more rings. So what we have now:

1. if the polygon is a singlepart polygon and contains only one ring, we extract its vertices immediately
2. if the number of rings in polygon is greater than 1 — we loop through the rings and extract vertices from them
3. if the polygon is a multipart polygon — loop through the parts and for each for part perform checks starting from the step 1

So, this is a classical recursion and the final function will look like this:

def dump_geometry(geom, file):
    geomType = geom.GetGeometryType() # geometry type
    # singlepart polygon
    if geomType == ogr.wkbPolygon or geomType == ogr.wkbPolygon25D:
        nRings = geom.GetGeometryCount() # number of rings
        print "Rings", nRings
        if nRings == 1: # there are no "holes" in the polygons
            print "Exterior ring"
            ring = geom.GetGeometryRef(0)
            for i in range(ring.GetPointCount()):
                print i, ring.GetX(i), ring.GetY(i)
        else: # polygon with "holes"
            for i in range(nRings):
                ring = geom.GetGeometryRef(i)
                if i == 0:
                    print "Exterior ring"
                else:
                    print "Interior ring", i

                for j in range(ring.GetPointCount()):
                    print j, ring.GetX(j), ring.GetY(j)
    # multipart polygon
    elif geomType == ogr.wkbMultiPolygon or geomType == ogr.wkbMultiPolygon25D:
        for i in range(geom.GetGeometryCount()):
            subGeom = geom.GetGeometryRef(i) # for each polygon part
            dump_geometry(subGeom) # recursively call dump_geometry() function

The code is a bit bloated, but my aim was to demonstrate a possible way of approaching the task, not to come up with an ideal solution.