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.

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. Eh!

I am a beginner to Python scripting and I know that this question is almost spam, but I appreciate any advice or tips like how I can 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:

enter image description 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")
share|improve this question

closed as too broad by PolyGeo, Vince, Devdatta Tengshe, Paul, whuber Aug 5 at 12:15

There are either too many possible answers, or good answers would be too long for this format. Please add details to narrow the answer set or to isolate an issue that can be answered in a few paragraphs.If this question can be reworded to fit the rules in the help center, please edit the question.

2  
Something is the matter with the procedure used to compute the flow lines if it allows them to cross. Why not use a better procedure at the outset and avoid this post-processing altogether? –  whuber Jul 30 at 17:32
3  
You seem to be asking for a code review which I think makes your question too broad for our focussed Q&A format. I think you should identify any bottlenecks by using more timing code and then look to ask about alternatives to those as individual code snippets. –  PolyGeo Jul 30 at 21:59
    
We compute the lines from DTM for whole republic. The DTM was in 5x5m resolution and the GRASS was used for it. There were many inaccuracies in DTM, like places with no outflow and digital noise. So the function Fill was used to generalised DTM and fill zones with no outflow. Now we will try to compute lines again with better DTM, but I think this kind of errors (like crossing and so on - this topic is only about crossing, but I try to solve many others problems) will still appears. –  david_p Jul 31 at 8:22
1  
As Whuber sugested, try to get cleaner results. I am not sure how did you get your results above, but it does not seem like straight product of ie. Hydrology tools in Arc. Lines does not connects etc. There is nothing wrong with the script since it is doing the job, but some Python developers would say that FOR(For(for())) is not the best way to do things :) Did you consider raster operation? I believe that you area must be massive, so you can always split raster in to tiles and process them in parallel mode. –  Tomek Jul 31 at 8:30
1  
For cleaner results from DEM you may try to use something like blogs.esri.com/esri/arcgis/2013/03/05/… or lago-consulting.com/fill_sinks_plus.html. Never tested myself but you may give it a shout and share some experiences. –  Tomek Jul 31 at 8:57

4 Answers 4

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)
share|improve this answer
    
Thank you, I'll try to use arcpy.da cursors and add the indexes to my code. I hope that will make the script faster a bit at least. But one question.. If I use AddIndex to field eg. ID, then I use SearchCursor on filed original field ID or on IndexName? –  david_p Jul 31 at 7:30
1  
What ever field you are searching against apply the index to. –  artwork21 Jul 31 at 10:41

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!

share|improve this answer
1  
+1 I was going to mention the multiprocessing module. It wouldn't be hard to break the feature class into (more or less) equally sized partitions that you feed each to a pool of workers. This would provide the most speedup. It just requires more legwork because you'd need to merge the results together at the end. –  Roman Jul 30 at 18:08
    
Thanks for the advices, I should learn more about subprocesses and multiprocessing. Its quite new to me, but I think that would help. I have lot of computing potencial here (16GB RAM, 8 cores Intel i7 3,4 GHz), but when i run the script, there is lots of potencial still unused. So its possible with subprocesses or multiprocessing to "split" the tasks and use all of computing potencial? Also I can use WM on servers.. –  david_p Jul 31 at 7:28
1  
Exactly! Multiprocessing has two benefits. First, you can make more use of your machine's resources (8 cores is great so I recommend having a pool of seven workers - more might be OK but will depend on how intense you application is). The second benefit for long running processes is that there are sometimes memory leaks in compiled code that ArcPy calls. Some where notorious but seem to have been fixed, but I believe there are others which still exist. Each time you close a subprocess, the memory leak (if one exists) is cleared. For long-running processes like yours, this is good. –  MappaGnosis Jul 31 at 7:52
    
I've added a link to the python documentation for multiprocessing. –  MappaGnosis Jul 31 at 7:55
    
I have idea that I could manage the script with multiprocessing and use subprocess to - for calling functions from outside the script. I could make it faster too, no? If I get it right, it possible to define some stand-alone functions (like erasing, merging eg. like in my script) and call them from the "main" script? Now I try to learn how to apply multiprocessing, cause I never try it. –  david_p Jul 31 at 8:16

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)]
share|improve this answer
1  
Depending on the size of the search cursor, he or she may want to use generators instead of list comprehensions, that way the whole thing doesn't have to be loaded into memory at once. –  MakinFlippyFloppy Jul 30 at 15:56
    
Out of curiosity, would that be really faster? –  Tomek Jul 31 at 8:18
    
I have found them to be faster. –  dklassen Jul 31 at 14:59

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.

share|improve this answer
5  
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. ha ha ha WHAT? Where did you hear this? The old cursors were never "broken," they were just "slow" and "the APIs didn't feel right to a professional Python developer." We couldn't break backward compatibility for existing scripts so we left the old cursors untouched and offered the arcpy.da ones in addition to what was already there. –  Jason Scheirer Jul 30 at 14:43
1  
they broke in the fact that after applying service pack 2 my script went from taking 18 hours to process to over 12 days. It may not have generated an error but in my book that counts as broken. –  Lemur Jul 30 at 14:56
1  
Did you contact support? Was it on a specific dataset? Again, they went untouched from version to version. There must be something else going on in your case. –  Jason Scheirer Jul 30 at 14:57
3  
I did not contact support for that particular case because at the time I was having issues getting them to fix other things they broke and decided to call "features". Even getting our sales rep to help did not work. Those "features" screwed up all of our other production and I was wasting time trying to fix our production workflows. It was easier to keep a single instance of sp1 running just for this process rather than trying to deal with ESRI support. –  Lemur Jul 30 at 15:05

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