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

• Hornbydd: Sounds like you are asking to start from a cell and travel in a downstream direction to the river which you have masked out with nodata? Get the mean of you coefficients and assign that to the cell then move onto the next cell and repeat, have I understood this correctly? If so sounds like you want to create a d/s flow path that could be used to select the cells, have a look at this Q&A.
• Maddalena: The problem is that the output of the Downslope flowpath tool is a raster where every cell of a flow path has the same value (first source -> first flow path -> cells value=1, 2nd source -> 2nd flow path -> cells value=2). That works pefectly till you have separate flow paths, but in my case, as I have as many sources as the number of cells in the raster, I would have too many paths and an useless output. Or, I should be able to do the process for every cell one-by-one (start from the first, calculate flowpath, get the mean of the coefficient, than go to the next cell and repeat). But how?

Answer 1

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

enter image description here

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

enter image description here

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?

• Hornbydd: FelixIP, been looking at your code above, very interesting, I think this was what was going on in my head but I have a simple question about 1 line. I'm not mega hot on python so when I saw this line (tot,m,nR,nC=S,1,nRow,nCol) I am not understanding it. Is this some clever way of initializing variables in python?
• FelixIP: @Hornbydd, yes it is. It is equivalent of multiple assignment statements, i.e. tot = S; m = 1 etc. Apparently if you need to swap values of 2 variables this (x, y)=(y, x) works! In general being able to trace cell by cell is mighty thing. I've used it before through excel vb and binary versions of fdir, it is very fast. The key line here is flipping values in dRow. This is due to difference between numpy array (rows changing from top down) and grid. Modify it slightly to get many useful things, like elev. Above discharge point. Tested on 2000×2000. Fine
• Maddalena: FelixIP, thanks a lot! It seems it works perfectly. It's a run-off coefficient, to calculate a final River Pollution index. If you search for "Estimating river pollution from diffuse sources in the Viterbo province using the potential non-point pollution index" you'll find the descriptive article.
• FelixIP: @Maddalena, glad it works. Do you mind to mark it answered. They say it helps others with searches...

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)

• Maddalena: Hi Pau. If I understand your answer correctly, I don't have a raster for every coefficient, I have one single coefficient raster with a value for every cell (except the ones overlapping the stream network). So I first have to locate the right cells (that constitute the flow path from a source cell to the river network) than get the mean of their coefficient and assign this value to the source cell.