How to extract interior polygon coordinates using Shapely?

Question

I am new to Shapely (but enthusiastic about it), and recently I've discovered a bit of a road bump.

I have a polygon shapefile that I am reading in via Fiona. This shapefile contains BOTH polygon and multipolygon items and I need to build an array for each feature of all the coordinates within it (i.e. both exterior and/or interior). Notably, two of the polygon items have interior rings (and they are valid).

I seem to have no problem accessing the exterior coordinates of the polygon(s)/multipolygon(s) ... but I am not pulling anything for the interior coordinates.

Do I need to take a new approach here (i.e. LinearRings)...?

def convert_polygons(inFile):

    for polys in fiona.open(inFile):
        myShape = shape(polys['geometry'])
        exterior_poly = 0
        interior_poly = 0
        if isinstance(myShape, Polygon):
            print "yes, I am a polygon"
            # count how many points for each interior polygon
            try:
                interior_poly += len(myShape.interior.coords)
            except:
                pass
            # count how many points for each exterior polygon
            exterior_poly += len(myShape.exterior.coords)
            geomArray = asarray(myShape.exterior)
            print geomArray
            print "number of interior points in polygon " + str(interior_poly)
            print "number of exterior points in polygon " + str(exterior_poly)
        elif isinstance(myShape, MultiPolygon):
            print "yes, I am a MultiPolygon"
            # count how many points for each interior polygon
            try:
                interior_poly += len(myShape.interior.coords)
            except:
                pass
            try:
                # count how many points for each exterior polygon
                exterior_poly += len(myShape.exterior.coords)
            except:
                pass
            try:
                geomArray = asarray(myShape.interior)
            except:
                pass
            try:
                geomArray = asarray(myShape.exterior)
            except:
                pass
            print geomArray
            print "number of interior points in polygon " + str(interior_poly)
            print "number of exterior points in polygon " + str(exterior_poly)

Wende

Answer 1

Interior and exterior rings are structured differently. For any polygon, there is always 1 exterior ring with zero or more interior rings.

So looking at the structure of a geometry, exterior is a LinearRing object, and interiors is a list of zero or more LinearRing objects. Any LinearRing object will have coords, which you can slice to see a list of the coordinates with coords[:].

The following is a function that returns a dict of lists of exterior and interior coordinates:

def extract_poly_coords(geom):
    if geom.type == 'Polygon':
        exterior_coords = geom.exterior.coords[:]
        interior_coords = []
        for interior in geom.interiors:
            interior_coords += interior.coords[:]
    elif geom.type == 'MultiPolygon':
        exterior_coords = []
        interior_coords = []
        for part in geom:
            epc = extract_poly_coords(part)  # Recursive call
            exterior_coords += epc['exterior_coords']
            interior_coords += epc['interior_coords']
    else:
        raise ValueError('Unhandled geometry type: ' + repr(geom.type))
    return {'exterior_coords': exterior_coords,
            'interior_coords': interior_coords}

E.g.:

extract_poly_coords(myShape)

Lightheaded

Okay...well, I got passed that block...but I can't delete a coordinate (i.e. item) from the list of .coords? I keep getting: >>TypeError: 'CoordinateSequence' object doesn't support item deletion on >> for i in sorted(removeList, reverse=True): del ring.coords[i] – Wende