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

I'm using ArcGIS 10 with Spatial Analyst.

I have a set of rasters (8 in total) that only every contain a 1 or a 0 for each cell. Each raster represents a different years worth of data. For arguments sake year 1 to year 8.

I can add all the rasters together which will give me a final grid with values ranging from 0 to 8. An 8 indicating that the cell was constantly 1 for the set of rasters (all years).

I would like to find out for each cell the longest consecutive number of 1's.

So for example the total grid may record for a single cell a value of say 5 but over the 8 grids that cell has the largest consecutive number of 1's equal to 3. Or another way of expressing this is for 3 years that cell was a 1 then it started to oscillate between zeros and ones.

My raster processing skills is not as hot as my vector processing skill and I have had a good look at the ESRI help file but I can't figure out how one would achieve this using off the shelf geo-processing tools?

Any ideas?

• That's actually a pretty cool analysis. As usual, there's more than one way to do what you're trying to do. I think you're going to need to do some programming to loop through all the combinations though. Jun 14, 2012 at 18:32
• A general comment (inspired by this remark by @MLowry): please upvote questions whenever they look interesting or are clearly articulated. Good questions drive everything on our site; please let's do what we can to reward those who ask them! Jun 14, 2012 at 18:44

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

• A huge +1: This is a very good idea. (The only limitation is that when more than 31 grids are involved, you will run out of bits to use.) I have taken the liberty to flesh out your idea a little so that others can see how easy it is to implement. Jun 15, 2012 at 11:59

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.

• There might be a typo in the type code block: count=count=1 should probably be count=count+1 Jun 14, 2012 at 19:25
• @ML Thanks (good eyes!): it's fixed now. It's hard to make pseudocode absolutely correct; human review is a real asset in finding errors. Also, although I did not test the solutions in ArcGIS, I did implement the first solution in R, so I have some assurance the approach is correct. Jun 14, 2012 at 20:33
• "whuber" yet again you are THE man who knows! God help the rest of us if you ever get run over by a bus! Your initial Python approach was the direction I was thinking but I know with rasters one can often do everything so much slicker which you have proven. If you find yourself in Britain it would be a honour to buy you a pint of our best room temperature flat beer! :) Jun 15, 2012 at 9:46
• Thanks, Duncan: but check out Andy Harfoot's excellent solution! Jun 15, 2012 at 12:19

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

• This is precisely what @Andy Harfoot suggested a few hours earlier, Ryan. :-) Jun 15, 2012 at 12:18
• Thanks and sorry. I am read this on my phone on vacation. Jun 15, 2012 at 13:16