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 shapefile consisting of around 4000 polygons. The polygons are created from buffering around a similar number of islands somewhere in the Norwegian archipelago. Due to an error in the dissolve functionality using arcpy dissolve fails. I am assuming this is because arcpy.Dissolve_management can not handle too many overlapping features.

To try to work around this I plan on splitting my initial file into ten separate files consisting of around 400 features which are then dissolved separately. These ten files are then merged together and a final dissolve is performed.

How can I do such a selection? The arcpy.select_analysis allows for a where clause. How can I use this to select the n first features, then the n+1 to n*2 features, then the n*2+1 to n*3 features etc?

My final goal is to make a small arcpy function which could be a preliminary fix for us mere mortals relying on ESRIs bug fixing.

share|improve this question
    
What exactly is the error you encounter? Have you considered simplifying the geometry? –  blah238 Nov 26 '12 at 6:11
    
Simplified the geometry. Didnt help. It seems like a certain number of overlapping features leads to this error. Norwegian customer contact at geodata.no referred to it as: "NIM079373: Running a large number of features through the Dissolve or Buffer with dissolve option, hangs during process" Supposed to be fixed for 10.1 sp1, but they failed to do so. –  ragnvald Nov 26 '12 at 14:52
    
I have failed to find the referred to error code while googling it. Puzzled by this. It could be an ESRI internal reference. –  ragnvald Nov 26 '12 at 15:02

5 Answers 5

up vote 5 down vote accepted

I ended up making the function suggested below. In part inspired by inputs from other contributors on this question. The object handling is coarse and it creates a lot of temporary files.

I am sure this function could be made a lot better. If properly polished it could be a decent preliminary fix for the dissolve bug by ESRI. I am making it a community wiki.

# Import modules
import arcpy, os

def shape_dissolve(file_in, file_out,group_by):  
    count           =  0
    current_min     =  0
    current_max     =  group_by
    joinstring      =  []

    features_total    =  int(arcpy.GetCount_management(file_in).getOutput(0))

    while (current_max < (features_total+group_by)):

        resulting_file  = "%sbuffer_result_%s.shp" % (path_maps_basis,count)

        where_clause    = '"FID"> %s AND "FID" <= %s' %(current_min,current_max)

        arcpy.Select_analysis(file_in, resulting_file, where_clause)

        resulting_file_d = "%sbuffer_result_d_%s.shp" % (path_maps_basis,count)

        arcpy.Dissolve_management(resulting_file, resulting_file_d,"","","SINGLE_PART","")

        # delete temporary files
        arcpy.Delete_management(resulting_file)

        joinstring.append(resulting_file_d)

        current_min = group_by*count
        current_max = current_min + group_by

        count +=1

    resultbuffer = "%sresultbuffer.shp" % (path_maps_basis)

    arcpy.Merge_management(joinstring, resultbuffer,"")

    # delete temporary files
    for shapefile in joinstring:

        arcpy.Delete_management(shapefile)

    arcpy.Dissolve_management(resultbuffer, file_out,"","","SINGLE_PART","")

    # delete temporary files
    arcpy.Delete_management(resultbuffer)

#Calling the function with these sample values
group_by = 700
path_maps_basis     = "C:/"
multi_feature_file  = "%sbuffer_complex.shp"    % (path_maps_basis)
resulting_file      = "%sbuffer_dissolved.shp"  % (path_maps_basis)


shape_dissolve(multi_feature_file,resulting_file,group_by)
share|improve this answer
1  
Just a note, this sequential grouping logic cannot be easily adapted to feature classes in a geodatabase because their OBJECTID values are not necessarily sequential (they can have gaps, and the highest OBJECTID can be greater than the feature count). The FID field in a shapefile is guaranteed to be sequential because it is generated on the fly and not stored in the shapefile: support.esri.com/es/knowledgebase/techarticles/detail/37480 Hence the slightly more complex partitioning logic used in my answer, but if you are only concerned with shapefiles then this should work fine. –  blah238 Nov 26 '12 at 18:36
    
I am doing this one working with shapefiles, so it works for me. Not the ideal thing to do considering the vast options available in PostGIS. But quick and dirty does it this time :-) –  ragnvald Nov 26 '12 at 20:22

You can adapt the partitioning methodology used in this script for your purposes: Port “Producing Building Shadows” Avenue code to ArcGIS 10. Once you've got it working you could even experiment with multiprocessing by adapting other parts of the script (its purpose was to improve performance using multiple CPU cores).

The relevant function:

def getOidRanges(inputFC, oidfield, count):
    oidranges = []
    if procfeaturelimit > 0:
        message("Partitioning row ID ranges ...")
        rows = arcpy.SearchCursor(inputFC, "", "", oidfield, "%s A" % oidfield)
        minoid = -1
        maxoid = -1
        for r, row in enumerate(rows):
            interval = r % procfeaturelimit
            if minoid < 0 and (interval == 0 or r == count - 1):
                minoid = row.getValue(oidfield)
            if maxoid < 0 and (interval == procfeaturelimit - 1 or r == count - 1):
                maxoid = row.getValue(oidfield)
            if minoid >= 0 and maxoid >= 0:
                oidranges.append([minoid, maxoid])
                minoid = -1
                maxoid = -1
        del row, rows
    return oidranges

Note: procfeaturelimit is a global variable that specifies the maximum number of features each "partition" should be assigned.

Which is used like this:

count = int(arcpy.GetCount_management(inputFC).getOutput(0))
# Get the appropriately delmited field name for the OID field
oidfielddelimited = arcpy.AddFieldDelimiters(inputFC, oidfield)
oidranges = getOidRanges(inputFC, oidfield, count)
for o, oidrange in enumerate(oidranges):
    # Build a where clause for the given OID range
    whereclause = "%s >= %d AND %s <= %d" % (oidfielddelimited, oidrange[0], oidfielddelimited, oidrange[1])
    # Do your Select_analysis and whatever else here with the given whereclause
share|improve this answer

Use the looping and variable value as shown in Aragon's answer with the FID field (this is the zero-based object ID field for shapefiles) as the select field as shown in L_Holcombe's answer to generate the where clause, and all should be good. To program in the total number of features use the Get Count tool and divide by 10 assigned to the variable. Will also need Make Feature Layer. The lines with the '#' character need to be adjusted to suit your data and the '#' removed when ready to run.

#lyr = arcpy.MakeFeatureLayer_management(specify input file here)
#getcount = int(arcpy.GetCount_management(lyr).getOutput(0))
getcount = 4000     # just an example - remove line and use get count from above
numfiles = 10
numfeatures = getcount / numfiles

for i in range(0, getcount, numfeatures):
    select_features = specify output file here
    #dissolve_features = specify output file here
    j = i + numfeatures
    if i != 0:
        where_clause = '"FID" > ' + str(i) + ' AND "FID" <= ' + str(j)
    else:
        where_clause = '"FID" <= ' + str(j)
    print where_clause
    #arcpy.Select_analysis(lyr, select_features, where_clause)
    #arcpy.Dissolve_management(select_features, dissolve_features, {dissolve_field})
share|improve this answer

Perhaps something simple would work? Assuming "gid" is a column with consecutive integers:

arcpy.Select_analysis("in_features", "out_features", '"gid" < 400')
arcpy.Select_analysis("in_features", "out_features", '"gid" < 800 AND "gid" >= 400')

etc.

I haven't tried this, but it should work.

ArcGIS SQL expression reference: http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00s500000033000000

share|improve this answer

You can use enumarete function in python.

myList = []

n = 99

for i, v in enumarete(myList):
    if i == n:
        print v 'first 100 item' #append your list
    if i > n and i < (n+2)*2:
        print v 

Or you can use this script from SO:

def split_list(alist, wanted_parts=1):
    length = len(alist)
    return [ alist[i*length // wanted_parts: (i+1)*length // wanted_parts] 
             for i in range(wanted_parts) ]

A = [0,1,2,3,4,5,6,7,8,9]

print split_list(A, wanted_parts=1)
print split_list(A, wanted_parts=2)
print split_list(A, wanted_parts=8)

i hope it helps you...

share|improve this answer
    
Thankyou - your code does the handling of the logic nicely. But the retrieval of the n features using arcpy still remains as a challenge. –  ragnvald Nov 25 '12 at 20:35
    
rows = arcpy.UpdateCursor(full_buffer_temp) can be used to iterate through the features. So I need to push the features into a separate file. Close it and start anew. –  ragnvald Nov 25 '12 at 20:40

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.