1

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

  • Multiple rasters? How many? – FelixIP Mar 17 '16 at 0:37
  • I have about 80 rasters (i.e. 80 years). – SilvanD Mar 17 '16 at 13:23
  • 3
    EqualToFrequency, applied once for each possible value, solves (1), while HighestPosition when applied to indicators of each value will solve (2). You will want to program the creation of the indicator rasters. – whuber Mar 18 '16 at 19:25
  • I ended up using equaltofrequency for objective (1), as opposed to the answer below, which was proving to be too long (it is taking over 3 hours to process my entire dataset, and has not yet finished). However, for objective (2) highestposition doesn't seem to do the trick for the following reason: assuming a particular cell has the same value for multiple years, highestposition returns the position of the first raster, and not the last, which is what i'm looking for. Any other ideas on how to address that second objective? – SilvanD Mar 20 '16 at 18:29
  • To expand upon my previous comment-- to cite the help from the tool HighestPosition, If two or more input rasters contain the maximum value for a particular cell location, the position of the first one is returned on the output raster. I need the position of the last one to be returned... Any thoughts? – SilvanD Mar 20 '16 at 19:10
2

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:

enter image description here

Histograms:

enter image description here

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()        
  • Thanks for the reply, FelixIP. I would like to do this using a script because I'm going to have to run it quite a number of times. I've taken a look at the link, and am still not sure how to proceed once I've converted the raster to a numpy array. Any further elaboration/clarification would be greatly appreciated. – SilvanD Mar 17 '16 at 13:25
  • See updated answer, you'll have to modify in and out folder names – FelixIP Mar 17 '16 at 21:33
  • Thank you very much for expanding upon your answer. ArcMap crashes everytime I try your code. Not sure how to get around that issue. I will keep on trying though! – SilvanD Mar 20 '16 at 13:55
  • FANTASTIC. After waiting some 20 minutes, it finally worked, and yields exactly what I was hoping for. Beautiful, elegant answer. Thank you. My only question is what does "myRaster" correspond to? The generated "totals" and "latest" files address objectives 1 and 2, respectively, but it is not clear what the resulting myRaster refers to. – SilvanD Mar 20 '16 at 15:39
  • Myraster=... line Is just the way to save totals and latest to disk – FelixIP Mar 20 '16 at 17:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.