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-26-2016

xander_bakker Thank you. I ran the code but following error massage i am getting in line no 39

Message File Name Line Position 
Traceback 
 <module> <module1> 65 
 main <module1> 39 
ValueError: The truth value of a raster is ambiguous. Invalid use of raster with Boolean operator or function. Check the use of parentheses where applicable. ‍‍‍‍‍

XanderBakker · ‎08-26-2016

Oops, you probably have to change:

if all([ras_APAR, ras_TSCALAR, ras_WSCALAR]):

... by:
if all([not ras_APAR is None, not ras_TSCALAR is None, not ras_WSCALAR is None]):

ShouvikJha · ‎08-26-2016

Xander Bakker‌. Thank you very much. I changed the line. Now code working perfectly. I am grateful to you for your help.

DanPatterson_Retired · ‎08-26-2016

Xander... boolean slicing ... maybe just as bad as ... is not None, Not is None, that's not the None etc...

a = [None, 'a', None]

sum([[1,0][i is None] for i in a])
1

# then check for > 0 etc
‍‍‍‍‍‍‍‍‍‍‍

XanderBakker · ‎08-26-2016

very true Dan, at least it didn't feel as the "best" solution.

Your example looks good, and in this example line 6 should check == 3, since all three must exist to perform the calculation.

ShouvikJha · ‎09-01-2016

Xander Bakker‌, In just trying to modify above code, instead of multiplying 0.98, i just want to replace it by constant raster , below is my updated code. but its not working,

i tried both line to replace 0.98 by constant raster.

 
 ras_LUE = getRaster(D:\MOD-REF\NDVI\LUE-CASA.img)
 ras_LUE = arcpy.raster (D:\MOD-REF\NDVI\LUE-CASA.img)‍‍‍‍‍‍‍‍‍‍‍‍

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)
        #ras_LUE 
        # verify if raster are valid
        if all([not ras_APAR is None, not ras_TSCALAR is None, not ras_WSCALAR is None]):

        
        
            # 0.98 will replace by constant raster (ras_LUE) 
            ras_NPP = (ras_APAR * ras_TSCALAR * ras_WSCALAR * 0.98)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

XanderBakker · ‎09-02-2016

To access a raster as a Raster object, you will need to use the path as a string in a proper format, like this:

ras_LUE = arcpy.raster(r'D:\MOD-REF\NDVI\LUE-CASA.img')‍

The "r" is used to specify a raw format and avoid the slashes to be interpreted differently (avoid "\N" to be converted to a new line). You can get a good explanation in a blog post by Dan_Patterson here: /blogs/dan_patterson/2016/08/14/filenames-and-file-paths-in-python

Then you can use it in the calculation and you may replace the constant value for the raster object:

ras_NPP = (ras_APAR * ras_TSCALAR * ras_WSCALAR * ras_LUE)‍

ShouvikJha · ‎09-02-2016

xander_bakker‌, Thank you. Now code perfectly running

XanderBakker · ‎09-02-2016

Glad to hear that!

ShouvikJha · ‎09-21-2016

Hi Xander, Please kindly Could you help me out, I am totally not understanding that how i can use the logic to take 8 days time spam raster and do calculation for multi-parameters ,

I am trying to calculate some parameters (parameters mentioned in my code ) using same time span from different band using ArcPy.

I have list of raster, raster name like

MOD09A1.MRTWEB.A2012001.005.sur_refl_b01.tif

MOD09A1.MRTWEB.A2012001.005.sur_refl_b02.tif

MOD09A1.MRTWEB.A2012001.005.sur_refl_b03.tif

MOD09A1.MRTWEB.A2012001.005.sur_refl_b04.tif

MOD09A1.MRTWEB.A2012001.005.sur_refl_b06.tif

MOD09A1.MRTWEB.A2012009.005.sur_refl_b01.tif

MOD09A1.MRTWEB.A2012009.005.sur_refl_b02.tif

MOD09A1.MRTWEB.A2012009.005.sur_refl_b03.tif

MOD09A1.MRTWEB.A2012009.005.sur_refl_b04.tif

MOD09A1.MRTWEB.A2012009.005.sur_refl_b06.tif

MOD09A1.MRTWEB.A2012017.005.sur_refl_b01.tif

MOD09A1.MRTWEB.A2012017.005.sur_refl_b02.tif

MOD09A1.MRTWEB.A2012017.005.sur_refl_b03.tif

MOD09A1.MRTWEB.A2012017.005.sur_refl_b04.tif

MOD09A1.MRTWEB.A2012017.005.sur_refl_b06.tif

In raster name different bands and 8 days interval time spam is there (e.g. A2012001, A2012009, A2012017……………A2012361)

I trying to calculate some parameters for each time spam using different band and will loop same function on all time span (time spam 8 days interval: A2012001, A2012009, A2012017……………A2012361).

For instance, whatever raster bands in 2012001 time span , will calculate different parameters , and save it 001_EVI, 001_LSWI, 001_PSCALAR, after calculation program will take 2012009 time span all bands and do same calculation and save it 002_EVI, 002_LSWI, 002_PSCALAR.

I have written some code. I think its not enough to do it. And I am not understanding how to loop 8 days interval time span to all raster band.

I have attached below, reference raster and raster name list.

def main():
    import arcpy
    from arcpy import env
    from arcpy.sa import *
    import os
    arcpy.env.overwriteOutput = True

    # Checkout extension
    arcpy.CheckOutExtension("Spatial")
    arcpy.env.overwriteOutput = True
    ws_in_mean = r'D:\NPP_GUJARAT\MODIS_DATA\MODIS_NDVI_2012\MASKED_2012C'
    ws_out_EVI= r'D:\MODIS_SURFACE_RF2012\EVI'
    ws_out_LSWI= r'D:\MODIS_SURFACE_RF2012\LSWI'
    ws_out_pscalar = r'D:\MODIS_SURFACE_RF2012\PSCALAR'

    # list "mean" rasters
    arcpy.env.workspace = ws_in_mean
    lst_ras_mean = arcpy.ListRasters()
    print "lst_ras_mean", lst_ras_mean

    for ras_name in lst_ras_mean:
        ras_num = ras_name[:3]

        ras_mean = arcpy.Raster(os.path.join(ws_in_mean, ras_name))
        print "ras_mean", ras_name

        # calculate EVI 2.5 * (band 01 - band 02) / ((band 02 + (6 * band 01) - (7.5 * band 03) + 1) (Min Raster Value )
        ras_EVI = 2.5 * (b01 - b02) / ((b02 + (6 * b01) - (7.5 * b03) + 1)

        # save the output raster for EVI (for 01, 001_EVI, 009_EVI so on.......for 361, 361 _EVI)
        out_name_EVI = os.path.join(ws_out_EVI, 'r{0}_EVI.TIF'.format(ras_num))
        ras_EVI.save(out_name_EVI)

        # Next calculate LSWI from BAND 2 and band 6
        ras_LSWI = (b02 - b06 ) / (b02 + b06)

        # save the output raster for LSWI (for 01, 001_LSWI, 009_LSWI so on.......for 361, 361 _LSWI)
        out_name_LSWI = os.path.join(ws_out_LSWI, 'r{0}_LSWI.TIF'.format(ras_num))
        ras_LSWI.save(out_name_LSWI)

        # Whatever raster we get from step 2, we just take each raster from step 2 and execute below euation
        ras_pscalar = 1 + ras_LSWI / 2

        # save the output raster for pscalar (for 01, 001_PSCALAR, 009_PSCALAR so on.......for 361, 361 _LSWI)
        out_name_pscalar = os.path.join(ws_out_pscalar, 'r{0}_PSCALAR.TIF'.format(ras_num))
        ras_pscalar.save(out_name_pscalar)

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