[Numpy-discussion] Help making better use of numpy array functions

[mdekauwe]

Hi I have written some code and I would appreciate any suggestions to make
better use of the numpy arrays functions to make it a bit more efficient and
less of a port from C. Any tricks are thoughts would be much appreciated.

The code reads in a series of images, collects a running total if the value
is valid and averages everything every 8 images.

Code below...

Thanks in advance...

#!/usr/bin/env python

"""
Average the daily LST values to make an 8 day product...

Martin De Kauwe
18th November 2009
"""
import sys, os, glob
import numpy as np
import matplotlib.pyplot as plt


def averageEightDays(files, lst, count, numrows, numcols):
	
    day_count = 0
    # starting day - tag for output file
    doy = 122
    
    # loop over all the images and average all the information
    # every 8 days...
    for fname in glob.glob(os.path.join(path, files)):
        current_file = fname.split('/')[8]
        year = current_file.split('_')[2][0:4]
		
        try:
            f = open(fname, 'rb')
        except IOError:
            print "Cannot open outfile for read", fname
            sys.exit(1)
	
        data = np.fromfile(f, dtype=np.float32).reshape(numrows, numcols)
        f.close()
		
        # if it is day 1 then we need to fill up the holding array...
        # copy the first file into the array...
        if day_count == 0:
            lst = data
            for i in xrange(numrows):
                 for j in xrange(numcols):
                     # increase the pixel count if we got vaild data
(undefined = -999.0
                     if lst[i,j] > -900.:
                         count[i,j] = 1	
                         day_count += 1
		
        # keep a running total of valid lst values in an 8 day sequence
        elif day_count > 0 and day_count <= 7:	
            for i in xrange(numrows):
                for j in xrange(numcols):
                    # lst valid pixel?
                    if data[i,j] > -900.:
                        # was the existing pixel valid?
                        if lst[i,j] < -900.:
                            lst[i,j] = data[i,j]
                            count[i,j] = 1
                        else:
                            lst[i,j] += data[i,j]
                            count[i,j] += 1					
                day_count += 1
		
        # need to average previous 8 days and write to a file...
        if day_count == 8:
            for i in xrange(numrows):
                for j in xrange(numcols):
                    if count[i,j] > 0:
                        lst[i,j] = lst[i,j] / count[i,j]
                    else:
                        lst[i,j] = -999.0
			
            day_count = 0
            doy += 8
			
            lst, count = write_outputs(lst, count, year, doy)
			
    # it is possible we didn't have enough slices to do the last 8day avg...
    # but if we have more than one day we shall use it	
    # need to average previous 8 days and write to a file...
    if day_count > 1:
        for i in xrange(numrows):
            for j in xrange(numcols):       
                if count[i,j] > 0:
                    lst[i,j] = lst[i,j] / count[i,j]
                else:
                    lst[i,j] = -999.0
			
        day_count = 0
        doy += 8
		
        lst, count = write_outputs(lst, count, year, doy)
	
def write_outputs(lst, count, year, doy):
    path = "/users/eow/mgdk/research/HOFF_plots/LST/8dayLST"
    outfile = "lst_8day1030am_" + str(year) + str(doy) + ".gra"
	
    try:
        of = open(os.path.join(path, outfile), 'wb')
    except IOError:
        print "Cannot open outfile for write", outfile
        sys.exit(1)
	
    outfile2 = "pixelcount_8day1030am_" + str(year) + str(doy) + ".gra"
    try:
        of2 = open(os.path.join(path, outfile2), 'wb')
    except IOError:
        print "Cannot open outfile for write", outfile2
        sys.exit(1)
	
    # empty stuff and zero 8day counts
    lst.tofile(of)
    count.tofile(of2)	
    of.close()
    of2.close()
    lst = 0.0
    count = 0.0
	
    return lst, count
	
	
if __name__ == "__main__":
	
    numrows = 332
    numcols = 667
	
    path = "/users/eow/mgdk/research/HOFF_plots/LST/gridded_03/"
    lst = np.zeros((numrows, numcols), dtype=np.float32)
    count = np.zeros((numrows, numcols), dtype=np.int)
    averageEightDays('lst_scr_2006*.gra', lst, count, numrows, numcols)		
	
    lst = 0.0
    count = 0.0
    averageEightDays('lst_scr_2007*.gra', lst, count, numrows, numcols)	
		
		
 
-- 
View this message in context: http://old.nabble.com/Help-making-better-use-of-numpy-array-functions-tp26503657p26503657.html
Sent from the Numpy-discussion mailing list archive at Nabble.com.