Join "Lines" result from OD Cost Matrix to point feature class for interpolation

[HeidiErdle]

Hi,

My overall goal is to use network analyst (OD CostMatrix tool) to determine the haulcost (expressed as minutes) from each polygon (forest stand) to each mill in New Brunswick. I have the road network required and a points feature class designed to represent my polygon "origins" or starting points.

1) I need to create 77 cost matricies (one matrix for the travel time of each polygon in the province to 1 mill... there are 77 mills, so therefore 77 matricies to be created).

2)From there, I want to join each "lines" output from the OD CostMatrix to my point feature class (same as my "origins") and interpolate the travel time from each point to the travel time for each stand by creating a raster.

3)I will then use spatial zoning to apply the average raster value to my polygon (forest stand) layer.

I am automating the process since I need to complete the above steps for each of 77 mills. I am able to generate an OD cost matrix layer for each mill using python.

Here is my problem. I am unable to join "lines" to my point feature class; i.e. I cant' relate my cost matrix outputs back to the points of origin to create my raster layer.

How can I create a new feature class or select out my "lines" in the od cost matix layer so that I can join them to my point feature class?

I have tried arcpy.AddJoin_management - but it fails to exacute because my "join table" aka my "lines" output is not in an appropriate data format.

Can anyone help? I have attached the script I've created below... it works up until line 76 (where the OD cost matrix is solved). Some of the following commented out lines are my various attempts to turn my OD cost matrix outputs into something I can apply the AddJoin tool to. Advice would be much appreciated! Thanks!

- Heidi

# haulcosting.py
#
# Get haul time from each polygon to each mill in NB
#
#

import arcpy
import sys
from arcpy import env
import os

#Check Network Extension License
#
arcpy.CheckOutExtension("Network")
arcpy.CheckOutExtension ("spatial")

#Set environment
#
env.workspace = r"C:/Users/he1999/Desktop/FMP_Costs/Haul/NBROADNET.gdb"

#Set variables for cost matrix
#
inNetworkDataset = "NBNET/NBNETWORK"
outNALayerName_base = "TimeTo"
impedenceAttribute = "Minutes"
accumulateAttributeName = ["Minutes","Meters"]
inOrigins = "ThirdJunctions"
inDestinations = "Mills/OperatingMills2012"
where_clause_base = '"DESTCODE" ='
OD_CostMatrix = "OD_CostMatrix"

print where_clause_base
print "So far good till looping"

#create new od cost matrix layer. Find minutes and meters to each mill
#
# Create List for Destinations
#
MillList = [row[0] for row in arcpy.da.SearchCursor(inDestinations,["DESTCODE"])]

print "Creating Matrix layer"
outNALayer = arcpy.na.MakeODCostMatrixLayer(inNetworkDataset,"OD_CostMatrix","Minutes",accumulate_attribute_name = ["Minutes","Meters"],UTurn_policy = "NO_UTURNS",hierarchy= "USE_HIERARCHY",\
                                            hierarchy_settings="",output_path_shape="STRAIGHT_LINES")
outNALayer = outNALayer.getOutput(0)
subLayerNames = arcpy.na.GetNAClassNames(outNALayer)

#Store sublayer names for origins and destinations
#
originsLayername = subLayerNames["Origins"]
destinationsLayerName = subLayerNames["Destinations"]


#Load ThirdJunctions as origins and map the DESTCODE field from mills feature as MillName
#
print "loading origins"
arcpy.na.AddLocations (outNALayer, originsLayername, inOrigins,"","#",snap_to_position_along_network = "SNAP")

# Generate Mill looping for Cost Matrix
#
for Mill in MillList:
    print "Analyzing Mill" + Mill
    outNALayerName = outNALayerName_base + Mill
    outLayerFile = r"C:/Users/he1999/Desktop/FMP_Costs/Haul" + "/" + outNALayerName + ".lyr"
    where_clause = where_clause_base + "'" + Mill + "'"
    print where_clause
    arcpy.MakeFeatureLayer_management(inDestinations,"temp",where_clause)
    candiateFields = arcpy.ListFields("temp")
    fieldmappings = arcpy.na.NAClassFieldMappings(outNALayer,destinationsLayerName,True,candiateFields)
    arcpy.na.AddLocations (outNALayer, destinationsLayerName, "temp",fieldmappings,"#", append = "CLEAR")
  
#Solve OD cost matrix
#
    print "Starting Solve"
    arcpy.na.Solve (outNALayer, "SKIP")
    print "Done"
   
#Save solved matrix layer as layer file on disk with relative paths
#
    arcpy.MakeFeatureLayer_management("OD_CostMatrix","C:/Users/he1999/Desktop/FMP_Costs/Haul/NBROADNET.gdb/Timeto" + Mill)
    print "OD Saved as feature layer"
   
    #arcpy.management.SaveToLayerFile(outNALayer,outLayerFile,"RELATIVE")
    #print "Script completed successfully"
   
    #arcpy.CopyFeatures_management(outLayerFile,"C:/Users/he1999/Desktop/FMP_Costs/Haul/NBROADNET.gdb")
    #print "Saved to Feature Class for" + Mill
   
#Join Lines from OD cost matrix to Third Junctions
    #inFeatures = "NBROADNET.gdb"
   
    layerName = inOrigins
    joinTable = "OD_CostMatrix/Lines"
    inField = "OBJECTID"
    joinField = "OriginID"
    outFeature = "NBROADNET.gdb/Join" + Mill

#Select Data for join
    arcpy.SelectData_management(outNALayer, "Lines")
    print "Lines Selected"

    #Add Join
    arcpy.AddJoin_management(layerName, inField, joinTable, joinField, "KEEP_COMMON")
    print "Join Complete"
   
    #Copy layer to a new permanent feature class"
    arcpy.CopyFeatures_management(layerName, outFeature)
    print "PointHaul join saved as feature for" + Mill

    inPointFeatures = outFeature
    zField = "Minutes"
    cellSize = 100
    power = 2
    searchRadius = RadiusVariable(12,1000000)

    outIDW = Idw(inPointFeatures, zField, cellSize, power, searchRadius)

    outIDW.save("C:Users/he1999/Desktop/FMP_Costs/Haul/NBROADNET/InterpTo" + Mill)

    print "IDW Complete for" + Mill
    break

arcpy.CheckInExtension("Network")
arcpy.CheckInExtension ("spatial")

print "We got er done Cowgirl!"

[MelindaMorang]

First of all, you do not need to calculate 77 different OD matrices.  Just load the forest stands as origins and all the mills as destinations.  You will get one OD Line from each origin to each destination.  That said, depending on how many forest stands there are, it could be a very large OD matrix.  If you start running out of memory on Solve, you'll want to chunk up the problem into smaller pieces, similar to what you're already doing.

You can do the following to get the output Lines table so you can use it in other tools like AddJoin:

# ODLayer is the NA Layer object returned by getOutput(0)
ODLayer = arcpy.MakeODCostMatrixLayer_na(inNetworkDataset, outNALayer_OD,
                                impedanceAttribute, BufferSize, "",
                                accumulate, uturns, restrictions,
                                hierarchy, "", PathShape).getOutput(0)

# To refer to the OD sublayers, get the sublayer names.
if ArcVersion in ["10.1", "10.2"]:
    naSubLayerNames = arcpy.na.GetNAClassNames(ODLayer)
elif ArcVersion == "10.0":
    naSubLayerNames = dict((sublayer.datasetName, sublayer.name) for sublayer in  arcpy.mapping.ListLayers(ODLayer)[1:])
points = naSubLayerNames["Origins"]
stops = naSubLayerNames["Destinations"]
lines = naSubLayerNames["ODLines"]

# Make layer objects for each sublayer
linesSubLayer = arcpy.mapping.ListLayers(ODLayer, lines)[0]
pointsSubLayer = arcpy.mapping.ListLayers(ODLayer, points)[0]
stopsSubLayer = arcpy.mapping.ListLayers(ODLayer, stops)[0]

# ... use sublayers in other tools

Use the Origin_ID and Destination_ID fields in the output Lines to join back to the ObjectID field of Origins and Destinations