Source Code for Module nictools.makemedmask

  1  #! /usr/bin/env python 
  2  # 
  3  # Author: Dave Grumm (based on work by Tomas Dahlen and Eddie Bergeron) 
  4  # Program: makemedmask.py 
  5  # Purpose: routine to create median mask for 'Finesky' 
  6  # History: 03/04/08 - first version 
  7     
  8  import numpy as N 
  9  import pyfits, sys, string, time 
 10  from convolve import boxcar 
 11  from optparse import OptionParser 
 12  import fsutil, opusutil 
 13   
 14  __version__ = "0.1 (2008 Mar 4)" 
 15   
 16  ASIZE = 256 # length of cal array 
 17  ERROR_RETURN = 2 
 18   
 19 -class Makemedmask: 
 20      """ Create and output a median mask from cal files and blt files 
 21   
 22      example:      
 23         m_mask = makemedmask.Makemedmask( medfile='medout2.fits', callist='/hal/data2/dev/nicmos_ped/inlist1.lst', 
 24                  thresh = 0.7,verbosity = 1) 
 25         m_mask.makemask() 
 26      """ 
 27   
 28 -    def __init__( self, thresh=None, medfile=None, callist=None, verbosity=0 ): 
 29          """constructor 
 30   
 31          @param thresh: threshold used in making mask 
 32          @type thresh: real 
 33          @param medfile: name of output masked median image 
 34          @type medfile: string 
 35          @param callist: name of text file containing cal file names 
 36          @type callist: string         
 37          @param verbosity: verbosity level (0 for quiet, 1 verbose, 2 very verbose) 
 38          @type verbosity: string 
 39          """ 
 40   
 41          if (thresh == None):  thresh = fsutil.thresh 
 42          if (medfile == None):  medfile = fsutil.medfile 
 43          if (callist == None):  callist = fsutil.callist 
 44          if (type(thresh) == str) : # ensure that thresh is a float 
 45             thresh = string.atof(thresh) 
 46           
 47          self.thresh = thresh 
 48          self.medfile = medfile 
 49          self.callist = callist 
 50          self.verbosity = verbosity 
 51          self.mm_version = __version__ 
 52          self.mm_run = time.asctime() 
 53   
 54          if (self.verbosity >=1): 
 55              print ' Makemedmask run on ',self.mm_run, ', using code version: ', self.mm_version 
 56   
 57   
 58 -    def makemask( self ): 
 59        """ Make and output mask 
 60        """ 
 61        thresh = self.thresh  
 62        medfile = self.medfile  
 63        callist = self.callist 
 64        verbosity = self.verbosity 
 65            
 66        blot1 = N.zeros((ASIZE,ASIZE), dtype=N.float64)    
 67   
 68        # open callist and read file names 
 69        cfile = open(callist, 'r') 
 70        calfiles = [] 
 71        num_files = 0 
 72        while 1: # first count number of files in list and generate file list 
 73            line = cfile.readline() 
 74            if not line: break 
 75            num_files += 1 
 76            calfiles.append( line )  
 77        cfile.close() 
 78   
 79        bltfiles = [] # list of blot files 
 80         
 81        if (verbosity >=1 ):  print 'There are' ,num_files,'cal files. They are : ' 
 82        for ii in range(num_files): 
 83           calfiles[ii].lstrip().rstrip() # strip leading and trailing whitespace 
 84           calfile_prefix = calfiles[ii].split('_')[0] 
 85           if (verbosity >=1 ): print '  calfiles[',ii,'] = ',calfiles[ii] 
 86   
 87        #  associate blt files with cal files 
 88           bltfile =  calfile_prefix+str("_cal_sci1_blt.fits") 
 89           bltfile.lstrip().rstrip() # strip leading and trailing whitespace       
 90           bltfiles.append( bltfile )  
 91           bltfiles[ii] = bltfile  
 92   
 93        im_cube = N.zeros((ASIZE, ASIZE, num_files), dtype=N.float64)  
 94        blot_cube = N.zeros((ASIZE, ASIZE, num_files), dtype=N.float64)    
 95   
 96        for kk in range(num_files): 
 97           fh_cal = pyfits.open(calfiles[ kk ]) 
 98           fh_blot = pyfits.open(bltfiles[ kk ]) 
 99           im_cube[:,:,kk] = fh_cal[1].data 
100           blot_cube[:,:,kk] = fh_blot[0].data                 
101   
102     # make mask from blotted images 
103        mask_cube = N.zeros((ASIZE, ASIZE, num_files), dtype=N.float64)     
104   
105        for ii in range(num_files): 
106           mm = N.zeros((ASIZE, ASIZE), dtype=N.float64)    
107           dif_0 = blot_cube[:,:,ii] 
108           dif = N.reshape( dif_0,((ASIZE,ASIZE)))  
109           ur =  dif > thresh 
110           mm[ ur ] = 1 
111   
112     # expand the mask 
113           mm = boxcar( mm,(3,3))  # smooth over 3x3 ; this will differ from IDL's "smooth" which ... 
114           #  ... leaves boundary values unchanged, which is not an option in convolve's boxcar 
115   
116           ur =  mm <> 0.0 
117           mm = N.zeros((ASIZE, ASIZE), dtype=N.float64)    
118           mm[ ur ] = 1 
119           mask_cube[:,:,ii] = mm 
120   
121     ## make the masked median image 
122        if (verbosity >=1 ):  print ' Making the masked median image ... ' 
123   
124        maskall= N.zeros((ASIZE, ASIZE), dtype=N.float64)    
125   
126        for jj in range(ASIZE):  
127          for kk in range(ASIZE): 
128             uu = mask_cube[ kk,jj,:] <> 1     
129             im_sub =  im_cube[kk,jj,uu]       
130             im_sub_size = im_sub.size         
131             im_1d = N.reshape( im_sub, im_sub.size)   
132             if ( im_sub_size  > 0 ):  maskall[ kk,jj ]= N.median(im_1d)     
133   
134     # get primary header of 1st cal file to copy to output 
135        fh_cal0 = pyfits.open(calfiles[ 0 ]) 
136        pr_hdr = fh_cal0[0].header 
137   
138        write_to_file(maskall, medfile, pr_hdr, verbosity)  
139   
140        if (verbosity >=1 ):  print 'DONE' 
141   
142 -    def print_pars(self): 
143          """ Print parameters. 
144          """ 
145          print 'The input parameters are :' 
146          print '  thresh:  ' , self.thresh 
147          print '  medfile:  ' , self.medfile 
148          print '  callist:  ' , self.callist 
149   
150      
151   
152 -def write_to_file(data, filename, hdr, verbosity): 
153      """ Write data to specified filename with specified header 
154   
155      @param data: numpy array 
156      @type data: float 
157      @param filename: name of output file 
158      @type filename: string 
159      @param hdr: header for output file 
160      @type hdr: pyfits header object 
161      @param verbosity: verbosity level (0 for quiet, 1 verbose, 2 very verbose) 
162      @type verbosity: string 
163      """ 
164      fimg = pyfits.HDUList() 
165      fimghdu = pyfits.PrimaryHDU( header = hdr) 
166      fimghdu.data = data 
167      fimg.append(fimghdu) 
168      fimg.writeto(filename) 
169      if (verbosity >=1 ): print '...wrote masked median image to: ',filename 
170   
171   
172  if __name__=="__main__": 
173   
174       """Get input file and other arguments, and call CalTempFromBias. 
175          The command-line options are: 
176              -q (quiet) 
177              -v (very verbose) 
178   
179          @param cmdline: command-line arguments 
180          @type cmdline: list of strings 
181       """ 
182   
183       usage = "usage:  %prog [options]" 
184       parser = OptionParser( usage) 
185              
186      # add options and set defaults for parameters 
187       parser.set_defaults( verbosity = fsutil.QUIET) 
188       parser.add_option( "-q", "--quiet", action = "store_const", 
189                          const = fsutil.QUIET, dest = "verbosity",default=None, 
190                          help = "quiet, print nothing") 
191       parser.add_option( "-v", "--verbose", action="store_const", 
192                          const = fsutil.VERY_VERBOSE, dest="verbosity",default=None, 
193                          help="very verbose, print lots of information") 
194       parser.add_option( "-t", "--thresh", dest = "thresh",default = fsutil.thresh, 
195                          help = "threshold for making masked median image.") 
196       parser.add_option( "-m", "--medfile", dest = "medfile",default = fsutil.medfile, 
197                          help = "name of output masked median file.") 
198       parser.add_option( "-c", "--callist", dest = "callist",default = fsutil.callist, 
199                          help = "name of file containing list of cal files.") 
200        
201       (options, args) = parser.parse_args() 
202   
203       fsutil.setVerbosity( options.verbosity) 
204       verbosity = options.verbosity              
205                       
206       fsutil.setThresh(options.thresh ) 
207       if options.thresh!=None: thresh = options.thresh 
208   
209       fsutil.setMedfile(options.medfile ) 
210       if options.medfile!=None: medfile = options.medfile 
211        
212       fsutil.setCallist(options.callist ) 
213       if options.callist!=None: callist = options.callist 
214   
215       try:             
216         m_mask = Makemedmask( thresh=thresh, medfile=medfile, callist=callist, verbosity=verbosity )   
217          
218         if (verbosity >=1 ):  m_mask.print_pars() 
219   
220         m_mask.makemask() 
221   
222         del m_mask 
223   
224       except Exception, errmess: 
225         opusutil.PrintMsg("F","FATAL ERROR "+ str(errmess)) 
226         sys.exit( ERROR_RETURN) 
227