calculate mean of selected raster chosen from multiple raster using ArcPy

ShouvikJha · ‎08-24-2016

I am trying to calculate mean of selected raster chosen from multiple raster. I have multiple raster file in different month folder, i am trying to select only two raster (e.g For JANUARY : 001_Max_Temper.tif, 001_Min_Temper.tif) from different month folder (JANUARY, FEBRUARY....DECEMBER) and calculate the mean of those selected raster and save it as on same month folder (e.g output name, JANUARY: 001_Mean_Temp).

I have written a code to do this task but i am getting error massage while i am running this code. Below i have attached my code.

import arcpy, os, calendar
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
arcpy.env.parallelProcessingFactor = "100%"
topWorkspace = r'D:\SWAT-WEATHER-DATA2'
arcpy.env.workspace = topWorkspace

# Get dict of months and month-number (i.e. January = 001, March = 003 etc.)
months = {calendar.month_name[i].upper(): str(i).zfill(3) for i in range(1, 13)} # Get dict of months and month number (i.e. January = 001, March = 003 etc.)

# Step through list of all folders
for folderPath in arcpy.ListWorkspaces():
    baseName = os.path.basename(folderPath).upper()
    if baseName in months: # Test that subfolder is a month name
        monthNumber = months[baseName] # Get month-number for use in output filename
        arcpy.env.workspace = folderPath
        rasterList = arcpy.ListRasters('*Max_Temper.tif', '*Min_Temper.tif')
        # Execute CellStatistics
        outCellStatistics = CellStatistics(rasterList, "MEAN", "NODATA")
        #Save the output
        outRasterName = outCellStatistics, "Mean_Temp_{}.tif".format(monthNumber)
        outCellStatistics.save(outRasterName)
        

#print done
print 'done'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

ShouvikJha · ‎08-25-2016

Xander Bakker‌. The expected output for month of October is 1 and we got the same. Through this calculation i am going to do research on monthly temperature effect on vegetation productivity, In India month of October is peak period for vegetation growth therefore October month temperature considering as a optimum temperature (ras_oct) for vegetation growth and its keep constant for all month's calculation.

If Temperature more or less from Optimum Temperature than productivity of vegetation also less therefore we kept constant the October month temperature.

Month of October the final output of TScalar is 1 and except October month TScalar vary between 0 to 1 because temperature can't meet with optimum temperature .

I am very thankful to you for your guidance.

XanderBakker · ‎08-25-2016

Thanks for the feedback! Interesting!

ShouvikJha · ‎08-26-2016

Xander Bakker‌. I am taking one more step to do calculation from all those raster from different subfolder. I am trying to execute this code for the next step calculation but i am getting error, below is code. Please cooperate me. The error code also attached

Message File Name Line Position
Traceback
 <module> <module2> 53
 main <module2> 32
UnboundLocalError: local variable 'ras_APAR' referenced before assignment‍‍‍‍‍

def main():
 import arcpy
 import os
 arcpy.env.overwriteOutput = True
 # Checkout extension
 arcpy.CheckOutExtension("Spatial")
 arcpy.env.overwriteOutput = True
# Workspace Folder
 ws_in_apar = r'D:\MODIS-NDVI-2012\APAR'
 ws_in_tscalar = r'D:\MODIS-NDVI-2012\TScalar'
 ws_in_wscalar = r'D:\MODIS-NDVI-2012\WScalar'
 ws_out_NPP = r'D:\MODIS-NDVI-2012\NPP'

# list "apar" rasters (e.g r001_APAR, r002_APAR so on  )
 arcpy.env.workspace = ws_in_apar
 lst_ras_APAR = arcpy.ListRasters()
 print "lst_ras_APAR", lst_ras_APAR

# list "tscalar" rasters (e.g r001_TSCALAR, r002_TSCALAR so on )
 arcpy.env.workspace = ws_in_tscalar
 lst_ras_TScalar = arcpy.ListRasters()
 print "lst_ras_TScalar", lst_ras_TScalar

 # list "wscalar" rasters (e.g r001_WSCALAR, r002_WSCALAR so on)
 arcpy.env.workspace = ws_in_wscalar
 lst_ras_WScalar = arcpy.ListRasters()
 print "lst_ras_WScalar", lst_ras_WScalar

for ras_name in lst_ras_APAR:
 ras_APAR = arcpy.Raster(os.path.join(ws_in_apar, ras_APAR))
 print "ras_apar", ras_apar

for ras_TSCALAR in lst_ras_TSCALAR:
 ras_TSCALAR = arcpy.Raster(os.path.join(ws_in_tscalar, ras_TSCALAR))
 print "ras_TSCALAR", ras_TSCALAR

 for ras_WSCALAR in lst_ras_WSCALAR:
 ras_WSCALAR = arcpy.Raster(os.path.join(ws_in_wscalar, ras_WSCALAR))
 print "ras_WSCALAR", ras_WSCALAR

# calculate 
 ras_NPP = (ras_APAR * ras_TSCALAR * ras_WSCALAR * 0.98)

#Process rectricted within mask file : Process will occur only on location that fall within the mask
 #else it will be assigned to nodata inthe output
 mask="D:\MOD-REF\NDVI\CROP-L-OCT-15.img"
 NPP_Mask = ExtractByMask(ras_NPP, mask)

 # save raster
 ras_num = ras_name[:3]
 out_name_NPP = os.path.join(ws_out_NPP, 'r{0}_NPP.TIF'.format(ras_num))
 NPP_Mask.save(out_name_NPP)

if __name__ == '__main__':
 main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

DanPatterson_Retired · ‎08-26-2016

ras_APAR = arcpy.Raster(os.path.join(ws_in_apar, ras_APAR))
print "ras_apar", ras_apar <<=====

ras_APAR I suspect

ShouvikJha · ‎08-26-2016

I have changed that line, but still error.

XanderBakker · ‎08-26-2016

shouvik jha , as Dan Patterson correctly mentions keep in mind that variable names are case sensitive.

There are however as lot more things going on. I assume that for each month you want to perform the calculation of the NPP raster. If that is true, you need to look through the month numbers and construct the names for the APAR, TSCALAR and WSCALAR rasters for that month, define the complete name (including the path) for each raster and construct the raster object. Only then you can calculate the NPP raster.

Within the same loop you will have to include the NPP mask raster creation.

ShouvikJha · ‎08-26-2016

xander_bakker‌. Thank you. My humble request Can i get a small example of that modification

ShouvikJha · ‎08-26-2016

While i am running the program i am getting the print of raster file (Upto 27 line), when program execute next line then only error occurring, please give suggestion

Printing file

lst_ras_APAR [u'r001_APAR.TIF', u'r002_APAR.TIF', u'r003_APAR.TIF', u'r004_APAR.TIF', u'r005_APAR.TIF', u'r006_APAR.TIF', u'r007_APAR.TIF', u'r008_APAR.TIF', u'r009_APAR.TIF', u'r010_APAR.TIF', u'r011_APAR.TIF', u'r012_APAR.TIF']
lst_ras_TScalar [u'r001_TSCALAR.TIF', u'r002_TSCALAR.TIF', u'r003_TSCALAR.TIF', u'r004_TSCALAR.TIF', u'r005_TSCALAR.TIF', u'r006_TSCALAR.TIF', u'r007_TSCALAR.TIF', u'r008_TSCALAR.TIF', u'r009_TSCALAR.TIF', u'r010_TSCALAR.TIF', u'r011_TSCALAR.TIF', u'r012_TSCALAR.TIF', u'T_SCALAR_001.tif', u'T_SCALAR_002.tif', u'T_SCALAR_003.tif', u'T_SCALAR_004.tif', u'T_SCALAR_005.tif', u'T_SCALAR_006.tif', u'T_SCALAR_007.tif', u'T_SCALAR_008.tif', u'T_SCALAR_009.tif', u'T_SCALAR_010.tif', u'T_SCALAR_011.tif', u'T_SCALAR_012.tif']
lst_ras_WScalar [u'r001_WSCALAR.TIF', u'r002_WSCALAR.TIF', u'r003_WSCALAR.TIF', u'r004_WSCALAR.TIF', u'r005_WSCALAR.TIF', u'r006_WSCALAR.TIF', u'r007_WSCALAR.TIF', u'r008_WSCALAR.TIF', u'r009_WSCALAR.TIF', u'r010_WSCALAR.TIF', u'r011_WSCALAR.TIF', u'r012_WSCALAR.TIF']
Traceback (most recent call last):
 File "<module2>", line 59, in <module>
 File "<module2>", line 32, in main
UnboundLocalError: local variable 'ras_apar' referenced before assignment
>>>‍‍‍‍‍‍‍‍

ShouvikJha · ‎08-26-2016

Xander Bakker‌. I have tried in my best and and corrected the code suggested by you and Dan Patterson‌. the code is successfully run But i am not getting the desire output from the program, the output name of NPP should be

001_NPP.TIF

but i am getting like this

rr001_WSCALAR.TIF_NPP.TIF

Please cooperation me to solve this issue . Below i attached my updated code

def main():
 import arcpy
 import os
 from arcpy import env
 from arcpy.sa import *
 arcpy.env.overwriteOutput = True
 # Checkout extension
 arcpy.CheckOutExtension("Spatial")
 arcpy.env.overwriteOutput = True
 ws_in_apar = r'D:\MODIS-NDVI-2012\APAR'
 ws_in_tscalar = r'D:\MODIS-NDVI-2012\TScalar'
 ws_in_wscalar = r'D:\MODIS-NDVI-2012\WScalar'
 ws_out_NPP = r'D:\MODIS-NDVI-2012\NPP'

# list "apar" rasters (e.g r001_APAR)
 arcpy.env.workspace = ws_in_apar
 lst_ras_APAR = arcpy.ListRasters()
 print "lst_ras_APAR", lst_ras_APAR

# list "tscalar" rasters (e.g r012_TSCALAR)
 arcpy.env.workspace = ws_in_tscalar
 lst_ras_TSCALAR = arcpy.ListRasters()
 print "lst_ras_TSCALAR", lst_ras_TSCALAR

 # list "wscalar" rasters (e.g r001_WSCALAR)
 arcpy.env.workspace = ws_in_wscalar
 lst_ras_WSCALAR = arcpy.ListRasters()
 print "lst_ras_WSCALAR", lst_ras_WSCALAR

for ras_name in lst_ras_APAR:
 ras_APAR = arcpy.Raster(os.path.join(ws_in_apar, ras_name))
 print "ras_APAR", ras_APAR

for ras_name in lst_ras_TSCALAR:
 ras_TSCALAR = arcpy.Raster(os.path.join(ws_in_tscalar, ras_name))
 print "ras_TSCALAR", ras_TSCALAR

 for ras_name in lst_ras_WSCALAR:
 ras_WSCALAR = arcpy.Raster(os.path.join(ws_in_wscalar, ras_name))
 print "ras_WSCALAR", ras_WSCALAR

# calculate
 ras_NPP = (ras_APAR * ras_TSCALAR * ras_WSCALAR * 0.98)

#Process rectricted within mask file : Process will occur only on location that fall within the mask
 #else it will be assigned to nodata inthe output
 mask="D:\MOD-REF\NDVI\CROP-L-OCT-15.img"
 NPP_Mask = ExtractByMask (ras_NPP, mask)

 # save raster
 out_name_NPP = os.path.join(ws_out_NPP, 'r{0}_NPP.TIF'.format(ras_name))
 NPP_Mask.save(out_name_NPP)
if __name__ == '__main__':
 main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

XanderBakker · ‎08-26-2016

You could try something like this:

import arcpy
import os

def main():
    arcpy.env.overwriteOutput = True

    # Checkout extension
    arcpy.CheckOutExtension("Spatial")
    arcpy.env.overwriteOutput = True
    ws_in_apar = r'D:\MODIS-NDVI-2012\APAR'
    ws_in_tscalar = r'D:\MODIS-NDVI-2012\TScalar'
    ws_in_wscalar = r'D:\MODIS-NDVI-2012\WScalar'
    ws_out_NPP = r'D:\MODIS-NDVI-2012\NPP'

    mask = "D:\MOD-REF\NDVI\CROP-L-OCT-15.img"

    # list "apar" rasters (e.g r001_APAR)
    arcpy.env.workspace = ws_in_apar
    lst_ras_APAR = arcpy.ListRasters()
    print "lst_ras_APAR", lst_ras_APAR

    # list "tscalar" rasters (e.g r012_TSCALAR)
    arcpy.env.workspace = ws_in_tscalar
    lst_ras_TSCALAR = arcpy.ListRasters()
    print "lst_ras_TSCALAR", lst_ras_TSCALAR

    # list "wscalar" rasters (e.g r001_WSCALAR)
    arcpy.env.workspace = ws_in_wscalar
    lst_ras_WSCALAR = arcpy.ListRasters()
    print "lst_ras_WSCALAR", lst_ras_WSCALAR

    for ras_num in [str(i).zfill(3) for i in range(1, 13)]:
        # get required rasters for current month number
        ras_APAR = getRasterFromList(ras_num, lst_ras_APAR, ws_in_apar)
        ras_TSCALAR = getRasterFromList(ras_num, lst_ras_TSCALAR, ws_in_tscalar)
        ras_WSCALAR = getRasterFromList(ras_num, lst_ras_WSCALAR, ws_in_wscalar)

        # verify if raster are valid
        if all([ras_APAR, ras_TSCALAR, ras_WSCALAR]):

            # calculate
            ras_NPP = (ras_APAR * ras_TSCALAR * ras_WSCALAR * 0.98)

            #Process rectricted within mask file : Process will occur only on location that fall within the mask
            #else it will be assigned to nodata inthe output
            NPP_Mask = arcpy.sa.ExtractByMask (ras_NPP, mask)

            # save raster
            out_name_NPP = os.path.join(ws_out_NPP, 'r{0}_NPP.TIF'.format(ras_num))
            NPP_Mask.save(out_name_NPP)


def getRasterFromList(ras_num, lst, ws):
    ras = None
    for ras_name in lst:
        if ras_num in ras_name:
            ras = arcpy.Raster(os.path.join(ws, ras_name))
            break
    if ras is None:
        print ras_num, "not found in lst", lst
    return ras


if __name__ == '__main__':
    main()‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

The idea behind the code is to loop through the months and create the ras_num, then check in each raster list is there is a name that contains that ras_num (see the function getRasterFromList).

If all 3 rasters are defined the calculation is performed.

Once again, I did not test this code.