topic Re: Adding Together Rasters in Geoprocessing Questions

Adding Together Rasters

Anonymous User — Tue, 20 Jul 2010 12:41:51 GMT

Original User: Kenyon07

Hi,
I am trying to add together a bunch or rasters (about 300), and I know of no way other than using the raster calculater to do this. (Ex. [RASHFT0.tif] + [RASHFT0.tif] + [RASHFT0.tif] )

But this seems not to work. How can I truely add together all of these rasters?

Attached is a text file that I put the values into the raster calculator to do the calculation (Adding.txt) and also a picture of my screen(Calculation.jpg). The calculation shows values in the TOC but returns NoData.

Re: Adding Together Rasters

KenyonAbbott — Tue, 20 Jul 2010 20:07:45 GMT

Bump. Anybody?

Re: Adding Together Rasters

Anonymous User — Tue, 20 Jul 2010 20:50:41 GMT

Original User: csny490

For the addition to work (or pretty much any other raster math), your rasters must have the same spatial extent. So basically (I assume these are the light poles) for each raster you need to create a new raster where the NoData within your study area is 0 (or some other non NoData value). Here's some Python code that might give you some ideas:

# Name: make_fast_mosaic_v93.py

#
# Description
# -----------
# This script will mosaic a binch of grids, but using a secret method that is
# mucho faster than the standard v9.3.1 SP1 mosaic tool.
#
# Written By: Chris Snyder, WA DNR, 03/01/2010, chris.snyder(at)wadnr.gov
#
# Written For: Python 2.5.1 and ArcGIS v9.3.1 SP1
#
# UPDATES:
#
# Notes on input parameters (for the toolbox):
# VARIABLE             PAREMETER_INDEX     PARAMETER_DATA_TYPE
# -------------------------------------------------------------------
# gridList             0                   Grids
# outputGrid           1                   Grid
# tempWorkspace        2                   Workspace
# mosaicMethod         3                   String (MAX, MIN, or MEAN)
# tempNoDataValue      4                   Long Integer

try:
    #Process: Import some modules
    import os, string, sys, time, traceback, arcgisscripting

    #Process: Create the gp object
    gp = arcgisscripting.create(9.3)

    #Process: Defines some functions used for getting messages from the gp and python
    def showGpMessage():
        gp.AddMessage(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpWarning():
        gp.AddWarning(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpError():
        gp.AddError(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showPyLog(): #just print to the log file!
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
    def showPyMessage():
        gp.AddMessage(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyWarning():
        gp.AddWarning(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyError():
        gp.AddError(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message

    #Specifies the root directory variable, defines the logFile variable, and does some minor error checking...
    dateTimeString = str(time.strftime('%Y%m%d%H%M%S'))
    scriptName = os.path.split(sys.argv[0])[-1].split(".")[0]
    userName = string.lower(os.environ.get("USERNAME")).replace(" ","_").replace(".","_")
    tempPathDir = os.environ["TEMP"]
    logFileDirectory = r"\\snarf\am\div_lm\ds\gis\tools\log_files"
    if os.path.exists(logFileDirectory) == True:
        logFile = os.path.join(logFileDirectory, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
        try:
            print >> open(logFile, 'a'), "Write test successfull!"
        except:
            logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
    else:
        logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
    if os.path.exists(logFile)== True:
        os.remove(logFile)
        message = "Created log file " + logFile; showPyMessage()
    message = "Running " + sys.argv[0]; showPyMessage()

    #Process: Attempts to check out a Spatial Analyst license
    try:
        gp.CheckOutExtension("Spatial")
    except:
        message = "ERROR: Spatial analyst license is unavailable... Exiting script!"; showPyError(); sys.exit()

    #Process: Check out the highest license available
    try:
        if gp.CheckProduct("ArcView") == "Available":
            gp.SetProduct("ArcView")
        elif gp.CheckProduct("ArcEditor") == "Available":
            gp.SetProduct("ArcEditor")
        elif gp.CheckProduct("ArcInfo") == "Available":
            gp.SetProduct("ArcInfo")
    except:
        message = "ERROR: Could not select an ArcGIS license level! Exiting script..."; showPyError(); sys.exit()
    message = "Selected an " + gp.ProductInfo() + " license"; showPyMessage()

    #Process: Sets some gp environment variables
    gp.overwriteoutput = True
    gp.pyramid = "NONE"
    gp.rasterStatistics = "NONE"

    #Process: Collect the input parameters
##    gridList = r"D:\csny490\usgs_dem10m\dem11746;D:\csny490\usgs_dem10m\dem11747"
##    outputGrid = r"D:\csny490\temp\temp\howaboutthis"
##    processingWorkspace = r"D:\csny490\temp\temp"
##    mosaicMethod = "MAX"
##    tempNoDataValue = "-32000"
    gridList = gp.GetParameterAsText(0)
    outputGrid = gp.GetParameterAsText(1)
    processingWorkspace = gp.GetParameterAsText(2)
    mosaicMethod = gp.GetParameterAsText(3)
    tempNoDataValue = gp.GetParameterAsText(4)

    #Process: Print out the input parameters
    message = "INPUT PARAMETERS"; showPyMessage()
    message = "----------------"; showPyMessage()
    message = "Input Grids          = " + gridList; showPyMessage()
    message = "Output Grid          = " + outputGrid; showPyMessage()
    message = "Processing Workspace = " + processingWorkspace; showPyMessage()
    message = "Mosaic Method        = " + mosaicMethod; showPyMessage()
    message = "Temp NoData Value    = " + tempNoDataValue + "\n"; showPyMessage()

    #Process: Figure out the maximum extent of all the input grids
    message = "Calculating maximum extent of input rasters..."; showPyMessage()
    i = 0
    gridsToProcessList = gridList.split(";")
    for grid in gridsToProcessList:
        i = i + 1
        dsc = gp.describe(grid)
        if i == 1:
            xMin = dsc.extent.xmin
            yMin = dsc.extent.ymin
            xMax = dsc.extent.xmax
            yMax = dsc.extent.ymax
        else:
            if dsc.extent.xmin < xMin:
                xMin = dsc.extent.xmin
            if dsc.extent.xmax > xMax:
                xMax = dsc.extent.xmax
            if dsc.extent.ymin < yMin:
                yMin = dsc.extent.ymin
            if dsc.extent.ymax > yMax:
                yMax = dsc.extent.ymax

    gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax)
    message = "Maximum extent of input rasters is " + str(gp.extent) + "..."; showPyMessage()

    #Process: Now build some new grids that have the max extent
    maxString = ""
    tempGridList = []
    i = 0
    tempGridDict = {}
    for grid in gridsToProcessList:
        i = i + 1
        message = "Normalizing extent of " + grid + "..."; showPyMessage()
        dsc = gp.describe(grid)
        if dsc.format != "GRID": #if it aint a grid, make it one
            tempConvertGrid = gp.createscratchname("","","RASTER",processingWorkspace)
            gp.extent = ""
            gp.CopyRaster_management(grid, tempConvertGrid, "", "", "", "", "", "")
            tempGridList.append(tempConvertGrid)
            somaExp = "con(isnull(" + tempConvertGrid + "), " + str(tempNoDataValue) + ", " + tempConvertGrid + ")"
        else:
            somaExp = "con(isnull(" + grid + "), " + str(tempNoDataValue) + ", " + grid + ")"
        tempGrid = gp.createscratchname("","","RASTER", processingWorkspace)
        tempGridList.append(tempGrid)
        maxString = maxString + tempGrid + ","
        gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax)
        gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid)

    #Process: Max, then put back to null
    message = "Mosaicing - This will take a while..."; showPyMessage()
    tempGrid = gp.createscratchname("","","RASTER",processingWorkspace)
    tempGridList.append(tempGrid)
    if mosaicMethod == "MIN":
        somaExp = "min(" + maxString + ")"
    if mosaicMethod == "MAX":
        somaExp = "max(" + maxString + ")"
    if mosaicMethod == "MEAN":
        somaExp = "mean(" + maxString + ")"
    gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid)
    somaExp = somaExp = "con(" + tempGrid + " == " + str(tempNoDataValue) + ", setnull(1), " + tempGrid + ")"
    if outputGrid.split(".")[-1] in (".bmp",".gif",".img",".jpg",".jp2",".png",".tif") or\
       ".gdb" in outputGrid or ".mdb" in outputGrid or ".sde" in outputGrid:
        tempMosaicGrid = gp.createscratchname("","","RASTER", processingWorkspace)
        tempGridList.append(tempMosaicGrid)
        gp.SingleOutputMapAlgebra_sa(somaExp, tempMosaicGrid)
        gp.CopyRaster_management(tempMosaicGrid, outputGrid, "", "", "", "", "", "")
    else:
        gp.SingleOutputMapAlgebra_sa(somaExp, outputGrid)

    #Process: Delete all the temp grids
    message = "Removing temporay processing grids..."; showPyMessage()
    for tempGrid in tempGridList:
        try:
            gp.Delete_management(tempGrid, "")
        except:
            message = "WARNING: Failed to delete processing grid " + str(tempGrid) + "..."; showPyWarning()

    message = "All done!"; showPyMessage()

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

Re: Adding Together Rasters

ChrisSnyder — Sat, 11 Dec 2021 06:24:58 GMT

For the addition to work (or pretty much any other raster math), your rasters must have the same spatial extent. So basically (I assume these are the light poles) for each raster you need to create a new raster where the NoData within your study area is 0 (or some other non NoData value). Here's some Python code that might give you some ideas:

# Name: make_fast_mosaic_v93.py
#
# Description
# -----------
# This script will mosaic a bunch of grids, but using a secret method that is
# mucho faster than the standard v9.3.1 SP1 mosaic tool.
#
# Written By: Chris Snyder, WA DNR, 03/01/2010, chris.snyder(at)wadnr.gov
#
# Written For: Python 2.5.1 and ArcGIS v9.3.1 SP1
#
# UPDATES:
#
# Notes on input parameters (for the toolbox):
# VARIABLE             PAREMETER_INDEX     PARAMETER_DATA_TYPE
# -------------------------------------------------------------------
# gridList             0                   Grids
# outputGrid           1                   Grid
# tempWorkspace        2                   Workspace
# mosaicMethod         3                   String (MAX, MIN, or MEAN)
# tempNoDataValue      4                   Long Integer

try:
    #Process: Import some modules
    import os, string, sys, time, traceback, arcgisscripting

    #Process: Create the gp object
    gp = arcgisscripting.create(9.3)

    #Process: Defines some functions used for getting messages from the gp and python
    def showGpMessage():
        gp.AddMessage(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpWarning():
        gp.AddWarning(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showGpError():
        gp.AddError(gp.GetMessages())
        print >> open(logFile, 'a'), gp.GetMessages()
        print gp.GetMessages()
    def showPyLog(): #just print to the log file!
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
    def showPyMessage():
        gp.AddMessage(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyWarning():
        gp.AddWarning(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message
    def showPyError():
        gp.AddError(str(time.ctime()) + " - " + message)
        print >> open(logFile, 'a'), str(time.ctime()) + " - " + message
        print str(time.ctime()) + " - " + message

    #Specifies the root directory variable, defines the logFile variable, and does some minor error checking...
    dateTimeString = str(time.strftime('%Y%m%d%H%M%S'))
    scriptName = os.path.split(sys.argv[0])[-1].split(".")[0]
    userName = string.lower(os.environ.get("USERNAME")).replace(" ","_").replace(".","_")
    tempPathDir = os.environ["TEMP"]
    logFileDirectory = r"\\snarf\am\div_lm\ds\gis\tools\log_files"
    if os.path.exists(logFileDirectory) == True:
        logFile = os.path.join(logFileDirectory, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
        try:
            print >> open(logFile, 'a'), "Write test successfull!"
        except:
            logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")  
    else:
        logFile = os.path.join(tempPathDir, scriptName + "_" + userName + "_" + dateTimeString + ".txt")
    if os.path.exists(logFile)== True:
        os.remove(logFile)
        message = "Created log file " + logFile; showPyMessage()
    message = "Running " + sys.argv[0]; showPyMessage()

    #Process: Attempts to check out a Spatial Analyst license
    try:
        gp.CheckOutExtension("Spatial")
    except:
        message = "ERROR: Spatial analyst license is unavailable... Exiting script!"; showPyError(); sys.exit()
    
    #Process: Check out the highest license available
    try:
        if gp.CheckProduct("ArcView") == "Available":
            gp.SetProduct("ArcView")
        elif gp.CheckProduct("ArcEditor") == "Available":
            gp.SetProduct("ArcEditor")
        elif gp.CheckProduct("ArcInfo") == "Available":
            gp.SetProduct("ArcInfo")
    except:
        message = "ERROR: Could not select an ArcGIS license level! Exiting script..."; showPyError(); sys.exit()
    message =  "Selected an " + gp.ProductInfo() + " license"; showPyMessage()

    #Process: Sets some gp environment variables
    gp.overwriteoutput = True
    gp.pyramid = "NONE"
    gp.rasterStatistics = "NONE"

    #Process: Collect the input parameters
##    gridList = r"D:\csny490\usgs_dem10m\dem11746;D:\csny490\usgs_dem10m\dem11747"
##    outputGrid = r"D:\csny490\temp\temp\howaboutthis"
##    processingWorkspace = r"D:\csny490\temp\temp"
##    mosaicMethod = "MAX"
##    tempNoDataValue = "-32000"
    gridList = gp.GetParameterAsText(0)
    outputGrid = gp.GetParameterAsText(1)
    processingWorkspace = gp.GetParameterAsText(2)
    mosaicMethod = gp.GetParameterAsText(3)
    tempNoDataValue = gp.GetParameterAsText(4)

    #Process: Print out the input parameters
    message = "INPUT PARAMETERS"; showPyMessage()
    message = "----------------"; showPyMessage()
    message = "Input Grids          = " + gridList; showPyMessage()
    message = "Output Grid          = " + outputGrid; showPyMessage()
    message = "Processing Workspace = " + processingWorkspace; showPyMessage()
    message = "Mosaic Method        = " + mosaicMethod; showPyMessage()
    message = "Temp NoData Value    = " + tempNoDataValue + "\n"; showPyMessage()
    
    #Process: Figure out the maximum extent of all the input grids
    message = "Calculating maximum extent of input rasters..."; showPyMessage()
    i = 0
    gridsToProcessList = gridList.split(";") 
    for grid in gridsToProcessList:
        i = i + 1
        dsc = gp.describe(grid)
        if i == 1:
            xMin = dsc.extent.xmin
            yMin = dsc.extent.ymin
            xMax = dsc.extent.xmax
            yMax = dsc.extent.ymax
        else:
            if dsc.extent.xmin < xMin:
                xMin = dsc.extent.xmin
            if dsc.extent.xmax > xMax:
                xMax = dsc.extent.xmax
            if dsc.extent.ymin < yMin:
                yMin = dsc.extent.ymin
            if dsc.extent.ymax > yMax:
                yMax = dsc.extent.ymax

    gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax)
    message = "Maximum extent of input rasters is " + str(gp.extent) + "..."; showPyMessage()

    #Process: Now build some new grids that have the max extent
    maxString = ""
    tempGridList = []
    i = 0
    tempGridDict = {}
    for grid in gridsToProcessList:
        i = i + 1
        message = "Normalizing extent of " + grid + "..."; showPyMessage()
        dsc = gp.describe(grid)
        if dsc.format != "GRID": #if it aint a grid, make it one
            tempConvertGrid = gp.createscratchname("","","RASTER",processingWorkspace)
            gp.extent = ""
            gp.CopyRaster_management(grid, tempConvertGrid, "", "", "", "", "", "")
            tempGridList.append(tempConvertGrid)
            somaExp = "con(isnull(" + tempConvertGrid + "), " + str(tempNoDataValue) + ", " + tempConvertGrid + ")"
        else:
            somaExp = "con(isnull(" + grid + "), " + str(tempNoDataValue) + ", " + grid + ")"
        tempGrid = gp.createscratchname("","","RASTER", processingWorkspace)
        tempGridList.append(tempGrid)
        maxString = maxString + tempGrid + ","
        gp.extent = str(xMin) + " " + str(yMin) + " " + str(xMax) + " " + str(yMax)
        gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid)
         
    #Process: Max, then put back to null
    message = "Mosaicing - This will take a while..."; showPyMessage()
    tempGrid = gp.createscratchname("","","RASTER",processingWorkspace)
    tempGridList.append(tempGrid)
    if mosaicMethod == "MIN":
        somaExp = "min(" + maxString + ")"
    if mosaicMethod == "MAX":
        somaExp = "max(" + maxString + ")"
    if mosaicMethod == "MEAN":
        somaExp = "mean(" + maxString + ")"
    gp.SingleOutputMapAlgebra_sa(somaExp, tempGrid)
    somaExp = somaExp = "con(" + tempGrid + " == " + str(tempNoDataValue) + ", setnull(1), " + tempGrid + ")"
    if outputGrid.split(".")[-1] in (".bmp",".gif",".img",".jpg",".jp2",".png",".tif") or\
       ".gdb" in outputGrid or ".mdb" in outputGrid or ".sde" in outputGrid:
        tempMosaicGrid = gp.createscratchname("","","RASTER", processingWorkspace)
        tempGridList.append(tempMosaicGrid)
        gp.SingleOutputMapAlgebra_sa(somaExp, tempMosaicGrid)
        gp.CopyRaster_management(tempMosaicGrid, outputGrid, "", "", "", "", "", "")
    else:
        gp.SingleOutputMapAlgebra_sa(somaExp, outputGrid)

    #Process: Delete all the temp grids
    message = "Removing temporay processing grids..."; showPyMessage()
    for tempGrid in tempGridList:
        try:
            gp.Delete_management(tempGrid, "")
        except:
            message = "WARNING: Failed to delete processing grid " + str(tempGrid) + "..."; showPyWarning() 
    
    message = "All done!"; showPyMessage()  
    
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()

Re: Adding Together Rasters

Anonymous User — Tue, 19 Jul 2011 14:27:42 GMT

Original User: susanwitherly

Have you tried adding just a few of your rasters together using the raster calculator? Might give you a sense of whether the problem is an extent issue or a format problem? There might be a limit to how much text the raster calculator can handle. Try using a script such as:

import arcpy
intCount = 0

MainWS = "C:\GISProj\RasterData"
from arcpy import env
from arcpy.sa import *
env.workspace = MainWS
arcpy.CheckOutExtension("Spatial")
print intCount
rasterList = arcpy.ListRasters()
for raster in rasterList:
      print raster
      if intCount == 0:
            output = raster
            intCount = intCount + 1
            print output
      else:
            intCount = intCount + 1
            output = arcpy.sa.Plus(output, raster)
output.save("C:\GISProj\Workspace\Processing.gdb\output")
print "Total Number of Rasters: " + str(intCount)

Note the loss of indentation in the python code!