Source Code for Module ComboCode.cc.tools.readers.FitsReader

  1  # -*- coding: utf-8 -*- 
  2   
  3  """ 
  4  A module for reading and managing fits-based line profile data files.  
  5   
  6  Author: R. Lombaert 
  7   
  8  """ 
  9   
 10  import os 
 11  from scipy import array,arange 
 12  #import pyfits 
 13  from astropy.io import fits as pyfits 
 14   
 15  import cc.path 
 16  from cc.tools.readers.LPDataReader import LPDataReader 
 17  from cc.tools.io import DataIO 
 18   
 19   
 20   
 21 -def fixFits(fn): 
 22   
 23      ''' 
 24      Fix fits files in case opening them with FitsReader complains about trailing 
 25      spaces and nonstandard values.  
 26       
 27      @param fn: The filename and path to the fitsfile 
 28      @type fn: str 
 29       
 30      ''' 
 31       
 32      hdulist = pyfits.open(fn,mode='update') 
 33      hdulist.flush() 
 34      hdulist.close() 
 35       
 36   
 37   
 38 -def changeFitsHeader(fn_old,fn_new,key,value,comment=None,**kwargs): 
 39   
 40      ''' 
 41      Change a key in the header of a fits file and save to a new file. 
 42   
 43      Additionali just keywords for the writeto pyfits function can be passed along 
 44      with kwargs. eg output_verify='ignore'. 
 45   
 46      @param fn_old: The old filename of the fits file 
 47      @type fn_old: str 
 48      @param fn_new: The new filename of the fits file 
 49      @type fn_new: str 
 50      @param key: The header key to be changed. Has to be present already 
 51                  (for now) 
 52      @type key: str 
 53      @param value: The new value 
 54      @type value: str, int, float 
 55   
 56      @keyword comment: If desired, the comment for the keyword can be changed 
 57                        here. By default, no change is made. 
 58   
 59                        (default: None) 
 60      @type comment: str 
 61   
 62      ''' 
 63       
 64      hdulist = pyfits.open(fn_old, mode = 'update') 
 65      hdulist[0].header.update(key,value,comment=comment) 
 66      hdulist.flush() 
 67      hdulist.close() 
 68   
 69   
 70 -class FitsReader(LPDataReader): 
 71   
 72      ''' 
 73      A FITS file reader for line profile data. 
 74   
 75      ''' 
 76   
 77 -    def __init__(self,fn,star_name=None,*args,**kwargs): 
 78   
 79          ''' 
 80          A FITS file reader for line profile data. 
 81   
 82          Filename of the FITS file is passed to the object upon creation. 
 83           
 84          Additional args/kwargs are used for the dict creation of the parent of  
 85          Reader. 
 86           
 87          @param fn: The FITS filename, including filepath. 
 88          @type fn: string 
 89   
 90          @keyword star_name: The star name if the filename doesn't follow naming 
 91                              conventions. None otherwise. 
 92   
 93                              (default: None) 
 94          @type star_name: str 
 95   
 96          ''' 
 97   
 98          super(FitsReader, self).__init__(fn=fn,star_name=star_name,\ 
 99                                           *args,**kwargs) 
100          self.type = 'fits'  
101          self.readFits() 
102   
103   
104   
105 -    def readFits(self): 
106   
107          ''' 
108          Read the FITS file. 
109   
110          Assumes Tmb flux values in K, with respect to velocity. 
111   
112          The source velocity is taken from the fits file if available. Otherwise 
113          it is taken from Star.dat. 
114   
115          ''' 
116   
117          #- get header, number of grid points and stepsize 
118          hdr = pyfits.getheader(self.fn) 
119          self['contents']['hdr'] = hdr 
120          n_points = hdr.get('naxis1') 
121          if n_points is None or n_points == 1: 
122              raw_input('Could not find array size of fits file.') 
123          crpix1 = hdr.get('CRPIX1') 
124   
125          #- load the data 
126          #try: 
127          lp = pyfits.getdata(self.fn) 
128          #except MemoryError: 
129              #print 'WARNING! Reading %s results in a '%self.fn + \ 
130                    #'MemoryError. Ignoring datafile for now.' 
131              #self['contents']['flux'] = [] 
132              #self['contents']['velocity'] = [] 
133              #return 
134          if len(lp.shape) == 4 and lp.shape[0] == lp.shape[1] == lp.shape[2] == 1: 
135              self['contents']['flux'] = lp[0][0][0] 
136          elif len(lp.shape) == 3 and lp.shape[0] == lp.shape[1] == 1: 
137              self['contents']['flux'] = lp[0][0] 
138          elif lp.shape[0] != 1: 
139              self['contents']['flux'] = lp 
140          else: 
141              raise IOError('Unknown fits file formate for the FitsReader.'+\ 
142                            'Check dimensions fits data. (Then contact Robin)') 
143   
144          #- Check vlsr 
145          vlsr_keys = ['VELO-LSR','VCORR','VLSR'] 
146          vlsr_fits = None 
147          while vlsr_fits is None and not vlsr_keys == []: 
148              vlsr_fits = hdr.get(vlsr_keys.pop(0)) 
149   
150          #- Create velocity OR frequency grid depending on fits file vlsr 
151          #- If vlsr_fits if 0.0, it is possible the velocity grid is already 
152          #- scaled to the right vlsr. In that case, take vlsr from an inputfile. 
153          if vlsr_fits is None: 
154              #- use input vlsr 
155              self.checkVlsr() 
156              #- if None, vlsr not in fits file, so take guess from Star.dat 
157              #- if 0.0, it IS in fits file, and needs to be used as such 
158              crval1 = hdr.get('CRVAL1') 
159              cdelt1 = hdr.get('CDELT1') 
160              ctype1 = hdr.get('CTYPE1') 
161              if ctype1.strip() == 'VRAD': 
162                  #-- If VRAD type on axis 1 use this: Does not need vlsr! 
163                  vel_grid = crval1 + (arange(n_points) - crpix1)*cdelt1 
164              else: 
165                  restfreq = hdr.get('RESTFREQ') 
166                  freq_grid = restfreq + crval1 + (arange(n_points) - crpix1)*cdelt1 
167                  #- c in cm/s, vlsr in km/s, conversion factor of 10**5, final: km/s 
168                  vel_grid = self.getVlsr() - (freq_grid-restfreq)/restfreq*self.c/100000. 
169          else: 
170              #- use vlsr from fits file 
171              #- if it is 0.0, save the initial guess from Star.dat, 
172              #if abs(vlsr_fits) < 1e03: self.checkVlsr() 
173              #- else save value from fits file. 
174              self['contents']['vlsr'] = vlsr_fits/1000. 
175              self.checkVlsr() 
176              deltav = hdr.get('deltav') 
177              #- vlsr_fits in m/s, we want km/s, conversion factor of 1000 
178              #- use vlsr_fits, regardless of it being 0 or not. 
179              vel_grid = (vlsr_fits + (arange(n_points) - crpix1)*deltav)/1000. 
180   
181          #- Make sure the velocity grid is ascending. 
182          if vel_grid[0] > vel_grid [-1]: 
183              vel_grid = vel_grid[::-1] 
184              self['contents']['flux'] = self['contents']['flux'][::-1] 
185          self['contents']['velocity'] = vel_grid 
186   
187          #- Get the date of observation 
188          self['contents']['date_obs'] = hdr.get('DATE-OBS') 
189