Hi mdekauwe,
as long as your data size is small enough to fit a 8x array in memory
and use it, why not just skip the running total and average 8 data
arrays each 8day period? And skip the x and y loops; these are real
performance killers. As a bonus, your code gets a lot smaller and more
readeable... :)
Please correct me is I'm wrong: the ultimate goal is to have a average
of the valid pixels 8 days of data, each 8 days, correct?
I'd do it this way (more or less pythonic pseudocode, untested, but your
should get the idea):
# read filenames
filenames = glob.glob.....
filenames.sort() # to be sure they are in proper order
filenamesChunked = [filenames[n:n+8] for n in range(0, len(filenames, 8)]
# this will produce a list of lists, with each sublist containing the
# next 8 filenames
nodatavalue = -999
for fchunk in filenamesChunked:
data8days = numpy.ones((8, ncols, nrows)) * -999 # fill with nodata
values, in case there are less than 8 days
for di in range(len(fchunk)):
f = fchunk[di]
data8days[di] = <read f>
weights = data8days > nodatavalue
# this way all nodata values get a weight of 0
avg8days = numpy.average(data8days, axis=0, weights=weights)
<save avg8days>
uhm... well, that's it, really :)
best regards,
Vincent.
mdekauwe wrote:
> 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)
>
>
>
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
