i tried to follow this procedure Spatial Analysis – Calculate a cells elevation (height) above the nearest stream cell but i have some problem to understand everything. I mean the 2th answer in the link.
I don't understand the 3th and 4th points. I have flowdirection, stream network and dem. I would to understand what about flow elevation difference and what i need to do to get the hand model.
[GIS] Explanation of HAND (height above nearest drainage) model procedure
arcgis-10.2arcpyfloodhandhydrology
Related Solutions
Have you consider to use GRASS GIS analysis? I have expirience that GRASS algorithms have very good accurance on hydrology analysis. For example, I want to generate something like drainage network on DTM with resolution 5x5m. I had compared tools from ArcMap (including ArcHydro Tools) and you can view the result on first picture (red lines). Then I tried to use GRASS GIS function 'r.stream.extract' and I had got result shown on picture 2 (red lines). Both drainage lines are generated with cathement area 3 hectares.
It is really different, and it has pretty accurance in comparsion with real streams (picture 3, real streams are blue). And GRASS GIS has many hydrological tools, i.e. for generate catchement area too.
![Drainage lines using ArcMap]](https://i.stack.imgur.com/65J2t.jpg)
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
Best Answer
I will assume you have your DEM, have been able to calculate the flow direction raster (the ArcGis help has a good explanation). You also need to have defined a raster of the drainage which has values of 0 for each cell along a drainage (height above drainage is 0 by definition) and has nulls for all other cells.
You need to know which cell is downslope and the elevation drop for each cell in your DEM, except for the cells that are along the drainage. To do this you use focal statistics with a mask for each possible flow direction.
The directions when you calculate a flow direction raster are defined by powers of 2 so the designation for the directions from cell X are:
32..64..128
16...X....1
8....4....2
So we calculate the elevation of the cell in each direction by using a masks that defines a 3x3 square neighborhood and says only use one cell in that neighborhood for the focal statistics. You need 8 mask files for this. For the cell to the east with value 1 the mask file looks like:
3 3
0 0 0
0 0 1
0 0 0
For the cell to the sw with value 8 it is:
3 3
0 0 0
0 0 0
1 0 0
I call these files DirCode1.txt, DirCode2.txt, Dircode4.txt...DirCode128.txt. Put them all in the same directory.
The function
returns the minimum elevation in the cells around each location where the mask has the value 1. Since there is only one such cell in the mask it doesn't really matter if you use minimum, maximum or whatever and you just subtract that value from the DEM value for the original cell to get the difference. But you only do that for the cell in the direction of flow. Do this by iterating through each flow direction mask file and use the Con function to return a value if the mask file direction is the same as the flow direction for each cell:
# This block only needs to be calculated once for the area
For some reason I forget, you can get cells where the flow direction value is not in the 8 we are using. Define the elevation difference for those bu using the minimum value in the 3x3 region around the cell.
Calculating the HAND is similar. Start by copying your drainage raster to the output raster. For each focal mask file, check each null cell in the output raster to see if that mask matches the flow direction, then add the value from the cell in that direction. At first that will only give a value for cells next to the drainage. When that's done, you have HAND values for more cells. Repeat the process to calculate values for null cells upslope from those cells. Loop, working your way upslope, until you have no more null cells.