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 have a layer with catchment polygons. For those that are smaller than a certain size, I want to merge with an adjoining polygon (namely the one downstream). In the attribute table I have a GridID, which i use to identify different polygons, and I can identify the downstream polygon by that attribute.

How can I, in this case inside a SearchCursor, merge these two polygons together?

I have tried with an UpdateCursor and arcpy.Dissolve_management() without getting it to work. Since Dissolve requires two tables with a value in common, I suppose I could just add a new field and make sure that this is the case for each iteration. But I can't use the same file as input and output (obviously), and I would like to get the file updated right away. If not, I might merge polygons that are of sufficient size since it has been merged with another earlier in the process, and after that meet the criterion. In the case of an UpdateCursor, I think this is mainly for updating field values and not the geometry itself?

Is there a (more sufficient) method for this?

# inpgs = Catchment polygons shapefile

where1 = '"' + 'Shape_Area' + '" < 5000'
polyrows = arcpy.SearchCursor(inpgs, where1)

for prow in polyrows:
    grd1 = prow.GridID

    ...

    # grd1 = GridID on the current polygon
    # grd2 = GridID on the polygon downstream

    where4 = '"' + 'GridID' + '" = ' + str(grd2)
    polyrows2 = arcpy.UpdateCursor(inpgs, where4)

    for prow2 in polyrows2:
        arcpy.Dissolve_management(inpgs, outpgs, ???)
        inpgs = outpgs # ?

I use ArcInfo 10.0, Python 2.6 on Windows 7.

EDIT: I have now implemented the changes dmahr suggested, and made some more myself. That resulted in the following code:

# Create dictionary for all geometries
polyrows = arcpy.SearchCursor(inpgs)
geometryDictionary = {}

for prow in polyrows:
    geometryDictionary[prow.GridID] = prow.Shape
del prow, polyrows

# Create dictionary (nDD) which contains the GridIDs of the downstream catchment area 
# as value and the current catchment area's GridID as key.

nextrows = arcpy.SearchCursor(indrln)
nextDownDictionary = {}
nDD = {}

for nrow in nextrows:
    nextd = nrow.NextDownID
    nextDownDictionary[nrow.GridId] = nrow.NextDownID

for k in dict.keys((nextDownDictionary)):
    clause = '"' + 'HydroID' + '" = ' + str(nextDownDictionary[k])
    nr = arcpy.SearchCursor(indrln, clause)
    for n in nr:
        nDD[k] = n.GridID
del n, nr

where1 = '"' + 'Shape_Area' + '" < 5000'
polyrows = arcpy.SearchCursor(inpgs, where1)
arcpy.Copy_management(inpgs, outpgs)

# Merge the small polygon into the downstream polygon
for prow in polyrows:
    grd1 = prow.GridID
    grd2 = nDD.get(grd1)

    if grd2:
        # grd1 = GridID on the current polygon
        # grd2 = GridID on the polygon downstream

        geometry1 = prow.Shape
        geometry2 = geometryDictionary[grd2]

        arcpy.Merge_management([geometry1, geometry2], tempMerged)
        arcpy.Dissolve_management(tempMerged, tempPgs)
        arcpy.Append_management(tempPgs, outpgs, "NO_TEST")
        grids.append(prow.GridID)


# After appending the dissolved polygons to the output file, the row for the merged polygon is being deleted
for g in grids:
    whereU = '"GridID" = ' + str(g)
    upd = arcpy.UpdateCursor(outpgs, whereU)
    for up in upd:
        upd.deleteRow(up)


del prow, polyrows
del up, upd
share|improve this question
    
Do you have to use a python script to accomplish this? What about a manual method using selections/definition queries/field calculator? –  Baltok Dec 11 '12 at 15:14
    
Yes, since I want to execute this as a script tool in Arcmap. The tool is not for me personally, but to be used by people who know nothing about programming and want this as an automated task. If that was not the case, I could do as you suggests. –  Martin Dec 11 '12 at 15:26
add comment

1 Answer

up vote 8 down vote accepted

There are a few reasons why this isn't working for you:

  1. Dissolve only works with a single input. You need to use a tool like Merge in order to dissolve two polygons.

  2. ArcGIS geoprocessing tools like Dissolve can run using arcpy geometries as inputs. However, your code does not pass arcpy geometry objects to the dissolve tool, it is entire layers/shapefiles. Remember that cursors are not definition queries applied to datasets that carry through to other tools, they are just an object used for accessing rows. Rows are just what they sound like: a row in the attribute table. The geometry is one attribute of the row, usually with the field name SHAPE (though you can get the field name using from the ShapeFieldName property of the arcpy.Describe(input) object).

  3. You should avoid using nested cursors because cursors, especially updateCursor, lock the input file to prevent other editing.

Here's how I would code it:

# inpgs = Catchment polygons shapefile
# tempMerged = temporary merged polygon shapefile (better to have this in-memory)
where1 = '"' + 'Shape_Area' + '" < 5000'
shapeField = arcpy.Describe(inpgs).ShapeFieldName

#Read polygon geometry into dictionary; key = GridID, value = geometry
geometryDictionary = {}
polyrows = arcpy.SearchCursor(inpgs)
for prow in polyrows:
    geometryDictionary[prow.GridID] = prow.shapeField
del prow
del polyrows

#Loop through catchments that are to be dissolved
polyrows = arcpy.SearchCursor(inpgs, where1)
for prow in polyrows:

    #Get geometries of two polygons
    # grd1 = GridID on the current polygon
    # grd2 = GridID on the polygon downstream
    geometry1 = prow.shapeField
    geometry2 = geometryDictionary[grd2]

    #Merge datasets
    arcpy.Merge_management([geometry1, geometry2], tempMerged)

    #Dissolve
    arcpy.Dissolve_management(tempMerged, outpgs)

A part I left out is how you'll have to delete the polygons that were used in merges. There are a number of different ways to code this.

Also, keep in mind that this code does not handle the edge case possibility that a polygon with area < 5000 is upstream of another polygon with area < 5000. If you happened to merged the downstream one first, you are going to get strange results.

EDIT: One other possibility is if you can get an upgrade to ArcGIS 10.1. Geometry objects in 10.1 (like Polygon) have a Union method that does a true geometric Union of two geomtries (shown below in blue). This is not to be confused with what the Union geoprocessing tool does. ArcGIS 10.1 union method

share|improve this answer
    
Thank you dmahr, i will try it out! –  Martin Dec 11 '12 at 15:27
    
+1 Nice tip on the union of geometries in 10.1 –  om_henners Dec 12 '12 at 0:52
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.