Identifying internal polygons using Arcpy

This is just a small python script that will identify any polygons that are within other polygons in a shapefile. It has come handy for me for a project I was working on where I had to identify the “internal” polygons. The main trick here is to compare each polygon with the outer boundary of every other

# -*- coding: utf-8 -*-
#!/usr/bin/env python
import arcpy
#-------------------------------------------------------------------------------
# Name:        flag_int_polygons
# Purpose:     Identifies interior polygons
#
#-------------------------------------------------------------------------------
def main():

    output_folder = u"C:\temp"
    inFC=u"C:\temp\polys.shp"

    # Check if the compare shapefile already exists and if so, recreate it
    if arcpy.Exists(output_folder+"/tmpCompare.shp")==True:
        arcpy.Delete_management(output_folder+"/tmpCompare.shp")
    arcpy.CopyFeatures_management(inFC, output_folder+"/tmpCompare.shp")

    try:
        arcpy.AddField_management(output_folder+"/tmpCompare.shp", "LABEL", "TEXT", "", "", 2, "", "", "REQUIRED")
    except:
        pass

     # Get shape column
    shapeName = arcpy.Describe(inFC).shapeFieldName
    shapeName1 = arcpy.Describe(output_folder+"/tmpCompare.shp").shapeFieldName

     # Loop through input features
    rows=arcpy.SearchCursor(inFC)

    # Start loop through input fc
    r=1
    for row in rows:
        partnum = 0
        parcel=row.getValue(shapeName)
        exBoundary=getPolyBoundary(parcel)
        fid=row.getValue("FID")
        rows1=arcpy.UpdateCursor(output_folder+"/tmpCompare.shp")
        for row1 in rows1:
            fid1=row1.getValue("FID")
            parcel1=row1.getValue(shapeName1)
            if fid1 <> fid:
                if parcel1.within(exBoundary)==True:
                    print "Feature " + str(fid1) + " is within feature " + str(fid)
            rows1.updateRow(row1)

        del row1
        del rows1
    del row
    del rows
#-------------------------------------------------------------------------------
# Name:        getPolyBoundary
# Purpose:     Returns the external polygon boundary
#
#-------------------------------------------------------------------------------
def getPolyBoundary(geom):
    partno=0
    pntObj = arcpy.Point()
    arrayObj = arcpy.Array()

    for part in geom:
        pntCounter=0;
        # Loop through points in part
        for pnt in geom.getPart(partno):
            if pnt:
                pntObj.X=pnt.X
                pntObj.Y=pnt.Y
                arrayObj.add(pntObj)
            else:
                break
        pntCounter=pntCounter+1
    bound=arcpy.Polygon(arrayObj)
    #print str(bound)
    return bound
if __name__ == '__main__':
    main()

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s