Looping through Raster to compute Mean in Cell Statistics

[KyleHartfield]

I am attempting to loop through a set of MODIS data to compute averages by period.  Since 2000 is missing periods 1 - 3 I am just computing the averages of periods 1 -3 for 2001 through 2012.  I am then computing the averages for periods 4 - 23 with data from 2000 through 2012.  I have written a script to do this, but keep getting an error 000732 about my data not existing or not being supported.  I opened one of the grids in Arc and it was fine.  I am not sure why I am getting this error and have included the code below.

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

# Set environment settings
env.workspace = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\"

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Path names for input and output data
path = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or"
output_path = "B:\\MODIS_NDVI\\Avg_per_period\\avg"

for period in range (1,3):
#MODIS does not have data for the first three periods of 2000.  This will calculate an average for the first three periods ignoring the year 2000.  It will also place a 0 in front of the single digit to match the file name.
    if period <= 3:
        period = ("0"+str(period))
        grids01 = Raster('"' + path + "01" + period + "ndvius" + '"')
        grids02 = Raster('"' + path + "02" + period + "ndvius" + '"')
        grids03 = Raster('"' + path + "03" + period + "ndvius" + '"')
        grids04 = Raster('"' + path + "04" + period + "ndvius" + '"')
        grids05 = Raster('"' + path + "05" + period + "ndvius" + '"')
        grids06 = Raster('"' + path + "06" + period + "ndvius" + '"')
        grids07 = Raster('"' + path + "07" + period + "ndvius" + '"')
        grids08 = Raster('"' + path + "08" + period + "ndvius" + '"')
        grids09 = Raster('"' + path + "09" + period + "ndvius" + '"')
        grids10 = Raster('"' + path + "10" + period + "ndvius" + '"')
        grids11 = Raster('"' + path + "11" + period + "ndvius" + '"')
        grids12 = Raster('"' + path + "12" + period + "ndvius" + '"')
        #grids = grids01 + ";" + grids02 + ";" + grids03 + ";" + grids04 + ";" + grids05 + ";" + grids06 + ";" + grids07 + ";" + grids08 + ";" + grids09 + ";" + grids10 + ";" + grids11 + ";" + grids12   
        avg_temp = ('"' + output_path + period + "temp" + '"')
        avg = (output_path + period)
        outCellStatistics = CellStatistics([grids01, grids02, grids03, grids04, grids05, grids06, grids07, grids08, grids09, grids10, grids11, grids12], "MEAN", "DATA")
        outCellStatistics.save(avg_temp)
#This is will place a 0 in front of the single digit to match the file name.
    elif period < 10:
        period = ("0"+str(period))
        grids00 = ('"' + path + "00" + period + "ndvius" + '"')
        grids01 = ('"' + path + "01" + period + "ndvius" + '"')
        grids02 = ('"' + path + "02" + period + "ndvius" + '"')
        grids03 = ('"' + path + "03" + period + "ndvius" + '"')
        grids04 = ('"' + path + "04" + period + "ndvius" + '"')
        grids05 = ('"' + path + "05" + period + "ndvius" + '"')
        grids06 = ('"' + path + "06" + period + "ndvius" + '"')
        grids07 = ('"' + path + "07" + period + "ndvius" + '"')
        grids08 = ('"' + path + "08" + period + "ndvius" + '"')
        grids09 = ('"' + path + "09" + period + "ndvius" + '"')
        grids10 = ('"' + path + "10" + period + "ndvius" + '"')
        grids11 = ('"' + path + "11" + period + "ndvius" + '"')
        grids12 = ('"' + path + "12" + period + "ndvius" + '"')
        grids = "[" + grids00 + "," + grids01 + "," + grids02 + "," + grids03 + "," + grids04 + "," + grids05 + "," + grids06 + "," + grids07 + "," + grids08 + "," + grids09 + "," + grids10 + "," + grids11 + "," + grids12 + "]"    
        avg_temp = (output_path + period + "temp")
        avg = (output_path + period)
        arcpy.gp.CellStatistics_sa(grids,avg_temp,"MEAN", "DATA")
#These are double digit periods and will match the file name without correction.
    else:
        period = (str(period))
        grids00 = ('"' + path + "00" + period + "ndvius" + '"')
        grids01 = ('"' + path + "01" + period + "ndvius" + '"')
        grids02 = ('"' + path + "02" + period + "ndvius" + '"')
        grids03 = ('"' + path + "03" + period + "ndvius" + '"')
        grids04 = ('"' + path + "04" + period + "ndvius" + '"')
        grids05 = ('"' + path + "05" + period + "ndvius" + '"')
        grids06 = ('"' + path + "06" + period + "ndvius" + '"')
        grids07 = ('"' + path + "07" + period + "ndvius" + '"')
        grids08 = ('"' + path + "08" + period + "ndvius" + '"')
        grids09 = ('"' + path + "09" + period + "ndvius" + '"')
        grids10 = ('"' + path + "10" + period + "ndvius" + '"')
        grids11 = ('"' + path + "11" + period + "ndvius" + '"')
        grids12 = ('"' + path + "12" + period + "ndvius" + '"')
        grids = "[" + grids00 + "," + grids01 + "," + grids02 + "," + grids03 + "," + grids04 + "," + grids05 + "," + grids06 + "," + grids07 + "," + grids08 + "," + grids09 + "," + grids10 + "," + grids11 + "," + grids12 + "]" 
        avg_temp = (output_path + period + "temp")
        avg = (output_path + period)
        arcpy.gp.CellStatistics_sa(grids,avg_temp,"MEAN", "DATA")

[RhettZufelt]

Do you actually have a "raster" data set named S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or0101ndvius ?

If so, what format is it, and you might try to remove the "'" from both ends of your grids## = assignments.  If I type them in as you have them, I get " ' S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or0101ndvius ' " and it might be the extra quotes causing the issue.  I don't have the data, so can't really test.

If you don't have a raster with that name, that is what it is looking for so you might check to make sure the resultant path actually matches an existing raster.

Running under the 64 bit IDE can give this error as well.  Normally it is when you are trying to access something through an ArcCatalog connection though (since that only supports 32 bit, the 64 bit one can't find the ArcCatalog/connectionto/ files).

R_

[KyleHartfield]

Do you actually have a "raster" data set named S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or0101ndvius ?

If so, what format is it, and you might try to remove the "'" from both ends of your grids## = assignments.  If I type them in as you have them, I get " ' S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or0101ndvius ' " and it might be the extra quotes causing the issue.  I don't have the data, so can't really test.

If you don't have a raster with that name, that is what it is looking for so you might check to make sure the resultant path actually matches an existing raster.

Running under the 64 bit IDE can give this error as well.  Normally it is when you are trying to access something through an ArcCatalog connection though (since that only supports 32 bit, the 64 bit one can't find the ArcCatalog/connectionto/ files).

R_

Rzufelt,

That is the name of my raster.  I have attempted various formats of naming and was trying to copy how the command looks when I export it from model builder.  I removed the additional quotes and received the same error.  The machine I am running is a 64 bit machine, but I am not accessing the data through ArcCatalog.  Is there a way to work around the 64 bit problem.  I have included the piece of my adjusted script.

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

# Set environment settings
env.workspace = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\"

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Path names for input and output data
path = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or"
output_path = "B:\\MODIS_NDVI\\Avg_per_period\\avg"

for period in range (1,3):
#MODIS does not have data for the first three periods of 2000.  This will calculate an average for the first three periods ignoring the year 2000.  It will also place a 0 in front of the single digit to match the file name.
    if period <= 3:
        period = ("0"+str(period))
        grids01 = Raster(path + "01" + period + "ndvius")
        grids02 = Raster(path + "02" + period + "ndvius")
        grids03 = Raster(path + "03" + period + "ndvius")
        grids04 = Raster(path + "04" + period + "ndvius")
        grids05 = Raster(path + "05" + period + "ndvius")
        grids06 = Raster(path + "06" + period + "ndvius")
        grids07 = Raster(path + "07" + period + "ndvius")
        grids08 = Raster(path + "08" + period + "ndvius")
        grids09 = Raster(path + "09" + period + "ndvius")
        grids10 = Raster(path + "10" + period + "ndvius")
        grids11 = Raster(path + "11" + period + "ndvius")
        grids12 = Raster(path + "12" + period + "ndvius")
        #grids = grids01 + ";" + grids02 + ";" + grids03 + ";" + grids04 + ";" + grids05 + ";" + grids06 + ";" + grids07 + ";" + grids08 + ";" + grids09 + ";" + grids10 + ";" + grids11 + ";" + grids12   
        avg_temp = ('"' + output_path + period + "temp" + '"')
        avg = (output_path + period)
        outCellStatistics = CellStatistics([grids01, grids02, grids03, grids04, grids05, grids06, grids07, grids08, grids09, grids10, grids11, grids12], "MEAN", "DATA")
        outCellStatistics.save(avg_temp)

[curtvprice]

   The machine I am running is a 64 bit machine, but I am not accessing the data through ArcCatalog.  Is there a way to work around the 64 bit problem.  I have included the piece of my adjusted script.

I don't think 64 bits is the issue.

Let me try an approach that may work for you better. Both in formatting your pathnames, and using a list of arcpy raster objects instead of putting monster strings together:

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

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Path names for input and output data
path = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or"
output_path = "B:\\MODIS_NDVI\\Avg_per_period"
env.workspace = env.scratchWorkspace = output_path

rasList = []
for period in range (1,3):  # [1,2]
    for m in range(1,13): # 1..12
        rasPath = "{p}{month}{per:02d}ndvius".format(
            p=path, month=m, period)
        if arcpy.Exists(rasPath):
            rasList.append(Raster(rasPath))
        else:
            raise Exception("Raster {0} not found".format(rasPath))
    outCellStatistics = CellStatistics(rasList, "MEAN", "DATA")
    outName =  "avg{0:02d}".format(period)
    outCellStatistics.save(outName)

[RhettZufelt]

Here is how I access the 32 bit (being that I installed ArcGIS desktop (32 bit) and AGS Server (64 bit) on same computer.

If I am running a scheduled script or one spawned from another, I use this code in a batch file:

cd C:\Python27\ArcGIS10.1

Python C:\arcgisserver\python\buildings.py

Exit

However, to run in a 32 bit IDE, I made a shortcut with the follwing info [ATTACH=CONFIG]24732[/ATTACH], double click on it, and opens IDE in 32 bit (as you can see in the info at the top of the IDE).
I have gotten very few scripts to actually run correctly in the 64 bit window.

Also, I am finding many bugs introduced in the 10.1 arcpy so, if you have the ability, you might try your code in the 10.0 version and see if you've found another bug....

R_

Also, not sure in this case, but a lot of ArcGIS tools use the raster extention to determine the file type.  without an extention (as in your data), and not in a FGDB, most ArcTools will assume that is a GRID dataset.  If not, perhaps this is causing an issue?

Just a thought

[curtvprice]

  I have gotten very few scripts to actually run correctly in the 64 bit window.

There are unsupported tools and data types that must be kept in mind when developing for 64 bit arcpy.

Desktop 10.1 help: Background Geoprocessing (64-bit)

[KyleHartfield]

Thanks to all those who chimed in to help.  I have attached the code I finally got to work.  At this point in my python education this code makes the most sense in my head.  I understand that there is a more efficient way to do this and I am open to learning how if you can point me in the right direction.  I explored the .format tool and could not get it to loop the way I wanted to. Thanks again for all the help.

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

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Path names for input and output data
path = "S:\\MODIS\\DATA_v05\\modis016_0250_ndvi\\us\\or\\or"
output_path = "B:\\MODIS_NDVI\\Avg_per_period\\avg"

for period in range (1,24):
#MODIS does not have data for the first three periods of 2000.  This will calculate an average for the first three periods ignoring the year 2000.  It will also place a 0 in front of the single digit to match the file name.
    if period <= 3:
        period = ("0"+str(period))
        grids01 = Raster(path + "01" + period + "ndvius")
        grids02 = Raster(path + "02" + period + "ndvius")
        grids03 = Raster(path + "03" + period + "ndvius")
        grids04 = Raster(path + "04" + period + "ndvius")
        grids05 = Raster(path + "05" + period + "ndvius")
        grids06 = Raster(path + "06" + period + "ndvius")
        grids07 = Raster(path + "07" + period + "ndvius")
        grids08 = Raster(path + "08" + period + "ndvius")
        grids09 = Raster(path + "09" + period + "ndvius")
        grids10 = Raster(path + "10" + period + "ndvius")
        grids11 = Raster(path + "11" + period + "ndvius")
        grids12 = Raster(path + "12" + period + "ndvius")
        grids = [grids01, grids02, grids03, grids04, grids05, grids06, grids07, grids08, grids09, grids10, grids11, grids12]   
        avg_temp = (output_path + period + "temp")
        avg = (output_path + period)
        arcpy.gp.CellStatistics_sa(grids, avg_temp,"MEAN", "DATA")
        print ("Cell Stats complete for" + period)
        arcpy.gp.Int_sa(avg_temp, avg)
        print ("Int complete for" + period)
#This is will place a 0 in front of the single digit to match the file name.
    elif period < 10:
        period = ("0"+str(period))
        grids00 = Raster(path + "00" + period + "ndvius")
        grids01 = Raster(path + "01" + period + "ndvius")
        grids02 = Raster(path + "02" + period + "ndvius")
        grids03 = Raster(path + "03" + period + "ndvius")
        grids04 = Raster(path + "04" + period + "ndvius")
        grids05 = Raster(path + "05" + period + "ndvius")
        grids06 = Raster(path + "06" + period + "ndvius")
        grids07 = Raster(path + "07" + period + "ndvius")
        grids08 = Raster(path + "08" + period + "ndvius")
        grids09 = Raster(path + "09" + period + "ndvius")
        grids10 = Raster(path + "10" + period + "ndvius")
        grids11 = Raster(path + "11" + period + "ndvius")
        grids12 = Raster(path + "12" + period + "ndvius")
        grids = [grids00, grids01, grids02, grids03, grids04, grids05, grids06, grids07, grids08, grids09, grids10, grids11, grids12]   
        avg_temp = (output_path + period + "temp")
        avg = (output_path + period)
        arcpy.gp.CellStatistics_sa(grids, avg_temp,"MEAN", "DATA")
        print ("Cell Stats complete for" + period)
        arcpy.gp.Int_sa(avg_temp, avg)
        print ("Int complete for" + period)
#These are double digit periods and will match the file name without correction.
    else:
        period = (str(period))
        grids00 = Raster(path + "00" + period + "ndvius")
        grids01 = Raster(path + "01" + period + "ndvius")
        grids02 = Raster(path + "02" + period + "ndvius")
        grids03 = Raster(path + "03" + period + "ndvius")
        grids04 = Raster(path + "04" + period + "ndvius")
        grids05 = Raster(path + "05" + period + "ndvius")
        grids06 = Raster(path + "06" + period + "ndvius")
        grids07 = Raster(path + "07" + period + "ndvius")
        grids08 = Raster(path + "08" + period + "ndvius")
        grids09 = Raster(path + "09" + period + "ndvius")
        grids10 = Raster(path + "10" + period + "ndvius")
        grids11 = Raster(path + "11" + period + "ndvius")
        grids12 = Raster(path + "12" + period + "ndvius")
        grids = [grids00, grids01, grids02, grids03, grids04, grids05, grids06, grids07, grids08, grids09, grids10, grids11, grids12]   
        avg_temp = (output_path + period + "temp")
        avg = (output_path + period)
        arcpy.gp.CellStatistics_sa(grids, avg_temp,"MEAN", "DATA")
        print ("Cell Stats complete for" + period)
        arcpy.gp.Int_sa(avg_temp, avg)
        print ("Int complete for" + period)