I am working with the Flow Accumulation tool in ArcMap 10's spatial analyst hydrology tools. This tool calculates the total number of raster cells that flow into a given cell. However, for my purposes, I would like to restrict how far uphill the accumulation originates.
IE: I would like to see how much flow there is into a cell within 1 km instead of all the way to the ridge-line.
Answer
You can approximate this using standard Spatial Analyst tools. (To do better, you will need to write the code from scratch in Python.)
(The approximation is that the flow accumulation can be limited to a fixed number of steps, where each step is a horizontal, vertical, or diagonal flow across a single cell. For a cellsize of c and n steps, then, you are guaranteed that the total flow distance is at least n * c and at most sqrt(2) * n * c, or about 41% greater. Because the "D8" flow direction methods used in Spatial Analyst are already crude approximations to flow directions and distances, the error of this approximation should be tolerable in many applications.)
The idea is to accumulate the flow one step at a time, then iterate. For instance, on a 30 m cellsize grid, 1 kilometer = 1000 meters = 33 steps, so you would iterate the accumulation 33 times.
I will show how to accumulate flow this way by starting with the illustration from the Spatial Analyst help page on "how flow accumulation works."
The right map is the alleged result of a flow accumulation calculation (assuming a unit amount of water falls into each cell). The left map indicates flow directions both with arrows and colors. It means, for instance, that water falling in the upper left corner will flow into the cell at its lower right. If you were to peer at the actual grid data underlying this (using the wonderful new "Pixel Inspector" for instance) they would look like
For each cell into which water flows, we need a way to find the cells which contribute water, read off the amount of water in them, and add these contributions: that's one step of flow accumulation. The built-in methods to read values within local neighborhoods are called "Focal" statistics. The trick is to use eight separate neighborhoods, each indicating exactly one neighboring cell. Let's look at a picture of how this works:
The first column, for instance, shows direction code "1", corresponding to a nominal easterly flow. At any given cell, the cell flowing eastwards into it must lie due west. This cell is picked out by the three-by-three neighborhood matrix shown in the second row. Thus, when we use this matrix to compute a focal sum (or focal mean), what we're really doing is just shifting all values in the grid one cell to the right. The bottom images highlight which cells have contributions from the given directions. For instance, at the very right the two cells flowing in direction 128 (northeasterly) will affect the two cells that have been shaded orange.
These focal operations can do two things for us: not only can they find how much is contributed by each neighboring cell, they can also identify which neighboring cells really are contributors. This is done by checking that the shifted direction codes match the direction of the shift.
For instance, we identify the shaded light green cells in the first column by comparing the focal sum of the direction grid to the direction code of 1; we identify the shaded dark green cells in the second column by comparing the focal sum of the direction grid to the direction code of 2; and so on.
Finally, we use this information to sum the shifted flows. Because the result of a logical comparison is 1 for true and 0 for false, all we have to do is multiply the results of our shifted-direction comparisons by the shifted flows and add them all up.
How much work is this? Well, we only need to shift the direction grid once: that's eight focal sums. At each incremental step, we need to shift all flow values: that's eight focal sums per step plus a sum of all those focal sums (to obtain the total contribution at each cell). Finally, we need to sum all these incremental grids. It's a lot of work--about 16 grid operations per step--but even for huge grids each step is very fast and the processes are repetitive, making them simple to code.
Let's watch the process in action using the ESRI example. As a check of the work, I will expand the grid's extent to let any flows spilling off the edge accumulate just outside. When we do this, no water is ever lost. A simple check, then, is to sum the total amount of water after each step and verify it has not changed.
At the outset, let's drop one unit of water on each cell. On this six by six grid, that's 36 total units of water. Here's how those units flow at each step using the process described above:
The locations of the 36 water units after the first step are shown in the first column, and then they are shown after the second, third, fourth, and fifth steps. The "incremental flow" row shows the grids computed by obtaining the focal sums and adding them, while the "total flow" row shows the cumulative sums of those incremental flows. For instance, the total flow in the second column (total flow after two steps) is obtained by adding the grid to its left (total flow after one step) and the grid above it (incremental flow between steps 1 and 2).
It's easy to check that the total amount of water in each incremental flow grid stays constant at 36: at least that much is correct!
The borders of these grids are all "sinks": they simply collect any flows that fall into them. Only two cells actually receive any flow: they appear along the bottom. At the end, they have received 26 and 42 units of water. (These values ordinarily would not be shown in a flow accumulation calculation.)
There is no change within the original grid extent after step 5: that's as far as we need to go. After all, six steps in any direction will always take us off the original six by six grid. So the (unrestricted) flow accumulation result appears in the lower right corner of the total flow diagrams. You might notice that it disagrees with the ESRI example. (For instance, the example has a flow accumulation of 35 in a cell where I compute only 27 and, next to it, the example has 2 where I have only 1. It sure is hard to see how that cell could ever receive more than one unit of water!) I think the ESRI example is wrong, but it's instructive to go through the steps yourself to see whether or not that is the case.
To obtain the restricted flow, you don't take this calculation to its end: instead, you stop after the desired number of steps. In the example, if we wanted only the flow accumulation over two cellsizes, we would stop after step 2 and use the result in column "2" of the "Total flow" row.
No comments:
Post a Comment