4

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.

Usually I would use raster calculator, 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. enter image description here

2

You could do this in ModelBuilder, but it would tricky and probably require nesting models to get the iteration by month to work properly. With python, it's pretty straightforward. Here's some code to get you started:

import os, collections, re
import arcpy

arcpy.CheckOutExtension("Spatial")
arcpy.env.workspace = "<path to rasters directory>"
outdir = "<path to output dir"  

months = ["jan", "feb", "mar", "apr", "may", "jun",
          "jul", "aug", "sep", "oct", "nov", "dec"]
pat = re.compile("|".join(months))
ddict = collections.defaultdict(list)

# Create groupings of the files by month.
for raster in arcpy.ListRasters():
    res = re.findall(pat, raster)
    if res != []:
        ddict[res[0]].append(raster)    

for month,items in ddict.items():
    output = os.path.join(outdir, month + ".png")
    # Create raster objects to perform map algebra
    rasterObjs = [arcpy.Raster(it) for it in items]
    summed = sum(rasterObjs)
    rasOut = summed/summed
    # Save to disk
    arcpy.CopyRaster_management(rasOut, output)        

    # Cleanup
    try:
        arcpy.Delete_management(summed)
        arcpy.Delete_management(rasOut)
        arcpy.Delete_management(rasterObjs)
    except:
        pass
  • Hi Paul. Very helpful so far. However, I will need to maintain the raster expression as it is, rather than using the sum function. How would this be possible? Especially as some months have more weekly rasters than others. I imagine I will have a list for each month, with each containing 4-5 items/weekly rasters. I will need some code to form the correct expression regardless of number of items. I.e. if the JAN list has 4 items, then (raster1+raster2+raster3+raster4)/(raster1+raster2+raster3+raster4), but if the DEC list has 2 items, then: (raster1+raster2)/(raster1+raster2) – G. Gip Aug 8 '16 at 11:17
  • @G.Gip, if you can't use sum, then you could use a for loop and +=, or a list comprehension to build a string that is passed to eval() (evil!). I'd prefer the former, though it'll create some temporary outputs. – Paul Aug 8 '16 at 19:08
0

I don't know about a model, but I have written some code that may give you some insight into working towards your own script. I do not have the Raster Calculator license to test it.

The idea here bing you want to iterate through the filenames (strings) of files inside your directory path. Then you want to grab the month of any file that has the nomenclature of ##(3 letter month)##. If that string equals the specified month then you add that filename to a list. I created a function to make the calculations look a little better and again I haven't tested it, but it should work something like that. run the created function to basically concatenate all the strings you stored in each array you have created per month. There is more than likely an easier way to do this, but I am not adept at python. The only thing I can think of that may become an issue is that when you enter in the parameters for the raster calculator it may only take layers vs actual .png files.

`

import arcpy
import os
mJan = []
# have to create one for each month mFeb, mMar
for filename in os.listdir(r'DirectoryPath'):
    month = filename[-9:-6]
    if month == "jan":
        mJan.append(filename)
    #continue elif statesments for each month
    #elif month == "feb":
        #mFeb.append(filename)


def rasterCalc(month, array):
    arcpy.gp.RasterCalculator_sa(" & ".join(array),r'DirectoryPath'+ '\\' + month + '.png')

rasterCalc("jan", mJan)
#continue placing statments`

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.