Route Layer with ArcPy and OSM Data - very slow

Question

This is my first attempt at using arcpy: I have two CSVs which have columns:

CSV1: join|id_a|lat_a|lon_a
CSV2: join|id_b|lat_b|lon_b

Each has got 100 rows of data. I would like to write a script that gives me the drive-time from row[join] in csv1 to row[join] in csv2.

First I convert the CSVs into layer files:

env.workspace = ".../Documents/ArcGIS/Test"
for x in ["a", "b"]:
    try:
        in_Table = "test_%s.csv" % x
        x_coords = "lon_%s" % x
        y_coords = "lat_%s" % x
        out_Layer = "test_%s_layer" % x
        saved_Layer = "test_%s.lyr" % x
        spRef = "Coordinate Systems/Geographic Coordinate Systems/World/WGS 1984.prj"
        arcpy.MakeXYEventLayer_management(in_Table, x_coords, y_coords, out_Layer, spRef)
        print arcpy.GetCount_management(out_Layer)
        arcpy.SaveToLayerFile_management(out_Layer, saved_Layer)
    except:
        # If an error occurred print the message to the screen
        print "Something went wrong"
        print arcpy.GetMessages()

Then I try to do the route layer:

#http://desktop.arcgis.com/en/arcmap/10.3/tools/network-analyst-toolbox/make-route-layer.htm
import arcpy
import time
from arcpy import env
try:
    # Check out the Network Analyst extension license
    arcpy.CheckOutExtension("Network")
    # Set environment settings
    env.workspace = ".../ArcGIS/Default.gdb"
    env.overwriteOutput = True
    # Set local variables
    inNetworkDataset = ".../ArcGIS/osm.gdb/france_osm/france_osm_nd"
    outNALayerName = "BestRoute"
    impedanceAttribute = "DriveTime"
    outLayerFile = ".../ArcGIS/Test/outRoutes"
    # Locations(should have common field: join)
    route_a = ".../ArcGIS/Test/test_a.lyr"
    route_b = ".../ArcGIS/Test/test_b.lyr"
    # My original excel had columns: join|id_a|lat_a|lon_a
    # Create a new Route layer
    outNALayer = arcpy.na.MakeRouteLayer(inNetworkDataset, outNALayerName, impedanceAttribute)
    # Get the layer object from the result object
    outNALayer = outNALayer.getOutput(0)
    # Get the names of all the sub-layers within the route layer.
    subLayerNames = arcpy.na.GetNAClassNames(outNALayer)
    # Stores the layer names that we will use later
    stopsLayerName = subLayerNames["Stops"]
    routesLayerName = subLayerNames["Routes"]
    fieldMappings = arcpy.na.NAClassFieldMappings(outNALayer, stopsLayerName)
    fieldMappings["RouteName"].mappedFieldName = "join"
    # Add my FROM location
    arcpy.na.AddLocations(outNALayer, stopsLayerName, route_a,
                          fieldMappings, "",
                          exclude_restricted_elements="EXCLUDE")
    # Add my TO location
    arcpy.na.AddLocations(outNALayer, stopsLayerName, route_b,
                          fieldMappings, "",
                          exclude_restricted_elements="EXCLUDE")
    # Solve
    print "Solving..."
    stime = time.time()
    arcpy.na.Solve(outNALayer, "SKIP")
    print "Finished in %.0f seconds" % (time.time() - stime)
    # Save the solved route layer as a layer file on disk with relative paths
    RoutesSubLayer = arcpy.mapping.ListLayers(outNALayer, routesLayerName)[0]
    arcpy.management.CopyFeatures(RoutesSubLayer, outLayerFile)
    print "Script completed successfully"
except Exception as e:
    import traceback, sys
    tb = sys.exc_info()[2]
    print "An error occured on line %i" % tb.tb_lineno
    print str(e)

However, this is extremely slow (and only for 100 locations?) taking more than 10 minutes. Is this because the OSM network database is broken? Or just unoptimised.

Edit: 100 rows ended up taking 428 seconds!!

It works but super slow. I hope my network dataset is not corrupt as it took a week to get it from OSM into ArcCatalog

And my network dataset looks like this:

Can I check it's health somehow?

---

So this solves the speed issue but leaves me even more confused. Instead of crossing all my stores and imposing a cutoff of 50km with scipy.spatial.cKDTree to generate around 20mill O-D pairings within 50km -> proxy for being within 30 minutes, which I then wanted to solve with the route layer (as above) I decided to just try running an OD cost matrix with a cut-off of 30 minutes.

This went through all 36,000 origins (and = destinations) and generated 5.7 million drive-times in just 55 minutes!

So (36253 origins * 36253 destinations), 30min cut-off, ND = France OSM:

        stime = time.time()
        arcpy.CheckOutExtension("Network")
        env.workspace = wpath
        env.overwriteOutput = True
        # Set local variables
        inNetworkDataset = "france_osm/france_osm_nd"
        outNALayerName = "ODCost_30min_36k"
        impedanceAttribute = "DriveTime"
        # Search tolerance: limits the distance from the location point
        # to network lines for network analyst to consider them valid
        searchTolerance = "1000 Meters"
        inBoth = "test_b.lyr"
        outLayerFile = outNALayerName + ".lyr"
        #Create a new OD Cost matrix layer
        print "Creating OD layer"
        outNALayer = arcpy.na.MakeODCostMatrixLayer(inNetworkDataset, outNALayerName,
                                                    impedanceAttribute, 30, "",
                                                    output_path_shape="NO_LINES")
        #Get the layer object from the result 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"]
        #Load the warehouse locations as origins using a default field mappings and
        #a search tolerance of 1000 Meters.
        print "Loading Origin"
        arcpy.na.AddLocations(outNALayer, originsLayerName, inBoth, "",
                              searchTolerance)
        #Load the store locations as destinations
        fieldMappings = arcpy.na.NAClassFieldMappings(outNALayer, destinationsLayerName)
        fieldMappings["Name"].mappedFieldName = "codeinsee"
        print "Loading Destination"
        arcpy.na.AddLocations(outNALayer, destinationsLayerName, inBoth,
                              fieldMappings, searchTolerance)
        #Solve the OD cost matrix layer
        print "Solving..."
        arcpy.na.Solve(outNALayer)
        print "Finished in %.0f" % (time.time() - stime)
        #Save the solved OD cost matrix layer as a layer file on disk with relative
        #paths
        arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"RELATIVE")
        print "Script completed successfully

How is possible to go from 4 seconds for a distance calculation to nearly 2k a second? Is it because the raw C code vectorises the calculation so 1 origin to 1 destination roughly takes the same as 1 origin to 10,000 destination? Is the OD cost matrix just faster than a route layer?

Maybe since no cut-off was imposed for the route layer it was harmed by small breaks in the network (as it would have to solve a maze through all of France), whereas here it just cuts off after 30minutes. Also I specified not to generate lines which I forgot with the route layer (I guess that speeds things up too).

Can anyone give me some average times for OD cost matrix and route layer so that I can orientate myself (e.g. seconds, or minutes, or hours ...)?

Answer 1

I wrote a script that uses OSRM for routing and I can get 1000 requests a second (not 4 seconds per request as stated somewhere in the comments).