topic Re: Sum 6,000+ rasters - Python script in Python Questions

Sum 6,000+ rasters - Python script

MathieuLacorde1 — Tue, 02 Feb 2016 10:05:24 GMT

Hi there,

I am trying to sum close to 6,000 ESRI grids and found a piece of script on this forum that I have tried to adapt (attached). The script runs fine and prompts every grid but eventually crashes when it reaches the last raster. No output is created.

Would you have any idea as to why it crashes? The script seems to work fine when I test it on 5-10 grids. Could the reason be the 5,000 ESRI grid limit that can be stored in a single workspace directory (Esri Grid format—Help | ArcGIS for Desktop) ? I did try on less than 5,000 grids but had the same outcome.

I am using ArcGIS 10.2.2 and Python 2.7.5.

Thanks for your help!

Mathieu

Re: Sum 6,000+ rasters - Python script

DanPatterson_Retired — Tue, 02 Feb 2016 11:10:31 GMT

the download has a ...py.zip filename, but the zip file contains an addin. do you have the script?

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Tue, 02 Feb 2016 11:28:49 GMT

Hi Dan,

The zip file contains the *.py that I edited and ran within IDLE. Here is the script:

import os,sys

import arcpy

from arcpy.sa import *

sPath = sys.path[0]

dataPath = 'E:\BOM\daily-temperature\daily-maximum-temperature\max_temp_1971_2000_grids' # ADD your workspace path here

outPath = 'E:\BOM\daily-temperature\daily-maximum-temperature'

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 = out1

i += 1

else:

out2 = out2 + out1

i += 1

#save final output

out2.save(os.path.join(outPath,'sumRas2'))

Thanks,

Mathieu

Re: Sum 6,000+ rasters - Python script

XanderBakker — Sun, 12 Dec 2021 00:02:12 GMT

Just to make the thread more readable I have included the script you attached:

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

sPath = sys.path[0]
dataPath = 'E:\BOM\daily-temperature\daily-maximum-temperature\max_temp_1971_2000_grids' # ADD your workspace path here
outPath = 'E:\BOM\daily-temperature\daily-maximum-temperature'

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 = out1
        i += 1
    else:
        out2 = out2 + out1
        i += 1

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

The first thing I notice is the long path with multiple minus signs. This can cause problems when working with grid files. The number of files in a single folder can be the cause of problems too. I also notice that for each raster you first set the NoData cells to 0 and then add it to the sum. Have you considered using Cell Statistics—Help | ArcGIS for Desktop with the SUM statistics type and the DATA option. This allows you to sum multiple rasters and ignore the NoData cells (they will be considered 0 in the SUM operation). I don't think you will be able to use 6000 rasters as inputs but you can chunk the list and perform it on smaller batches and afterwards on the intermediate results.

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Tue, 02 Feb 2016 11:49:03 GMT

Hi Xander,

I did try without setting the NoData cells to Zero and it didn't change the result. The Cell Statistics tool only allows 1,000 rasters if I am not mistaken so I could give it a go in several batches. The thing is I have to run this script on another 10,000 files and I was hoping to get it to work on the whole workspace.

It looks like I will have to subdivide my data though. I'll try to shorten the path too.

Thanks

Re: Sum 6,000+ rasters - Python script

XanderBakker — Sun, 12 Dec 2021 00:02:15 GMT

When I mentioned using chunk, I was referring to let Python divide the list of rasters into chunks of less than 1000 rasters, perform the sum and store the result temporarily and at the end of the chunks, combine the temporal results into one.

An example of how to "chunk":

def main():
    my_list = ["Raster{0}".format(i) for i in range(100)]

    for my_list_chunk in chunks(my_list, 10):
        print my_list_chunk

def chunks(l, n):
    """ Yield successive n-sized chunks from l."""
    for i in xrange(0, len(l), n):
        yield l

if __name__ == '__main__':
    main()

This will yield 10 lists with 10 elements each:

['Raster0', 'Raster1', 'Raster2', 'Raster3', 'Raster4', 'Raster5', 'Raster6', 'Raster7', 'Raster8', 'Raster9']
['Raster10', 'Raster11', 'Raster12', 'Raster13', 'Raster14', 'Raster15', 'Raster16', 'Raster17', 'Raster18', 'Raster19']
['Raster20', 'Raster21', 'Raster22', 'Raster23', 'Raster24', 'Raster25', 'Raster26', 'Raster27', 'Raster28', 'Raster29']
['Raster30', 'Raster31', 'Raster32', 'Raster33', 'Raster34', 'Raster35', 'Raster36', 'Raster37', 'Raster38', 'Raster39']
['Raster40', 'Raster41', 'Raster42', 'Raster43', 'Raster44', 'Raster45', 'Raster46', 'Raster47', 'Raster48', 'Raster49']
['Raster50', 'Raster51', 'Raster52', 'Raster53', 'Raster54', 'Raster55', 'Raster56', 'Raster57', 'Raster58', 'Raster59']
['Raster60', 'Raster61', 'Raster62', 'Raster63', 'Raster64', 'Raster65', 'Raster66', 'Raster67', 'Raster68', 'Raster69']
['Raster70', 'Raster71', 'Raster72', 'Raster73', 'Raster74', 'Raster75', 'Raster76', 'Raster77', 'Raster78', 'Raster79']
['Raster80', 'Raster81', 'Raster82', 'Raster83', 'Raster84', 'Raster85', 'Raster86', 'Raster87', 'Raster88', 'Raster89']
['Raster90', 'Raster91', 'Raster92', 'Raster93', 'Raster94', 'Raster95', 'Raster96', 'Raster97', 'Raster98', 'Raster99']

One question, I don't think your goal is to get the total sum of maximum temperatures, right? You are probably looking for an average (dividing the sum by the number of rasters)? The Cell Statistics also comes with the MEAN statistics type and other perhaps interesting statistics types.

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Tue, 02 Feb 2016 12:09:34 GMT

Thanks for that, I wouldn't have thought about it!

You're right, my goal is not to get the total sum but it is a first step into a more complex calculation. Once I get this script running smoothly, I will fine tune the calculation.

Cheers

Re: Sum 6,000+ rasters - Python script

XanderBakker — Tue, 02 Feb 2016 12:18:47 GMT

I noticed another thing... the way you specify your paths. This will probably yield problems. There is a nice post by Dan Patterson that explains in more detail about paths in Python: Filenames and file paths in Python .

I usually prefer the raw notation (put an r in front of the path, outside the quotes).

Re: Sum 6,000+ rasters - Python script

XanderBakker — Sun, 12 Dec 2021 00:02:17 GMT

To give you a start (code below is not tested!):

def main():
    import os
    import arcpy

    data_path = r'E:\BOM\daily-temperature\daily-maximum-temperature\max_temp_1971_2000_grids'
    # out_path = r'E:\BOM\daily-temperature\daily-maximum-temperature'
    res_ws = r'E:\BOM\daily-temperature\daily-maximum-temperature\yourGDB.gdb' # I prefer gdb over folder

    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension('Spatial')
    arcpy.env.workspace = data_path

    # create fgdb res_ws
    path, name = os.path.split(res_ws)
    arcpy.CreateFileGDB_management(path, name)

    #create a list of rasters in the workspace
    rasters = arcpy.ListRasters()

    i = 0
    lst_cellstat = []
    for rasters_chunk in chunks(rasters, 500):
        i += 1
        cellstat = arcpy.sa.CellStatistics(rasters_chunk, "SUM", "DATA")
        outname = os.path.join(res_ws, "cellstat{0}".format(i))
        cellstat.save(outname)
        lst_cellstat.append(outname)

    # create the final raster
    cellstat = arcpy.sa.CellStatistics(lst_cellstat, "SUM", "DATA")
    outname = os.path.join(res_ws, "total_sum")
    cellstat.save(outname)

    # eliminate intermediate results
    for ras in lst_cellstat:
        arcpy.Delete_management(ras)

def chunks(l, n):
    """ Yield successive n-sized chunks from l."""
    for i in xrange(0, len(l), n):
        yield l

if __name__ == '__main__':
    main()

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Tue, 02 Feb 2016 13:53:16 GMT

I'll try it as soon as I can!

Re: Sum 6,000+ rasters - Python script

JamesCrandall — Tue, 02 Feb 2016 13:58:48 GMT

res_ws = r'E:\BOM\daily-temperature\daily-maximum-temperature\yourGDB.gdb' # I prefer gdb over folder

I agree. Not sure how many times it has occurred but simply moving things into a FGDB resolves various issues I've run into. Pretty much I've just accepted the fact that there's probably a darn good reason why ESRI engineers have fully integrated the fgdb (Default.gdb) into the ArcGIS product, so why not follow?

Re: Sum 6,000+ rasters - Python script

ChrisDonohue__GISP — Tue, 02 Feb 2016 16:05:02 GMT

I agree with Xander Bakker that the dashes in the file path may be an issue. Maybe try using a different folder structure that does not have dashes?

A general suggestion. With Grid files, shorter is usually better in terms of paths and filenames. Oftentimes if the combined length of filename and path is long, there can be issues. Yet the same process will often work fine if the path is simplified to be very short. In theory it shouldn't matter, but there can be a difference. So try a shorter path, that may resolve it.

ArcGIS Help (10.2, 10.2.1, and 10.2.2) - ESRI Grid format

Chris Donohue, GISP

Re: Sum 6,000+ rasters - Python script

TonySoeller — Thu, 04 Feb 2016 18:56:32 GMT

Your computer might not have enough memory. As Xander Bakker states, try smaller chunks of rasters.

Re: Sum 6,000+ rasters - Python script

AndreasSutter2 — Fri, 05 Feb 2016 07:07:57 GMT

Hi Mathieu

I would also recommend trying the FGDB ideas given above. And if Xander Bakker's good looking code doesn't do the trick, I have some further options that might be worthy to check out. We ran into similar problems at work while processing a large number of (raster)data in Python scripts. The following helped us solving similar issues as you are currently experiencing:

Monitor your computers memory while the tool runs using Task Manager. Should there be a memory issue, or if the above ideas using the FGDB does not work try the following:

Write a .BAT file that accesses the Python.exe and then opens your script. Insert all the necessary arguments in the command line script. You can prepare the batches using an Editor or Excel.
This way you'll be able to process a large number of iterations of your script outside the ArcCatalog environment, which (in case there are any in the first place) would solve the memory issues, because running a BAT-script means that the Python environment is closed after each iteration and loaded again for the next.

Hope this helps.

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Sun, 06 Mar 2016 09:09:46 GMT

Hi,

Sorry to only get back to you now, I just had a chance to look at the script today... I removed the dashes and tried to use shorter pathnames. I ran the script and got an error message:

File "E:/Software/new_script_geonet.py", line 44, in <module>

main()

File "E:/Software/new_script_geonet.py", line 26, in main

cellstat.save(outname)

RuntimeError: ERROR 010240: Could not save raster dataset to E:\BOM\daily_temp\daily_max_temp\MaxTemp.gdb\cellstat1 with output format FGDBR.

Any ideas on what I could modify in the script to get it right?

Thanks!

Mathieu

Re: Sum 6,000+ rasters - Python script

DanPatterson_Retired — Sun, 06 Mar 2016 11:54:55 GMT

try a different format... just add .tif to the output name

Re: Sum 6,000+ rasters - Python script

MathieuLacorde1 — Sun, 06 Mar 2016 12:55:22 GMT

Hi there,

Struggling a bit as to where I should add the extension ".tif". I would say:

outname = os.path.join(res_ws, "cellstat{0}".format(i).tif)

However, I still get an error message, different one though:

File "E:\Software\new_script_geonet_revB.py", line 45, in <module>

main()

File "E:\Software\new_script_geonet_revB.py", line 26, in main

outname = os.path.join(res_ws, "cellstat{0}".format(i).tif)

AttributeError: 'str' object has no attribute 'tif'

Re: Sum 6,000+ rasters - Python script

DanPatterson_Retired — Sun, 06 Mar 2016 14:11:39 GMT

Well according to this Cell Statistics—Help | ArcGIS for Desktop

you should be able to save to an Esri grid format (with no file extension given)

But the only example that they give there is one with an img or grid format, but saved to a folder rather than a geodatabase. So, I don't know why it won't save to the geodatabase but I do know that esri grids don't like anything in their paths other than letters, the underscore and maybe the odd number (not beginning any folder segment), although Xander's example suggests otherwise. So you can try saving to a folder until Xander Bakker can comment

Re: Sum 6,000+ rasters - Python script

XanderBakker — Sun, 06 Mar 2016 19:11:13 GMT

The error refers to not being able to write to the output File Geodatabase. The Help states:

Check that a raster with the same name and format does not already exist in the output location. Also, check the Help for the technical specifications of raster dataset formats to make sure that the expected range of values in the output is compatible with the specified format

Just to be sure "MaxTemp.gdb" is a file geodatabase and you have privileges to write to that location?

And at what moment does the error occur,when writing the first raster to that location or when part of the data has been processed?

You could include "arcpy.env.overwriteOutput = True" at the beging of the scipt (after "import arcpy")

Re: Sum 6,000+ rasters - Python script

XanderBakker — Sun, 06 Mar 2016 19:13:28 GMT

The line you included doesn't add the ".tif" correctly. Change:

outname = os.path.join(res_ws, "cellstat{0}".format(i).tif)

into

outname = os.path.join(res_ws, "cellstat{0}.tif".format(i))

and be sure that the res_ws points to a folder and not a File Geodatabase.