topic Consecutive values in a set of rasters in Python Questions

Consecutive values in a set of rasters

anTonialcaraz — Sat, 11 Dec 2021 11:46:01 GMT

I'm running ArcGIS 10.2.1 Advanced with SA and Python 2.7

I'm trying to replicate in Python something done in excel:

I have a set of 12 rasters with values 0 and 1, each raster representing a month of the year. I need to create an output raster that will count the maximum number of consecutive months occurring. In the excel example above the result is 4 as it starts in December and accumulates through until March. So it is looping through all months.

I have the column called "Rating Count" above sorted. The problem I encounter is with Jan in the "Consec Mo" column as this value is not fixed because it is fed from the December value in the same column.

I have tried this:

# Set monthly values to 0 and 1


jan = Con(IsNull(Rating_jan), 0, 1)
feb = Con(IsNull(Rating_feb), 0, 1)
mar = Con(IsNull(Rating_mar), 0, 1)
apr = Con(IsNull(Rating_apr), 0, 1)
may = Con(IsNull(Rating_may), 0, 1)
jun = Con(IsNull(Rating_jun), 0, 1)
jul = Con(IsNull(Rating_jul), 0, 1)
aug = Con(IsNull(Rating_aug), 0, 1)
sep = Con(IsNull(Rating_sep), 0, 1)
oct = Con(IsNull(Rating_jan), 0, 1)
nov = Con(IsNull(Rating_jan), 0, 1)
dec = Con(IsNull(Rating_jan), 0, 1)




# Calculate consecutive and MAX value


jan1 = jan
feb2 = Con((feb) == 1, (jan1) + 1, 0)
mar3 = Con((mar) == 1, (feb2) + 1, 0)
apr4 = Con((apr) == 1, (mar3) + 1, 0)
may5 = Con((may) == 1, (apr4) + 1, 0)
jun6 = Con((jun) == 1, (may5) + 1, 0)
jul7 = Con((jul) == 1, (jun6) + 1, 0)
aug8 = Con((aug) == 1, (jul7) + 1, 0)
sep9 = Con((sep) == 1, (aug8) + 1, 0)
oct10 = Con((oct) == 1, (sep9) + 1, 0)
nov11 = Con((nov) == 1, (oct10) + 1, 0)
dec12 = Con((dec) == 1, (nov11) + 1, 0)




Final_Output = CellStatistics([(jan1), (feb2), (mar3), (apr4), (may5), (jun6), (jul7), (aug8), (sep9), (oct10), (nov11), (dec12)]
, "MAXIMUM", "NODATA")


Final_Output.save(OUTWorkspace + "\\" + "Final_Output")

But as I said the "jan" value (row 23) is not fixed as in the script so it is not quite what I'm after.

I wonder if someone could give me some hints about how to solve this issue.

Thanks in advanced

Tony

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

NeilAyres — Thu, 22 Oct 2015 13:25:21 GMT

Sorry bit confused here.

I can see that Dec14, Jan15, Feb15, Mar15 are consecutive.

But how do you know that Dec is the previous year and not 8 month later?

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Thu, 22 Oct 2015 15:45:51 GMT

Thanks Neil. Sorry, I probably didn't explain myself properly.

I'm not dealing with years here; just months of any year. I need to calculate cumulative number of months that hold data (in this case value 1). The calculation should loop through the 12 months with no fixed starting point.

Hope it is a bit more clear now. Thanks again!

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

NeilAyres — Thu, 22 Oct 2015 16:56:25 GMT

That any consecutive in a loop of 12 will be quite tricky.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DarrenWiens2 — Thu, 22 Oct 2015 17:33:42 GMT

If I understand correctly, rather than hard-coding month names into your variable names, I believe you should try naming based on something like "number of months before now (or given date)". So, if you ran the code now (October) the value for September would be held in the variable "months_ago_1", August in "months_ago_2", etc. When you run the code in November, October would move into "months_ago_1" etc.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Fri, 23 Oct 2015 08:06:04 GMT

Thanks Neil/Darren,

Darren, WHEN running the code it is not relevant here. I'll try to graphically explain again:

I have a stack of 12 cells each one representing presence or absence of data (1/0) for 12 months.

I need to produce an output raster cell that tells me the maximum number of cumulative months with data.

In example 1 above the output is 5 because there is data for sep, oct, nov, dec and jan

In example 2 above the output is 3 because there is data for apr, may and jun

Cheers

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Fri, 23 Oct 2015 08:11:11 GMT

It is something that I have in excel but it seems not that straight forward to code in Python (or maybe it is...)

A function in python that calculates the "last occurrence"; in this case the last occurrence of 0 so it can start counting from there, I guess it would do the job. However, I cannot find in the documentation anything related to that.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DanPatterson_Retired — Sat, 11 Dec 2021 11:46:04 GMT

Numpy format would be best

>>> #----------------------------------------
>>> data = [0,0,1,1,1,1,0,0,0,1,1,1,1,1,0,0,1,1,1,1]
>>> a = np.array(data,dtype=bool)  # boolean data
>>> chg = (a[:-1] ^ a[1:])        # change
>>> pos = np.arange(len(chg))[chg] # positions
>>> diff = np.diff(pos)            # differences
>>> max_seq = np.max(diff)        # max in sequences
>>>
>>> chg
array([False,  True, False, False, False,  True, False, False,  True,
      False, False, False, False,  True, False,  True, False, False, False], dtype=bool)
>>> pos
array([ 1,  5,  8, 13, 15])
>>> diff
array([4, 3, 5, 2])
>>> max_seq
5

I am sure it can also be done with cursors and other stuff, if you can't get your data into a numpy array using numpyarraytotable or similar

EDIT

to prevent wraparound you need to pad the array at the beginning and the end. It can be done int this fashion, in this example yielding the same result

>>> # to prevent wrap around sequence starts at zero
>>> an = np.hstack(([0],a,[0]))  # a is from above
>>> diff2 = np.diff(an)
>>> begin = np.where(diff2 > 0)
>>> end = np.where(diff2 < 0)
>>> # tuples are returned to get the 0th result
>>> end[0] - begin[0]
array([4, 5, 4], dtype=int64)
>>>

Arrays are powerful tools

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DanPatterson_Retired — Fri, 23 Oct 2015 19:52:40 GMT

I posted a numpy arcpy example on sequences on my blog

If the data can be brought out into tabular form then this example may help. You can use RasterToNumpyArray to get rasters to arrays or TableToNumPyArray to get tabular data into array format.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DarrenWiens2 — Sat, 11 Dec 2021 11:46:07 GMT

I see. I think you can basically use your existing code, wrapping around the year one more time to pick up the dec->jan connection. Then, use the second "year" of rasters to get the max.

jan1 = jan  
feb2 = Con((feb) == 1, (jan1) + 1, 0)  
mar3 = Con((mar) == 1, (feb2) + 1, 0)  
apr4 = Con((apr) == 1, (mar3) + 1, 0)  
may5 = Con((may) == 1, (apr4) + 1, 0)  
jun6 = Con((jun) == 1, (may5) + 1, 0)  
jul7 = Con((jul) == 1, (jun6) + 1, 0)  
aug8 = Con((aug) == 1, (jul7) + 1, 0)  
sep9 = Con((sep) == 1, (aug8) + 1, 0)  
oct10 = Con((oct) == 1, (sep9) + 1, 0)  
nov11 = Con((nov) == 1, (oct10) + 1, 0)  
dec12 = Con((dec) == 1, (nov11) + 1, 0) 

# wrap back around
jan1_2 = Con((jan) == 1, (dec12) + 1, 0) 
feb2_2 = Con((feb) == 1, (jan1_2) + 1, 0) 
mar3_2 = Con((mar) == 1, (feb2_2) + 1, 0) 
apr4_2 = Con((apr) == 1, (mar3_2) + 1, 0) 
may5_2 = Con((may) == 1, (apr4_2) + 1, 0) 
jun6_2 = Con((jun) == 1, (may5_2) + 1, 0) 
jul7_2 = Con((jul) == 1, (jun6_2) + 1, 0) 
aug8_2 = Con((aug) == 1, (jul7_2) + 1, 0) 
sep9_2 = Con((sep) == 1, (aug8_2) + 1, 0) 
oct10_2 = Con((oct) == 1, (sep9_2) + 1, 0) 
nov11_2 = Con((nov) == 1, (oct10_2) + 1, 0) 
dec12_2 = Con((dec) == 1, (nov11_2) + 1, 0)

Final_Output = CellStatistics([(jan1_2), (feb2_2), (mar3_2), (apr4_2), (may5_2), (jun6_2), (jul7_2), (aug8_2), (sep9_2), (oct10_2), (nov11_2), (dec12_2)]  , "MAXIMUM", "NODATA")

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Mon, 26 Oct 2015 10:25:00 GMT

Thanks Darren,

That works nicely.

The only thing is that the Max value you get is 24 so that needs to be set back to 12 in the final output.

Other than that, that's good.

Thanks again

Tony

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Mon, 26 Oct 2015 10:28:09 GMT

Thanks Darren,

That seems a very elegant way of doing it.

I'm trying to adapt your script to my set of 12 rasters using the RasterToNumpyArray.

Thanks again!

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DanPatterson_Retired — Mon, 26 Oct 2015 11:48:48 GMT

I will post a 3d example sometime this week. The only issue I have seen so far is array size vs memory. You may need to use an generator to pass your arrays in, instead of holding them in all at once.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

anTonialcaraz — Mon, 26 Oct 2015 12:32:08 GMT

(I meant "Thanks Dan" on my previous reply)

I'll keep an eye on your 3rd example.

Cheers

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DanPatterson_Retired — Sat, 11 Dec 2021 11:46:09 GMT

Busy now, will finish later... consider these huge binary arrays (aka 3x3) forming a time stack for a location. Just an idea, I will let you play with the other axes should you be interested in other differences. As indicated, it may be too much to load a large sequence in memory, but the whole process could be parsed. I took out some of the extra print lines between the tiles

>>> a
array([[0, 1, 0],
      [1, 1, 1],
      [0, 1, 0]])
>>> b
array([[1, 0, 1],
      [0, 0, 0],
      [1, 0, 1]])
>>> c
array([[  1,    0,  10],
      [  0,  100,    0],
      [1000,    0,    1]])

>>> stack
array([[[  0,    1,    0],
        [  1,    1,    1],
        [  0,    1,    0]],

      [[  0,    1,    0],
        [  1,    1,    1],
        [  0,    1,    0]],

      [[  0,    1,    0],
        [  1,    1,    1],
        [  0,    1,    0]],

      [[  1,    0,    1],
        [  0,    0,    0],
        [  1,    0,    1]],

      [[  1,    0,    1],
        [  0,    0,    0],
        [  1,    0,    1]],

      [[  1,    0,  10],
        [  0,  100,    0],
        [1000,    0,    1]]])
>>> # stack = np.array((a,a,a,b,b,c))

... np.diff(stack,axis=0)
...
array([[[  0,  0,  0],
        [  0,  0,  0],
        [  0,  0,  0]],

      [[  0,  0,  0],
        [  0,  0,  0],
        [  0,  0,  0]],

      [[  1,  -1,  1],
        [ -1,  -1,  -1],
        [  1,  -1,  1]],

      [[  0,  0,  0],
        [  0,  0,  0],
        [  0,  0,  0]],

      [[  0,  0,  9],
        [  0, 100,  0],
        [999,  0,  0]]])

later

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

curtvprice — Sat, 11 Dec 2021 11:46:12 GMT

Speaking of elegant, Darren, we can do better than that. 0 is false and anything else is true so:

jan1 = jan
feb2 = Con(feb, jan1 + 1, 0)   
mar3 = Con(mar, feb2 + 1, 0)   
apr4 = Con(apr, mar3 + 1, 0)   
may5 = Con(may, apr4 + 1, 0)   
jun6 = Con(jun, may5 + 1, 0)   
jul7 = Con(jul, jun6 + 1, 0)   
aug8 = Con(aug, jul7 + 1, 0)   
sep9 = Con(sep, aug8 + 1, 0)   
oct10 = Con(oct, sep9 + 1, 0)   
nov11 = Con(nov, oct10 + 1, 0)   
dec12 = Con(dec, nov11 + 1, 0) 

# wrap back around 
jan1_2 = Con(jan, dec12 + 1, 0)
feb2_2 = Con(feb, jan1_2 + 1, 0) 
mar3_2 = Con(mar, feb2_2 + 1, 0) 
apr4_2 = Con(apr, mar3_2 + 1, 0) 
may5_2 = Con(may, apr4_2 + 1, 0) 
jun6_2 = Con(jun, may5_2 + 1, 0) 
jul7_2 = Con(jul, jun6_2 + 1, 0) 
aug8_2 = Con(aug, jul7_2 + 1, 0) 
sep9_2 = Con(sep, aug8_2 + 1, 0) 
oct10_2 = Con(oct, sep9_2 + 1, 0) 
nov11_2 = Con(nov, oct10_2 + 1, 0) 
dec12_2 = Con(dec, nov11_2 + 1, 0) 
counts = CellStatistics(
    [jan1_2, feb2_2, mar3_2, apr4_2, may5_2, jun6_2,
    jul7_2, aug8_2, sep9_2, oct10_2, nov11_2, dec12_2],
  "MAXIMUM", "NODATA")   
# "wraps" may give us more than 12 (actually 24 if no gaps)
final_raster = Con(counts > 12, 12, counts)

Note these are all local operators so with deferred raster calculations this may be extremely efficient, and you don't have to load all the rasters into numpy arrays - so this approach is more likely to be successful with large rasters.

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DarrenWiens2 — Sat, 11 Dec 2021 11:46:15 GMT

I think Dan Patterson had it 95% right the first time (the chg and pos arrays are great), and then bypassed the problem and took it to the moon. Really, you just need to append the first array (12 items) to itself (24 items) to capture the wraparound.

>>> import numpy as np
... year1 = np.array([1,0,0,1,1,1,1,0,1,1,1,1],dtype=bool) # year data
... print year1
... doubleYear1 = np.append(year1,year1) # append year data to itself
... print doubleYear1
... chg = (doubleYear1[:-1] ^ doubleYear1[1:]) # change array. Ingenious, really.
... print chg
... pos = np.arange(len(chg))[chg] # indices of change
... print pos
... diff = np.diff(pos) # difference between n+1 and n
... print diff
... max_seq = np.max(diff) # max
... print max_seq
... if max_seq > 12:
...     max_seq = 12

[ True False False  True  True  True  True False  True  True  True  True]

[ True False False  True  True  True  True False  True  True  True  True
  True False False  True  True  True  True False  True  True  True  True]
[ True False  True False False False  True  True False False False False
  True False  True False False False  True  True False False False]
[ 0  2  6  7 12 14 18 19]

[2 4 1 5 2 4 1]

5

Re: CONSECUTIVE VALUES IN A SET OF RASTERS

DanPatterson_Retired — Sat, 11 Dec 2021 11:46:18 GMT

Darren Wiens yes, definitely, but it sometimes depends on perspective...lets try rolling/padding by 1 position/year.

In this example the first year and the last year are the same.

>>> a = np.array([1,0,0,1,1,1,1,0,1,1,1,1],dtype=bool) # with first and last month the same
>>> b = np.roll(a,-1)
>>> chg = (b ^ a)
>>> pos = np.arange(len(chg))[chg]
>>> diff = np.diff(pos)
>>> max_seq = np.max(diff)
>>> max_seq
4
>>> a.astype(int)
array([1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1])
>>> b.astype(int)
array([0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1])
>>> chg.astype(int)
array([1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0])
>>> pos
array([0, 2, 6, 7])
>>> diff
array([2, 4, 1])
>>>

So the question becomes how far do you roll (12 years as in your case)?

But what if the, first year and the last year are different?

>>> a = np.array([1,0,0,1,1,1,1,0,1,1,1,0],dtype=bool) # with last month different than first
>>> b = np.roll(a,-1)
>>> chg = (b ^ a)
>>> pos = np.arange(len(chg))[chg]
>>> diff = np.diff(pos)
>>> max_seq = np.max(diff)
>>> # results........
>>> a.astype(int)
array([1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0])
>>> b.astype(int)
array([0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1])
>>> chg.astype(int)
array([1, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1])
>>> pos
array([ 0,  2,  6,  7, 10, 11])
>>> diff
array([2, 4, 1, 3, 1])
>>> max_seq
4

The same results are yielded, largely by coincidence since the first and last year are different

It is interesting that you can roll forward or backward as long as you keep track of your positions.

You can also use padding with different pad widths and types AND direction. In this example the edges were padded with their neighbor.

>>> a = np.array([1,0,0,1,1,1,1,0,1,1,1,1],dtype=bool) # with first and last month the same
>>> b = np.pad(a,1,mode="edge")
>>> b.astype(int)
array([1, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1])
>>> chg = (b[:-1] ^ b[1:]) 
>>> pos = np.arange(len(chg))[chg]
>>> diff = np.diff(pos)
>>> max_seq = np.max(diff)
>>> max_seq
4

It just becomes more interesting when you are trying to shift, difference and pad sequences of 2D+ arrays with or without padding or replication. As I indicated, I will try to summarize in an upcoming blog because raster size can play a role. Also, if you are only interested in limited positions within the raster, reducing the input raster to a sparse array will reduce memory needs and speed up processing time.