I think that the original poster may want to calculate for each cell the height above the first stream cell that would be reached by water flowing from the cell. So the 'nearest stream' is calculated along the downslope flow path, not euclidian distance. The references for this Height Above Nearest Drainage (HAND) are:
Rennó, C. D., Nobre, A. D., Cuartas, L. A., Soares, J. V., Hodnett, M.
G., Tomasella, J. and Waterloo, M. J. (2008) HAND, a new terrain
descriptor using SRTM-DEM: Mapping terra-firme rainforest environments
in Amazonia. Remote Sensing of Environment 112, 3469-3481.
Nobre, A. D., Cuartas, L. A., Hodnett, M., Rennó, C. D., Rodrigues,
G., Silveira, A., Waterloo, M. and Saleska, S. (2011) Height Above the
Nearest Drainage - a hydrologically relevant new terrain model. J.
Hydrol. 404, 13-29.
My rather kludgy implementation was:
1) Create a flow direction raster and flow accumulation raster from the DEM. The help files can walk you through this.
2) Create a stream raster from a flow accumulation raster by setting a threshold value for what is considered a stream (in my case 250 m^2). Merge with a sink raster because Australia is full of hydrologically relevant closed depressions. This is the drainage raster so set all values to a HAND of 0.
3) Calculate a raster of height above the next downslope cell for the entire area (flow elevation difference). This will be used in lots of iterations but only needs to be calculated once. You need to store text files for each irregular neighbourhood direction, e.g. for flow direction 64:
DirCode64.txt -
3 3
0 1 0
0 0 0
0 0 0
The calculation is performed once for each determinate value in the flow direction raster (powers of 2). For the odd-ball values that are not powers of 2, I used the height above the minimum elevation cell in the surrounding 3x3 neighbourhood.
4) Calculate the height above drainage iterating out from the drainage lines as:
flow elevation difference + HAND of next downslope cell.
This is performed once for each flow direction, added to the HAND raster, then iterating until the number of null cells in the HAND raster stops changing. I included an escape if the number of iterations got too high and save the output periodically so I can restart if/when the process crashes. The saving seems to be a slow step so I don't do it every iteration.
Hope this is clear enough. I'm sure the code could be cleaned up but I stopped when it was working. Thanks for help from other threads that pointed me in the right direction. Here is my code:
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
env.workspace = "C:/Data/GIS_Data/DEM"
# Starting HAND raster with 0 for streams/sinks
outHandRaster = Raster("hnd20strsnk")
# use below for restarting iterations
#outHandRaster = Raster("hand20gh2")
inElevRaster = Raster("dtm20m")
inFDirRaster = Raster("fdir20m")
lstDirection = [1,2,4,8,16,32,64,128]
# Count the number of null cells in the initial stream raster
nullOutRaster = Con(IsNull(outHandRaster),1)
nullOutRaster.save("handNull")
cursor = arcpy.da.SearchCursor("handNull","Count")
nullCount = cursor.next()[0] # Count of null cells in outHandRaster
print "nullCount = ", nullCount
nullDif = 1 # anything but 0
# Calculation of floweldif raster – contains the elevation difference
# between each cell and the cell in the downslope direction.
# This block only needs to be calculated once for the area
for idx in lstDirection:
focalMaskFile = "C:/Data/GIS_Data/DEM/FocalStatNeighbor/" + "DirCode" + str(idx) + ".txt"
outElDifRaster = Con(inFDirRaster == idx, inElevRaster - FocalStatistics(inElevRaster, NbrIrregular(focalMaskFile), "MINIMUM"))
else:
# Calculate values for indeterminate flow direction cells
outElDifRaster = Con(inFDirRaster,inElevRaster - FocalStatistics(inElevRaster, NbrRectangle(3,3), "MINIMUM"))
outElDifRaster.save("floweldif")
# Iterative calculation of HAND raster
inElDiffRaster = Raster("floweldif")
maxIter = 100
i = 0 # iteration counter limits number of loops for testing
while nullDif != 0:
for idx in lstDirection:
focalMaskFile = "C:/Data/GIS_Data/DEM/FocalStatNeighbor/" + "DirCode" + str(idx) + ".txt"
outHandRaster = Con(inFDirRaster == idx,Con(IsNull(outHandRaster),inElDiffRaster + FocalStatistics(outHandRaster, NbrIrregular(focalMaskFile), "MAXIMUM"),outHandRaster), outHandRaster)
else:
# Calculate values for indeterminate flow direction cells
outHandRaster = Con(IsNull(outHandRaster),inElDiffRaster + FocalStatistics(outHandRaster, NbrRectangle(3,3), "MINIMUM"),outHandRaster)
i += 1
print(str(i) + " iterations complete")
if i % 5 == 0:
outHandRaster.save("hand20")
nullOutRaster = Con(IsNull(outHandRaster),1)
nullOutRaster.save("handNull")
cursor = arcpy.da.SearchCursor("handNull","Count")
newCount = cursor.next()[0]
print "nullcount = ", nullCount, "newCount = ", newCount
nullDif = nullCount - newCount
nullCount = newCount
print "nullDif = ", nullDif
cursor.reset()
if i >= maxIter:
outHandRaster.save("hand20gh2") #restart file
break
I found a way! I used a couple of different sources and manipulated the script to fit my needs (mostly from gotanuki on StackExchange). A variation of this is already on the site, but the of data that I was using is different. In case someone else runs into this issue, this is the solution I used:
# Define orientation (start, mid, end) and draw Matchlines function
# This function is defined before it is called
def orientMatchlines(feat,ix,iy):
# If the line is horizontal or vertical, the slope and
# negative reciprocal calculations will fail, so do this instead
if starty==endy or startx==endx:
if starty == endy:
y1 = iy + matchDist
y2 = iy - matchDist
x1 = ix
x2 = ix
if startx == endx:
y1 = iy
y2 = iy
x1 = ix + matchDist
x2 = ix - matchDist
else:
# Get slope of line
m = float((starty - endy)/(startx - endx))
# Get negative reciprocal
negativereciprocal = -1*((startx - endx)/(starty - endy))
# For all values of slope, calculate perpendicular line
# with length = matchDist
if m > 0:
if m >= 1:
y1 = negativereciprocal*(matchDist)+ iy
y2 = negativereciprocal*(-matchDist) + iy
x1 = ix + matchDist
x2 = ix - matchDist
if m < 1:
y1 = iy + matchDist
y2 = iy - matchDist
x1 = (matchDist/negativereciprocal) + ix
x2 = (-matchDist/negativereciprocal)+ ix
if m < 0:
if m >= -1:
y1 = iy + matchDist
y2 = iy - matchDist
x1 = (matchDist/negativereciprocal) + ix
x2 = (-matchDist/negativereciprocal)+ ix
if m < -1:
y1 = negativereciprocal*(matchDist)+ iy
y2 = negativereciprocal*(-matchDist) + iy
x1 = ix + matchDist
x2 = ix - matchDist
point1.X = x1
point1.Y = y1
point2.X = x2
point2.Y = y2
lineArray.add(point1)
lineArray.add(point2)
del x1
del x2
del y1
del y2
# Create search cursor in inLine
rows = arcpy.SearchCursor(inLine)
# Get number of records in inLine
numRecords = int(arcpy.GetCount_management(inLine).getOutput(0))
# Create new point files and array to collect values
point1 = arcpy.Point()
point2 = arcpy.Point()
lineArray = arcpy.Array()
# Define counter
counter = 0
# Loop over rows in outLine
for row in rows:
# Create the geometry object
feat = row.Shape
# Get coordinate values as lists
firstpoint = feat.firstPoint
lastpoint = feat.lastPoint
midpoint = feat.centroid
# Get X and Y values for each point
startx = firstpoint.X
starty = firstpoint.Y
endx = lastpoint.X
endy = lastpoint.Y
midx = midpoint.X
midy = midpoint.Y
m = ((starty - endy)/(startx - endx+.0001))
# For all points besides the last one
if counter < numRecords - 1:
orientMatchlines(feat,startx,starty)
# For the last point
else:
orientMatchlines(feat,endx,endy)
#Create insert cursor
cur = arcpy.InsertCursor(matchlines)
#Insert new row from array
feat = cur.newRow()
#feat.slope = m
feat.shape = lineArray
cur.insertRow(feat)
lineArray.removeAll()
del cur
# Increase counter by 1 and start again
counter = counter + 1
del row
del rows
Best Answer
It would be easier to convert your line directly to raster, and then use spatial analyst tools to modify the DEM