Extend XY Line to Z Value from DEM

[AndrewMellick]

Hi All,

Hopefully I can get some help....

I have a series of lines and a single raster DEM.

I'd like to extend my line segments until their Z-component has changed by at least 20 metres. My initial thoughts on this were:

1.Write Z_Min and Z_Max for the line

2.Calculate Z_Diff = Z_Max - Z_Min, for a check

3.If not, extend the line 30m (random number), check, if not, extend, check, etc.

I've managed to get some code that will run through the process but only stops once ALL features are satisfied (not each individually). I'd like to iterate through each feature and do the check (or check all, extend those that don't satisfy, etc.)

My attempts at querying the datasets via using a Where statement within UpdateCursor function statement were unsucessful - I'm fairly certain I was messing with the created tuples of coordinates (the tuple of OID@ and SHAPE@XY) but I'm not certain.

Current Code - loops through just once and despite my best efforts - more or less resembles some sample script already posted on another website. If anybody has suggestions/samples, that would VERY much be appreciated.

Is this a simple "While" Loop?

#Code from <http://gis.stackexchange.com/questions/71645/a-tool-or-way-to-extend-line-by-specified-distance> by Paul.

from math import hypot
import collections
from operator import add
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("3D")

#Local variables:
EnfDEM_Bird = "EnfDEM_Bird"

#Need to set these as imported Script for ArcGIS
layer = arcpy.GetParameterAsText(0)
distance = float(arcpy.GetParameterAsText(1))

# Process: Add Surface Information
arcpy.AddSurfaceInformation_3d(layer, EnfDEM_Bird, "Z_MIN;Z_MAX;SURFACE_LENGTH", "BILINEAR", "", "1", "0", "NO_FILTER")
# Process: Calculate Field
arcpy.CalculateField_management(layer, "Z_Diff", "[Z_Max] - [Z_Min]", "VB", "")

#Computes new coordinates x3,y3 at a specified distance
#along the prolongation of the line from x1,y1 to x2,y2
def newcoord(coords, dist):
    (x1,y1),(x2,y2) = coords
    dx = x2 - x1
    dy = y2 - y1
    linelen = hypot(dx, dy)

    x3 = x2 + dx/linelen * dist
    y3 = y2 + dy/linelen * dist    
    return x3, y3

#accumulate([1,2,3,4,5]) --> 1 3 6 10 15
#Equivalent to itertools.accumulate() which isn't present in Python 2.7
def accumulate(iterable):    
    it = iter(iterable)
    total = next(it)
    yield total
    for element in it:
        total = add(total, element)
        yield total

#OID is needed to determine how to break up flat list of data by feature.
coordinates = [[row[0], row[1]] for row in
               arcpy.da.SearchCursor(layer, ["OID@", "SHAPE@XY"], explode_to_points=True)]

oid,vert = zip(*coordinates)

#Construct list of numbers that mark the start of a new feature class.
#This is created by counting OIDS and then accumulating the values.
vertcounts = list(accumulate(collections.Counter(oid).values()))

#Grab the last two vertices of each feature
lastpoint = [point for x,point in enumerate(vert) if x+1 in vertcounts or x+2 in vertcounts]

#Convert flat list of tuples to list of lists of tuples.
#Obtain list of tuples of new end coordinates.
newvert = [newcoord(y, distance) for y in zip(*[iter(lastpoint)]*2)]    

j = 0

with arcpy.da.UpdateCursor(layer, ["SHAPE@XY"], explode_to_points=True) as rows:
    for i,row in enumerate(rows):
        if i+1 in vertcounts:    
            row[0] = newvert
            j+=1
            rows.updateRow(row)
            # Process: Add Surface Information
            arcpy.AddSurfaceInformation_3d(layer, EnfDEM_Bird, "Z_MIN;Z_MAX;SURFACE_LENGTH", "BILINEAR", "", "1", "0", "NO_FILTER")
            # Process: Calculate Field
            arcpy.CalculateField_management(layer, "Z_Diff", "[Z_Max] - [Z_Min]", "VB", "")