0

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

enter image description here

And my network dataset looks like this:

enter image description here

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 ...)?

Improve this question
13
  • 1
    I think your problem may be at the beginning here: arcpy.SaveToLayerFile_management(out_Layer, saved_Layer). You should create an actual GIS layer rather than a .lyr from the csv data. A .lyr is just a pointer file to an actual GIS layer (e.g. shp or gdb/fc). Try using feature class to feature class arcpy method to create the layer from the make xy event layer – artwork21 Feb 12 '16 at 16:31
  • edit, ...rather than a .lyr from the event layer. The event layer is stored in memory and may not exist depending on how your script is written, hence the .lyr may not be able to source the layer. – artwork21 Feb 12 '16 at 16:38
  • @artwork21 Thanks! However, I'm not sure what you mean exactly (maybe because I edited the question). The solver works for 100 calculations but ends up taking 428 seconds to do it. – mptevsion Feb 12 '16 at 17:02
  • @artwork21 I tried adding this line arcpy.CopyFeatures_management(saved_Layer,out_shape) and then using the .shp file, however it hasn't really made any difference – mptevsion Feb 12 '16 at 17:08
  • 1
    4-8 million seconds. – Vince Feb 12 '16 at 21:14
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).

Improve this answer
1
  • Good you solved this. Looking at the link you specified, this seems to indicate what I thought: there is talk about the need to use osrm-prepare to pre-process the data. Also looking at this help.openstreetmap.org/questions/30272/… link and Richard's response referring to "Contraction hierarchies", it seems OSRM uses preprocessing of the data to create subsequent highly efficient routing requests. It may be that ArcGIS's network dataset currently is a bit outdated in this respect, and uses older algorithms to compute routes. – Marco_B Feb 29 '16 at 15:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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