I have searched for some ideas on how to formulate this analysis, but have come up short. I imagine there is a pretty straight forward answer; any help would be greatly appreciates.

Data: DEM, 1/9th arc second (~10-meter), 239,891,470 cells, projected; Streams, National Hydrology Data set, high resolutions, 9,470 features vector; Streams, Same as above, but converted to a raster

Problem: Need to create a raster of the same resolution and extent as the DEM where each cells's value is its height (in meters) above the nearest stream (line or cell).

I am working with ArcGIS 10.1 with spatial analyst. I would prefer to run this analysis in Python, but can do it in R if it is easier.

Edit to clarify: each cell in my stream raster has the value of the elevation at that point. All cells that are not within a stream have a value of zero. For every cell in the extent of the DEM, I am looking to find the nearest stream cell and subtract the value of that cell (the elevation) from the elevation of the cell being analyzed. Calculate the difference in elevation between each cell and the nearest cell representing a stream.

Thank you for your consideration!

  • Are you asking to extract from the DEM just the cells that are the cells in your rasterized version of your rivers? If so the CON tool in Spatial Analyst toolbox can do that. – Hornbydd Nov 12 '13 at 22:02
  • Thanks Hornbydd. Not exactly, the resulting raster will cover the entire extent and each cell will contain a value representing its height above the nearest stream. The cells that overlay the rasterized streams will have a value of zero since they coincide with the nearest stream cell and have the same elevation. – Mr.ecos Nov 12 '13 at 22:06
  • Just throwing in a pathological case...If your cell is at the watershed divide and equidistant from 2 streams and one stream is much lower than the other (because one side is steeper) which would you go for? – Hornbydd Nov 12 '13 at 22:19
  • great question. In the case of a tie, I would go with the minimum elevation difference. – Mr.ecos Nov 12 '13 at 22:22
  • If "nearest" means in terms of distance on the map, then just subtract the Euclidean allocation of the stream elevations from the DEM. The procedure is illustrated in an answer to a related question. – whuber Nov 12 '13 at 23:09

whuber has it right.

First, get the euclidean allocation. Second, use the minus operation to get the difference between the two.

Can be easily done in python.

  • Not really: it sounds like OP is looking to calculate the HAND (height above nearest drainage). Therefore the Euclidian distance allocation is inappropriate. – Richard Law Jun 23 at 23:02

I am late here, but this question was my entry into the forum as I was researching the same question.

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 *
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)
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"))
    # Calculate values for indeterminate flow direction cells
    outElDifRaster = Con(inFDirRaster,inElevRaster - FocalStatistics(inElevRaster, NbrRectangle(3,3), "MINIMUM"))

# 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)
        # 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:
        nullOutRaster = Con(IsNull(outHandRaster),1)
        cursor = arcpy.da.SearchCursor("handNull","Count")
        newCount = cursor.next()[0]
        print "nullcount = ", nullCount, "newCount = ", newCount
        nullDif = nullCount - newCount
        nullCount = newCount
        print "nullDif = ", nullDif
    if i >= maxIter:
        outHandRaster.save("hand20gh2") #restart file
  • Thank you for the late but thorough response. Yes, this approach was what I was thinking, but I did make it work with Euclidean dist. When I do it again, I will work through this approach. Thanks again! – Mr.ecos Nov 30 '14 at 23:40

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.