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

[RyanDeBruyn]

Magnus,

What version are you using?

In model builder you can create variables which will iterate through rasters in a workspace which would help you automate the process for Con().

Also using a python script would provide you with an easy way to create a loop through rasters in a workspace do the Con() and sum the rasters together as you go. 

If you have version 10 this is even easier to do and will only result in your final output being created not 1000's of itermediate rasters.

-R

[MagnusLund]

Thanks! Yes I am using version 10. I am not used to Python however, if anyone has a possibility to providing me with a general script I would be very grateful. I do know other programming languages so I should be able to adapt to my own needs given I know the basic structure.

Best regards,
Magnus

[RyanDeBruyn]

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

[RyanDeBruyn]

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

[SimonAndersson2]

rasters = arcpy.ListRasters('','')
rastersum = rasters[0]
for raster in rasters[1:]:
    rastersum+=arcpy.Raster(raster)‍‍‍‍‍‍‍‍‍

I know this is a bit late but I came across this from google and thought I'd clean up your for loop.

[XanderBakker]

A few comments on your code:

• you probably need to change line 2 to: rastersum = arcpy.Raster(rasters[0]) in order to use the "+" operator with rasters
• The code does not account for NoData cells, NoData cells in any of the raster will propagate to NoData in the end result. To avoid this it would be better to use Cell Statistics—Help | ArcGIS Desktop with the SUM statistics type and the NODATA option
• In order to speed the process up, you could use "chunks". See an example here: python - How do you split a list into evenly sized chunks? - Stack Overflow 

[MagnusLund]

Thanks Ryan! I have tried it and get an error message that I cannot really understand:

Runtime error <type 'exceptions.TypeError'>: expected a raster or layer name

The routine does start well, the print command executes and out1 is created correctly. However, out2 is never created so I guess that's where the routine crashes. But out1 should be a raster name? Please see adapted code below:

Best regards,
Magnus


import os,sys
import arcpy
from arcpy.sa import *

sPath = sys.path[0]
dataPath = 'C:/footprintz/py_indata/results.gdb' # ADD your workspace path here
outPath = 'C:/footprintz/py_outdata/'

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

[MagnusLund]

Thanks Ryan! It worked great.

Best regards,
Magnus

[fabefabe]

Ryan,

I also have to add the rasters. I tried to work on your posted code but I am having a error message. It says:

Traceback (most recent call last):
  File "D:\Netcdf\sumrasters.py", line 31, in <module>
    out2 = arcpy.Raster(out1)
TypeError: expected a raster or layer name

when I tried to print the name of out1, it printed out:
D:/Netcdf/AS13894/rasterout\ifthe_ras

I am thinking that I have to trimmed out the file path and give the name of the file only. Am I right? But I have only limited knowledge on Python. So, I would really appreciate if you could help me on this.

Thanks
Sami