2

I am working on a task to calculate Height Above Nearest Drainage (HAND) - that is the height above the nearest point along the flow path that is on a drainage - not the geographically nearest drainage point. I have created a raster where the value is the elevation of the next downstream point, based on a flow direction raster. I plan to populate my result raster by (unless I can figure out anything better):

Start with a raster where all stream cells are assigned a height of 0 above drainage and everything else is null. Use focal statistics to find the height above drainage for each cell's downslope cell (initially the values only along the stream) then calculate the HAND for each cell from:

HAND = cell elevation - downslope cell elevation + downslope cell HAND

This will only populate cells where downslope cell HAND is not null (i.e. adjacent to the stream initially). Then I will add those cells to the HAND raster and iterate until I have calculated all cells in the area.

How do I know when I'm done iterating?

My current thought is to either make a copy of the HAND raster before each iteration and compare it to the raster after to see if all cells have the same value or if any cell has a value where the old raster has a null, or to count the number of non-zero cells in each raster and stop when the number is the same. I don't know how to do either or if there is a better way.

  • Perhaps create a binary raster using IsNull and as you modify the cells set from 0 (not null) to 1 (is null) you will be finished when the binary raster is all 1. What programming language are you using for this? – Michael Stimson Nov 13 '14 at 2:01
  • 1
    This may or may not help but it was my first thought upon reading your question title. You may want to throw the two rasters at Raster Compare (Data Management) and see if the result is useful to you. – PolyGeo Nov 13 '14 at 2:04
  • I looked at raster compare but am unclear how I would use it in my iteration. So far I have been using a manual step-by-step process using the toolbox but it looks like I will need to use the python window or an external python script (my python skill are pretty rudimentary). – haresfur Nov 13 '14 at 21:26
  • I tried Michael Miles-Stimson's method of creating a binary raster. I then used a cursor to retrieve the Count field for that raster. The idea is that that count field won't change when the iteration is done. My cursor code works from the command prompt but fails when I include it in the loop - I think because file locks won't let me use the cursor. So, I'm still stuck. Thanks all. – haresfur Nov 23 '14 at 21:20
1

I believe this code solves the problem, although it is a bit ugly. It creates a binary raster with 1 values for each null cell in my HAND raster. I then used a cursor to retrieve the count field for the binary raster's table and compare that to the count from the previous iteration in a while loop. The loop ends when the count stops changing (it may be that there are still null values at the end of the process so I can't just let the count go to 0). The HAND raster is recalculated with focal statistics in an inner For loop. To use the cursor to retrieve the count for the binary raster, it seems you have to save it - you can't just use the variable name. I had quite a bit of difficulty with crashing leaving the rasters and my scratch geodatabase in broken states so there is probably some clean-up that could be done on the code. I left in a lot of print statements and I limit the maximum iterations. For some reason, when restarting the process, I have to save the HAND raster with a new name or it crashes. Thanks for the help, and any further comments are welcome. Hope this helps someone else.

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
env.workspace = "C:/Data/GIS_Data/DEM"
outHandRaster = Raster("HAND_250")
#for restarting iterations
#outHandRaster = Raster("hand20")
inElDiffRaster = Raster("floweldiff")
inFDirRaster = Raster("fdir20m")
lstDirection = [1,2,4,8,16,32,64,128]
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

maxIter = 50
i = 0 # iteration counter limits number of loops for testing

while nullDif != 0:
    outHandRaster.save("hand20baka")
    for idx in lstDirection:
        focalMaskFile = "C:/Data/GIS_Data/DEM/FocalStatNeighbor/" + "DirCode" + str(idx) + ".txt"
        outHandRaster = Con(inFDirRaster == idx,inElDiffRaster + FocalStatistics(outHandRaster, NbrIrregular(focalMaskFile), "MAXIMUM"),outHandRaster)
    outHandRaster.save("hand20") # change to new name if restarting iterations
    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()
    print(str(i+1) + " iterations complete")
    i += 1
    if i >= maxIter:
        break
0

This is not the answer for your question but may be for your problem. I would populate stream raster, for each point (cell) in streams populate watershad. Within populated watershads, iterate trough all points and calculate HAND something like: cell_value - min(watershad).

  • If I understand you correctly, I don't think this will work because each cell will be a member of watersheds for multiple stream cells down the drainage. That is quite a large number of water sheds to calculate and hard to sort out which watershed stream value to use for HAND. Thanks and please elaborate on your answer if I'm wrong. – haresfur Nov 22 '14 at 3:27
  • I do have problems typing my comment (3rd attempt). Yes it will create many miniwatershads, is that a problem? This proces is very the same, as creating watershads for stream gauges, just in much larger scale, as you treat each stream point as gauge. My recommendation is to go with raster. Oryginal way of delinetion watershads for gauges in ArcHydro would produce vector. Raster is way faster. – Tomek Nov 24 '14 at 7:52

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.