Calculate sum of 2334 raster layers

[MagnusLund]

Hello,

I have 2334 raster grids (Geodatabase raster dataset, 1 band, 16 bit signed integer) that I want to sum into one raster. These rasters all have different spatial extent. Due to the large number of grids, using Raster Calculator is not really an option.

I have tried using the Cell Statistics tool to combine these grids. I unchecked the "Ignore NoData in calculations", set the processing extent to "Union of inputs", however it doesn't seem to work. I only get the sum where grids overlap.

I have read in the forum that one can get around this problem by replacing NoData with zeros, ie by using Con( IsNull[gridname], 0, [gridname] ) in Raster Calculator. However, this will once again be a bit boring when having thousands of grids.

Is there any solution to this problem?

Best regards,
Magnus

[fabefabe]

I found the solution. Just incase if someone like me is wondering about this problem:

just modify the code posted by Ryan for line below:

if i == 0:  
        out2 = out1      
        i += 1
    else:
        out2 = out2 + out1
        i += 1

#save final output
out2.save(os.path.join(outPath,'sumRas'))

[MichelineSnively]

Thanks Ryan for posting the python code, I was able to modify and use it to add hundreds of rasters together.  The following is the modified code for my purposes, hopefully it endevours to help someone else out in the future.

import arcpy
import os
arcpy.CheckOutExtension('Spatial')
arcpy.env.scratchWorkspace = "datapath"
arcpy.env.workspace = "datapath"
#create a list of rasters in the workspace
rasters = arcpy.ListRasters("*", "GRID")

i = 0
#loop through rasters in list
for raster in rasters:
    print "processing raster: %s" %os.path.join("datapath",raster)
     #get first raster
    out1 = raster
     #sum rasters together
    if i == 0:
         out2 = arcpy.Raster(out1)
         i += 1
    else:
         out2 = out2 + out1
         i += 1
#save final output       
out2.save(os.path.join("outpath",'sumRas'))

- Micheline

[RyanDeBruyn]

Sami and Micheline, glad this worked out for you.

I just want to make note of the importance of the Raster object in this case. All Spatial Analyst tools/functions return a raster object as a result. You can use this raster (variable) as input to another tool or query its properties for use. This is what I did in my initial example for Magnus from the Con(), and Sami found the error in the syntax (my mistake as out1 is already a raster object and didn't need to cast it).

Micheline, in your example you are setting a variable "out1 = raster" to the name (string) of the raster in the list. This then is being cast as a raster object in the sum logic.

If you need a raster object you can cast one by using the Raster() method, this is required in some cases and here it is neccesary to intialize the first raster in the list so we can sum the Rasters together using the '+' operator in a map algebra expression.

I have modified the code a bit, to cast the raster directly instead of creating the extra variable. Perhaps this will help you in your future workflows.

    import arcpy
    import os
    arcpy.CheckOutExtension('Spatial')
    arcpy.env.scratchWorkspace = r"c:\temp\scratch.gdb"
    arcpy.env.workspace = r"c:\temp\work.gdb"
    #create a list of rasters in the workspace
    rasters = arcpy.ListRasters("*", "GRID")

    i = 0
    #loop through rasters in list
    for raster in rasters:
        print ("processing raster: %s" %os.path.join("datapath",raster))
        #sum rasters together
        if i == 0:
            outSUM = arcpy.Raster(raster)
            i += 1
        else:
            outSUM = outSUM + raster
            i += 1

    #save final output to the current workspace
    outSUM.save('sumRas')

Good luck
-R

[JerryGarcia1]

I am trying to do something similar.  Is there a script I can use for 9.3.1?

[SarahBurns]

I have reclassified a raster to be only 0's and 1's and I would like to sum the 1's together.
I have used the code from this thread but my output is not correct and python keeps shutting down each time I run it so I do not even get an error message.
Can someone please tell me what I am doing wrong.

Thanks

# Import system modules
import arcpy
from arcpy import env
from arcpy.sa import *
import os

# Set environment settings
arcpy.env.workspace = "D:/sumrasters/rasterimgs"
outputfolder="D:/sumrasters/output_sumRaster"
datapath="D:/DSE_work/reclass/sumrasters/rasterimgs"
arcpy.env.overwriteOutput = 1
arcpy.CheckOutExtension('Spatial')
arcpy.env.scratchWorkspace = r"c:\temp\scratch.gdb"


#create a list of rasters in the workspace
rasters = arcpy.ListRasters("*","IMG")

i = 0
#loop through rasters in list
for raster in rasters:
    print rasters
    #sum rasters together
    if i == 0:
        outSUM = arcpy.Raster(raster)
        i += 1
    else:
        outSUM = outSUM + raster
        i += 0

#save final output to the current workspace
    outSUM.save(os.path.join(outputfolder,"sumRas.img"))

print "end of processing"

[deleted-user-IR249IovB3CN]

The issue of adding rasters of different extents does not seem to be adressed fully here. There are many posts addressing this issue but none taht I have found completely solve the problem, or I am setting something incorrectly.  I have changed all NoData values to 0, set the processing extent to "Union of Inputs" but when adding the rasters, the output raster is only the extent of the intersection of all rasters.

How do I add rasters of different extents to "sum" all values for all rasters to the full extent of the largest raster?

[paulcripps]

Is there some limit on processing? The code works perfectly for smaller numbers of inputs but my main project has 2500 rasters. I've added some debug addMessage statements and it iterates fines through the calculation, seems to be falling over at the save statement. Any ideas?

[StephenTeet1]

I am getting a name error that I am not sure how to resolve:

Traceback (most recent call last):
  File "C:\Python26\ArcGIS10.0\Lib\site-packages\Pythonwin\pywin\framework\scriptutils.py", line 312, in RunScript
    exec codeObject in __main__.__dict__
  File "F:\Thesis\BME\S\BME\sumrasters2.py", line 26, in <module>
    out2.save(os.path.join("outPath",'sumRas'))
NameError: name 'out2' is not defined

is out2 not defined by "out2 = arcpy.Raster(out1)" and "out2 = out2 + out1"?

I am new to using python, any advice on where my problem may be?

-S

[MarthaSnyderwine]

Hello!  I can sum the rasters but often the results are +/ 3.4 e+038.  Which is the NoData value.  I tried resetting all the NoData values to 0, but then this effects the min, max, and average which is unsat.  I've inspected the input shapfiles and individual rasters.  I don't see any extreme values.  

Any other suggestions?

[PaulLeonard]

Hello.

Since ascii files are soooo slow.  I have a question.  Can you provide guidance for doing this exact same process with .gzip ascii files?  I was wondering if it would be possible without decompressing.

-Cheers

Hi Magnus,  here is abit of code that could help get you started. I didn't test it but give it a try (see attached file). 

good luck
-Ryan

arcpy.env.overwriteOutput = 1
arcpy.CheckOutExtension('Spatial')
arcpy.env.scratchWorkspace = outPath
arcpy.env.workspace = dataPath
#create a list of rasters in the workspace
rasters = arcpy.ListRasters('','')

i = 0
#loop through rasters in list
for raster in rasters:
    print "processing raster: %s" %os.path.join(dataPath,raster)

    #convert nodata to zero
    out1 = Con(IsNull(raster), 0, raster)

    #sum rasters together
    if i == 0:
        out2 = arcpy.Raster(out1)
        i += 1
    else:
        out2 = out2 + out1
        i += 1

#save final output
out.save(os.path.join(outPath,'sumRas'))