Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
    
What do you get if you just do print singlepolys instead of looping through the list? –  Chad Cooper Dec 13 '12 at 22:46
    
Something I discovered in writing a script that extracted individual extents from multi-part features was that very often, what we think of as single-part features can actually be deemed multi-part where there is only one real "good" part and there is a "bad" part with null geometry. I don't know what causes this, but I did an evaluation on over 700 datasets with 6 million individual features total, and there were many datasets that had nothing but false multi-part features with this issue. I'll dig around and see if I can find anything on it. –  egdetti Dec 13 '12 at 22:48
    
@ChadCooper: I've tried that. I get something that looks correct namely 26 arcpy.Polygon objects with different memory locations. They still have the same area, though. –  mrdr Dec 14 '12 at 22:54
add comment

2 Answers

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...
share|improve this answer
    
+1 Good catch, excellent insight and explanation. –  Chad Cooper Dec 14 '12 at 0:59
    
Thanks for your answer, @dmahr! I also believe it has something to do with the kind of object being passed around. Problem 1: Neither copy.copy() or copy.deepcopy() works for a Geometry/Polygon object. Problem 2: If I loop over a nested list (mutable objects as you mention) instead of the rows things work perfectly anyway. I don't think it's as simple as just a problem with mutable vs immutable objects. Thank you! –  mrdr Dec 14 '12 at 22:52
    
Ok, copy.copy() actually works but gives the same problem as I mentioned in my original question. –  mrdr Dec 17 '12 at 10:56
add comment

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.

share|improve this answer
    
Very strange, but good info. –  egdetti Dec 13 '12 at 23:24
    
@ChadCooper: I've constructed a work-around but haven't solved the original problem. Just like you, I've been able to extract things like the area to a list and that works. But I want the real Geometry objects because I want to input them to the arcpy.Buffer_analysis(...) tool (it's possible to send in individual ArcPy geometry objects to various tools instead of shapefiles). –  mrdr Dec 14 '12 at 22:58
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.