Are two rasters the same, including null values, in ArcGIS desktop?

Question

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.

Answer 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

Answer 2

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).