dagnu70

Aggregatin Daily Raster to Monthly

Discussion created by dagnu70 on May 30, 2011
Latest reply on Jun 4, 2011 by dagnu70
Dear All,

I am trying to write a script that would automate the conversion of daily raster datasets to monthly. My script seems to work but is very slow. I know it can be improved (efficiecny and esthetics). Can it be parallelized? Can any one help? Can a for ... loop be used instead of doing same thing ten times ...

Thanks
DD

 
#Daily to Monthly conversion of satellite data
import sys
import datetime
import arcpy
import os
from arcpy import env
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
print "Script is running, please wait ..."
#Set scratch workspace
arcpy.env.overwriteOutput = True
arcpy.env.scratchWorkspace = "c:\\mytempfile.gdb"
# Set environment settings


def doy2date(doy, yyyy): 
        day=doy
        year=yyyy
        if isinstance(year, str):
                iyear = int(year)
        else:
                iyear = year
        if isinstance(day, str):
                iday = int(day)
        else:
                iday = day
        iday-=1
        thedate = datetime.date(iyear,1,1) + datetime.timedelta(days=iday)
        return (thedate.month)

#try:
basedir="F:\\NILE2011\\Data\\Gabriel_Nile\\resampled\\RFE\\"

#In the basedir I have 10 subdirs (2001 - 2010) inside of which I have daily values
indirs = ['rfe_2001\\','rfe_2002\\', 'rfe_2003\\', 'rfe_2004\\', 'rfe_2005\\',  \
           'rfe_2006\\', 'rfe_2007\\', 'rfe_2008\\', 'rfe_2009\\', 'rfe_2010\\']
for indir in indirs:
    env.workspace = basedir + indir
    outdir="F:\\NILE2011\\Data\\Gabriel_Nile\\resampled\\monthly\\RFE\\"
    outpath=outdir + indir
    if not os.path.isdir(outpath):
       os.mkdir(outpath)
    i=[0,0,0,0, 0, 0, 0, 0, 0, 0, 0, 0]
    rasterList = arcpy.ListRasters("*", "GRID")
    for infile in rasterList:
        yyyy = int(infile[5:9])  #e.g. of filename: rain_20021
        doy = int(infile[9:])
        mm = doy2date(doy, yyyy)
        if (mm == 1):
            if i[0] == 0:
                rsum1 = arcpy.Raster(infile)
                i[0]+=1
            else:
                rsum1 = rsum1 + infile
                i[0] += 1
        elif (mm == 2):
            if i[1] == 0:
                rsum2 = arcpy.Raster(infile)
                i[1]+=1
            else:
                rsum2 = rsum2 + infile
                i[1] += 1
        elif (mm == 3):
            if i[2] == 0:
                rsum3 = arcpy.Raster(infile)
                i[2]+=1
            else:
                rsum3 = rsum3 + infile
                i[2] += 1

        elif (mm == 4):
            if i[3] == 0:
                rsum4 = arcpy.Raster(infile)
                i[3]+=1
            else:
                rsum4 = rsum4 + infile
                i[3] += 1
        elif (mm == 5):
            if i[4] == 0:
                rsum5 = arcpy.Raster(infile)
                i[4]+=1
            else:
                rsum5 = rsum5 + infile
                i[4] += 1
        elif (mm == 6):
            if i[5] == 0:
                rsum6 = arcpy.Raster(infile)
                i[5]+=1
            else:
                rsum6 = rsum6 + infile
                i[5] += 1
        elif (mm == 7):
            if i[6] == 0:
                rsum7 = arcpy.Raster(infile)
                i[6]+=1
            else:
                rsum7 = rsum7 + infile
                i[6] += 1
        elif (mm == 8):
            if i[7] == 0:
                rsum8 = arcpy.Raster(infile)
                i[7]+=1
            else:
                rsum8 = rsum8 + infile
                i[7] += 1
        elif (mm == 9):
            if i[8] == 0:
                rsum9 = arcpy.Raster(infile)
                i[8]+=1
            else:
                rsum9 = rsum9 + infile
                i[8] += 1
        elif (mm == 10):
            if i[9] == 0:
                rsum10 = arcpy.Raster(infile)
                i[9]+=1
            else:
                rsum10 = rsum10 + infile
                i[9] += 1
        elif (mm == 11):
            if i[10] == 0:
                rsum11 = arcpy.Raster(infile)
                i[10]+=1
            else:
                rsum11 = rsum11 + infile
                i[10] += 1
        else:
            if i[11] == 0:
                rsum12 = arcpy.Raster(infile)
                i[11]+=1
            else:
                rsum12 = rsum12 + infile
                i[11] += 1

    for ii in range(0,11):
        print i[ii]            
    rsum1 /=i[0]       
    minfile="etmm_" + str(yyyy) + str(1)
    rsum1.save(os.path.join(outpath,minfile))
    rsum2 /=i[1]        
    minfile="etmm_" + str(yyyy) + str(2)
    rsum2.save(os.path.join(outpath,minfile))
    rsum3 /=i[2]        
    minfile="etmm_" + str(yyyy) + str(3)
    rsum3.save(os.path.join(outpath,minfile))
    rsum4 /=i[3]        
    minfile="etmm_" + str(yyyy) + str(4)
    rsum4.save(os.path.join(outpath,minfile))
    rsum5 /=i[4]        
    minfile="etmm_" + str(yyyy) + str(5)
    rsum5.save(os.path.join(outpath,minfile))
    rsum6 /=i[5]        
    minfile="etmm_" + str(yyyy) + str(6)
    rsum6.save(os.path.join(outpath,minfile))
    rsum7 /=i[6]        
    minfile="etmm_" + str(yyyy) + str(7)
    rsum7.save(os.path.join(outpath,minfile))
    rsum8 /=i[7]        
    minfile="etmm_" + str(yyyy) + str(8)
    rsum8.save(os.path.join(outpath,minfile))
    rsum9 /=i[8]        
    minfile="etmm_" + str(yyyy) + str(9)
    rsum9.save(os.path.join(outpath,minfile))
    rsum10 /=i[9]        
    minfile="etmm_" + str(yyyy) + str(10)
    rsum10.save(os.path.join(outpath,minfile))
    rsum11 /=i[10]        
    minfile="etmm_" + str(yyyy) + str(11)
    rsum11.save(os.path.join(outpath,minfile))
    rsum12 /=i[11]        
    minfile="etmm_" + str(yyyy) + str(12)
    rsum12.save(os.path.join(outpath,minfile))
    print "In the name of the ..... " + str(yyyy)


#except:
#    print arcpy.GetMessages(2)

Outcomes