Building python script to compute monthly average from 8 days MODIS data

Question

Goal- Building python script to compute monthly average from 8 days MODIS data

I have MODIS data which is contain 8 days product. Some of month has 3 images and some of month has four imgaes. So how i can build a Arc GIS python scrip to calculate monthly average from this data

I have listed below number of images in each month.

Jan - 4, Feb- 4, Mar- 4, Apr- 3, May- 3, June-4, July- 4, Aug - 4, Sept - 4, Oct - 3, Nov - 4, Dec - 4,

Total 45 images are in same foldar with 8 days interval name

e.g. MOD15A2.MRTWEB.A2005337.005.Fpar_1km, MOD15A2.MRTWEB.A2005345.005.Fpar_1km

After calculation i just i want to save it in month wise name.

I have Build one code but i can take only defined interval e.g 3 days or 4 days. not 3 days and 4 days at a same time.

Below i have shown my python code

import arcpy, os, sys 
from arcpy import env 
from arcpy.sa import * 
arcpy.env.overwriteOutput = True 
arcpy.CheckOutExtension("Spatial") 
arcpy.env.extent = "MAXOF" 
env.workspace = 'D:\MODIS-NDVI\FPAR-2005\MASKED-FPAR-05B' 
rasters = arcpy.ListRasters() 
out_ws = 'D:\MODIS-NDVI\FPAR-2005\AVG-FPAR-05' 
rastersCount = len(rasters) 
counter = 0 
dek = 1 
while counter < rastersCount: 
 dekad = [] 
 for i in range(counter,(counter+2)): 
 print i 
 dekad.append(rasters[i]) 
 outCellStatistics = CellStatistics(dekad, "MEAN", "NODATA") 
 outRasterName = out_ws + 'tmax_dek_{:03d}.tif'.format (dek) 
 #outRasterName = "00" + str(dek) + ".tif" 
 outCellStatistics.save(outRasterName) 
 counter += 2 
 dek += 1 

Above can take only defined interval e.g 3 days or 4 days. not 3 days and 4 days at a same time.

Answer 1

It sounds like you just need to iterate over each month's folder(?), add each month's tif files to the list "dekad", perform the CellStatistics command and save result as the monthly average raster, correct?

Instead of the "for i in range... dekad.append[]..." method, you could do something like below using "glob"...

import arcpy, os, glob
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
rootdir = 'D:/MODIS-NDVI/FPAR-2005/MASKED-FPAR-05B'
dekad = []
for month in os.listdir(rootdir):
 if os.path.isdir(os.path.join(rootdir, item)):
 print month 
 work_dir = os.path.join(rootdir, month)
 final_raster = os.path.join(work_dir, month) + '_mean.tif'
 for raster in glob.glob(work_dir + "*.tif"):
 dekad.append(raster)
 print dekad
 print "Creating a MEAN raster!"
 outfile = CellStatistics(dekad, "MEAN", "NODATA")
 outfile.save(final_raster)
 dekad = [] 

• Shouvik Jha: Thank you. I have call Raster in same folder only. Just for better understanding purpose i mentioned month wise number. But all images in same folder only. I have shared image details with folder
• Paul: You should use os.path.join() to concatenate your paths.
• Shouvik Jha: Thank you Sir. We have all data in same folder only. I just mentioned no of images in different month for better understanding purpose. Below i have mention my all input file names MOD15A2.MRTWEB.A2005001.005.Fpar_1km MOD15A2.MRTWEB.A2005009.005.Fpar_1km MOD15A2.MRTWEB.A2005017.005.Fpar_1km MOD15A2.MRTWEB.A2005025.005.Fpar_1km MOD15A2.MRTWEB.A2005033.005.Fpar_1km MOD15A2.MRTWEB.A2005041.005.Fpar_1km So on.......

Answer 2

I'm not familiar with the MODIS file naming conventions, but maybe you could try the regular expression module (re) to add rasters matching a pattern in your root folder, add them to "dekad" based on the pattern, then perform CellStatistics. Below, I've assumed the month is the 3 digits before the "Fpar_1km" (??) You could adjust it to suit your needs.

import arcpy, os, glob, re
from arcpy import env
from arcpy.sa import *
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
rootdir = 'D:/MODIS-NDVI/FPAR-2005/MASKED-FPAR-05B/'
env.workspace = rootdir
dekad = []
mons = ['005', '006'] # <-- adjust as necessary
for mon in mons:
 patt = r"MOD15A2\.MRTWEB\.A2005\d{3}\." + mon + r"\.Fpar_1km\.tif"
 prog = re.compile(patt)
 for raster in glob.glob(rootdir + "*.tif"):
 rastername = os.path.basename(raster)
 if prog.match(rastername):
 work_file = os.path.join(rootdir, rastername)
 print 'Work file: ' + work_file
 dekad.append(rastername)
 print dekad
 print "Creating a MEAN raster!"
 final_raster = os.path.join(rootdir, mon) + '_mean.tif'
 outfile = CellStatistics(dekad, "MEAN", "NODATA")
 outfile.save(final_raster)
 dekad = [] 

• Shouvik Jha: Thank you Sir, which you assumed as month of images "3 digits before the "Fpar_1km" 005" its not month. it is common in every images. name format of images for different time is like A200500.005, A2005009.005, A2005017.005, A2005025.005, A2005033.005. name of different time period is stared after the year "2005" and its 8 days interval its increase . similarly i have the raster till A2005361 . First massage i have mentioned the full image name. Image name is changes 3 digit after the year (e.g 2005)

Answer 3

I don't know of any other way to do this unless you could access the image's EXIF data. Test "datetime" or "datetimeoriginal" and check the month the image was taken. The Python image library (PIL) module can do this like below.

import arcpy, os, glob, re
from arcpy import env
from arcpy.sa import *
from PIL import Image
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("Spatial")
rootdir = 'D:/MODIS-NDVI/FPAR-2005/MASKED-FPAR-05B/'
env.workspace = rootdir
dekad = []
mons = ['05', '06'] # <-- adjust as necessary
def get_date_taken(path):
 # returns datetimeoriginal like "2014:02:18 14:25:24"
 return Image.open(path)._getexif()[36867]
for mon in mons:
 # match datetimeoriginal on month 
 patt = r"\d{4}\:" + mon + r"\:\d{2}\s\d{2}\:\d{2}\:\d{2}"
 prog = re.compile(patt)
 for raster in glob.glob(rootdir + "*.jpg"):
 image_date = get_date_taken(raster)
 print image_date
 rastername = os.path.basename(raster)
 print rastername
 normrast = os.path.normpath(raster)
 if prog.match(image_date):
 print 'it matched!'
 work_file = os.path.join(rootdir, rastername)
 dekad.append(rastername)
 print dekad

• cm1: Glad to hear it. Regards - cm1
• cm1: Happy to help, and welcome to Stack Overflow. If this answer or any other one solved your issue, please mark it as accepted.