1

Following the question I asked in this thread, I wrote a python script designed to compute an OD cost matrix given a road network and origins/destinations locations. I am completely new to python, so I have no doubt my script would need some improvements.

Here are the steps included in the script:

  1. Load a csv file containing the filepaths to the different datasets
  2. For each dataset (a line in the csv), I want to compute the OD cost matrix (this allows batch processing)
  3. I check the size of the Origins and Destinations. If the number of solutions is higher than 10,000,000 I split the work by chopping up the Origins. If the number of solutions is small, then I go on computing the OD matrix.

The most delicate steps are inside the loop when I have a large amount of Origins and Destinations. In this second for loop, I follow these steps:

  1. I calculate the number of pieces I need in order not to have more than 5,000,000 solutions per loop
  2. I select a number of features in the Origins dataset that will be used as origins
  3. I add these Origins in the OD Matrix layer
  4. I solve the OD Matrix
  5. I export to a .csv file the columns of the "Lines" table in which I am interested

In details, here is the code:

# Import modules
import arcpy
import numpy
import multiprocessing
import time
import csv
start = time.time()
# Check out any necessary licenses
arcpy.CheckOutExtension("Network")
# Input arguments
print "Importing data"
Arguments_Array = numpy.loadtxt('C:\Users\ddeltenre\Desktop\data.csv',dtype='string',delimiter=',',skiprows=0)
processes = len(Arguments_Array)
for x in xrange(0,processes):
    # Script arguments
    network = Arguments_Array[x][0]
    Origins = Arguments_Array[x][1]
    Sort_Field_Origins = Arguments_Array[x][2]
    Destinations = Arguments_Array[x][3]
    Sort_Field_Destinations = Arguments_Array[x][4]
    Output_Feature = Arguments_Array[x][5]
    Output_Location = Arguments_Array[x][6] 
    ODMatrix = 'OD_Cost_Matrix' 
    print "Computing OD cost matrix for network "
    print network
    # Benchmarking
    print "Checking benchmarks"
    nOrigins = int(arcpy.GetCount_management(Origins).getOutput(0))
    nDestinations = int(arcpy.GetCount_management(Destinations).getOutput(0))
    if nOrigins*nDestinations > 100000000:
        print "Too many points! The dataset will be split"
        xrow = (10000000//nDestinations)+1
        loops = nOrigins/xrow
        # Process: Make OD Cost Matrix Layer
        arcpy.MakeODCostMatrixLayer_na(network, ODMatrix, "Length", "", "", "Length", "ALLOW_UTURNS", "", "NO_HIERARCHY", "", "NO_LINES", "")
        # Process: Add Destinations
        arcpy.AddLocations_na(ODMatrix, "Destinations", Destinations, "Name ID #", "5000 Meters", Sort_Field_Destinations, "", "MATCH_TO_CLOSEST", "CLEAR", "NO_SNAP", "5000 Meters", "INCLUDE", "")
        print "Destinations added: {} s".format(time.time()-start)
        tmpOrigins = "{}\\tmp_Origins.shp".format(Output_Location)
        print "Benchmarking complete: {} s".format(time.time()-start)
        for x in xrange(0, loops):
            loopstart = time.time()
            print "Processing... {} / {}".format(x+1,loops)
            tmpOutput_Feature = "{}_{}{}".format(x*xrow,(x+1)*xrow,Output_Feature)
            #pool = multiprocessing.Pool(4)
            # Select part of the origins
            query = '"OBJECTID" > {} AND "OBJECTID" <= {}'.format(x*xrow,(x+1)*xrow)
            arcpy.Select_analysis(Origins, tmpOrigins, query)
            # Process: Add Origins
            originsstart = time.time()
            arcpy.AddLocations_na(ODMatrix, "Origins", tmpOrigins, "Name ID #", "5000 Meters", Sort_Field_Origins, "", "MATCH_TO_CLOSEST", "CLEAR", "NO_SNAP", "5000 Meters", "INCLUDE", "")
            print "Origins loaded: {} s".format(time.time()-originsstart)
            # Process: Solve
            solvestart = time.time()
            arcpy.Solve_na(ODMatrix, "SKIP", "TERMINATE")
            print "Net solved: {} s".format(time.time()-solvestart)
            # Process: Feature Class to Feature Class
            exportstart = time.time()
            Selection = arcpy.SelectData_management(ODMatrix, "Lines")
            fields = arcpy.ListFields(Selection)
            field_names = [field.name for field in fields]
            field_names = [field_names[2], field_names[6]]
            with open("{}\\{}".format(Output_Location,tmpOutput_Feature),'wb') as f:
                dw = csv.DictWriter(f,field_names)
                #--write all field names to the output file
                dw.writeheader()
                #--now we make the search cursor that will iterate through the rows of the table
                with arcpy.da.SearchCursor(Selection,field_names) as cursor:
                    for row in cursor:
                        dw.writerow(dict(zip(field_names,row)))
            print "Data exported: {} s".format(time.time()-exportstart)
            #pool.close()
            #Delete useless select points
            arcpy.Delete_management(tmpOrigins)
            print "Time in loop: {} s".format(time.time()-loopstart)            
    else:
        # Process: Make OD Cost Matrix Layer
        arcpy.MakeODCostMatrixLayer_na(network, ODMatrix, "Length", "", "", "Length", "ALLOW_UTURNS", "", "NO_HIERARCHY", "", "NO_LINES", "")
        # Process: Add Origins
        arcpy.AddLocations_na(ODMatrix, "Origins", Origins, "Name ID #", "5000 Meters", Sort_Field_Origins, "", "MATCH_TO_CLOSEST", "CLEAR", "NO_SNAP", "5000 Meters", "INCLUDE", "")
        # Process: Add Destinations
        arcpy.AddLocations_na(ODMatrix, "Destinations", Destinations, "Name ID #", "5000 Meters", Sort_Field_Destinations, "", "MATCH_TO_CLOSEST", "CLEAR", "NO_SNAP", "5000 Meters", "INCLUDE", "")
        # Process: Solve
        arcpy.Solve_na(ODMatrix, "SKIP", "TERMINATE")
        # Process: Feature Class to Feature Class
        Selection = arcpy.SelectData_management(ODMatrix, "Lines")
        fields = arcpy.ListFields(Selection)
        field_names = [field.name for field in fields]
        field_names = [field_names[2], field_names[6]]
        with open("{}\\{}".format(Output_Location,Output_Feature),'wb') as f:
            dw = csv.DictWriter(f,field_names)
            #--write all field names to the output file
            dw.writeheader()
            #--now we make the search cursor that will iterate through the rows of the table
            with arcpy.da.SearchCursor(Selection,field_names) as cursor:
                for row in cursor:
                    dw.writerow(dict(zip(field_names,row)))
        print "Data exported: {} s".format(time.time()-exportstart)
        #pool.close()
        #Delete useless select points
        arcpy.Delete_management(tmpOrigins))
        print "Done"

This is the version I am working on, that explains the many print and time commands. Regarding what happens in the second for loop, I was able to see that the main time-consuming operations are:

  1. Adding the Origins arcpy.AddLocations_na()
  2. Solving the OD cost matrix arcpy.Solve_na()
  3. Writing the .csv on the disk

I am a bit worried by the fact that the amount of time needed for these operations is not constant:

enter image description here

Could you help me find out how I can optimize my script so it will be time- and memory-efficient?

Also, I was wondering which strategy is best: processing a large batch of data then write a big file in the disk, or processing and writing small parts one after the other?

Finally, I tried using the multiprocessing pool functions, but I could not make it work. Is there a place for multiprocessing in this script?

closed as too broad by PolyGeo Sep 7 '18 at 17:46

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

  • Not a time efficiency improvement, but you might find readability is better if you use str.format(), e.g "Processing... "+str(x+1)+"/"+str(loops) == "Processing... {0}/{1}".format(x+1, loops) – Paul Jun 17 '14 at 12:19
  • "MakeFeatureLayer" instead of "select" could help a bit. – radouxju Jun 17 '14 at 13:49
2

I have a similar problem but because my process is iterative, I can have it easier by reusing the path and updating the feature every iteration. I also coded the program to be a standalone python rather than trying to run it inside arcGIS. I've basically coded everything into its own function, so for example, the part where I solve and export the OD table into CSV is a function that is called in the loop. This way whenever it finishes whatever objects that is created will be erased automatically, which probably would help in memory management.

In terms of multiprocessing, there's surely a place for it in your script. Again this is if you are running this as a standalone script. In my case I simply do it manually by running three or four separate python console running different versions of the same script. In your case since the calculation is independent of each other (not iterative like mine), one easy way to manually parallel process is to split the source CSV and run multiple instances of the code. Regarding using pool or other MP tools I also didnt get it to work, but my issue was primarily with splitting the OD calculation into four chunks, merging the results together, and continue with the rest of the program.

This is how I did the solving and exporting. Not super clean but works:

import arcpy
from arcpy import env
arcpy.env.overwriteOutput = True

def NAtoCSV(inSpace, inGdb, inNetworkDataset, impedanceAttribute, accumulateAttributeName, inOrigins, inDestinations, outNALayerName, outFile):
try:
    #Check out the Network Analyst extension license
    if arcpy.CheckExtension("Network") == "Available":
        arcpy.CheckOutExtension("Network")
    else:
        # Raise a custom exception
        print "Network license unavailable, make sure you have network analyst extension installed."


    #Check out the Network Analyst extension license
    arcpy.CheckOutExtension("Network")

    #Set environment settings
    env.workspace = inSpace + inGdb
    env.overwriteOutput = True

    #Create a new OD Cost matrix layer.
    outNALayer = arcpy.na.MakeODCostMatrixLayer(inNetworkDataset, outNALayerName,
                                                impedanceAttribute, "#", "#",
                                                accumulateAttributeName,
                                                "ALLOW_UTURNS","Oneway","NO_HIERARCHY","#","NO_LINES","#")


    #Get the layer object from the result object. The OD cost matrix layer can 
    #now be referenced using the layer object.
    outNALayer = outNALayer.getOutput(0)

    #Get the names of all the sublayers within the OD cost matrix layer.
    subLayerNames = arcpy.na.GetNAClassNames(outNALayer)

    #Stores the layer names that we will use later
    originsLayerName = subLayerNames["Origins"]
    destinationsLayerName = subLayerNames["Destinations"]
    linesLayerName = subLayerNames["ODLines"]

    #Adjust field names
    #Exploit the fact that the detector feature is named hd_ML_snap,
    #change the field mapping of Name to id_stn
    oriField = "Name ID_TAZ12A #"
    oriSort = "ID_TAZ12A"
    destField = "Name ID_TAZ12A #"
    destSort = "ID_TAZ12A"
    if inOrigins[-4:] == "snap":
        oriField = "Name id_stn #"
        oriSort = "id_stn"

    if inDestinations[-4:] == "snap":
        destField = "Name id_stn #"
        destSort = "id_stn"

    #Add locations
    arcpy.AddLocations_na(outNALayer, originsLayerName, inOrigins,
                          oriField, sort_field = oriSort, append = "CLEAR")
    arcpy.AddLocations_na(outNALayer, destinationsLayerName, inDestinations,
                          destField, sort_field = destSort, append = "CLEAR")    

    #Solve the OD cost matrix layer
    print "Begin Solving"
    arcpy.na.Solve(outNALayer)
    print "Done Solving"


    # Extract lines layer, export to CSV
    # SOMEHOW THIS WORKS, SO LETS KEEP IT THAT WAY
    fields = ["OriginID", "DestinationID", "Name", "Total_Cost"]
    for lyr in arcpy.mapping.ListLayers(outNALayer):
        if lyr.name == linesLayerName:
            with open(outFile, 'w') as f:
                #f.write(','.join(fields)+'\n') # csv headers
                with arcpy.da.SearchCursor(lyr, fields) as cursor:
                    print "Successfully created lines searchCursor.  Exporting to " + outFile
                    for row in cursor:
                        f.write(','.join([str(r) for r in row])+'\n')

    # Deleteing using del outNALayer is not enough.  Need to delete within arcpy to release
    arcpy.Delete_management(outNALayer)

except Exception as e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "An error occurred in ModRealSolve line %i" % tb.tb_lineno
    print str(e)        


finally:
    #Check the network analyst extension license back in, regardless of errors.
    arcpy.CheckInExtension("Network")

and then you just call the program:

#TT COST CALCULATION STARTS HERE
    outNALayerName = "ODTT"
    inOrigins = "Trans/"+fcTAZ
    inDestinations = "Trans/"+fcTAZ
    outFile = inSpace+"CSV/TT.csv"
    NAtoCSV(inSpace, inGdb, inNetworkDataset, impedanceAttribute, accumulateAttributeName, inOrigins, inDestinations, outNALayerName, outFile)
print "TT Solved"
  • About loading the origins, if you have sorted your networks and features in your CSV according to their complexity, then it make sense that later on the reading will become slower... – SWOT Oct 14 '14 at 18:56

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