Getting cell value along Flow Direction using ArcPy?

Question

I have:

• a flow direction raster where I took off (assigned NoData value) the cells that overlayed the river network (from flow accumulation).
• a perfectly overlapping raster where I asigned a coefficient value to the cells (except the ones overlapping the river network)

I need to calculate an index for every cell of my raster that is calculated as the mean between the coefficient value of the cell and those of all the cells between the first and the river network, along the flow direction.

To explain better: Let's begin with a cell, who has its coefficient value (c1) and has a 8 in the flow direction raster. I want to pick the coefficient value and "move" to the south-west cell (that is the next one along the flow direction). Then I want to pick its coefficient value (c2), see in the FlowDirection raster which is the next cell and move there, to pick its coefficient value (c3)....and so on till I find a cell with NoValue (river network).

Then I need to calculate the mean beetween these values and assign the result, in a new raster, to the beginning cell

(c1+c2+....cn)/n

Answer 1

I've used filled elevation model (with streams removed) to replicate your coefficients:

Keeping in mind that water usually runs downhill :), it is reasonable to expect that all the values in output should be less or equal to the value at start point. I've checked results using difference [Filled]-[ScriptOutput] and found that script (see below) failed (produced slightly negative values) for 2 out of 63414 cells. Both points are sitting on very flat sites next to 'streams'. This is a nature of flow direction algorithm though, called most likely direction guess...

Script:

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

fdir=r'D:\Scratch\fdir1'
coeffs=r'D:\Scratch\fill1'

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    dirArray = arcpy.RasterToNumPyArray(fdir,"","","",-9999)
    coefArray = arcpy.RasterToNumPyArray(coeffs,"","","",-9999)
    blankArray=coefArray
    nRows,nCols=dirArray.shape
    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, cellsTotal)
    for nRow in range (nRows):
            for nCol in range (nCols):
                arcpy.SetProgressorPosition()
                S=coefArray[nRow,nCol]
                if S==-9999:continue
                tot,m,nR,nC=S,1,nRow,nCol
                while True:
                    direction=dirArray[nR,nC]
                    i=fDirs.index(direction)
                    dX=dCol[i];nC+=dX
                    if nC<0 or nC==nCols: break
                    dY=dRow[i];nR+=dY
                    if nR<0 or nR==nRows: break
                    S=coefArray[nR,nC]
                    if S==-9999:break
                    tot+=S;m+=1
                tot=tot/m
                blankArray[nRow,nCol]=tot
    myRaster = arcpy.NumPyArrayToRaster(blankArray,origin,cSize,cSize)
    oneGrid=Con(myRaster<>-9999,myRaster)
    oneGrid.save(r'D:\Rubbish\avers')
    del dirArray,coefArray,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()        

You'll have to replace 3 lines of code, pointing to input flow direction, coefficient and output raster names.

Also note that it is slow, it took 2 minutes to process 234*271 raster on my rather good machine. Script designed to run from Arcmap (User toolbox-Add-Script). If you want it to work stand alone, replace arcpy.AddMessage with print statement. Also progressor statements might not work. Just remove them. Script can be easily modified for all sort of studies, e.g. avalanche run out distance...

Any questions, let me know. Thanks for good question. What kind of coefficients, please? Roughness?

Answer 2

Try this:

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")


# Local variables:
c1 = "F:\\img1.tif"
c2 = "F:\\img2.tif"
c3 = "F:\\img3.tif"
output = "F:\\output.tif"

# Process: Raster Calculator
arcpy.gp.RasterCalculator_sa("(\"%c1%\"+\"%c2%\"+\"%c3%\")/3", output)