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])
Improve this question
7
  • Did you try the method explained here? gis.stackexchange.com/questions/26369/…
    – Richard Morgan
    Jun 12, 2019 at 17:19
  • I have updated the code in my sample to add more at the beginning, if that is helpful.
    – AttilaTheHun
    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.
    – AttilaTheHun
    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.
    – Richard Fairhurst
    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.
    – Richard Fairhurst
    Jun 12, 2019 at 18:34

1 Answer 1

Reset to default
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
Improve this answer

Your Answer

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

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