nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

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

Find more posts tagged with

Accepted answers

All comments

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

sumrasters.py.zip

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

sumrasters.py.zip

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:

/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

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

ChristopherStrother

I am having the same problem as Sarah. Every method I have attempted results in an output raster with values of all zeroes. I think this has something to do with the 1 bit depth of the binary input rasters. Anyone have any ideas? I've been working on this for several weeks with no success!

SimonAndersson2

rasters <SPAN class="operator token">=</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>ListRasters<SPAN class="punctuation token">(</SPAN><SPAN class="string token">''</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="string token">''</SPAN><SPAN class="punctuation token">)</SPAN>
rastersum <SPAN class="operator token">=</SPAN> rasters<SPAN class="punctuation token">[</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">]</SPAN>
<SPAN class="keyword token">for</SPAN> raster <SPAN class="keyword token">in</SPAN> rasters<SPAN class="punctuation token">[</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">:</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">:</SPAN>
&nbsp;&nbsp;&nbsp; rastersum<SPAN class="operator token">+=</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>Raster<SPAN class="punctuation token">(</SPAN>raster<SPAN class="punctuation token">)</SPAN>‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

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

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.