Model Builder/Python: batch raster calculation

[TomHawkins]

I wish to carry out a raster calculation on a batch of rasters.
I want to add rasters by month, as identified by their file name.

I.e. all rasters within the month of January should add together and output a single jan.png/raster file. Then all the rasters within Feb should add together and produce a single FEB.png.

Rather than manually write a new rast calc expression for each month, I would prefer to automate it via either python or modelbuilder.

My plan outline is as follows:

Input Files. 36 Weekly raster files

5jan15-11jan15.png

12jan15-18jan15.png

19jan15-25jan15.png

26jan15-31jan15.png

01feb15-07feb15.png....

Raster operation.

(raster1+raster2+raster3+raster4)/(raster1+raster2+raster3+raster4) = month.png

output. 12 rasters for each month.

jan.png

feb.png

march.png

is it possible to write a python expression or create a modelbuilder which will add rasters by month using their filename?

Current Attempts

I have built a model, however this only runs on the rasters I provide it, requiring me to write a new expression each time. In other words it's the same as manually running raster calculator and writing a new expression for each month, I wish to automate this via modelbuilder.

[DanPatterson_Retired]

you will have to read up on iterators

A quick tour of using iterators—Help | ArcGIS for Desktop

Examples of using iterators in ModelBuilder—Help | ArcGIS for Desktop

[TomHawkins]

Hi Dan, I've been reading up on iterators and trying to wrap my head around them.

In the examples I have seen it seems it will cause the raster Calc to run on EACH of the rasters. Whereas I only want it to run on GROUPS of the raster, which are grouped by month.

i.e. I want the following...

(janweek1+janweek2+janweek3+janweek4) /((janweek1+janweek2+janweek3+janweek4))= jan.png

(febweek1+febweek2+febweek3+febweek4)/(janweek1+janweek2+janweek3+janweek4) = feb.png

yet from my understanding the iterator will produce something akin to

janweek1+janweek2+janweek3+janweek4 + febweek1+febweek2+febweek3+febweek4....=month.png

[DarrenWiens2]

You can make a new table to hold the month names (row values = jan, feb, etc.). Then Iterate through those with Iterate Field Values. Finally, use the outputs of that simple submodel to feed into the wildcard parameter of Iterate Rasters.

Unfortunately, using multiple iterators is a relatively complex way to do something very straight forward in Python. Something like:

This is a list of month names.
For each month name, loop through rasters.
If the raster name contains the month name, add it and save it.

[TomHawkins]

I've read that this is something that would be easier in Python. If you could provide some advice as to a python-like solution that would be brilliant and perhaps less mind boggling.

I would like Python to run the calculator expression by months.

[DarrenWiens2]

Untested - there may be errors:

import calendar, os # import libraries
arcpy.env.workspace = r'C:\myfolder\anotherfolder' # your workspace
for month in calendar.month_abbr: # loop through 3-letter month abbreviations
     rasters = arcpy.ListRasters('*' + month.lower() + '*') # make list of matching rasters in workspace
     monthsum = sum(rasters) # add matching rasters together
     monthsum.save(os.path.join(r'C:\myfolder\anotherfolder',month.lower()+'sum.png' ) # save to new raster

[TomHawkins]

Thanks Darren.
Does line 5 represent the raster calculator expression?

Unfortunately my actual expression is a little more complicated than provided...

Float(Con(IsNull("w1") ,0,"w1") + Con(IsNull("w2"), 0,"w2")+Con(IsNull("w3"), 0,"w3")+Con(IsNull("w4"), 0,"w4") )/(Con(IsNull("w1"), 0,1) + Con(IsNull("w2"), 0,1)+Con(IsNull("w3"), 0,1)+Con(IsNull("w4"), 0,1))

That expression is only valid if there is 4 rasters for that month of course, there could only be three etc.

A crux of problem is trying to tell Arcmap that the variables 'week1, week2, week3, week4` are dynamic, meaning that the input rasters change for each month.

In other words I'm trying to create a "dynamic expression"

The basic expression is

Float(Con(IsNull("w1") ,0,"w1") + Con(IsNull("w2"), 0,"w2") /(Con(IsNull("w1"), 0,1) + Con(IsNull("w2"), 0,1)

I would like this expression to adapt to each month, automatically inputting the rasters of that group and returning a `month.png/raster`

[DarrenWiens2]

Here's what I would try. Again, untested (sorry, I'm too lazy to duplicate your raster folder setup). I'm not 100% sure that you can sum rasters using sum(), but it's worth a shot.

import calendar, os # import libraries  
arcpy.env.workspace = r'C:\myfolder\anotherfolder' # your workspace  
for month in calendar.month_abbr: # loop through 3-letter month abbreviations  
    rasters = arcpy.ListRasters('*' + month.lower() + '*') # make list of matching rasters in workspace  
    con_rasters = [] # empty list
    for raster in rasters: # loop through matching rasters
        con_rasters.append(Con(IsNull(raster),0,raster)) # run Con expression and add to list

    monthsum = sum(con_rasters) # add Con rasters
    monthsum.save(os.path.join(r'C:\myfolder\anotherfolder',month.lower()+'sum.png' ) # save to new raster 

[TomHawkins]

Sorry, but I'm a little bit confused, this is most likely due to my naivety in reading Python syntax.

When looking at lines 05-07. of code, I see

1. con_rasters = [] # empty list 
2.     for raster in rasters: # loop through matching rasters 
3.         con_rasters.append(Con(IsNull(raster),0,raster)) # run Con expression and add to list 

I feel that it is only running the first half of the expression, which is...

          Float(Con(IsNull("w1") ,0,"w1") + Con(IsNull("w2"), 0,"w2")+Con(IsNull("w3"), 0,"w3")+Con(IsNull("w4"), 0,"w4")

and leaving out the second half of the expression which is....

         / (Con(IsNull("w1"), 0,1) + Con(IsNull("w2"), 0,1)+Con(IsNull("w3"), 0,1)+Con(IsNull("w4"), 0,1))

The expression is essentially the number of times a cell had a value and divide by the number of times it was possible for a value to be detected. The numerator in this case is how many times the pixel had the value 1. The denominator is now many times the pixel did not have No Data (i.e. it had the value 1 or 0).

Expression:

number of times a cell had a value / number of times the cell could have had a value

Float(Con(IsNull(A),0,A) + Con(IsNull(B),0,B) + Con(IsNull(C),0,C)) / (~IsNull(A) + ~IsNull(B) + ~IsNull(C))

I hope that clarifies what I'm trying to achieve and apologies if I've misunderstood your code!

[DarrenWiens2]

I need to hop off here for a while - just wanted you to know that yes, I completely missed the division part.