Perform Raster Calculation from Multiple Sub-folder Using ArcPy

[ShouvikJha]

I have five subfolders with Rasters name with the same parameter. I am trying to figure out a way to Loop average of all  raster and save output raster in outfolder . For instance 

Subfolder1   Subfolder2   Subfolder3   Subfolder4    Subfolder5   Outfolder
r001_NPP.tif r001_npp.tif r001_NPP.tif r_001_NPP.tif r001_NPP.tif r001_MEAN_NPP.tif‍‍‍‍‍‍

‍‍‍‍‍‍‍

I have tried to write code for this, Below i post my code what i have so far

import arcpy, os, glob
arcpy.CheckOutExtension("Spatial")

ws1 = glob.glob(r'D:\subfolder1\*.tif')
ws2 = glob.glob(r'D:\subfolder2\*.tif')
ws3 = glob.glob(r'D:\subfolder3\*.tif')
ws4 = glob.glob(r'D:\subfolder4\*.tif')
ws5 = glob.glob(r'D:\subfolder5\*.tif')
outws = r'D:\outfolder'

for r in ws1, ws2, ws3, ws4:
    basename = os.path.basename(r).split("_")[0]
    r1 = arcpy.sa.Raster(os.path.join(ws1, basename + "r_NPP.tif"))
    r2 = arcpy.sa.Raster(os.path.join(ws2, basename + "r_NPP.tif"))
    r3 = arcpy.sa.Raster(os.path.join(ws3, basename + "r_NPP.tif"))
    r4 = arcpy.sa.Raster(os.path.join(ws4, basename + "r_NPP.tif"))
    r5 = arcpy.sa.Raster(os.path.join(ws5, basename + "r_NPP.tif"))
    

    result = (r1 + r2 + r3 + r4 + r5) / 5
    outname = basename + "r_Mean_NPP.tif"

    result.save(os.path.join(outws, outname))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[ShouvikJha]

xander_bakker‌, We need your cooperation to solve above problem. Kindly cooperate  with us. 

[DanPatterson_Retired]

Here is how you test things... when you get an error, narrow down the code and test it, like in the following example to ensure that you output file names are correct..

>>> import os
>>> out_workplace = 'D:/Test/MEAN_NPP'  # note I am using os independent path specs
>>> num = 1   # obviously if you had this in a loop 'num' would be changing
>>> p = 'r{0}_MEAN_NPP.TIF'.format(num)   # this is the path I am trying to make
>>> x = os.path.join(out_workplace, p)  # joining them together
>>> x
'D:/Test/MEAN_NPP/r1_MEAN_NPP.TIF'      # this is what we get
‍‍‍‍‍‍‍

Now... is that right? 

[XanderBakker]

I'm looking at the code and trying to figure what you want to do exactly, since the indentation is not correct.

First of all I assume you want to perform this calculation for each month (001, .. , 012). To do so, you will need increase the range to (1, 13) (see line 9 below).

I also assume that you do not want to have the same output results for each month and that the list of rasters should be filtered by the month value. Therefore I included line 14.

Please be aware that in you original code you create the raster r001_MEAN_NPP.TIF and since that raster is stored in a subfolder of the input workspace, it will be included in the result for the second month. I don't think that is what you want. Also remember your closing brackets on line 19 (there were two missing in your code)

import arcpy
arcpy.CheckOutExtension("spatial")
from arcpy.sa import *

# generate list of rasters
workspace = r"D:\TEST"
out_workplace = "D:\TEST\MEAN_NPP"

for ras_num in [str(i).zfill(3) for i in range(1, 13)]:
    rasters = []
    walk = arcpy.da.Walk(workspace, datatype="RasterDataset", type="TIF")
    for dirpath, dirnames, filenames in walk:
        for filename in filenames:
            if ras_num in filename:
                rasters.append(os.path.join(dirpath, filename))

    # calculate mean
    ras_mean = CellStatistics(rasters, "MEAN")
    ras_mean.save(os.path.join(out_workspace, 'r{0}_MEAN_NPP.TIF'.format(ras_num)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍))‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[DanPatterson_Retired]

Xander, thanks.... line 7, perhaps...

now before people get sloppy and forget to use 'raw' or unix based folder formatting... lets get those little 'r' s in there, especially if people get sloppy and/or use lower case, or upper case gets translated to lower case.... sooo. a 'case' in point...

>>> w = 'D:\Test\Folder'
>>> print(w)
D:\Test\Folder
>>>
>>> # now let's get sloppy
>>>
>>> w = 'D:\test\folder'
>>> print(w)
D:     est 
older‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

[XanderBakker]

Absolutely Dan, you are right!

[ShouvikJha]

Thank you very much  all of you for your kind cooperation.

Xander Bakker  I ran the above code , but i am getting error. I think not having problem with code. Might be we need more explanation what i am trying to do, please follow the below 

Message     File Name     Line     Position     
Traceback                    
    <module>     <module1>     20          
    CellStatistics     C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py     2951          
    swapper     C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Utils.py     53          
    Wrapper     C:\Program Files\ArcGIS\Desktop10.3\ArcPy\arcpy\sa\Functions.py     2947          
RuntimeError: ERROR 999998: Unexpected Error.                    
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

I am giving here  what I am trying to do, I have three sub-folder (Name : NPP_2008, NPP_2009, NPP_2010),

NPP_2008      NPP_2009       NPP_2010       MEAN_NPP (output Folder)
r001_NPP.tif  r001_npp.tif   r001_NPP.tif   r001_MEAN_NPP.tif (mean of all r001_NPP)
r002_NPP.tif r002_npp.tif    r002_NPP.tif   r002_MEAN_NPP.tif (mean of all r002_NPP)
So on ..............

Every sub-folder hold monthly raster data for different year with the name of e.g r001_NPP, r002_NPP……..r012_NPP.  So I want calculate mean of those raster which is r001_NPP series in every sub-folder, similarly program will loop through till r012_NPP series of all raster to calculate mean and save it another subfolder (Name : Mean_NPP) with output name 001_MEAN_NPP(which is mean of all r001_NPP), 002_MEAN NPP…..012_MEAN_NPP.

Here with i have attached reference raster which i am using for my raster calculation

[XanderBakker · Accepted solution]

I can't reproduce the error you mentioned. I did encounter some other errors that should be corrected like:

• include import os at beginning of script and
• change "out_workplace" into "out_workspace"

import arcpy
import os
arcpy.CheckOutExtension("spatial")

# generate list of rasters
##workspace = r"D:\TEST"
##out_workplace = r"D:\TEST\MEAN_NPP"
workspace = r"C:\GeoNet\MeanRasters\Case20161009"
out_workspace = r"C:\GeoNet\MeanRasters\Case20161009\MEAN_NPP"

for ras_num in [str(i).zfill(3) for i in range(1, 13)]:
    print ras_num
    rasters = []
    walk = arcpy.da.Walk(workspace, datatype="RasterDataset", type="TIF")
    for dirpath, dirnames, filenames in walk:
        for filename in filenames:
            if ras_num in filename:
                rasters.append(os.path.join(dirpath, filename))

    # calculate mean
    if len(rasters) > 0:
        print "Calculating mean of rasters: ", ', '.join(rasters)
        ras_mean = arcpy.sa.CellStatistics(rasters, "MEAN")
        ras_mean.save(os.path.join(out_workspace, 'r{0}_MEAN_NPP.TIF'.format(ras_num)))
    else:
        print "no rasters for ras_num", ras_num‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

This created the mean raster for ras_num 001 to 005, see attached ZIP with resulting rasters

[ShouvikJha]

xander_bakker‌. Thank you very much. Above Code running perfectly.

Regarding  error issue which you got different error massage, because i corrected that line before running script.  

[XanderBakker]

I'm glad it's working now...

[ShouvikJha]

Thank you. I am not getting your points. i am trying plot two variable of raster like this ...