Finding shortest route between many points

[Jades1]

Hi all,

I'm trying to solve a problem that is not very complicated, but I can't seem to find a good solution in ESRI. I have different medi-cal clients and each one sees a different provider. In the dataset, I have client geolocation and provider geolocation, and I want to find the distance (shortest in time/or distance) for each client-provider pair. This is not an OD-cost matrix problem. It would seem to be a shortest path problem; however, it seems that shortest path is not really built for finding the shortest path between many different pairs. Instead, it is setup for finding the shortest path between many different points on a route. Essentially, if I wanted to find Euclidean distance, in geopandas it would be easy with geopandas.distance() command; however, I want to incorporate roadways.

I have to use ArcGIS pro, since there is protected health information involved. I mostly code in Python when using ArcGIS. What is the best way of doing this?

Thanks much!

James

[DanPatterson]

travel along roads requires the Network Analyst Extension, which is a separate package to ArcGIS Pro

Some of its functionality

Network Analyst tutorials—ArcGIS Pro | Documentation

[LindsayRaabe_FPCWA]

Assuming you have access to the Network Analyst extension that Dan has linked, you setup your Start and Finish points with Route Names. You need a start point to go with each corresponding finish point. Think of it like this - if you have Point A and you want to calculate routes to Points M, N & O, you need 3 copies of Point A, each with a unique Route Name that is also referenced by each finish point (i.e. Route A-M, A-N & A-O). When you run the Solve tool to calculate routes, it recognises those route names and calcs the route from A to N, then another for A-M and again for A-O. The same dataset could also then contain start points B, C, etc and other finish points, and as long as you have those start points copied and named to match each required finish location, it will only calculate those routes. It seems like a simple concept, but it took me a while to figure this out so, thought I'd share it. 

Below is an extract from a python code that generates the Start and Finish points based on 2 input datasets (the start points actually being polygons so you can see the dissolve and conversion steps added in to deal with this). The code goes on to solve the network too, but I haven't added that as I haven't got it running properly in an automated manner (which is how it is supposed to run). I have a model that does this too for users to use in ArcGIS Pro (what I used to build the below code in the first place). 

# Variables - update these lines
SDPRoutes = "Output feature service" # ArcGIS Online feature service containing routes to be updated
origins = r"\\path\to\input\origin\features"  # Input feature class or table (can be .csv/.xls/.xlsx/etc.) or service URL of feature service in ArcGIS Online
destinations = "\\path\to\input\destination\features"    # Input feature class or table (can be .csv/.xls/.xlsx/etc.) or service URL of feature service in ArcGIS Online
nd_path = r"\\path\to\network\dataset"
travelmode = "RAV 3"

import arcpy, time, os
from arcgis.gis import GIS
from arcgis.features import FeatureLayer

# Start Timer
startTime = time.time()

# Overwrite Output
arcpy.env.overwriteOutput = True

try:
    arcpy.CheckOutExtension("Network")
    print ("Network Analyst license checked-out")

    # Create GIS object
    print("Connecting to AGOL")
    #...login to ArcGIS Online section omitted...

    # Create temporary FGDB
    sdproutesgdb = arcpy.CreateFileGDB_management(arcpy.env.scratchFolder, "sdproutes")
    print ("Temp FGDB created")

    arcpy.env.workspace = arcpy.env.scratchFolder + "\\sdproutes.gdb"
    print ("Workspace set to temp FGDB")
    
    #...checks of routes in existing dataset omitted...
    
    # Prep Origin data
    origins_temp = arcpy.conversion.ExportFeatures(
        in_features = origins, 
        out_features = arcpy.env.workspace + "\\Origins_Temp",
        where_clause = "PatchClass = 1 AND Map_Symbol = 'Softwood'",
        use_field_alias_as_name = "NOT_USE_ALIAS",
        field_mapping = ,
        sort_field = None
    )
    print("Origins exported")
    
    origins_dissolved = arcpy.management.Dissolve(
        in_features = origins_temp,
        out_feature_class = arcpy.env.workspace + "\\Origins_Dissolved",
        dissolve_field = "Ops_Code",
        statistics_fields = None,
        multi_part = "MULTI_PART",
        unsplit_lines = "DISSOLVE_LINES",
        concatenation_separator = ""
    )
    print("Origins dissolved by Ops_Code")

    origin_points = arcpy.management.FeatureToPoint(
        in_features = origins_dissolved,
        out_feature_class = arcpy.env.workspace + "\\Origin_Points",
        point_location = "CENTROID"
    )
    print("Dissolved origins converted to points")

    origin_points = arcpy.management.AddFields(
        in_table = origin_points,
        field_description = "Name TEXT # 500 # #;RouteName TEXT # 1024 # #;Sequence LONG # # # #",
        template = None
    )
    print("Route fields added to Origin points")

    origin_points = arcpy.management.CalculateField(
        in_table = origin_points,
        field = "Name",
        expression = "!Ops_Code!",
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print("Name field calculated")

    origin_points = arcpy.management.CalculateField(
        in_table = origin_points,
        field = "Sequence",
        expression = "1",
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print("Sequence field calculated")

    static_origins = arcpy.analysis.Select(
        in_features = origin_points,
        out_feature_class = arcpy.env.workspace + "\\Static_Origins",
        where_clause = ""
    )
    print("Static Origin points generated")

    # Prep Destination data
    fl = FeatureLayer(destinations)
    fs = fl.query(where = "Delivery_Location_ID IN ('1015DA1', '1015DA2', '1043DA1', '1043DA6', '344D13', '410D10', '410DA1', '410DA2', '410DA8', '426DA1', '503DA1', '5097DA2', '5097DA7', '5102DA1')") # Query to only return required features
    destinations_temp = fs.save(arcpy.env.workspace, "destinations_temp")
    print("Required destinations exported from feature service")

    destination_points = arcpy.analysis.Select(
        in_features = destinations_temp,
        out_feature_class = arcpy.env.workspace + "\\destination_points",
        where_clause = ""
    )
    print("Destination Points generated with ObjectID's reset from 1")
    
    destination_points = arcpy.management.AddFields(
        in_table = destination_points,
        field_description = "Name TEXT # 500 # #;RouteName TEXT # 1024 # #;Sequence LONG # # # #",
        template = None
    )
    print("Route fields added to Destination points")

    destination_points = arcpy.management.CalculateField(
        in_table = destination_points,
        field = "Name",
        expression = "!Delivery_Location_ID!",
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print("Name field calculated")

    destination_points = arcpy.management.CalculateField(
        in_table = destination_points,
        field = "Sequence",
        expression = "2",
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print("Sequence field calculated")
_rows = static_destinations
    )[0])
    print("Static Destination count: " + str(DestinationCount))

    TargetCount = OriginCount*DestinationCount
    print ("Target Count = " + str(TargetCount))

    # Iterate through Origins to replicate points for each Destination
    print ("Start iterating Origin points to create 1 set of all Origin points for each Destination point")
    while OriginCount < TargetCount:
        print ("Origin count (" + str(OriginCount) + ") less than Target Count (" + str(TargetCount) + ") - Append data...")
        
        arcpy.management.Append(
        inputs = static_origins,
        target = origin_points,
        schema_type = "TEST",
        field_mapping = None,
        subtype = "",
        expression = "",
        match_fields = None,
        update_geometry = "NOT_UPDATE_GEOMETRY"
        )

        OriginCount = int(arcpy.management.GetCount(
        in_rows = origin_points
        )[0])
    print ("Origins Replicated")

    sorted_origins = arcpy.management.Sort(
        in_dataset = origin_points,
        out_dataset = arcpy.env.workspace + "\\sorted_origins",
        sort_field = "Name ASCENDING",
        spatial_sort_method = "UR"
    )
    print ("Origins sorted by Name")

    # Iterate through Destinations to replicate points for each Origin
    print ("Start iterating Destination points to create 1 set of all Destination points for each Origin point")
    while DestinationCount < TargetCount:
        print ("Destination count (" + str(DestinationCount) + ") less than Target Count (" + str(TargetCount) + ") - Append data...")
        
        arcpy.management.Append(
        inputs = static_destinations,
        target = destination_points,
        schema_type = "NO_TEST",
        field_mapping = r'Name "Name" true true false 500 Text 0 0,First,#,static_destinations,Name,0,500,destination_points,Name,0,500;RouteName "RouteName" true true false 1024 Text 0 0,First,#,static_destinations,RouteName,0,1024,destination_points,RouteName,0,1024;Sequence "Sequence" true true false 4 Long 0 0,First,#,static_destinations,Sequence,-1,-1,destination_points,Sequence,-1,-1',
        subtype = "",
        expression = "",
        match_fields = None,
        update_geometry = "NOT_UPDATE_GEOMETRY"
        )

        DestinationCount = int(arcpy.management.GetCount(
        in_rows = destination_points
        )[0])
    print ("Destinations Replicated")

    # Join Origins and Destinations to calculate route details
    sorted_origins = arcpy.management.JoinField(
        in_data = sorted_origins,
        in_field = "OBJECTID",
        join_table = destination_points,
        join_field = "OBJECTID",
        fields = "Name",
        fm_option = "NOT_USE_FM",
        field_mapping = None
    )
    print ("Destination names added to Origin points")

    sorted_origins = arcpy.management.CalculateField(
        in_table = sorted_origins,
        field = "RouteName",
        expression = '!Name!+"_"+!Name_1!',
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print ("Origin Route Names calculated")

    arcpy.management.JoinField(
        in_data = destination_points,
        in_field = "OBJECTID",
        join_table = sorted_origins,
        join_field = "OBJECTID",
        fields = "Name",
        fm_option = "NOT_USE_FM",
        field_mapping = None
    )
    print ("Origin names added to Destination points")
    
    destination_points = arcpy.management.CalculateField(
        in_table = destination_points,
        field = "RouteName",
        expression = '!Name_1!+"_"+!Name!',
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print ("Destination Route Names calculated")
    
    merged_locations = arcpy.management.Merge(
        inputs = [sorted_origins, destination_points],
        output = arcpy.env.workspace + "\\MergedLocations",
        field_mappings = 'Ops_Code "Ops_Code" true true false 8000 Text 0 0,First,#,sorted_origins,Ops_Code,0,8000;ORIG_FID "ORIG_FID" true true false 4 Long 0 0,First,#,sorted_origins,ORIG_FID,-1,-1;Name "Name" true true false 500 Text 0 0,First,#,sorted_origins,Name,0,500,destination_points,Name,0,500;RouteName "RouteName" true true false 1024 Text 0 0,First,#,sorted_origins,RouteName,0,1024,destination_points,RouteName,0,1024;Sequence "Sequence" true true false 4 Long 0 0,First,#,sorted_origins,Sequence,-1,-1,destination_points,Sequence,-1,-1;ORIG_FID_1 "ORIG_FID_1" true true false 4 Long 0 0,First,#,sorted_origins,ORIG_FID_1,-1,-1;Name_1 "Name" true true false 500 Text 0 0,First,#,sorted_origins,Name_1,0,500,destination_points,Name_1,0,500;Delivery_Location_ID "Delivery Location ID" true true false 255 Text 0 0,First,#,destination_points,Delivery_Location_ID,0,255',
        add_source = "NO_SOURCE_INFO"
    )
    print ("Origins and Destinations merged into single dataset")
    static_destinations = arcpy.analysis.Select(
        in_features = destination_points,
        out_feature_class = arcpy.env.workspace + "\\Static_Destinations",
        where_clause = ""
    )
    print("Static Destination points generated")

    # Get counts of input locations
    OriginCount = int(arcpy.management.GetCount(
        in_rows = static_origins
    )[0])
    print("Static Origin count: " + str(OriginCount))

    DestinationCount = int(arcpy.management.GetCount(
        in_rows = static_destinations
    )[0])
    print("Static Destination count: " + str(DestinationCount))

    TargetCount = OriginCount*DestinationCount
    print ("Target Count = " + str(TargetCount))

    # Iterate through Origins to replicate points for each Destination
    print ("Start iterating Origin points to create 1 set of all Origin points for each Destination point")
    while OriginCount < TargetCount:
        print ("Origin count (" + str(OriginCount) + ") less than Target Count (" + str(TargetCount) + ") - Append data...")
        
        arcpy.management.Append(
        inputs = static_origins,
        target = origin_points,
        schema_type = "TEST",
        field_mapping = None,
        subtype = "",
        expression = "",
        match_fields = None,
        update_geometry = "NOT_UPDATE_GEOMETRY"
        )

        OriginCount = int(arcpy.management.GetCount(
        in_rows = origin_points
        )[0])
    print ("Origins Replicated")

    sorted_origins = arcpy.management.Sort(
        in_dataset = origin_points,
        out_dataset = arcpy.env.workspace + "\\sorted_origins",
        sort_field = "Name ASCENDING",
        spatial_sort_method = "UR"
    )
    print ("Origins sorted by Name")

    # Iterate through Destinations to replicate points for each Origin
    print ("Start iterating Destination points to create 1 set of all Destination points for each Origin point")
    while DestinationCount < TargetCount:
        print ("Destination count (" + str(DestinationCount) + ") less than Target Count (" + str(TargetCount) + ") - Append data...")
        
        arcpy.management.Append(
        inputs = static_destinations,
        target = destination_points,
        schema_type = "NO_TEST",
        field_mapping = r'Name "Name" true true false 500 Text 0 0,First,#,static_destinations,Name,0,500,destination_points,Name,0,500;RouteName "RouteName" true true false 1024 Text 0 0,First,#,static_destinations,RouteName,0,1024,destination_points,RouteName,0,1024;Sequence "Sequence" true true false 4 Long 0 0,First,#,static_destinations,Sequence,-1,-1,destination_points,Sequence,-1,-1',
        subtype = "",
        expression = "",
        match_fields = None,
        update_geometry = "NOT_UPDATE_GEOMETRY"
        )

        DestinationCount = int(arcpy.management.GetCount(
        in_rows = destination_points
        )[0])
    print ("Destinations Replicated")

    # Join Origins and Destinations to calculate route details
    sorted_origins = arcpy.management.JoinField(
        in_data = sorted_origins,
        in_field = "OBJECTID",
        join_table = destination_points,
        join_field = "OBJECTID",
        fields = "Name",
        fm_option = "NOT_USE_FM",
        field_mapping = None
    )
    print ("Destination names added to Origin points")

    sorted_origins = arcpy.management.CalculateField(
        in_table = sorted_origins,
        field = "RouteName",
        expression = '!Name!+"_"+!Name_1!',
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print ("Origin Route Names calculated")

    arcpy.management.JoinField(
        in_data = destination_points,
        in_field = "OBJECTID",
        join_table = sorted_origins,
        join_field = "OBJECTID",
        fields = "Name",
        fm_option = "NOT_USE_FM",
        field_mapping = None
    )
    print ("Origin names added to Destination points")
    
    destination_points = arcpy.management.CalculateField(
        in_table = destination_points,
        field = "RouteName",
        expression = '!Name_1!+"_"+!Name!',
        expression_type = "PYTHON3",
        code_block = "",
        field_type = "TEXT",
        enforce_domains = "NO_ENFORCE_DOMAINS"
    )
    print ("Destination Route Names calculated")
    
    merged_locations = arcpy.management.Merge(
        inputs = [sorted_origins, destination_points],
        output = arcpy.env.workspace + "\\MergedLocations",
        field_mappings = 'Ops_Code "Ops_Code" true true false 8000 Text 0 0,First,#,sorted_origins,Ops_Code,0,8000;ORIG_FID "ORIG_FID" true true false 4 Long 0 0,First,#,sorted_origins,ORIG_FID,-1,-1;Name "Name" true true false 500 Text 0 0,First,#,sorted_origins,Name,0,500,destination_points,Name,0,500;RouteName "RouteName" true true false 1024 Text 0 0,First,#,sorted_origins,RouteName,0,1024,destination_points,RouteName,0,1024;Sequence "Sequence" true true false 4 Long 0 0,First,#,sorted_origins,Sequence,-1,-1,destination_points,Sequence,-1,-1;ORIG_FID_1 "ORIG_FID_1" true true false 4 Long 0 0,First,#,sorted_origins,ORIG_FID_1,-1,-1;Name_1 "Name" true true false 500 Text 0 0,First,#,sorted_origins,Name_1,0,500,destination_points,Name_1,0,500;Delivery_Location_ID "Delivery Location ID" true true false 255 Text 0 0,First,#,destination_points,Delivery_Location_ID,0,255',
        add_source = "NO_SOURCE_INFO"
    )
    print ("Origins and Destinations merged into single dataset")

[Jades1]

Awesome, thanks so much for this! I work for UCSDH, so I checked licenses, and it does look like I have a license for Network Analyst (is this the extension, or is there something else called Network Analyst Extension?)

Yeah, so currently every row has both a point A and a point B, so that should be good to go.

I'm just using point data, so a bit easier than polygons.

I'll look through this in more depth and let you know if I have questions. I've mostly worked through AGOL, so Pro is a bit new to me. Are you using Visual Studio Code with Arc vs their Jupyter Notebook extension?

[LindsayRaabe_FPCWA · Accepted solution]

Yep - it's just the Network Analyst extension as you've shown. I do most of my python work using IDLE or within Pro itself. I've never used Jupyter Notebooks or Visual Studio Code. 

[Dale_Honeycutt]

Once you have the Network Analyst Extension, take a look at Closest Facility Analysis.  Your provider would be the facilities and the clients would be the incidents.  Set the problem up so that clients are traveling to the facilities and you set the number of facilities to find to 1 (one).  

[Jades1]

Yeah, I've used finding nearest point in AGOL, but I thought it was really solving the problem of one location to many (since all providers are at different locations). I will look through closest facility in more depth and let you know if I have questions. Thanks!

[Jades1]

Wouldn't this just be finding the closest provider though? Vs the actual distance to each client's provider (which may not be the closest provider?

[Dale_Honeycutt]

My bad -- I didn't catch that you already had the provider-client pairs.  You can find the paths between them using the Route Analysis solver.  Each stop (the point to visit along the route) can be assigned its own unique RouteName.  The solver will create a shortest path for each unique RouteName.  So, just give each of your provider-client pairs a unique RouteName.  Set the client point as the route origin (sequence = 0) and the provider as route destination (sequence=1).  Of course this means you'll have to duplicate points for the providers, one point for each client assigned to that provider.  It seems you know Python pretty well, so maybe write a script that takes the provider-client pairs and their locations and outputs a point feature class with two features for each pair along with the RouteName/Sequence. (This would be useful code to share)

See https://pro.arcgis.com/en/pro-app/latest/help/analysis/networks/route-analysis-layer.htm for more information on the different settings for route analysis.