Combining information from multiple rasters into aggregate raster and populating fields accordingly

Question

I am stuck on how to best combine the information from multiple rasters into one raster. Specifically, I am working with a time series of raster files, whereby each raster corresponds to a different year (they are all the same extent). The pixels of each raster range from values 1-4 (excluding nodata). I am trying to use ArcGIS to combine the information from each year into one aggregate raster whereby each cell will have the following information: 1) the total count of each value (e.g. Total occurrence of "one", total occurrence of "two", total occurrence of "three", etc.); and 2) the last instance of each value (e.g. Most recent year of instance of "one", most recent year of instance of "two", etc.). I have been toying around with field calculator and raster calculator, but am not quite sure where to start on this problem (I'm a newbie in ArcGIS and Python).

Answer 1

Nice way to do so is using script, where you can convert each raster to numpy array and have access to individual cells, e.g.. It is much easier than it looks like. Note it work only on len(rasters) % 20 =0

No programming, brutal force approach:

1. Convert one of the rasters into points
2. Use Extract multivalues to point
3. Analyse resulting massive table in Excel, e.g. COUNTIF(range,=1)
4. Bring totals back to GIS by using pointID field

UPDATED ANSWER

Script below assumes that all rasters stored in 1 folder and they are identical in terms of extents and cell size:

import arcpy, os, traceback, sys, numpy
from arcpy import env
env.overwriteOutput = True
from arcpy.sa import *

inputFolder="D:/AERIALS/Clipped/"
outFolder="D:/AERIALS/Mosaics/"

env.workspace = inputFolder
inRasters = arcpy.ListRasters()
layoutR=arcpy.Raster(inRasters[0])
d=arcpy.Describe(layoutR)
origin=d.extent.lowerLeft
cSize=layoutR.meanCellHeight
zeroArray = arcpy.RasterToNumPyArray(layoutR,"","","",-9999)
nRows,nCols=zeroArray.shape

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    for testValue in range(1,5):
        arcpy.AddMessage("Testing %i" %testValue)
        totalArray=numpy.zeros((nRows, nCols),dtype=numpy.int)
        lateArray=numpy.zeros((nRows, nCols),dtype=numpy.int)
        arcpy.SetProgressor("step", "", 0, len(inRasters),1)
        for j,R in enumerate(inRasters):
            raster2test=arcpy.Raster(R)
            array2test = arcpy.RasterToNumPyArray(raster2test,"","","",-9999)
            for nRow in range (nRows):
                    for nCol in range (nCols):
                        v=array2test[nRow,nCol]
                        if v!=testValue:continue
                        totalArray[nRow,nCol]+=1
                        lateArray[nRow,nCol]=j+1
            arcpy.SetProgressorPosition()
        myRaster = arcpy.NumPyArrayToRaster(totalArray,origin,cSize,cSize)
        myRaster.save('%sTOTAL_%i' %(outFolder,testValue))

        myRaster = arcpy.NumPyArrayToRaster(lateArray,origin,cSize,cSize)
        myRaster.save('%sLATEST_%i' %(outFolder,testValue))
    del totalArray,lateArray
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()        

In order to test it I generated 100 random rasters using:

for i in range(100):
    outRandomRaster = CreateRandomRaster(i)
    outInt = Int(outRandomRaster*5)
    outF="D:/AERIALS/Clipped/Rand_%s" %str(i).zfill(2)
    outInt.save(outF)

Due to nature of generated rasters it is reasonable to expect:

1. for each of tested values (1,2,3,4) mean value of the rasters storing totals will be closed to 20, i.e. 1/5th of 100.
2. Distribution of values inside individual raster close to normal with mode 20

RESULTS:

Means:

Histograms:

UPDATE #2

This second approach is using Spatial Analyst Combine tool. Before testing it use GUI to combine 1st 20 rasters. Inside script Combine is used 4*4 times, thus if test taking too long forget it.

import arcpy, os, traceback, sys, numpy
from arcpy import env
env.overwriteOutput = True
from arcpy.sa import *

inputFolder="D:/AERIALS/Clipped/"
outFolder="D:/AERIALS/Mosaics/"
lookupField="LATEST"

env.workspace = inputFolder
inRasters = arcpy.ListRasters()

try:
    def showPyMessage():
        arcpy.AddMessage(str(time.ctime()) + " - " + message)

    for testValue in range(1,5):
        arcpy.AddMessage("Testing value %i group" %testValue)
        resultS=[]
        for j in range(4):
            arcpy.AddMessage("Combining %i group" %(j+1))
            list2combine=inRasters[j*20:j*20+20]
            outCombine =Combine(list2combine)
            outCombine.save('%sCOMBO_%i' %(outFolder,j))
            arcpy.MakeTableView_management('%sCOMBO_%i' %(outFolder,j), "COMBO")
            arcpy.AddField_management("COMBO", lookupField, "LONG")
            fNames=[f.name for f in arcpy.ListFields(outCombine)]
            reqFields = filter(lambda x: x not in ['Rowid','VALUE','COUNT'], fNames)
            with arcpy.da.UpdateCursor("COMBO",reqFields) as cursor:
                for row in cursor:
                    revList=row[:-1]
                    revList.reverse()
                    try: row[20]=(j+1)*20-revList.index(testValue)-1
                    except:row[20]=-1
                    cursor.updateRow(row)
            outRaster = Lookup(outCombine, lookupField)
            outRaster.save('%sRECLAS_%i' %(outFolder,j))
            resultS.append('%sRECLAS_%i' %(outFolder,j))
            arcpy.Delete_management("COMBO")
        outCellStatistics = CellStatistics(resultS, "MAXIMUM", "NODATA")
        outCellStatistics.save('%sLATEST_%i' %(outFolder,testValue))
        for j in range(4):
            arcpy.Delete_management('%sCOMBO_%i' %(outFolder,j))
            arcpy.Delete_management('%sRECLAS_%i' %(outFolder,j))

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