Get specific vertex coordinates from polyline geometry using indexing with Python?

[JoeBryant1]

Hi, I'm trying create 2 new line feature classes based on just the first vertex to vertex line segment and last vertex to vertex segment of an existing line feature class. This is for calculating upstream and downstream gravity main horizontal angles for manhole turbulence modelling. Some of the existing lines are a single straight segment (just 2 vertices), but most have multiple vertices. I've already densified the true curves and none of the lines are multipart. If the new upstream and downstream line segments are identical for straight lines, that is ok. I will be joining the calculated horizontal angle back to new upstream and downstream horizontal angle fields in the original gravity main line feature class.

How do I return the coordinate geometry of the first and second vertex in one pass, and then the last and second to last vertex in another pass using Python in the ArcPy pane in ArcGIS Pro (3.1). Everything I've seen only returns the first and last, or the centroid coordinates of the feature. Can I use some sort of indexing of the vertices (i.e: [0] and [1], then [-1] and [-2])?

Here is the code I have so far that returns ALL of the vertex coordinates for each row and part. I need to create the upstream and downstream array objects from the 2 coordinates at each end.

# Python script to create upstream and downstream line segments
# based on the first and last segment of lines with multiple vertices

import arcpy

# Set your workspace and feature class paths
workspace = r"C:\path\to\your\geodatabase.gdb"
input_feature_class = "your_input_feature_class"
output_upstream_feature_class = "upstream_lines"
output_downstream_feature_class = "downstream_lines"

# Set Unique ID field name from input feature class
ID = "CSDFeatureID"

# Get the spatial reference of the input feature class
spatial_reference = arcpy.Describe(input_feature_class).spatialReference

# Create feature classes with the same spatial reference
arcpy.CreateFeatureclass_management(workspace, output_upstream_feature_class, "POLYLINE", spatial_reference=spatial_reference)
arcpy.CreateFeatureclass_management(workspace, output_downstream_feature_class, "POLYLINE", spatial_reference=spatial_reference)

# Add Facility ID field to the result feature classes
arcpy.AddField_management(output_upstream_feature_class, ID, "TEXT")
arcpy.AddField_management(output_downstream_feature_class, ID, "TEXT")

# Start an edit session
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()

# Open a search cursor to iterate through the input feature class
with arcpy.da.SearchCursor(input_feature_class, [ID, "SHAPE@"]) as cursor:
    for row in cursor:
        print("Line Feature {}:".format(row[0]))
        for part in row[1]:
            for pnt in part:
                print("Vertex coordinates: {}, {}".format(pnt.X, pnt.Y))
                
                # Get the first two vertex coordinates and store in upstream_array
                # Get the last two vertex coordinates and store in downstream_array
        
        # Create new features from the upstream and downstream array coordinates
        with arcpy.da.InsertCursor(output_upstream_feature_class, ["SHAPE@", ID]) as up_cursor:
            up_line = arcpy.Polyline(upstream_array)
            up_cursor.insertRow((up_line, line_id))

        with arcpy.da.InsertCursor(output_downstream_feature_class, ["SHAPE@", ID]) as down_cursor:
            down_line = arcpy.Polyline(downstream_array)
            down_cursor.insertRow((down_line, line_id))

# Stop the edit session
edit.stopOperation()
edit.stopEditing(True)

print("Script completed.")

[DanPatterson]

a = np.array([[0., 0.], [2., 10.], [4., 10.], [2., 0.], [0., 0.]])

a[[0, 1]] 
array([[  0.00,   0.00],
       [  2.00,  10.00]])

a[[-2, -1]] 
array([[  2.00,   0.00],
       [  0.00,   0.00]])

# -- two point array/line
b = np.array([[0., 0.], [2., 10.]])

b[[0, 1]]
array([[  0.00,   0.00],
       [  2.00,  10.00]])

a[[-2, -1]] 
array([[  2.00,   0.00],
       [  0.00,   0.00]])

See if "part" supports fancy indexing, so rather than slicing one at a time, you can get both

Or you could take a different tact and convert the searchcursor row to an array

def _fc_as_narray_(in_fc, with_id=False):
    """Return geometry from a featureclass using `as_narray`."""
    flds = ["SHAPE@X", "SHAPE@Y"]
    if with_id:
        flds = ["OID@", "SHAPE@X", "SHAPE@Y"]
    with arcpy.da.SearchCursor(in_fc, flds, explode_to_points=True) as cursor:
        a = cursor._as_narray()
    del cursor
    return a

then 'a' will be your polyline as an array

[DuncanHornby]

Have you tried the indexing which you suggest as that is what I would have done?

[JoeBryant1]

I have not been able to figure out the syntax for retrieving the vertex points in a polyline in the proper format to index them and get their XY coordinates. I will try to adapt Dan's example above when I'm able. It sounds like converting the Search Cursor row to an array may be what I'm looking for? Then I have to index not only the "part", but the vertex? I almost need a flowchart of the hierarchy I'm trying to follow down to the coordinate pair level.

[DanPatterson]

if your polylines are singlepart (but can contain multiple segments) then what I suggested is easiest.

Don't confuse singlepart polylines are polylines that contain a single or multiple segments.

A multipart polyline is a polyline with disjoint geometry represented by a single record in the attribute table.

The slicing example works for on a singlepart polyline with one or many segments

[DuncanHornby]

This should get you going, this is how I would have done it.  Code assumes single part and a line with a minimum of 4 vertices:

import arcpy
lyrLine = "lines"
lstRes = list()
with arcpy.da.SearchCursor(lyrLine, "SHAPE@") as cursor:
    for row in cursor:
        geom = row[0]  # A polyline
        if geom.pointCount >= 4:
            arr =geom.getPart(0)
            p1 = arr[0] # first
            p2 = arr[1] # second
            p3 = arr[-2] # second from last
            p4 = arr[-1] # last
            s1 = str(p1.X) + "," + str(p1.Y)
            s2 = str(p2.X) + "," + str(p2.Y)
            s3 = str(p3.X) + "," + str(p3.Y)
            s4 = str(p4.X) + "," + str(p4.Y)
        lstRes.append([s1,s2,s3,s4])
print(lstRes)

[DanPatterson]

as in my example, slicing using 0, 1 and -2, -1 only requires 2 points

[JoeBryant1 · Accepted solution]

# Python script to create upstream and downstream line segments
# based on the first and last segment of lines with multiple vertices

import arcpy

# Set your workspace and feature class paths
workspace = r"C:\path\to\your\geodatabase.gdb"
input_feature_class = "your_input_feature_class"
output_upstream_feature_class = "upstream_lines"
output_downstream_feature_class = "downstream_lines"

# Set Unique ID field name from input feature class
ID = "CSDFeatureID"

# Get the spatial reference of the input feature class
spatial_reference = arcpy.Describe(input_feature_class).spatialReference

# Create feature classes with the same spatial reference
# arcpy.CreateFeatureclass_management(workspace, output_upstream_feature_class, "POLYLINE", spatial_reference=spatial_reference)
# arcpy.CreateFeatureclass_management(workspace, output_downstream_feature_class, "POLYLINE", spatial_reference=spatial_reference)

# Add Facility ID field to the result feature classes
# arcpy.AddField_management(output_upstream_feature_class, ID, "TEXT")
# arcpy.AddField_management(output_downstream_feature_class, ID, "TEXT")

# Start an edit session
edit = arcpy.da.Editor(workspace)
edit.startEditing(False, True)
edit.startOperation()

# Open a search cursor to iterate through the input feature class
with arcpy.da.SearchCursor(input_feature_class, ["SHAPE@", ID]) as cursor:
    for row in cursor:
        line_id = row[1] # ID of line
        geom = row[0]  # A polyline geometry object
        arr =geom.getPart(0) # the vertices as an array of arrays of points
        p1 = arr[0] # first vertex as an array of points
        p2 = arr[1] # second
        p3 = arr[-2] # second from last
        p4 = arr[-1] # last
        
        # Create arrays to hold vertex coordinates
        upstream_array = arcpy.Array([p1, p2])
        downstream_array = arcpy.Array([p3, p4])
        
        # Create new features from the upstream and downstream arrays
        with arcpy.da.InsertCursor(output_upstream_feature_class, ["SHAPE@", ID]) as up_cursor:
            up_line = arcpy.Polyline(upstream_array)
            up_cursor.insertRow((up_line, line_id))

        with arcpy.da.InsertCursor(output_downstream_feature_class, ["SHAPE@", ID]) as down_cursor:
            down_line = arcpy.Polyline(downstream_array)
            down_cursor.insertRow((down_line, line_id))

# Stop the edit session
edit.stopOperation()
edit.stopEditing(True)

print("Script completed.")

Thanks for the input @DanPatterson and @DuncanHornby  I was able to achieve success today. To answer my original question, indexing from the end does work with arcpy.Array.

Instead of using NumPy as Dan suggests, I used the arcpy.Array class as Duncan shows, but without the 4 point requirement. I still had trouble with which elements needed to be arrays versus lists (I could still benefit from some general object hierarchy education on this, as I indicated above). There was no need to convert to string or create a list. My main takeaway is that an arcpy geometry object for a polyline is an array of arrays of point coordinate pairs. And you need an array to create the new polylines, not a list.

Here is the script I got to work:

[DanPatterson]

In case someone wants to use numpy

d0_  # -- a shape converted to an array as in my previous post
array([[  0.000,   5.000],
       [  3.000,   7.000],
       [  5.000,  10.000],
       [  7.000,   7.000],
       [ 10.000,   5.000],
       [  7.000,   3.000],
       [  5.000,   0.000],
       [  3.000,   3.000],
       [  0.000,   5.000]])

# -- slice what you want from the shape and add to a list
subs = [d0_[:2], d0_[-2:]]

# -- some arcpy imports
from arcpy import Array, Point, Polyline

# -- construct the polyline
aa = []
for sub in subs:
    aa.append([Point(*pairs) for pairs in sub])
poly = Polyline(Array(aa), SR)

poly[0]  # -- view the points in the array
<Array [<Point (0.0, 5.0, #, #)>, <Point (3.0, 7.0, #, #)>]>

now snag the output in Spyder as an svg