Polygon inner rings/holes issue

Question

I posted the same question on GIS SE 8 days ago, but have had no response and have decided to broaden the knowledge base.

I am working on a script to smooth polygons. Essentially it, separates the inner and outer rings/polygons and send the coordinates to a smoothing def/function. The "smoothed" coordinates are then appended and written to a shapefile using pyshp.

The issue I am facing is as follows: For some polygons, after smoothing, not all the inner rings/holes are drawn/inserted. However, for similar complex polygons the holes are drawn/inserted. The reason for this I am not entirely sure.

What I have tried thus far to identify were the issue might be was to print the coordinates going into the smoothing def/function and it seems that the missing inner rings are, for some reason, not sent to the smoothing def/function.

Below is the script and example images. Please note I am fairly new to python and scripting in general.

Any suggestions to solve this issue???

EDIT: I can confirm that the missing inner rings in the smoothed shapefile is due to it not being sent to the smoothing function.

import os
import sys
import ogr
import fileinput
import shapefile
import Tkinter
import tkSimpleDialog
import gdalnumeric
import tkFileDialog
import numpy as np
from scipy.ndimage import gaussian_filter
import shapely
import time
import shutil


root = Tkinter.Tk()

#Smoothing formula for outer rings
def smoothingouter(pts):
    for i in pts:

        array = np.array(pts)
        x, y = array.T
        t = np.linspace(0, 1, len(x))
        t2 = np.linspace(0, 1, 100)

        x2 = np.interp(t2, t, x)
        y2 = np.interp(t2, t, y)

        sigma = 1.25
        x3 = gaussian_filter(x2, sigma)
        y3 = gaussian_filter(y2, sigma)

        x4 = np.interp(t, t2, x3)
        y4 = np.interp(t, t2, y3)

        zipped = zip(x4,y4)
        return zipped

#Smoothing formula for inner rings       
def smoothinginner(pts):
    for i in pts:

        array = np.array(pts)
        x, y = array.T
        t = np.linspace(0, 1, len(x))
        t2 = np.linspace(0, 1, 100)

        x2 = np.interp(t2, t, x)
        y2 = np.interp(t2, t, y)

        sigma = 0.2
        x3 = gaussian_filter(x2, sigma)
        y3 = gaussian_filter(y2, sigma)

        x4 = np.interp(t, t2, x3)
        y4 = np.interp(t, t2, y3)

        zipped = zip(x4,y4)
        return zipped

#Select input shapefile
inshp = tkFileDialog.askopenfilename(title='Select shapefile to smooth:')
mdir, mfilename = os.path.split(inshp)
mnam, mext = os.path.splitext(mfilename)

sf = shapefile.Reader(inshp)
shapeout = mdir + '/' + mnam + '_smoothed.shp'
w = shapefile.Writer(shapefile.POLYGON)
w.autoBalance = 1
w.field("ID")
id=8 

shape = sf.shapes()

# a: feature index tracker
a = 0

#Steps and cnt for the progress bar
steps = len(shape)/10
cnt = 0

#for each feature in layer 
for i in shape:
    b = shape[a].points
    # add feature points to numpy array
    d = np.array(b)
    exterior = True
    #q is the index of the first point for each ring (includes both inner and outer rings)
    q = 0
    #v is the index tracker of each point in a feature regardless whether it is part of inner or outer ring
    v = 0
    #create a list to add smoothed points to
    wparts = []
    #loops through each point of a polygon
    for c in b:
        #get the first point (this is the first point of the exterior polygon since all points of the outer ring are first in the index, all inner rings thereafter)
        first = b[q]
        #compares each point to the first point. Every ring will have identical coords for its first and last points. From this we can seperate each rings coords      
        if c == first:
            # q must be smaller than v. E.g For the first coord, q = 0 and v = 0 therefore d[q] and d[v] are the same point and no range will be found
            if (q < v):
                #All points from the exterior ring of a polygon gets processed using the smoothingouter function
                if exterior == True:
                    ss = smoothingouter(d[q:v])
                #Points from each inner ring of a polygon gets processed using the smoothinginner function 
                else:
                    ss = smoothinginner(d[q:v])
                #Add the resulting coordinates to the list
                wparts.append(ss)

                if (exterior == True):
                    w.record(id)
                #Update the starting index of the next exterior polygon
                q = v+1
                exterior = False



        v+=1
    w.poly(wparts)

    a+=1

#Create new smoothed poly shapefile
w.save(shapeout)

prjcopy = shapeout[:-3] + 'prj'

shutil.copy(inshp[:-3] + 'prj', prjcopy)

print 'done'

root.destroy()

Unsmoothed polygon with inner rings:

Smoothed polygon with missing inner rings:

Answer 1

I would recommend using shapely to iterate over your exterior and interior polygons.

When you have a polygon which contains all your exterior and interior polygons:

let's say your polygon is called R,

R.exterior.coords , gives you an object with your exterior coordinates

R.interiors, gives you an object with an interior ring sequence with all the interior rings, you can slice it or iterate over it.

therefore if you want the first interior ring,

r.interiors[0] if you want the second interior ring, r.interiors[1] if you want the coordinates r.interiors[1].coords

In this way you can include every interior polygon in your operations .