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'‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

XanderBakker · ‎09-22-2016

Hi shouvik jha

See below a snippet which shows how you can use that list you posted to obtain the bands for your calculations:

def main():
    import os
    import arcpy

    # define all your workspaces here


    # get the list of rasters (bands)
##    ws_bands = r'your path to bands'
##    arcpy.env.workspace = ws_bands
##    lst_ras = arcpy.ListRasters()

    # tmp list with bands to show flow
    lst_ras = ['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']

    dct = {}
    for ras_name in lst_ras:
        date_id =  ras_name[16:23]
        if date_id in dct:
            lst_bands = dct[date_id]
            lst_bands.append(ras_name)
            dct[date_id] = lst_bands
        else:
            lst_bands = [ras_name]
            dct[date_id] = lst_bands

    cnt = 0
    for date_id, lst_bands in sorted(dct.items()):
        cnt += 1
        ras_num = "%03d" % (cnt,)
        print ras_num, date_id, lst_bands
        b01_name = getBand(lst_bands, 1)
        b02_name = getBand(lst_bands, 2)
        b03_name = getBand(lst_bands, 3)
        b04_name = getBand(lst_bands, 4)
        b06_name = getBand(lst_bands, 6)

        b01 = os.path.join(ws_bands, b01_name)
        b02 = os.path.join(ws_bands, b02_name)
        b03 = os.path.join(ws_bands, b03_name)
        b04 = os.path.join(ws_bands, b04_name)
        b06 = os.path.join(ws_bands, b06_name)

        if sum([[1,0][i is None] for i in [b01, b02, b03, b04, b06]]) == 5:
            # do your calculations with the bands here
            print "Bands found, proceed with calculations..."
            pass


def getBand(lst_bands, band_number):
    band_id = "b{0}.tif".format("%02d" % (band_number,))
    band_found = None
    for band in lst_bands:
        if band.lower().endswith(band_id):
            band_found = band
            break
    return band_found


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

If you run the code it will print this:

001 2012001 ['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']
Bands are OK, proceed with calculations
002 2012009 ['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']
Bands are OK, proceed with calculations
003 2012017 ['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']
Bands are OK, proceed with calculations
‍‍‍‍‍‍

There is no need to convert the date in any way, you just simply perform the calculations for each "date" on the band of that date.

ShouvikJha · ‎09-23-2016

Xander Bakker‌, Hearty gratitude for your guidance, I have written the code as per your suggestions above, But line no 70 shows error massage

"SyntaxError: invalid syntax". I tried many times to solve but i could not able to solve.

And whatever raster name we mentioned in line no 17 that is only for 3 time span, similarly 8 days interval i have till A2012001, A2012009, A2012017.............A2012361 raster. So All raster name i need to mention?

Line no 12 : ws_bands - I have to give folder path where my raster data are stored ?

Below is updated code,

def main():
    import os
    import arcpy

    # define all your workspaces here
    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'

    # get the list of rasters (bands)
    ws_bands = r'D:\NPP_GUJARAT\MODIS_DATA\MODIS_NDVI_2012\MASKED_2012C' # is it correct to give folder path for all bands?
    arcpy.env.workspace = ws_bands
    lst_ras = arcpy.ListRasters()

    # tmp list with bands to show flow (Here i have to mention all raster name?)
    lst_ras = ['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']

    dct = {}
    for ras_name in lst_ras:
        date_id =  ras_name[16:23]
        if date_id in dct:
            lst_bands = dct[date_id]
            lst_bands.append(ras_name)
            dct[date_id] = lst_bands
        else:
            lst_bands = [ras_name]
            dct[date_id] = lst_bands

    cnt = 0
    for date_id, lst_bands in sorted(dct.items()):
        cnt += 1
        ras_num = "%03d" % (cnt,)
        print ras_num, date_id, lst_bands
        b01_name = getBand(lst_bands, 1)
        b02_name = getBand(lst_bands, 2)
        b03_name = getBand(lst_bands, 3)
        b04_name = getBand(lst_bands, 4)
        b06_name = getBand(lst_bands, 6)

        b01 = os.path.join(ws_bands, b01_name)
        b02 = os.path.join(ws_bands, b02_name)
        b03 = os.path.join(ws_bands, b03_name)
        b04 = os.path.join(ws_bands, b04_name)
        b06 = os.path.join(ws_bands, b06_name)

        if sum([[1,0][i is None] for i in [b01, b02, b03, b04, b06]]) == 5:

                          # do your calculations with the bands here#


        # calculate EVI 2.5 * (band 01 - band 02) / ((band 02 + (6 * band 01) - (7.5 * band 03) + 1)
         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)


         print "Bands found, proceed with calculations..."
         pass

def getBand(lst_bands, band_number):
    band_id = "b{0}.tif".format("%02d" % (band_number,))
    band_found = None
    for band in lst_bands:
        if band.lower().endswith(band_id):
            band_found = band
            break
    return band_found


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

XanderBakker · ‎09-25-2016

Hi shouvik jha ,

The first calculation had an unpaired parenthesis. I added a closing bracket at the end and made some additional changes.

Please try the code below:

def main():
    import os
    import arcpy

    arcpy.env.overwriteOutput = True
    arcpy.CheckOutExtension("Spatial")

    # define all your workspaces here
    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'

    # get the list of rasters (bands)
    ws_bands = r'D:\NPP_GUJARAT\MODIS_DATA\MODIS_NDVI_2012\MASKED_2012C'
    arcpy.env.workspace = ws_bands
    lst_ras = arcpy.ListRasters()

    dct = {}
    for ras_name in lst_ras:
        date_id =  ras_name[16:23]
        if date_id in dct:
            lst_bands = dct[date_id]
            lst_bands.append(ras_name)
            dct[date_id] = lst_bands
        else:
            lst_bands = [ras_name]
            dct[date_id] = lst_bands

    cnt = 0
    for date_id, lst_bands in sorted(dct.items()):
        cnt += 1
        ras_num = "%03d" % (cnt,)
        print ras_num, date_id, lst_bands
        b01_name = getBand(lst_bands, 1)
        b02_name = getBand(lst_bands, 2)
        b03_name = getBand(lst_bands, 3)
        b04_name = getBand(lst_bands, 4)
        b06_name = getBand(lst_bands, 6)

        if sum([[1,0][i is None] for i in [b01_name, b02_name, b03_name, b04_name, b06_name]]) == 5:

            b01 = arcpy.Raster(os.path.join(ws_bands, b01_name))
            b02 = arcpy.Raster(os.path.join(ws_bands, b02_name))
            b03 = arcpy.Raster(os.path.join(ws_bands, b03_name))
            b04 = arcpy.Raster(os.path.join(ws_bands, b04_name))
            b06 = arcpy.Raster(os.path.join(ws_bands, b06_name))

            # calculate EVI 2.5 * (band 01 - band 02) / ((band 02 + (6 * band 01) - (7.5 * band 03) + 1)
            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.0 + ras_LSWI / 2.0

            # 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)


def getBand(lst_bands, band_number):
    band_id = "b{0}.tif".format("%02d" % (band_number,))
    band_found = None
    for band in lst_bands:
        if band.lower().endswith(band_id):
            band_found = band
            break
    return band_found


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

XanderBakker · ‎09-25-2016

For additional questions, you can better start a new question. This one is getting very long and slow to edit.

DanPatterson_Retired · ‎09-25-2016

Xander... coolest line in the whole script

... sum([[1,0][i is None] for i in [b01_name, b02_name, b03_name, b04_name, b06_name]]) == 5 ...

XanderBakker · ‎09-26-2016

I guess it is the coolest line, since it is your line, right? 😉

Thanks again Dan!

ShouvikJha · ‎09-26-2016

Xander Bakker‌. Thank you so much. An excellent Code written by you. Code running perfectly. Once again my Hearty gratitude to you.

Regarding additional question definitely next time i will start my query as a new question.

XanderBakker · ‎09-26-2016

You are welcome!

ShouvikJha · ‎09-27-2016

Xander Bakker‌, I ran the code, code running successfully.

But when I compared the output between from our above program generated and by manually Arc GIS raster calculator generated output of EVI parameter of single time span (A2012001, A2012009, A201217) using Reference raster which i had attached earlier. I am not getting same output from both program. Manually generated output are differ from program generated output. We assumed manually generated output would be correct. kindly give suggestion to resolve the issue .

As you mentioned to start a new question for additional question.My humble request you, After this topic we will start new query as a new question. If slow issue make more trouble to you then i will re-post the code as a new question. Sorry for inconvenience cause.

XanderBakker · ‎09-27-2016

As I mentioned before your line:

ras_EVI = 2.5 * (b01 - b02) / ((b02 + (6 * b01) - (7.5 * b03) + 1)

... is missing a closing bracket. I simply added it to the end, but perhaps that was not correct. Could you verify that? And also to be sure change the "6" into 6.0 and the "1" into 1.0 to avoid intermediate rounding?