How to iterate traces and feature selections

833
10
Jump to solution
11-05-2020 12:27 PM
TimWhiteaker
Occasional Contributor II

I have a GP model that performs a trace and then assigns a long integer stored in the HydroID variable to a field in the edge layer.  Calculate Field has the Updated Trace Network as a precondition, so it waits until the trace is finished and there is a selection.  This works fine when I run the model.  The model is shown below, and is named "_Single Junction To Edges".

GP model to trace and assign values to results

The problem arises when I encapsulate this model in another model with an iterator.  I want to iterate through a selected set of junctions, and for each one, perform the trace and calculate the field.  However, calculate field appears to be running on the entire edge layer rather than the selection.  As the model runs, I can see the selections being made in the map, but the tool ignores them.  The model with the iterator, which includes "_Single Junction To Edges" as one of the tools, is shown below.

Iterator GP model

What do I need to do to get the iterator version to work like the standalone version with respect to using feature selections when it comes time to calculate the field?

I do acknowledge that outputting a group layer with the selection from the Trace tool in Pro 2.7, as discussed in this thread on working with trace results, may be the solution.  But it seems like the single model was working correctly, so I wonder if I'm just missing something that will make it work in the iterator.

I'm on Pro 2.6.

Tags (1)
0 Kudos
1 Solution

Accepted Solutions
TimWhiteaker
Occasional Contributor II

I figured out how to prevent the downstream edge from being returned in the upstream trace.  My mistake was simply copying my selected junction to a separate flag feature class.  To make this work, I also have to add and populate a couple of fields that you'll notice in the TN_Temp_Starting_Points feature class in your default geodatabase that stores the flags you create interactively (via Data tab on ribbon), namely, SOURCEID and FEATUREGLOBALID.  The SOURCEID refers to the identifier given to the junction layer in the trace network.  The FEATUREGLOBALID is the GLOBALID of the junction feature where you want the starting point placed.  I wish this information had been in a white paper that I had access to.

Here's code that demonstrates how to trace upstream given a junction layer.

import arcpy


def get_source_id(network, network_layer):
    d = arcpy.Describe(network_layer)
    fc_name = d.featureClass.name
    
    d = arcpy.Describe(network)
    for s in d.sources:
        if s.name == fc_name:
            return s.sourceID

    arcpy.AddError(
        'Could not determine network source ID for ' + network_layer.name)
    raise arcpy.ExecuteError


def point_to_flag(shape_xy, spatial_ref, source_id, global_id):
    flag_fc = arcpy.CreateFeatureclass_management(
        out_path='memory',
        out_name='flag',
        geometry_type='POINT',
        has_m='ENABLED',
        has_z='ENABLED',
        spatial_reference=spatial_ref)[0]
    flag_fc = arcpy.management.AddField(flag_fc, 'SOURCEID', 'LONG')[0]
    flag_fc = arcpy.management.AddField(flag_fc, 'FEATUREGLOBALID', 'GUID')[0]
    fields = ['SHAPE@XY', 'SOURCEID', 'FEATUREGLOBALID']
    with arcpy.da.InsertCursor(flag_fc, fields) as cursor:
        cursor.insertRow((shape_xy, source_id, global_id))
    return flag_fc


arcpy.env.overwriteOutput = True

network = arcpy.GetParameter(0)
junction_layer = arcpy.GetParameter(1)

spatial_ref = arcpy.Describe(junction_layer).spatialReference
source_id = get_source_id(network, junction_layer)

with arcpy.da.SearchCursor(junction_layer, ['SHAPE@XY', 'GLOBALID']) as cursor:
    for i, row in enumerate(cursor):        
        flag = point_to_flag(row[0], spatial_ref, source_id, row[1])

        Updated_Trace_Network = arcpy.tn.Trace(
            in_trace_network=network,
            trace_type='UPSTREAM',
            starting_points=flag)[0]

arcpy.SetParameter(2, junction_layer)

 

So, until Pro 2.7 is released in which the trace tools return selection layers, the keys to iterating traces are:

  • Create a point feature class for your starting point(s), and populate SOURCEID and  FEATUREGLOBALID attributes in that feature class.
  • Return aggregated geometry, and use that geometry to select features from your network layer(s).
  • If you're updating a network layer's attributes while iterating, when you call arcpy.tn.Trace remember to set validate_consistency='DO_NOT_VALIDATE_CONSISTENCY' so the tool executes even as you update features.
  • Run arcpy.ValidateNetworkTopology_tn when you need to validate the network after updating it.

I've implemented this via Python. It seems like it should also work via ModelBuilder iteration, but ModelBuilder iteration confuses me, so I haven't verified it works there.  My id_to_edges.py script, and the utils.py script it relies on, are attached in case folks find it helpful.

View solution in original post

0 Kudos
10 Replies
DavidPike
MVP Frequent Contributor

What's the output destination set to in _SingleJunctionoEdges in the iterated tool?  Have you used inline variable substitution?

If not, I would guess your output is just being overwritten each time.

Or script that bad boy...

0 Kudos
TimWhiteaker
Occasional Contributor II

I don't think there is an output destination in _SingleJunctionoEdges.  The input to that tool is a layer, and the tool calculates a field on the layer, so the output labeled Output Feature Class is actually the input layer.  Or maybe I'm misunderstanding.

I put the geodatabase and toolbox here in case you want to have a look:

NetworkTest.zip - Box 

To test:

1. Add the Hydro network in the geodatabase to a map.

2. Select a single WRProxyMZ point (or more, but it just takes longer).

3. In the toolbox, run Junctions To Edges Iterator, with the Hydro network, WRProxyMZ, and NHDFlowline as inputs.

The tool's job is to trace upstream from each selected point, using all other WRProxyMZ points as barriers.  For each upstream edge found between the starting point and the upstream points, assign the JunctionID attribute of the edges to be equal to the HydroID of the starting point.

The freshly unzipped geodatabase has all JunctionIDs as -1.  When the tool finishes, the edges in the trace results for the final iteration are selected in the map.  However, the calculate field task ran on the entire layer and not just the selection, as evidenced by all NHDFlowlines having the same JunctionID, e.g., 1392.

0 Kudos
TimWhiteaker
Occasional Contributor II

I attempted the same task in a script tool and found the same result: The input edge layer does not have the trace results as a selection when arcpy.tn.Trace completes.  I'm thinking the only recourse until 2.7 is released, is to output an aggregated feature and then use a spatial filter to select the appropriate edges, which will be very slow.  I'm also worried about geometries not matching exactly since exported junctions don't seem to return the same results as an interactive starting point.

import arcpy


arcpy.env.overwriteOutput = True

network = arcpy.GetParameter(0)
junction_layer = arcpy.GetParameter(1)
edge_layer = arcpy.GetParameter(2)


# Trace network overwrites selections, so grab HydroIDs before tracing
hydro_ids = []
with arcpy.da.SearchCursor(junction_layer, ['HydroID']) as cursor:
for row in cursor:
hydro_ids.append(row[0])

for hydro_id in hydro_ids:
# Using junction layer even with filter applied results in all points used,
# so send the junction to a temporary feature class for use in tracing
start_pt = arcpy.conversion.FeatureClassToFeatureClass(
in_features=junction_layer,
out_path='memory',
out_name='current_junction',
where_clause='HydroID = ' + str(hydro_id))[0]
Updated_Trace_Network = arcpy.tn.Trace(
in_trace_network=network,
trace_type='UPSTREAM',
starting_points=start_pt,
condition_barriers=[['HydroID', 'DOES_NOT_EQUAL', 'SPECIFIC_VALUE', hydro_id, 'OR']],
traversability_scope='JUNCTIONS_ONLY',)[0]
if Updated_Trace_Network:
with arcpy.da.SearchCursor(edge_layer, ['JunctionID']) as cursor:
count = 0
for row in cursor:
count += 1
# Should be a selected subset, but we're getting the entire feature class
arcpy.AddMessage('Edges: ' + str(count))

arcpy.SetParameter(3, edge_layer)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍
0 Kudos
DavidPike
MVP Frequent Contributor

Are you sure of the brackets for some of those?

where clause - is hydro_id a list you get the index 0 from?

condition_barriers - looks funky, brackets and list in a list, as well as the [0] also?

if condition barriers is wrong wouldn't that cause what you're experiencing?

0 Kudos
TimWhiteaker
Occasional Contributor II

hydro_id is a number, like 1392.  It's used to create a whereclause like "HydroID = 1392".  I use hydro_id to uniquely identify each feature in my feature class.  It's unique like ObjectID, but I have control over it.  In line 24, whereclause is one of the optional arguments in the FeatureClassToFeatureClass function that is called on line 20.  There is a closing paren on line 24 that completes the function.  Geoprocessing tools can return one or more outputs, so the [0] is grabbing the first output, which is the output feature class.

The same thing goes for the [0] bracket at the end of the arcpy.tn.Trace call that starts on line 25.  I think condition barriers is right.  I created that code by exporting it as a Python script from a working GP Model in ModelBuilder, and then tweaking it to run as a script tool.  Also, you can have multiple conditions, and each condition has a list of parameters, so the conditions argument is a list of conditions, i.e., a list of lists.  In my case, I only have one condition, so the outer list only has one element in it (the inner list), and I can see how that looks funny.

Even if those condition barriers were wrong, if it didn't cause an error, the worst it would do is to make all upstream edges be selected instead of just the ones between the starting point and the nearest upstream barriers.  However, even downstream edges are included in that search cursor on line 32, so it's like the layer isn't yet seeing the selection from the network trace. 

I observed that as the tool executes, I do see the selections from each trace happening in the map, and when the tool finishes, the selection from the last trace is still there.  But, the edge layer as seen by the code doesn't see that selection; ditto when run within an iterator in ModelBuilder.

0 Kudos
DavidPike
MVP Frequent Contributor

You're very patient Tim, I see now (so used referencing the variables before the tool).

I'm also still confused about what's going on! I'll leave this to the smart people.

Good luck!

0 Kudos
TimWhiteaker
Occasional Contributor II

I'm thinking the only recourse until 2.7 is released, is to output an aggregated feature and then use a spatial filter to select the appropriate edges

That's what I wound up doing.  I used Python rather than ModelBuilder because I find it easier to follow.  Some example code is below.  I also heard Arc Hydro is using trace networks now, so I will peak at their code later. 

For this to work, I have to deal with the fact that when I copy a given junction to a new feature class to use as a starting point, when I perform an upstream trace, the trace tends to also include the edge just downstream of that junction.  Before running the code below, I precompute length downstream to the network outlet for each junction, and also from the downstream end of each edge.  If an edge has a length downstream value less than the junction, then it is downstream from the junction.  I only want upstream edges.

Updated_Trace_Network = arcpy.tn.Trace(
in_trace_network=network,
trace_type='UPSTREAM',
starting_points=start_pt,
validate_consistency='DO_NOT_VALIDATE_CONSISTENCY', ## proceed even if we've updated other features
condition_barriers=[['HydroID', 'DOES_NOT_EQUAL', 'SPECIFIC_VALUE', hydro_id, 'OR']],
traversability_scope='JUNCTIONS_ONLY',
result_types=['AGGREGATED_GEOMETRY'],
trace_name=str(hydro_id),
aggregated_lines=agg_lines)[0]

if Updated_Trace_Network:
# Select the edges corresponding to the aggregated trace result
lyr = arcpy.SelectLayerByLocation_management(
edge_layer, 'SHARE_A_LINE_SEGMENT_WITH', agg_lines)[0]
edge_count = int(arcpy.GetCount_management(lyr)[0])
if edge_count > 0:
# Sometimes downstream edges are selected. Remove them.
arcpy.SelectLayerByAttribute_management(
lyr, 'SUBSET_SELECTION', 'LengthDown >= ' + str(length_down))
edge_count = int(arcpy.GetCount_management(lyr)[0])
if edge_count > 0:
arcpy.AddMessage('Upstream edge count: ' + str(edge_count))
arcpy.CalculateField_management(lyr, 'JunctionID', hydro_id)
else:
arcpy.AddMessage('No upstream edges found')
else:
arcpy.AddMessage('No upstream edges found')
0 Kudos
TimWhiteaker
Occasional Contributor II

I figured out how to prevent the downstream edge from being returned in the upstream trace.  My mistake was simply copying my selected junction to a separate flag feature class.  To make this work, I also have to add and populate a couple of fields that you'll notice in the TN_Temp_Starting_Points feature class in your default geodatabase that stores the flags you create interactively (via Data tab on ribbon), namely, SOURCEID and FEATUREGLOBALID.  The SOURCEID refers to the identifier given to the junction layer in the trace network.  The FEATUREGLOBALID is the GLOBALID of the junction feature where you want the starting point placed.  I wish this information had been in a white paper that I had access to.

Here's code that demonstrates how to trace upstream given a junction layer.

import arcpy


def get_source_id(network, network_layer):
    d = arcpy.Describe(network_layer)
    fc_name = d.featureClass.name
    
    d = arcpy.Describe(network)
    for s in d.sources:
        if s.name == fc_name:
            return s.sourceID

    arcpy.AddError(
        'Could not determine network source ID for ' + network_layer.name)
    raise arcpy.ExecuteError


def point_to_flag(shape_xy, spatial_ref, source_id, global_id):
    flag_fc = arcpy.CreateFeatureclass_management(
        out_path='memory',
        out_name='flag',
        geometry_type='POINT',
        has_m='ENABLED',
        has_z='ENABLED',
        spatial_reference=spatial_ref)[0]
    flag_fc = arcpy.management.AddField(flag_fc, 'SOURCEID', 'LONG')[0]
    flag_fc = arcpy.management.AddField(flag_fc, 'FEATUREGLOBALID', 'GUID')[0]
    fields = ['SHAPE@XY', 'SOURCEID', 'FEATUREGLOBALID']
    with arcpy.da.InsertCursor(flag_fc, fields) as cursor:
        cursor.insertRow((shape_xy, source_id, global_id))
    return flag_fc


arcpy.env.overwriteOutput = True

network = arcpy.GetParameter(0)
junction_layer = arcpy.GetParameter(1)

spatial_ref = arcpy.Describe(junction_layer).spatialReference
source_id = get_source_id(network, junction_layer)

with arcpy.da.SearchCursor(junction_layer, ['SHAPE@XY', 'GLOBALID']) as cursor:
    for i, row in enumerate(cursor):        
        flag = point_to_flag(row[0], spatial_ref, source_id, row[1])

        Updated_Trace_Network = arcpy.tn.Trace(
            in_trace_network=network,
            trace_type='UPSTREAM',
            starting_points=flag)[0]

arcpy.SetParameter(2, junction_layer)

 

So, until Pro 2.7 is released in which the trace tools return selection layers, the keys to iterating traces are:

  • Create a point feature class for your starting point(s), and populate SOURCEID and  FEATUREGLOBALID attributes in that feature class.
  • Return aggregated geometry, and use that geometry to select features from your network layer(s).
  • If you're updating a network layer's attributes while iterating, when you call arcpy.tn.Trace remember to set validate_consistency='DO_NOT_VALIDATE_CONSISTENCY' so the tool executes even as you update features.
  • Run arcpy.ValidateNetworkTopology_tn when you need to validate the network after updating it.

I've implemented this via Python. It seems like it should also work via ModelBuilder iteration, but ModelBuilder iteration confuses me, so I haven't verified it works there.  My id_to_edges.py script, and the utils.py script it relies on, are attached in case folks find it helpful.

View solution in original post

0 Kudos
GintautasKmieliauskas
Regular Contributor

Hi Tim,

I use ArcGIS Pro 2.7.2. In ArcGIS Pro Trace geoprocessing tool works as expected and selects some network elements. I have copied Python script and placed it in my application, but I don't know how to get selection from it. I have checked all results of arcpy.tn.Trace and my selection must be in element number 2.  My python script:

    result = arcpy.tn.Trace(traceNetworkDataset, "UPSTREAM", startPointsFC, None, "NO_DIRECTION", '', "EXCLUDE_BARRIERS",
                            "DO_NOT_VALIDATE_CONSISTENCY", "DO_NOT_IGNORE_BARRIERS_AT_STARTING_POINTS", "IGNORE_INDETERMINATE_FLOW",
                            None, None, "BOTH_JUNCTIONS_AND_EDGES", None, None, "SELECTION", "NEW_SELECTION", "CLEAR_ALL_PREVIOUS_TRACE_RESULTS", '',
                            "Trace_Results_Aggregated_Points", "Trace_Results_Aggregated_Lines", None, "DO_NOT_USE_TRACE_CONFIGURATION", '')

Your scripts don't use selection but maybe you know how to get it?

 

0 Kudos