How to filter out multipart polygons and then buffer?

Question

I've written a small script trying to accomplish the task described in the question title. It basically looks like:

in_fc = '...' # Some shapefile goes here

# Get the fields (column names) of the attribute table
desc = arcpy.Describe(in_fc)
shape_field = desc.ShapeFieldName

# Initialize an empty polygon list
singlepolys = []

# Open the attribute table of the shapefile
rows = arcpy.SearchCursor(in_fc)
for row in rows:
    poly = row.getValue(shape_field)

    # Check if single part (not multipart)
    if not poly.isMultipart:
        singlepolys.append(poly)

# Print areas of single polygons - just for debugging
for p in singlepolys:
    print p.area

# Buffer code would go here...

The final printing of the areas shows that only multiples of the last polygon have been copied to the list singlepolys. I really don't understand why.

Any help appreciated!

Edit

A clumsy work-around is possible (I've found a faster one but this is closer to the original thought). Here goes:

# ... Skipping initialization

# Open the attribute table of the shapefile
rows = arcpy.SearchCursor(in_fc)
for row in rows:
    poly = row.getValue(shape_field)

    # Check if single part (not multipart)
    if not poly.isMultipart:
        # Buffer this polygon only and add to list
        geom = arcpy.Geometry()
        result = arcpy.Buffer_analysis(poly, geom, '1000 Meters')
        singlepoly_bufs.append(result[0])

# Now union the buffered features
out_geom = arcpy.Geometry()
union = arcpy.Union_analysis(singlepoly_bufs, out_geom)

# ... Skipping the end

Answer 1

You have an issue with the poly variable pointing to the wrong object. I believe that this is called the evaluation strategy. Here's a quick summary of what's happening:

• The variable poly is assigned to the value of row.getValue(shape_field)

• A reference to the poly variable is added to the singlepolys list.

• In the next iteration of the for loop, the variable poly is reassigned to point to a new row.getValue(shape_field). But the first element of singlepolys is still just a reference to poly (not a copy of poly), so it now points to the new value.

This kind of behavior is sometimes called "pass-by-reference" since your list is being passed a reference to poly. This may seem contrary to your understanding of looping in Python, where passing values from within a loop to a list does not pass a reference. Well this is because when you are dealing with certain object types that are immutable (e.g. strings, numbers, tuples) then Python actually passes a copy and not a reference. But a Geometry object is a mutable type, so it passes a reference. This StackOverflow question goes into a lot more detail about Python's object model, which technically isn't pass-by-reference or pass-by-value. This blog post is also a good explanation.

Fortunately, there are ways to force Python to pass a copy of an object, using the copy module. Here's how you would do it:

import copy
in_fc = '...' # Some shapefile goes here

# Get the fields (column names) of the attribute table
desc = arcpy.Describe(in_fc)
shape_field = desc.ShapeFieldName

# Initialize an empty polygon list
singlepolys = []

# Open the attribute table of the shapefile
rows = arcpy.SearchCursor(in_fc)
for row in rows:
    poly = row.getValue(shape_field)

    # Check if single part (not multipart)
    if not poly.isMultipart:
        singlepolys.append(copy.deepcopy(poly))

# Print areas of single polygons - just for debugging
for p in singlepolys:
    print p.area

# Buffer code would go here...

Answer 2

Hmm...not too sure why, but it appears that it doesn't like it when you add the geometry to the singlepolys list, it just adds the last one repeatedly. This does work:

desc = arcpy.Describe('KOP_FoV')
shape_field = desc.ShapeFieldName
# Initialize an empty polygon list
singlepolys = []
# Open the attribute table of the shapefile
rows = arcpy.SearchCursor("KOP_FoV")
for row in rows:
    poly = row.getValue(shape_field)
    print poly.isMultipart
    # Check if single part (not multipart)
    if poly.isMultipart is False:         
        print poly.area
        singlepolys.append(poly.area)

# Print areas of single polygons - just for debugging
for p in singlepolys:
       print p 

Gives me this:

288629552.707
288629550.405
288629553.336
288629551.851
288629552.394
288629552.411

288629552.707
288629550.405
288629553.336
288629551.851
288629552.394
288629552.411

So your original code added 288629552.411 to singlepolys 6 times. Not sure what's up with that, but maybe this will help.