Remove unconnected polylines from a shapefile ArcPy

Question

How do you isolate a group of polylines that are not connected to the bigger polyline group.

Ideally, I would like to rank all the polyline groups based on how many features there are in each group.

The arrow in red shows the smaller polyline group that has to be isolated.

I have used,

arcpy.SelectLayerByLocation_management(temp, "BOUNDARY_TOUCHES", temp, "", "NEW_SELECTION", "INVERT")

But, this only works if the "group" has only one polyline.

I could also use the code below but it takes hours to run

arcpy.SelectLayerByAttribute_management("temp", "NEW_SELECTION", "FID = 100")
arcpy.SelectLayerByLocation_management("temp", "BOUNDARY_TOUCHES", temp, "", "ADD_TO_SELECTION")
arcpy.GetCount_management("temp")

Answer 1

I would solve the problem with NetworkX. It’s an Python package for the creation, manipulation, and study of the structure, dynamics, and functions of complex networks. It supports import and export of shapefile format. Install NetworkX in your ArcGIS Python environment: pip install networkx

The following code should give you one possibility to solve your problem: First, add an unique ID Field (e.g. "MY_FID") to your polyline shapefile. Than you can build a python dictionary, where the key is the ID field and the value is the number of the polyline group.

import arcpy
import networkx as nx

polyline_shapefile = r"D:\Test\polylines.shp"

# Add id and group fields
arcpy.AddField_management(polyline_shapefile, "GROUP", "LONG")
arcpy.AddField_management(polyline_shapefile, "MY_FID", "LONG")
oid_fieldname = arcpy.Describe(polyline_shapefile).OIDFieldName
arcpy.CalculateField_management(polyline_shapefile, "MY_FID", "!{oid}!".format(oid=oid_fieldname), "PYTHON_9.3")

polyline_group_fid_dict = {}
polyline_graph = nx.read_shp(polyline_shapefile, simplify=False).to_undirected()

# loop all connected components (polyline groups) of the polyline network graph 
for component_number, polyline_group in enumerate(nx.connected_component_subgraphs(polyline_graph)):
    for polyline in polyline_group.edges():
        fid = polyline_group.get_edge_data(polyline[0], polyline[1])["MY_FID"]
        polyline_group_fid_dict[fid] = component_number

# add the polyline group id to the "GROUP" field.
with arcpy.da.UpdateCursor(polyline_shapefile, ["MY_FID", "GROUP"]) as cursor:
    for row in cursor:
        row[1] = polyline_group_fid_dict[row[0]]
        cursor.updateRow(row)

Now you can count the lines per group or dissolve the polylines by the "GROUP" field and calculate the length of the polyline groups.

Answer 2

I have a function that comes really close to this. Instead of attributing with the count, it attributes with a unique "neighborhood" ID. It could easily be modified to report the count instead. Or, you could run this first and then get the counts with a simple cursor and a dictionary.

To directly attribute with counts, replace all instances of

row[0] = id_out

with:

row[0] = count_post

Also, the function below has a few additional parameters to fine-tune your groupings of features, but they may be superfluous for your needs.

I edited it to remove some dependencies, and I haven't tested the edits, but it should be really close.

import os
import datetime
import arcpy

def assign_neighborhood_id(feature_class, neighborhood_field, dissolve_fields=None, search_distance='0 FEET', where_clause=None):
    """
    (str, [str], [collection], [str]) -> str

    To each disjunct cluster of features, assign a unique neighborhood ID.

    Parameters:

    feature_class - The path/name of the feature to which you wish to assign
        neighborhood IDs

    neighborhood_field - The name of the field to which the IDs will be
        assigned.

    dissolve_fields - A field name or a collection of field names within which
        features must have identical values in order to be considered neighbors.

    search_distance - By default, neighbors will be identified by an
        intersection, however, if you wish to allow for gaps between features,
        you may enter an ArcGIS-compatible search distance here. By default or
        if None is passed, then '0 FEET' will be used.

    where_clause - Only features satisfying the where clause will be considered
        part of a neighborhood; all other features will have a Null value in the
        neighborhood field.
    """
    if not search_distance:
        search_distance = '0 FEET'

    if neighborhood_field not in [i.name for i in arcpy.ListFields(feature_class)]:
        arcpy.AddField_management(feature_class, neighborhood_field, 'LONG')

    already_layer = arcpy.Describe(feature_class).dataType == 'FeatureLayer'
    if (already_layer and where_clause) or not already_layer:
        # Ensure a unique, unused name for the in-memory layer
        lyr = os.path.join('in_memory', 'f{}'.format(datetime.datetime.now().strftime('%Y%m%d%H%M%S%f')))
        arcpy.MakeFeatureLayer_management(feature_class, lyr, where_clause)
    else:
        lyr = feature_class

    fn_oid = arcpy.Describe(feature_class).OIDFieldName
    vs_oid = set([row[0] for row in arcpy.da.SearchCursor(lyr, fn_oid)])

    select_type = 'INTERSECT'
    id_out = int()

    # If there are fields on which to dissolve
    if dissolve_fields:
        # Get a list of the names of string/date fields
        fns_string = [i.name for i in arcpy.ListFields(lyr) if i.type in ('String', 'Date')]
        # Compile a list of where clauses, one for each field
        wcs_components = list()
        for index, field in enumerate(dissolve_fields):
            if field in fns_string:
                wc_base = "{} = '{}'"
            else:
                wc_base = '{} = {}'
            wcs_components.append(wc_base.format(field, '{}'))
        # Compile the full where clause, but without values
        wc_diss = ' AND '.join(wcs_components)

        for v_oid in vs_oid:
            wc_seed = '{} = {}'.format(fn_oid, v_oid)
            # If the feature is already assigned to a neighborhood, skip it.
            if [x[0] for x in arcpy.da.SearchCursor(lyr, neighborhood_field, wc_seed)][0] is not None:
                continue
            # Create the where clause for querying suitable neighbors
            with arcpy.da.SearchCursor(lyr, dissolve_fields, wc_seed) as cur:
                for row in cur:
                    vs_seed = row
                    break
            if None in vs_seed:
                wcs_components_null = list()
                for index, v in enumerate(vs_seed):
                    if v is None:
                        wcs_components_null.append('{} IS NULL'.format(dissolve_fields[index]))
                    else:
                        wcs_components_null.append(wcs_components[index].format(v))
                wc_diss_use = ' AND '.join(wcs_components_null)
            else:
                wc_diss_use = wc_diss.format(*vs_seed)

            # Run the initial selection
            arcpy.SelectLayerByAttribute_management(lyr, None, wc_seed)
            count_initial = 1
            arcpy.SelectLayerByLocation_management(lyr, select_type, lyr, search_distance)
            arcpy.SelectLayerByAttribute_management(lyr, 'SUBSET_SELECTION', wc_diss_use)
            count_post = len(arcpy.Describe(lyr).FIDSet.split(';'))

            while count_initial < count_post:
                # Update the initial count
                count_initial = count_post
                # Select intersecting
                arcpy.SelectLayerByLocation_management(lyr, select_type, lyr, search_distance)
                # Remove those that don't match the attributes
                arcpy.SelectLayerByAttribute_management(lyr, 'SUBSET_SELECTION', wc_diss_use)
                # Update the current select count and the iterations count
                count_post = len(arcpy.Describe(lyr).FIDSet.split(';'))

            # Assign the IDs
            with arcpy.da.UpdateCursor(lyr, neighborhood_field) as cur:
                for row in cur:
                    row[0] = id_out
                    cur.updateRow(row)
            arcpy.SelectLayerByAttribute_management(lyr, 'CLEAR_SELECTION')
            # Increment the ID
            id_out += 1

    # If no dissolve fields are required
    else:
        for v_oid in vs_oid:
            wc_seed = '{} = {}'.format(fn_oid, v_oid)
            # If the feature is already assigned to a neighborhood, skip it.
            if [x[0] for x in arcpy.da.SearchCursor(lyr, neighborhood_field, wc_seed)][0] is not None:
                continue
            # Run the initial selection
            arcpy.SelectLayerByAttribute_management(lyr, None, wc_seed)
            count_initial = 1
            arcpy.SelectLayerByLocation_management(lyr, select_type, lyr, search_distance)
            count_post = len(arcpy.Describe(lyr).FIDSet.split(';'))

            while count_initial < count_post:
                # Update the initial count
                count_initial = count_post
                # Select intersecting
                arcpy.SelectLayerByLocation_management(lyr, select_type, lyr, search_distance)
                # Update the current select count and the iterations count
                count_post = len(arcpy.Describe(lyr).FIDSet.split(';'))

            # Assign the IDs and increment the ID
            with arcpy.da.UpdateCursor(lyr, neighborhood_field) as cur:
                for row in cur:
                    row[0] = id_out
                    cur.updateRow(row)
            arcpy.SelectLayerByAttribute_management(lyr, 'CLEAR_SELECTION')
            id_out += 1

    return neighborhood_field