Optimising and speeding up ArcPy code? [closed]

Question

Again my lines of surface runoff. I made a script, solving the issue when lines cross other lines. Something like generalization. Like on the picture below. The problem is that I have to solve 30 million. of these lines. I tested it on 20 000 lines and time needed to run was 2h 30 min. So approximately, if I will do it for whole 30 million lines, the time needed to run is 5 - 6 months.

I am a beginner to Python.

How can I adjust my script to run it faster?

I'm using Python 2.7.3, ArcGIS 10.2 and Windows 7. I put my script here:

'

#IMPORT MODULES, SETTINGS    
# ##############################################################################
import arcpy, datetime
start = datetime.datetime.now()
arcpy.env.overwriteOutput = True
arcpy.env.workspace = r"C:\Users\david\Desktop\GAEC\TEST\TEST1.gdb"

# DEFINITION OF VARIABLES
# ##############################################################################
OL = "OL"                             # LAYER WITH LINES
OL_FC = "in_memory" + "\\" + "OL_FC"
OL_buff = "in_memory" + "\\" + "OL_buff"
memOne = "in_memory" + "\\" + "OL_one"
memSelect = "in_memory" + "\\" + "memSelect"
memErase = "in_memory" + "\\" + "OL_erase"
memOut = "in_memory" + "\\" + "OL_out"
memOutTemp = "in_memory" + "\\" + "OL_outTemp"
memDissolve = "in_memory" + "\\" + "memDissolve"
memSplit = "in_memory" + "\\" + "memSplit"
memUnsplit = "in_memory" + "\\" + "memUnsplit"
OL_final = "OL_final"

# TESTING AND CREATING FIELDS IN ATRB TABLE, CALCULATING LENGTH
# ##############################################################################
if len(arcpy.ListFields(OL, "LENGTH")) > 0:
    print "Table already has a field LENGTH!..."
else:
    arcpy.AddField_management(OL, "LENGTH", "LONG")
    print "Missinf field LENGTH -> creating..."
arcpy.CalculateField_management(OL, "LENGTH",  "!shape.length!", "PYTHON")
print "Calculate length for lines..."

if len(arcpy.ListFields(OL, "EDIT")) > 0:
    print "Table already has a field EDIT!..."
    arcpy.DeleteField_management(OL, "EDIT")
    arcpy.AddField_management(OL, "EDIT", "SHORT")
else:
    arcpy.AddField_management(OL, "EDIT", "SHORT")
    print "Missinf field EDIT -> creating..."

arcpy.MakeFeatureLayer_management(OL, OL_FC)
rows = arcpy.SearchCursor(OL_FC)
sql = '"LENGTH" < ' + str(10)
lengths = list()
lengths0 = list()
i = 1
a = -1

# CYCLE WHICH MAKE LIST OF LENGHTS, USED FOR NUMBER OF LINES (CYCLES)
# ##############################################################################
for row in rows:
    x = row.getValue("LENGTH")
    lengths.append(x)

# WHOLE PROGRAMM CYCLE
# ##############################################################################
for q in (lengths):
    # COUNTING LENGTHS FOR LINES, NEED TO BE UPDATED CAUSE LENGHTS ARE CHANGING
    rows0 = arcpy.SearchCursor(OL_FC)
    for row0 in rows0:
        x0 = row0.getValue("LENGTH")
        lengths0.append(x0)
        lengths0.sort()
    # ##########################################################################
    leng = lengths0[a]
    # SELECTION OF ACTUAL LONGEST LINE, STARTS FROM END OF THE LIST AND DESCENDING
    sql1 = '"LENGTH" = ' + str(leng)
    arcpy.SelectLayerByAttribute_management(OL_FC, "NEW_SELECTION", sql1)
    # ##########################################################################
    curA = arcpy.SearchCursor(OL_FC)
    for rowA in curA:
        # TESTING IF LINE WAS EDITED BEFORE (IF YES, FIELD 'EDIT' = 1, THEN SKIP)
        test0 = rowA.getValue("EDIT")
        if test0 != 1:
        # ######################################################################
            arcpy.MakeFeatureLayer_management(OL_FC, memOne)
            print "Line exported into independently layer"
            arcpy.Buffer_analysis(memOne, OL_buff, "1.3 Meters")
            print "Created buffer around the line"
            arcpy.SelectLayerByLocation_management(OL_FC, "INTERSECT", OL_buff, "", "NEW_SELECTION")
            curAA = arcpy.UpdateCursor(OL_FC)
            for w2 in curAA:
                # FILLING FIELD EDIT WITH 1, SETTING FIELD LENGTH TO 0
                # ##############################################################
                w2.setValue("EDIT", 1)
                w2.setValue("LENGTH", 0)
                curAA.updateRow(w2)
            arcpy.MakeFeatureLayer_management(OL_FC, memSelect)
            print "Exporting affected lines"
            test = int(arcpy.GetCount_management(OL_FC)[0])
            # TESTING IF LINE TOUCHES AT LEAST ONE ANOTHER LINE, IF NO, SKIP
            # ##################################################################
            if test != 1:
                arcpy.Dissolve_management(memSelect, memDissolve)
                arcpy.Erase_analysis(memDissolve, OL_buff, memErase)
                arcpy.SplitLine_management(memErase, memSplit)
                arcpy.UnsplitLine_management(memSplit, memUnsplit)
                arcpy.Merge_management([memUnsplit, memJedna], memOut)
                print "Merge with original line (which buffer erases)"
                arcpy.Snap_edit(memOut, [[memOne, "EDGE", "1.4 Meters"]])
                print "LINES ARE SNAPPED TO ORIGINAL LINES"
                # DELETING LINES THAT ARE SHORTER THAN 10 M
                # ##############################################################
                arcpy.MakeFeatureLayer_management(memOut, memOutTemp)
                arcpy.CalculateField_management(memOutTemp, "LENGTH",  "!shape.length!", "PYTHON")
                arcpy.SelectLayerByAttribute_management(memOutTemp, "NEW_SELECTION", sql)
                arcpy.DeleteFeatures_management(memOutTemp)
                print "Lines shorter than 10 m deleted"
                # FIRST CYCLE CREATES OL_finan, AFTER LINES ARE APPENDING TO FINAL
                # ##############################################################
                if i == 1:
                    arcpy.CopyFeatures_management(memOutTemp, OL_final)
                else:
                    arcpy.Append_management(memOutTemp, OL_final)
                    print "New layer is append to layer before"
            # LINES WITHOUT CROSSING ANY OTHER LINE
            # ##################################################################
            else:
                for w2 in curAA:
                    w2.setValue("EDIT", 1)
                    w2.setValue("LENGTH", 0)
                    arcpy.Append_management(OL_FC, OL_final)
                    curAA.updateRow(w2)
                print "Standalone line appends to final layer"
            print "Imputed i, i = " + str(i) + "..."
    # DELETE RAM MEMORY
    # ##############################################################
    arcpy.Delete_management("in_memory")
    # IMPUTED TO NEXT CYCLE
    # ##########################################################################
    i+=1
    a-=1

# DELETE ASSISTANT FILEDS
# ##############################################################
arcpy.DeleteField_management(memOutTemp, "LENGTH")
arcpy.DeleteField_management(memOutTemp, "EDIT")

# ENDING OF TIMER
# ##############################################################################
end = datetime.datetime.now()
print "\nStart: " + start.strftime("%Y-%m-%d %H:%M:%S")
print "End of script: " + end.strftime("%Y-%m-%d %H:%M:%S")

Answer 1

In addition to using the new arcpy.da cursor, I would also suggest:

• You have many different search cursors on the same layer, see if you can eliminate some of those and pull your attributes from one or two
• Apply an Add Attribute Index on any column that you are querying against
• See if you can remove the select layer by attribute logic and apply that query in the where clause argument in the cursor (this will eliminate many records that you do not need to iterate through)

Answer 2

In addition to the answers about changing your cursor and reducing the number of them, you can also parallelise this a little by using subprocesses. Subprocesses in Python get around the limitation of the Global Interpreter Lock (GIL) - as opposed to multithreading in Python (a different thing to multiprocessing). Although ArcGIS 10 uses multiprocessing now, I believe you must still explicitly tell ArcPy to do so and create your subprocesses accordingly. The number you create is a balancing act though. Start up too many and you will actually slow things down. For heavyweight geoprocessing tasks I work on a rule of thumb of number of subprocesses = cores -1 (leaving one core for ArcGIS and the operating system and the rest can be hogged by your process).

You can also create a geoprocessing farm and press many machines into service. There is a way of doing this with a single copy of ArcGIS but it enters the 'grey area' of the software licence agreement. I have done it in the past with ESRI's knowledge and interest (burning up approximately 50 cores across about 20 machines) but you need to be careful about where you stand legally!.

Your alternative would be to recreate your process in FOSS4G software and press as many machines on your network as you can into service, farming the task out across them. Given the time scales you are talking about, you may need to consider this option if the changes suggested plus a few local subprocesses doesn't give you the speed you need.

EDIT For documentation on multiprocessing in Python see here. Just be careful in case there are any changes between v2.x and 3.x!

Answer 3

Switch over to arcpy.da for all of your cursors. The old arcpy cursor objects broke at 10.1 sp2 and the only way ESRI could fix it was to create the new arcpy.da cursor objects.

I had a process that ran in 10.1 sp1 that took 18 hours using the old cursor objects. After sp2 it took over 12 days before I killed it. The new arcpy.da cursors takes 4 hours.

Answer 4

In addition to the above suggestions, I would also suggest to use list comprehensions for simple loops. For example:

rows = arcpy.SearchCursor(OL_FC)
for row in rows:
    x = row.getValue("LENGTH")
    lengths.append(x)

Change to:

lengths = [row.getValue("LENGTH") for row in arcpy.SearchCursor(OL_FC)]