Polígono interno anillos/agujeros tema

Question

Estoy trabajando en un script para suavizar los polígonos. Esencialmente, separa los anillos interior y exterior/polígonos y enviar las coordenadas para un suavizado def/función. El "alisado" coordenadas, a continuación, se anexa y se escriben en un archivo de forma de usar pyshp.

El problema que me estoy enfrentando es la siguiente: Para algunos polígonos, después de suavizado, no todos los anillos interiores/agujeros son dibujados/inserta. Sin embargo, para similar polígonos complejos en los agujeros son dibujados/inserta. La razón de esto no estoy del todo seguro.

Lo que he intentado hasta ahora para identificar fueron el problema podría ser era para imprimir las coordenadas de entrar en el suavizado de def/función y parece que la falta de anillos interiores son, por alguna razón, no se envía a la suavizado def/función.

A continuación está el guión y las imágenes de ejemplo. Por favor nota que soy bastante nuevo en python y de secuencias de comandos en general.

Alguna sugerencia para solucionar este problema???

EDITAR: Puedo confirmar que la falta de anillos interiores en el alisado shapefile es debido a que no se envía a la función de suavizado.

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()

Sin alisamiento polígono con anillos interiores:

Alisado polígono con la ausencia de anillos interiores:

Answer 1

Después de algunas pruebas adicionales sobre diferentes conjuntos de datos, llegué a la conclusión de que el shapefile de datos parece ser el problema y no tanto de la secuencia de comandos (abierto para el debate). Para los conjuntos de datos sin ningún tipo de islas o complejo anillos interiores/agujeros de la suavizado script funciona muy bien. Sin embargo, al interior de los anillos/agujeros están presentes todo tipo de cosas interesantes suceden a la salida del shapefile, tales como las características que faltan, polígonos donde no debe haber agujeros y así sucesivamente. Los resultados también varían considerablemente entre los diferentes conjuntos de datos con polígonos complejos).