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

Question

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
        # arcpy.AddField_management(outFC, 'WS_Elev', 'FLOAT')

        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
    arcpy.AddZInformation_3d(outFC, 'Z')

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

print ('Finished!')

Answer 1

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.

Answer 2

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.