Python code for selecting and averaging raster layers within subfolders

[DeeMcElligott1]

Hi there,

Basically I have a folder with ten subfolders (one for each year of data...1992-2012), within each of these subfolders are 52 satellite images, one for each week of the corresponding year. I want to convert the weekly data into monthly data, thereby resulting in a single output for each month and twelve for each year. So I want to run a script that picks every four images consecutively, i.e tell it: the first four images = Jan, second four = Feb ............and so on, and get the average of these four images.


I can get the average by running the tool 'Cell Statistic' in Model builder and use 'Calculate Value' to allow me to enter the Python script to tell 'Cell Statistic' which four images to average.

I'm afraid my Python scripting capabilities are rather limited. Does anyone know a Python code to enable me to do this?

I hope it makes sense!

Thanks

Dee

[ThomasEmge]

Here is some pseudo code for 10.1:

arcpy.env.workspace = "folder for 1992"
listOfRastersByWeek = arcpy.ListRasters()

# repeat this inside a loop for all 12 months
week1 = arcpy.RasterToNumPyArray(listOfRastersByWeek[0])
week2 = arcpy.RasterToNumPyArray(listOfRastersByWeek[1])
week3 = arcpy.RasterToNumPyArray(listOfRastersByWeek[2])
week4 = arcpy.RasterToNumPyArray(listOfRastersByWeek[3])

month1 = (week1 + week2 + week3 + week4) / 4
raster_month1 = arcpy.NumPyArraryToRaster(month1)

raster_month1.save("unique name of the month")

[DeeMcElligott1]

Thanks very much for your reply Thomas.

Unfortunately it didn't work 😞 I'm using ArcMap 10, would this make a difference?

When I ran it, I got the following error:

Executing (Calculate Value): CalculateValue Raster "arcpy.env.workspace = "1998"\nlistOfRastersByWeek = arcpy.ListRasters()\n\n# repeat this inside a loop for all 12 months\nweek1 = arcpy.RasterToNumPyArray(listOfRastersByWeek[0])\nweek2 = arcpy.RasterToNumPyArray(listOfRastersByWeek[1])\nweek3 = arcpy.RasterToNumPyArray(listOfRastersByWeek[2])\nweek4 = arcpy.RasterToNumPyArray(listOfRastersByWeek[3])\n\nmonth1 = (week1 + week2 + week3 + week4) / 4\nraster_month1 = arcpy.NumPyArraryToRaster(month1)\n\nraster_month1.save("Jan")" Variant
Start Time: Thu Aug 09 21:16:46 2012

ERROR 000539: Runtime error <type 'exceptions.IndexError'>: list index out of range
Failed to execute (Calculate Value).

Am I making a blatant, Python amateur, mistake?

Cheers

Dee

[DanPatterson_Retired]

it doesn't look like you specified the full path in pythonic format for use by the arcpy.env variable

[DeeMcElligott1]

Thanks fro that Dan,

Still no joy though. I'm getting the same error:

Executing (Calculate Value): CalculateValue Raster "import arcpy\n\narcpy.env.workspace = "folder for 1992"\nlistOfRastersByWeek = arcpy.ListRasters()\n\n# repeat this inside a loop for all 12 months\nweek1 = arcpy.RasterToNumPyArray(listOfRastersByWeek[0])\nweek2 = arcpy.RasterToNumPyArray(listOfRastersByWeek[1])\nweek3 = arcpy.RasterToNumPyArray(listOfRastersByWeek[2])\nweek4 = arcpy.RasterToNumPyArray(listOfRastersByWeek[3])\n\nmonth1 = (week1 + week2 + week3 + week4) / 4\nraster_month1 = arcpy.NumPyArraryToRaster(month1)\n\nraster_month1.save("unique name of the month")" Variant
Start Time: Thu Aug 09 22:20:31 2012

ERROR 000539: Runtime error <type 'exceptions.IndexError'>: list index out of range
Failed to execute (Calculate Value).

Any other suggestions would be extremely welcome!!

Apologies for my obvious lack of Python programming!

[ThomasEmge]

for ArcGIS 10.0 the code should look somewhat like this

# we are using the python calendar module to help us out with the weeks
import calendar 

# define your current workspace as the arcpy listing function use this information
arcpy.env.workspace = r'd:\archive\1992'

# let's generate a dictionary with the names for the months
month_dict = dict((v,k) for k,v in enumerate(calendar.month_name))

# within the workspace let's get all the rasters
# the assumption here is that the naming of the individual rasters actually orders 
# them by weeks starting at the beginning of the year all the way to the end
# (this might or might not be true in your case)
listOfRastersByWeek = arcpy.ListRasters()

# now this getting a little funky and might be too much
# how do you distribute the 52 weeks over the year, or how many weeks are in each month?
year = []
weekCounter = 0
for index in range(1,13):
    # now month contains an arrangement of calendar days in weekly chunks
    month = calendar.monthcalendar(1992, index)
    # this array tells us how many days in a week for the specified month
    daysinweek = [len[1 for day in weeks if day] for weeks in month]
    # if a weekly chuck has more than 4 days we'll keep it and count it as a full week
    weeksInMonth = sum(1 for daycount in daysinweek if daycount > 3)
    # since the assumption of 4 days is somewhat random, let's doublecheck the number
    # of weeks
    weekCounter = weekCounter + weeksInMonth
    if (weekCounter > 52):
        year.append(52 - (weekCounter - weeksInMonth))
    else:
        year.append(weeksInMonth)    

# at this point the year array has the number of weeks per month
# with the understanding of how many weeks per month let's compute the average

startIndex = 0
for index in range(12):
    # initialize the monthly raster with the first week
    monthRaster = arcpy.Raster(listOfRastersByWeek[startIndex])
    for rasterIndex in range(1,year[index]):
        if startIndex < 52:
            # add each week into the monthly raster
            monthRaster = monthRaster + arcpy.Raster(listOfRastersByWeek[startIndex])
            # push the index to the next week
            startIndex = startIndex + 1

    # average the added weekly values
    monthRaster = monthRaster / year[index]
    # save the raster with the name of the month
    monthRaster.save(month_dict[index + 1])