topic Re: Error in multiple raster processing in Python Questions

Error in multiple raster processing

DuminduJayasekera — Sat, 11 Dec 2021 12:50:14 GMT

Hi,

I am trying to calculate the average using multiple rasters using years from 1981 to 2016. For example, average for week 1 for the month of 9 from years 1981 to 2016 is as follows:

Normal_Week_1_1981_2016_9.tif =average [Week_1_Sum1981_9.tif

Week_1_Sum2016_9.tif]

My following code outputs all the tif files at the "out" path but at the end it gives me the error below. What I am doing wrong here.? Just noticed the output raster values are also wrong.

Any help is appreciated.

Thanks in advance.

arcpy.CheckOutExtension("Spatial")
arcpy.OverwriteOutput = True

years = [1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013,2014,2015,2016]

out = r'H:\PRISM_800m_weekly_sum\Normal_out'          

path = r'H:\PRISM_800m_weekly_sum_Normal'
lstFiles = []  
for root, dirs, files in os.walk(path):  
    arcpy.env.workspace = root  
    files = arcpy.ListRasters("*", "TIF")  
    for file in files:  
        (filename, extension) = os.path.splitext(file)  
        if(extension == ".tif"):  
            lstFiles.append(os.path.join(root, file))  
print (lstFiles)  

new_list1=[]
new_list2=[]
new_list3=[]
new_list4=[]
new_list5=[]

months=[9,10,11]
for month in months:
     for file in lstFiles:
          if fnmatch.fnmatch(file, "H:\\PRISM_800m_weekly_sum_Normal" + "*" + "\\" + "Week_1" + "_*_" + str(month) + ".tif"):
               new_list1.append(file)
               new_list1.sort()
               print (new_list1)
               
     env.workspace = out
     finras1 = out + "\\" + "Normal_Week_1" + "_1981_2016" + "_" + str(month) + ".tif"
     calc1 = arcpy.sa.CellStatistics(new_list1, statistics_type = "MEAN")
     arcpy.CopyRaster_management(calc1,finras1,"","","","","","32_BIT_FLOAT")

     for file in lstFiles:          
          if fnmatch.fnmatch(file, "H:\\PRISM_800m_weekly_sum_Normal" + "*" + "\\" + "Week_2" + "_*_" + str(month) + ".tif"):
               new_list2.append(file)
               new_list2.sort()
               print (new_list2)
               
     env.workspace = out
     finras2 = out + "\\" + "Normal_Week_2" + "_1981_2016" + "_" + str(month) + ".tif"
     calc2 = arcpy.sa.CellStatistics(new_list2, statistics_type = "MEAN")
     arcpy.CopyRaster_management(calc2,finras2,"","","","","","32_BIT_FLOAT")               
                    
               
     for file in lstFiles:               
          if fnmatch.fnmatch(file, "H:\\PRISM_800m_weekly_sum_Normal" + "*" + "\\" + "Week_3" + "_*_" + str(month) + ".tif"):
               new_list3.append(file)
               new_list3.sort()
               print (new_list3)
               
     env.workspace = out
     finras3 = out + "\\" + "Normal_Week_3" + "_1981_2016" + "_" + str(month) + ".tif"
     calc3 = arcpy.sa.CellStatistics(new_list3, statistics_type = "MEAN")
     arcpy.CopyRaster_management(calc3,finras3,"","","","","","32_BIT_FLOAT")
     
     for file in lstFiles:          
          if fnmatch.fnmatch(file, "H:\\PRISM_800m_weekly_sum_Normal" + "*" + "\\" + "Week_4" + "_*_" + str(month) + ".tif"):
               new_list4.append(file)
               new_list4.sort()
               print (new_list4)
               
     env.workspace = out
     finras4 = out + "\\" + "Normal_Week_4" + "_1981_2016" + "_" + str(month) + ".tif"
     calc4 = arcpy.sa.CellStatistics(new_list4, statistics_type = "MEAN")
     arcpy.CopyRaster_management(calc4,finras4,"","","","","","32_BIT_FLOAT")
     
     
     for file in lstFiles:          
          if fnmatch.fnmatch(file, "H:\\PRISM_800m_weekly_sum_Normal" + "*" + "\\" + "Week_5" + "_*_" + str(month) + ".tif"):
               new_list5.append(file)
               new_list5.sort()
               print (new_list5)     

     env.workspace = out
     finras5 = out + "\\" + "Normal_Week_5" + "_1981_2016" + "_" + str(month) + ".tif"
     calc5 = arcpy.sa.CellStatistics(new_list5, statistics_type = "MEAN")
     arcpy.CopyRaster_management(calc5,finras5,"","","","","","32_BIT_FLOAT")‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Re: Error in multiple raster processing

DanPatterson_Retired — Thu, 14 Sep 2017 14:55:42 GMT

indentation is wrong in line 17

lines 18 - end should be dedented to line up with the for in line 15

throw in some print statements

try using

arcpy.Exists( some file )

to see if your files actually exist and/or can be found

A continuation of this thread which you might want to close

Re: Error in multiple raster processing

RandyBurton — Sat, 11 Dec 2021 12:50:17 GMT

Hi Dumindu,

I agree with Dan Patterson‌ that you have some indentation issues with your code.

I may not understand exactly what you are trying to accomplish, but it looks like you want to average the raster files for the first week (week 1) of September (month 9) for the years 1981 to 2016. Then move to the next month and eventually loop through all the weeks. This is suggested by the file you are saving to found in line 19 of your code.

finras1 = out + "\\" + "Normal_Week_" + str(week) + "_1981_2016" + "_" + str(month) + ".tif"
‍

If so, you may have your loops in the wrong order. You may want something like this:

for week in range(1, 5+1): # for weeks 1 to 5
    for month in range(9, 11+1): # for months 9 to 11
        for year in range(1981, 2016+1): # for years 1981 to 2016
            print "Week_{}_Sum{}_{}.tif".format(week, year, month)
'''
Result:
Week_1_Sum1981_9.tif
Week_1_Sum1982_9.tif
...
Week_1_Sum2015_9.tif
Week_1_Sum2016_9.tif
Week_1_Sum1981_10.tif
Week_1_Sum1982_10.tif
...
Week_1_Sum2015_10.tif
Week_1_Sum2016_10.tif
Week_1_Sum1981_11.tif
Week_1_Sum1982_11.tif
...
Week_1_Sum2015_11.tif
Week_1_Sum2016_11.tif
Week_2_Sum1981_9.tif
Week_2_Sum1982_9.tif
...
...
Week_5_Sum2015_11.tif
Week_5_Sum2016_11.tif
'''‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

If this is the case, you may want to change the directory structure of your rasters, with the weeks folders containing the months folders which contain the years folders. This would simplify what you are trying to do with arcpy.env.workspace and arcpy.ListRasters.

It would probably be helpful if you could describe the directory structure you are working with.

Re: Error in multiple raster processing

DuminduJayasekera — Fri, 15 Sep 2017 16:00:38 GMT

Thanks Randy for the reply. I have revised my code and posted above. There is no error now and producing outputs but the output values are NOT correct. I am doing something wrong here.

I have folders names "1981", "1982"....."2016" in "H:\PRISM_800m_weekly_sum_Normal" and each folder (year) has Weekly .tif files (for example Week_1_Sum1981_9.tif, Week_2_Sum1981_9.tif, ...., Week_5_Sum1981_9.tif,.....Week_5_Sum1981_9.tif) for the months from 9 to 11.

Can you help me to fix this?

Thanks in advance.

Re: Error in multiple raster processing

RandyBurton — Sat, 11 Dec 2021 12:50:20 GMT

Does the following code print a listing of your files in the proper grouping for analysis? (I'm still trying to understand your directory structure and the objective.)

groupNum = 1
for week in range(1, 5+1): # for weeks 1 to 5
    print "Files in group {}:".format(groupNum)
    groupNum += 1
    for month in range(9, 11+1): # for months 9 to 11
        if month > 9:
            print "Files in group {}:".format(groupNum)
            groupNum += 1
        for year in range(1981, 2016+1): # for years 1981 to 2016
            print "\tH:\\PRISM_800m_weekly_sum_Normal\\{}\\Week_{}\\Week_{}_Sum{}_{}.tif".format(year, week, week, year, month)
‍‍‍‍‍‍‍‍‍‍

If the listing is correct, then we can modify the code to do the analysis.

Re: Error in multiple raster processing

DuminduJayasekera — Fri, 15 Sep 2017 19:06:44 GMT

Randy, Thanks again. That is the output I want exactly but I would revise the last line like below. But, how can i modify the above code?

print "\tH:\\PRISM_800m_weekly_sum_Normal\\{}\\Week_{}_Sum{}_{}.tif".format(year, week, year, month)

Re: Error in multiple raster processing

RandyBurton — Sat, 11 Dec 2021 12:50:22 GMT

Glad that the output is showing the grouping that you want with your modification. I have added back in some of your earlier code, but I may have missed something.

import arcpy

# make sure this directory exists
out = r'H:\PRISM_800m_weekly_sum\Normal_out'

for week in range(1, 5+1): # for weeks 1 to 5

    for month in range(9, 11+1): # for months 9 to 11

        rasterlist = [] # initialize a list for the rasters
        
        for year in range(1981, 2016+1): # for years 1981 to 2016

            rasterfile = "H:\\PRISM_800m_weekly_sum_Normal\\{}\\Week_{}_Sum{}_{}.tif".format(year, week, year, month)
                
            # probably should add some error checking to see that file exists before adding it to rasterlist
            if arcpy.Exists(rasterfile):
                rasterlist.append(rasterfile) # if found, append file to the raster list
            else:
                print "Not found: {}".format(rasterfile)

        # same indentation as "for year in..."
        # rasterlist.sort() # you probably don't need to sort the list

        finras1 = out + "\\Normal_Week_{}_1981_2016_{}.tif".format(week, month)
        print "Processing file: {}".format(finras1)

        calc1 = arcpy.sa.CellStatistics(rasterlist, statistics_type = "MEAN", ignore_nodata="DATA")
        arcpy.CopyRaster_management(calc1,finras1,"","","","","","32_BIT_FLOAT")

print "Done."
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Also, keep in mind what Dan Patterson‌ said in his comments in another of your questions regarding working with rasters. Remember that indentation is important, and it is often good to have some print statements in your code to help you know where things are going right or wrong.

Hope this helps.

Re: Error in multiple raster processing

DuminduJayasekera — Fri, 15 Sep 2017 21:40:02 GMT

Thanks a lot. this helps.

Re: Error in multiple raster processing

RandyBurton — Fri, 15 Sep 2017 21:42:32 GMT

Glad it is working out.