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

[mdekauwe]

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?



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)	
> 		
> 		
>  
> 

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