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

[mdekauwe]

Thanks. Agreed it will break down under that scenario but that shouldn't be
encountered as I am simply checking the value is greater than what I have
set the undefined to be (-999.0). Here is the refjigged logic using numpy
functionality for anyone who it might help. Would appreciate any suggestions
how I can further improve this. The code is significantly more efficient now
(thankfully ;P)!

#!/usr/bin/env python

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

Martin De Kauwe
26th November 2009
"""
import sys, os, glob
import numpy as np

def averageEightDays(files, lst, count, numrows, numcols):
	
	day_count = 0
	doy = 122
	
	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()
		
		# copy the first file into the array...
		# increase the pixel count if we have some valid data at timestep 1
		if day_count == 0:
			lst = data
			count = np.where(lst > -900.0, 1, 0)
			day_count += 1
			
		# keep a running total of valid lst values in an 8 day sequence
		# need to check whether the lst pixel is valid
		elif day_count > 0 and day_count <= 7:	
			lst = np.where(data > -900.0, np.where(lst < -900.0, data, lst + data),
lst)
			count = np.where(data > -900.0, np.where(lst < -900.0, 1, count + 1),
count)
			day_count += 1
		
		# need to average previous 8 days and write to a file...
		if day_count == 8:
			print year, ':', doy
			lst = np.where(count > 0, lst / np.float32(count), -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:
		
		lst = np.where(count > 0, lst / np.float32(count), -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
	lst.tofile(of)
	count.tofile(of2)	
	of.close()
	of2.close()
	lst = np.zeros((numrows, numcols), dtype=np.float32)
	count = np.zeros((numrows, numcols), dtype=np.int)
	
	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 = np.zeros((numrows, numcols), dtype=np.float32)
	count = np.zeros((numrows, numcols), dtype=np.int)
	averageEightDays('lst_scr_2007*.gra', lst, count, numrows, numcols)	
		

Pierre GM-2 wrote:
> 
> On Nov 25, 2009, at 4:13 PM, mdekauwe wrote:
>> I tried redoing the internal logic for example by using the where
>> function
>> but I can't seem to work out how to match up the logic. For example (note
>> slightly different from above). If I change the main loop to
>> 
>> lst = np.where((data > -900.0) & (lst < -900.0), data, lst)
>> lst = np.where((data > -900.0) & (lst > -900.0), 5.0, lst)
>> 
>> If I then had a data array value of 10.0 and lst array value of -999.0.
>> In
>> this scenario the first statement would set the LST value to 10.0.
>> However
>> when you run the second statement, data is still bigger than the
>> undefined
>> (-999.0) but now so is the LST, hence the lst is set to 5.0
>> 
>> This doesn't match with the logic of            
>> 
>> if data[i,j] > -900.:
>>    if lst[i,j] < -900.:
>>        lst[i,j] = data[i,j]
>>    elif lst[i,j] > -900.:
>>        lst[i,j] = 5.0
>> 
>> as it would never get to the second branch under this logic.
>> 
>> Does anyone have any suggestions on how to match up the logic?
>> 
> 
> 
> np.where(data>-900,np.where(lst<-900,data,5.),lst)
> 
> Note that this should fail if lst=-900
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
> 
> 

-- 
View this message in context: http://old.nabble.com/Help-making-better-use-of-numpy-array-functions-tp26503657p26525315.html
Sent from the Numpy-discussion mailing list archive at Nabble.com.