topic Re: Aggregatin Daily Raster to Monthly in Python Questions

Aggregatin Daily Raster to Monthly

DagnachewLegesse — Mon, 30 May 2011 14:50:59 GMT

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)

Re: Aggregatin Daily Raster to Monthly

DagnachewLegesse — Sat, 11 Dec 2021 08:50:06 GMT

Given some more time and thinking, I amended my own code as follows: it seems to work but Once it finishes a subdirectory containing list of rasters, it 'walks' further into the GRID subfolder itself thus giving error. Is there any way that we can prevent it from looking into the GRID folder itself?

try:
    basedir="C:\\N2011\\Data\\NILE\\"
    subdir = "resampled\\"
    datadir=basedir + subdir

    for dirs, subdirs, files in os.walk(datadir):
        env.workspace = os.path.join(dirs)
        outtemp = env.workspace
        rasterList = arcpy.ListRasters("*", "GRID")
#        print rasterList
        if len(rasterList) > 0:
#            print outtemp
            outpath = outtemp.replace ('resampled', 'monthly')
            for i in range(1,13):
                counter = 0
                for infile in rasterList:
                    if not os.path.isdir(outpath):
                        os.makedirs(outpath)
                        print outpath
                    yyyy = int(infile[2:4])
                    mm = int(infile[4:6])            
                    if (mm == i):
                        if counter == 0:
                            rsum1 = arcpy.Raster(infile)
                            counter+=1
                        else:
                            rsum1 = rsum1 + infile
                            counter += 1
                            print infile
                yyyy += 2000
                moutfile="mmpet_" + str(yyyy) + str(i)
                rsum1.save(os.path.join(outpath,moutfile))
                print " ===> " + moutfile
except Exception as e:
    print e.message

Re: Aggregatin Daily Raster to Monthly

DagnachewLegesse — Sat, 11 Dec 2021 08:50:09 GMT

Still trying tom improve it. I really need suggestions though.

#Aggregating daily climate GRIDS to monthly
#DD June 2011

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="C:\\N2011\\Data\\"
    subdir = "resampled\\"
    datadir=basedir + subdir

    for dirs, subdirs, files in os.walk(datadir):
        env.workspace = os.path.join(dirs)
        outtemp = env.workspace
        rasterList = arcpy.ListRasters("*", "GRID")
#        print rasterList
        if len(rasterList) > 1:     #if subdirectory contains only one GRID file --> don't walk in
            print outtemp
            outpath = outtemp.replace ('resampled', 'monthly')
            for i in range(1,13):
                counter = 0
                #Create expression for Cell Statistics
                rstListExp = []
                for infile in rasterList:
                    if not os.path.isdir(outpath):
                        os.makedirs(outpath)
                        print outpath
                    try:
                        yyyy = int(infile[2:4]) + 2000  #if filename et010101
                        mm = int(infile[4:6])
                    except ValueError:
                        yyyy = int(infile[5:9])           #if filename == rain_2001100
                        doy = int(infile[9:])
                        mm = doy2date(doy, yyyy)

                    if (mm == i):
                        rstListExp.append(infile)
 #               print rstListExp
                rsum1 = CellStatistics(rstListExp, "SUM", "DATA")
                moutfile="mm_" + infile[:2] + str(yyyy) + str(i).zfill(2)
                rsum1.save(os.path.join(outpath,moutfile))
                print " ===> " + moutfile
except Exception as e:
    print e.message
    sys.exit(1)
finally:
     print 'Goodbye, world!'

Re: Aggregatin Daily Raster to Monthly

DagnachewLegesse — Sat, 04 Jun 2011 17:01:14 GMT

[PHP]#Aggregating daily climate GRIDS to monthly
#DD June 2011

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="C:\\N2011\\Data\\"
    subdir = "resampled\\"
    datadir=basedir + subdir

    for dirs, subdirs, files in os.walk(datadir):
        env.workspace = os.path.join(dirs)
        outtemp = env.workspace
        rasterList = arcpy.ListRasters("*", "GRID")
#        print rasterList
        if len(rasterList) > 1:     #if subdirectory contains only one GRID file --> don't walk in
            print outtemp
            outpath = outtemp.replace ('resampled', 'monthly')
            for i in range(1,13):
                counter = 0
                #Create expression for Cell Statistics
                rstListExp = []
                for infile in rasterList:
                    if not os.path.isdir(outpath):
                        os.makedirs(outpath)
                        print outpath
                    try:
                        yyyy = int(infile[2:4]) + 2000 #if filename et010101
                        mm = int(infile[4:6])
                    except ValueError:
                        yyyy = int(infile[5:9])           #if filename == rain_2001100
                        doy = int(infile[9:])
                        mm = doy2date(doy, yyyy)

                    if (mm == i):
                        rstListExp.append(infile)
#               print rstListExp
                rsum1 = CellStatistics(rstListExp, "SUM", "DATA")
                moutfile="mm_" + infile[:2] + str(yyyy) + str(i).zfill(2)
                rsum1.save(os.path.join(outpath,moutfile))
                print " ===> " + moutfile
except Exception as e:
    print e.message
    sys.exit(1)
finally:
     print 'Goodbye, world!'[/PHP]