5

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:

  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. 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[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]
            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", "")
1

Apologies for not understanding everything in your code, but I suggest a different method in general. Mostly it is to use GetCellValue_management. Just kind of winging this, so I'm sure there will be some kinks to work out. I was having a tough time with the formatting of the pseudo code, too...

I won't do the whole thing since that would take too much time on my end.

  1. get 2 xy pairs for each line. I'm guessing start and end nodes, but you might have curved lines or some other basis. You'll also need to define which direction you want to go. Could use a SearchCursor() to extract these values to a nice, simple array.
  2. figure out the increment you want to extend by (you said 30m) and turn that into adjustment increments for the "x" and "y" dimensions. I'm guessing you'll want to figure out the slope of the line defined by the two pairs of XY coords. Store in another array or list.
  3. use GetCellValue_management to find start and end elevations for all lines
  4. probably want to izip the elevations with the ending coordinates for iteration. let's say the fields are lineEnd_x, LineEnd_y, "start" elevation, "end" elevation. Could also include the adjustment factors in what you're zipping together for convenience. Let's call this lines.
  5. some pseudo code:

    for line in lines: 
        if line[2] - line[3] < threshold: # 20 meters or whatever, also might want abs()
          continue
        else: # enclosing your while loop.
            while (line[2] - line[3]) < threshold:
                line[0] = line[0] + xadj
                line[1] = line[1] + yadj
                line[3] = arcpy.GetCellValue_management(dem, arcpy.Point(line[0], line[1]))  
    
  6. iterate over updated lines array to extend lines. Could do this a number of ways, from creating an entirely new set of vectors to extending the old ones with something like ExtendLine_edit (would have to comparing new XY coords with the initial ending XY coords to get the offset distance).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.