# Create points along 3D stream where elevation changes by 1ft

Using ArcPy (Python 3), I am trying to create points along a line where the elevation changes by 1 ft between vertices. I can successfully convert my Stream Polyline feature class to 3D and create a Point feature class where the elevation changes every 1 ft. I now need to add the Z elevation from the vertex to the Point feature class. Can I do by using an Update Cursor within the Insert cursor? I'm still pretty new with Python, what I remember from class is that it's bad practice to create nested cursors? What other approach could I use to do this?

Update: Was able to figure it out. If I create point feature class that is Z aware the point vertices will have the elevations. Then I can use the AddZInformation tool to transfer those values to the Point's attribute table.

Update: Now I am trying to create points every 1000 feet between my existing points if no points currently exist. How can I do this? I'm thinking the Generate Points Along Line tool might work. https://pro.arcgis.com/en/pro-app/tool-reference/data-management/generate-points-along-lines.htm

``````import arcpy
import os
from arcpy import env

arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("3D")
arcpy.CheckOutExtension("Spatial")
env.workspace = r"C:\GIS\Tools\BFE_Tool\BFE_Tool.gdb"

streams = 'FY16_Streams_Diss' # Streams must be dissolved by stream name
outGDB = r'C:\GIS\Tools\BFE_Tool\BFE_Tool.gdb'

print ('Interpolating shape...')
stream3D = arcpy.InterpolateShape_3d(r'C:\GIS\Tools\BFE_Tool\Test_Data\wsg100x', streams, 'FY16_Streams3D')
print ('Interpolate shape done!')

print ('Creating BFE Points...')
try:

myLayer = stream3D
prev = -1
heightDiff = 1
pntList = list()
with arcpy.da.SearchCursor(myLayer,["SHAPE@"]) as cursor:
for row in cursor:
geom = row[0][0] # A whole polyline
for i,point in enumerate(geom):
if prev == -1:
prev = point.Z
else:
if abs(point.Z - prev) > heightDiff:
pntList.append(geom[i-1])
pntList.append(geom[i])
prev = point.Z

if len(pntList) > 0:
print('Creating Empty point feature class...')
outFC = arcpy.CreateFeatureclass_management(outGDB, 'BFE_Points', 'POINT','', '','ENABLED') # Create Empty Feature Class to store Points

# print ('Adding Elevation field to Points...') # Add field in Points to store Water Surface Grid Elevations

with arcpy.da.InsertCursor(outFC, ["SHAPE@"]) as cursor:
for point in pntList:
print(point.Z)
cursor.insertRow([point])
print ('Point Created')

print ('Adding Elevation to Attribute table...') # Add field in Points to store Water Surface Grid Elevations

except arcpy.ExecuteError:
# Geoprocessor threw an error
print(arcpy.GetMessages(2))

print ('Finished!')
``````
• "SHAPE@Z" is one z value so you cant iterate over it. You probably want SHAPE@
– Bera
Commented May 28, 2020 at 13:18
• Thanks! That worked, I was able to print the Z value for each vertex @BERA Commented May 28, 2020 at 13:21
• Now struggling to figure out how to compare each Z value to create a point if the elevation changes by 1 ft. Somehow flag the vertices where the point needs to be. Then run Feature Vertices to points using a SQL Expression to filter those flagged vertices Commented May 28, 2020 at 13:29
• You want to insert additional vertices? Or just create a Point feature class? How important is the 1 feet limit? For example can it be 0.8-1 feet?
– Bera
Commented May 28, 2020 at 13:32
• I would like to create a new point feature class where the difference is 1 ft or greater. So for example if the values go: 400, 400.5, 401, 403, 403.2 then a new point should be created at the 401 and 403 vertices Commented May 28, 2020 at 13:51

The following code will insert into an existing point featureclass the points where a Z change of greater than 1 unit occurs based upon the specifications you state in your comments:

``````import arcpy
try:

myLayer = "fcPolylineZ"
prev = -1
heightDiff = 1
pntList = list()
with arcpy.da.SearchCursor(myLayer,["SHAPE@"]) as cursor:
for row in cursor:
geom = row[0][0] # A whole polyline
for i,point in enumerate(geom):
if prev == -1:
prev = point.Z
else:
if (point.Z - prev) > heightDiff:
pntList.append(geom[i-1])
pntList.append(geom[i])
prev = point.Z

if len(pntList) > 0:
outFC = r"C:\Scratch\fGDB_Scratch.gdb\fcPoints"   # This must exist already
with arcpy.da.InsertCursor(outFC , ["SHAPE@"]) as cursor:
for point in pntList:
cursor.insertRow([point])
except arcpy.ExecuteError:
# Geoprocessor threw an error
print(arcpy.GetMessages(2))
``````

The assumption of this algorithm is that Z is increasing with direction of flow of line, which is of cause the exact opposite of water flowing downhill. So if you can't guarantee flow direction of your line then it would be an absolute difference of 1 you need to test for.

This script is designed to be run from the analysis window in ArcMap and makes a whole load of assumptions about your Z values. For example there are always a Z value for each vertex, so no `NaN`. But the code should be enough to get you going.

• Awesome! I can get this to run, but yes to your point, the Z elevations actually travel in the opposite direction of the flow line. I cant guarantee the flow direction so it will need to be absolute difference of 1 Commented May 28, 2020 at 15:53

SHAPE@Z —A double of the feature's z-coordinate.

So it is one z value which you cant iterate over. It's like trying to do:

``````for z in 123.456:
print(z)
builtins.TypeError: 'float' object is not iterable
``````

You probably want

SHAPE@ —A geometry object for the feature.

to later access z of each vertice.