Counting consecutive pixel values for set of rasters using ArcGIS Spatial Analyst?

Just chatting about this and wondering whether you could approach the problem by treating the input grids as a binary stream. This would allow you to combine them to give a unique summary integer for the sequence - ie 01110101 = 117. This value could then be reclassified to give the maximum number of consecutive 1s.

Here's an example showing one way to combine eight grids:

2*(2*(2*(2*(2*(2*(2*"g8" + "g7") + "g6") + "g5") + "g4") + "g3") + "g2") + "g1"

Bitwise operations could also be pressed into service for this step. Alternatively, you can use combine followed by a field calculation. (The field calculation will have an expression similar to the preceding one.)

The reclassification table has to provide max run lengths for all values between 00000000B = 0 and 11111111B = 255. In order, here they are:

0, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 6, 7, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 1, 1, 1, 2, 1, 1, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 1, 1, 1, 2, 1, 1, 2, 3, 1, 1, 1, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 5, 6, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 3, 4, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 3, 3, 4, 5, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 7, 8

This approach is limited to about 20 grids in ArcGIS: using more than this can create an unwieldy attribute table. (Combine is specifically limited to 20 grids.)

Because this is a local operation, let's work out how to do it for a single cell: Map Algebra will take care of the rest.

First note that the order of the rasters obviously matters. Therefore single-shot cell statistics, such as a cell sum, won't do it.

If we were to encounter a sequence such as 01110101 at a given cell, we would process this from start to end and

1. Begin with a count at zero.

2. Increment the count each time we encounter a 1.

3. Reset the count each time we encounter a 0, after saving the last count.

4. At the end, take the maximum saved count (including the final count).

Step 1 is implemented with a constant zero grid. Steps 2 and 3 depend on what we encounter: it is therefore a conditional operation. Step 4 clearly is a local maximum. We would code this, then, a little more formally as:

count = 0
result = 0
For each value:
    If (value==1):
        count=count+1
    else
        result = max(result, count)
        count=0
result = max(result, count)

That's best done with a Python script when you have many grids, but with eight it's not onerous to unroll the loop and write out the steps by hand. This reveals a slight problem: the result=max(longest,count) is sort of a "side-effect" that's hard to code with raster operations. (But it can be done, as shown in the second solution below.) It's also inefficient, because it adds an extra calculation at each step. We therefore modify the approach a little, with an aim to putting off the max operation until the end. This will require saving a separate count at each stage.

In going through this process I also found a shortcut for the first step. This leads to the following solution, which although a little long and RAM-intensive, is simple and involves quickly-executed steps:

result1 = "grid1"
result2 = con("grid2"==1, "result1"+1, 0)
result3 = con("grid3"==1, "result2"+1, 0)
result4 = con("grid4"==1, "result3"+1, 0)
result5 = con("grid5"==1, "result4"+1, 0)
result6 = con("grid6"==1, "result5"+1, 0)
result7 = con("grid7"==1, "result6"+1, 0)
result8 = con("grid8"==1, "result7"+1, 0)
CellStatistics(["result1", "result2", "result3", "result4", "result5", "result6", "result7" "result8"], "max")

The actual syntax varies with your version of ArcMap. (For instance, CellStatistics is new to version 10, I believe, but a local maximum operation has always been available.)

In the example with input 01110101, the sequence of "result*" grids would contain the values 0, 1, 2, 3, 0, 1, 0, 1, so at the end CellStatistics would return 3, the length of the longest string of 1's.

If RAM is scarce, the solution can be modified to re-use the intermediate results, at a cost of approximately doubling the execution time:

result = "grid1"
temp = con("grid2"==1, "result"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
temp = con("grid3"==1, "temp"+1, 0)
result = CellStatistics[["temp", "result"], "max"]
...
temp = con("grid8"==1, "temp"+1, 0)
CellStatistics[["temp", "result"], "max"]

In the example with input 01110101, the ("temp", "result") values would be (NoData, 0) after the first line and after each pair of ("Con", "CellStatistics") operations the values would be (1, 1), (2, 2), (3, 3), (0, 3), (1, 3), (0, 3), (1, 3). Once again the final value is 3.

The regular pattern of Map Algebra expressions in either solution indicates how to code the algorithm as a loop in a script, changing indexes as appropriate with each iteration.

Have you thought of changing the values from 0 and 1 to values with the power of 2 (1,2,4,8,16,32). When you combine the 8 grids you will get unique values for each cell which will give you consecutive info (ie: a value of 3 means year 1 and 2, where a value of 54 would be years 6 to 8).

Just a thought