1

I am trying to speed up some current code that I have for looping through lines, getting their segments and ultimately the vertex points in Python.

Can anyone recommend a faster method than what I currently have (perhaps using a map function or dictionary)?

My code is very slow, especially when dealing with large numbers of features.

    sup_list = [row for row in arcpy.da.SearchCursor(sup_lyr, fields)]

    # with arcpy.da.SearchCursor(sup_lyr, fields) as sup_cursor:
    for sup_row in pbar(sup_list, total=result, ascii=True):
        exp = str(tracking_field) + " = {0}".format(sup_row[3])
        sup_sel = arcpy.SelectLayerByAttribute_management(sup_lyr, "NEW_SELECTION", exp)
        inf_sel = arcpy.SelectLayerByLocation_management(inf_lyr, "WITHIN_A_DISTANCE_GEODESIC", sup_sel, "50 Meters"
                                                     , "NEW_SELECTION", "NOT_INVERT")
        sup_vertexList = []
        sup_proj_geom = sup_row[1].projectAs(calc_utm(sup_row[1].firstPoint))
        for part in sup_proj_geom:
            for pnt in part:
                sup_vertexList.append([pnt.Y, pnt.X])
        with arcpy.da.SearchCursor(inf_sel, fields) as inf_cursor:
            for inf_row in inf_cursor:
                inf_vertexList = []
                inf_proj_geom = inf_row[1].projectAs(calc_utm(inf_row[1].firstPoint))
                for part in inf_proj_geom:
                    for pnt in part:
                        inf_vertexList.append([pnt.Y, pnt.X])
7
  • Did you try the method explained here? gis.stackexchange.com/questions/26369/… Jun 12, 2019 at 17:19
  • I have updated the code in my sample to add more at the beginning, if that is helpful. Jun 12, 2019 at 17:35
  • Richard, thank you for the link. I have looked over the code, but I don't see that it would be much faster. Jun 12, 2019 at 17:37
  • 1
    A Spatial Join using the Many To One option would probably make the whole thing faster with the layer of the inner loop being the target layer and the layer with the attribute selection criteria being the join layer. Then there would be a common attribute selection on the Spatial Join features that would match your outer loop layer and you could use a dictionary to match them. The spatial join and attribute dictionary matching would run hundreds of times faster than large sets of spatial selections, and it becomes increasingly efficient as you add more features. Jun 12, 2019 at 17:49
  • 2
    Using Desktop, the easiest way to get a code snippet for a tool is to open the tool like you normally would (Analysis Tools.tbx->Overlay-Spatial Join), set up the parameters the way you want (Target layer, Join Layer, Output Name, JOIN_MANY_TO_ONE, optionally keep only matched features for performance, adjust field list to just include fields you need for performance, WITHIN_A_DISTANCE_GEODESIC, 50 Meters) and run it to test. Then you can open the Results tab (Geoprocessing menu->Results), right click the tool and choose Copy As Python Snippet. Paste into Python and customize as needed. Jun 12, 2019 at 18:34

1 Answer 1

0

Here is a version of a function I regularly use:

def get_point_coords_from_lines(shp_in, target_epsg, line_id_fieldname, output_location):
    """
    Get vertices for lines with NODE_ID order field
    :shp_in:            input polyline shapefile
    :target_epsg:       desired EPSG for point coordinates
    :line_id_fieldname: fieldname for line identifier; must be unique, like a primary key
    :output_location:   folder path for output files
    """

    import arcpy

    # copy input SHP just in case...
    temp1 = output_location + '\\' + 'temp_lines_in.shp'
    arcpy.CopyFeatures_management(shp_in, temp1, config_keyword = "")

    temp2 = output_location + '\\' + 'temp_lines_in_proj.shp'
    arcpy.Project_management(temp1, temp2, arcpy.SpatialReference(target_epsg))
    arcpy.Delete_management(temp1)

    my_pts = output_location + '\\' + 'temp_proj_as_Points.shp'
    arcpy.FeatureVerticesToPoints_management(temp2, my_pts, point_location = "ALL")
    arcpy.Delete_management(temp2)
    arcpy.AddXY_management(my_pts)

    ### this gives a "node_id" aka a point order so that you can rebuild the 
    ### lines with this as the sort field and 'line_id_fieldname' as the group field
    arcpy.AddField_management(my_pts, "NODE_ID", "SHORT", field_is_nullable="NULLABLE")
    line_id_previous = arcpy.da.SearchCursor(my_pts, (line_id_fieldname,)).next()[0]
    counter = 0
    with arcpy.da.UpdateCursor(my_pts, ['OID@', line_id_fieldname, 'NODE_ID']) as cursor:
        for row in cursor:
            line_id_current = row[1]
            if line_id_current == line_id_previous:
                counter += 1
            else:
                counter = 1
            row[2] = counter
            line_id_previous = line_id_current
            cursor.updateRow(row)

    return my_pts

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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