Correcting Errors in Height Above Nearest Drainage (HAND) model

Question

I have used this tread Spatial Analysis - Calculate a cells elevation (height) above the nearest stream cell to create a HAND map in ArcGIS. First I used the hydrology toolset to fill in sink in the SRTM DEM, the generated the flow direction, flow accumulation,and the drainage networks. I then computed the Zonal Statistics using the drainge network and the filled DEM as the data input. Then I used the Euclidean Allocation tool to calculates, for each cell, the zone of the closest source location in Euclidean distance. I subtracted the resulted map from the filled DEM and I obtained the the HAND (Height Above Nearest Drainage)map of the area.

the map is shown below. However, negative values appears in the legend which should no be the case. Can anyone suggest how I can rectify this?

Answer 1

Its possible the negatives are the result of real depressions, so it wouldn't be wrong to just reset any negative values to 0, IMO. Map up the location of the negatives and compare with your input DEM and a decent satellite image before you make a decision.

---

That aside, the key thing with HAND is that it traces a flowpath downhill to the nearest stream cell and then calculates the elevation difference. Avoiding a Euclidean calculation is the whole point.

It took me a while to track down the actual software that Renno et al used for HAND since they didn't bother quoting it in their paper ( >:-( ). Its a bit clunky, but it does the job. Start with the TerraHidro 0.5 console app from http://www.dpi.inpe.br/terrahidro/doku.php and follow the instructions to get d8 grid, contributing area, and flow accumulation rasters using a hydrologically corrected DEM, e.g.

th d8 input_dem.tif output_d8.tif
th d8ca output_d8.tif contareagrid.tif
th d8drainage contareagrid.tif drainagegrid_100.tif 100

then get Terraview 4.2.2 from http://www.dpi.inpe.br/terraview/php/dow.php?body=DowFiles and TerraHidro 0.4.2 from http://wiki.dpi.inpe.br/doku.php?id=download and install them sequentially. Now:

• Open TerraHidro (this opens TerraView with the plugin enabled)
• Go to File > Open Database and follow the prompts to create a new Access database (I know, I know, but its simple...)
• Go to File > Import Raster and follow the prompts to import the SRTM DEM-H, the D8 flow grid and the drainage grid
• Click on the Hydrological Tools button, which is in the upper right corner and has a tiny M on it
• Tab rightwards across the top to find HAND. Follow the prompts to get an output layer
• Outputs can be exported by right-clicking on the theme.

I'm recommending TerraHidro 0.5 for generating the inputs because it can handle larger datasets and is miles faster. It just doesn't have a HAND command.

You can get a very similar output to 'official' HAND by chaining r.stream.extract and r.stream.distance in GRASS 7. The bonus is that you can use dInf, not just d8 if you want, but it does produce some negatives, just as your procedure has.

Answer 2

Solution you refer to might work on very dense stream network. Obviously it is not your case.

I assume that what you call HAND is height above elevation of nearest drainage points downstream, i.e. first point downhill.

Use raster calculator to define your ‘streams’ from flow accumulation raster and assign elevation to it from Filled raster, e.g.:

Con("facc"=>350000,"filled")

Call output raster “flood”. Replace 3 lines in below script pointing to flow direction grid, “flood” and output grid name (at the bottom of the script)

Script assigns elevation of drainage point to cells within point catchment, and saves result to output raster. In order to obtain HAND you’ll need to subtract output raster from filled DEM.

import arcpy, os, traceback, sys, numpy
from arcpy import env
from arcpy.sa import *

fdir=r'D:\Scratch\fdir'
demG=r'D:\Scratch\flood'

d=arcpy.Describe(demG)
##SR=d.spatialReference
try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    dirArray = arcpy.RasterToNumPyArray(fdir,"","","",-9999)
    dem = arcpy.RasterToNumPyArray(demG,"","","",-9999)
    nRows,nCols=dirArray.shape
    blankArray=arcpy.RasterToNumPyArray(demG,"","","",-9999)
    cellsTotal=nCols*nRows
    d=arcpy.Describe(fdir)
    origin=d.extent.lowerLeft
    cSize=arcpy.Raster(fdir).meanCellHeight
    ##  directions to find cell neighbour
    fDirs=(1,2,4,8,16,32,64,128)
    dCol=(1,  1,  0, -1, -1,-1, 0,1)
    dRow=(0, -1, -1, -1,  0,  1, 1,1)
    ##  flipped 
    dRow=(0,  1,  1,  1,  0, -1, -1,-1)
    ##  main loop
    arcpy.SetProgressor("step", "", 0, nRows)
    for nRow in range (nRows):
            for nCol in range (nCols):
                S=dirArray[nRow,nCol]
                if S in (-1,-9999): continue
                nR,nC=nRow,nCol
                cells=[(nR,nC)]
                z=dem[nR,nC]
                while True:
                    direction=dirArray[nR,nC]
                    i=fDirs.index(direction)
                    dX=dCol[i];nC+=dX
                    ## if nC in (0,nCols): break
                    if nC not in range(nCols): break
                    dY=dRow[i];nR+=dY
                    ## if nR in (0,nRows): break
                    if nR not in range(nRows): break
                    cells.append((nR,nC))
                    S=dirArray[nR,nC]
                    z=dem[nR,nC]
                    if S<-9998 or z>-9998:
                        break
                    if S==-1:
                        z=blankArray[nR,nC]
                        break
                for nR,nC in cells:
                    blankArray[nR,nC]=z
                    dirArray[nR,nC]=-1
            arcpy.SetProgressorPosition()
    myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
    oneGrid=Con(myRaster<>-9999,myRaster)
    oneGrid.save(r'D:\Aerials\nztm\toriver')
    del dirArray,blankArray
except:
    message = "\n*** PYTHON ERRORS *** "; showPyMessage()
    message = "Python Traceback Info: " + traceback.format_tb(sys.exc_info()[2])[0]; showPyMessage()
    message = "Python Error Info: " +  str(sys.exc_type)+ ": " + str(sys.exc_value) + "\n"; showPyMessage()   

Note: It takes about 5 minutes to process 3000*3000 cells inputs.