Long time reader, first time poster. So here's my goal:
Working in ArcGIS, trying to work with Model Builder and Python Scripts.
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:
Write Z_Min and Z_Max for the line
Calculate Z_Diff = Z_Max - Z_Min, for a check
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. 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: DEM = "DEM" #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, row] 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 = newvert[j] 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", "")