casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables
task_importgmrt.py
Go to the documentation of this file.
00001 ########################################################################3
00002 #  task_importgmrt.py
00003 #
00004 #
00005 # Copyright (C) 2009
00006 # Associated Universities, Inc. Washington DC, USA.
00007 #
00008 # This script is free software; you can redistribute it and/or modify it
00009 # under the terms of the GNU Library General Public License as published by
00010 # the Free Software Foundation; either version 2 of the License, or (at your
00011 # option) any later version.
00012 #
00013 # This library is distributed in the hope that it will be useful, but WITHOUT
00014 # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00015 # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00016 # License for more details.
00017 #
00018 # You should have received a copy of the GNU Library General Public License
00019 # along with this library; if not, write to the Free Software Foundation,
00020 # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00021 #
00022 # Correspondence concerning AIPS++ should be adressed as follows:
00023 #        Internet email: aips2-request@nrao.edu.
00024 #        Postal address: AIPS++ Project Office
00025 #                        National Radio Astronomy Observatory
00026 #                        520 Edgemont Road
00027 #                        Charlottesville, VA 22903-2475 USA
00028 #
00029 # <author>
00030 # Shannon Jaeger (University of Calgary)
00031 # </author>
00032 #
00033 # <summary>
00034 # </summary>
00035 #
00036 # <reviewed reviwer="" date="" tests="" demos="">
00037 # </reviewed
00038 #
00039 # <etymology>
00040 # importgmrt stands for import GMRT UV FITS data.
00041 # </etymology>
00042 #
00043 # <synopsis>
00044 # task_importgmrt.py is a Python script providing an easy to use task
00045 # for importing GMRT FITS data files into the CASA system.  A measurement
00046 # set is produced from the GMRT FITS file.
00047 #
00048 # The given flag files are read and data is flagged according to
00049 # the information found in these files. It is expected that the
00050 # FLAG files are in the old AIPS TVFLAG file format.
00051 #
00052 # </synopsis> 
00053 #
00054 # <example>
00055 # <srcblock>
00056 
00057 # </srblock>
00058 # </example>
00059 #
00060 # <motivation>
00061 # To provide a user-friendly method for importing GMRT FITS files.
00062 # </motivation>
00063 #
00064 # <todo>
00065 # </todo>
00066 
00067 import numpy
00068 import os
00069 from importuvfits import *
00070 from taskinit import *
00071 #from fg import *
00072 
00073 def importgmrt( fitsfile, flagfile, vis ):
00074     retValue=False
00075 
00076     # First check if the output file exists. If it
00077     # does then we abort.  CASA policy prevents files
00078     # from being over-written.
00079     #
00080     # Also if the vis name given is the empty string
00081     # we use the default name, input fitsfile with
00082     # "fits" removed
00083     if ( len( vis ) < 1 ):
00084         vis = ''
00085         fits_list=fitsfile.split( '.' )
00086         last = len(fits_list)
00087         if ( fits_list[last-1].lower() == 'fits' or \
00088              fits_list[last-1].lower() == 'fit' ):
00089             last = last -1
00090         for i in range( last ):
00091             if ( len ( vis ) < 1 ):
00092                 vis = fits_list[i]
00093             else:
00094                 vis = vis + '.'+fits_list[i]
00095         vis=vis+'.MS'
00096         
00097         casalog.post( "The vis paramter is empty, consequently the" \
00098                       +" the default visibilty file\nname will be used, "\
00099                       + vis, 'WARN' )
00100     if ( len( vis ) > 0 and os.path.exists( vis ) ):
00101         casalog.post( 'Visibility file, '+vis+\
00102               ' exists. import GMRT can not proceed, please\n'+\
00103               'remove it or set vis to a different value.', 'SEVERE' )
00104         return retValue
00105     
00106 
00107     # Loop through the list of flag files and make sure each one
00108     # of them exists
00109     if isinstance( flagfile, str ):
00110         if ( len( flagfile ) > 0 ):
00111             flagfile=[flagfile]
00112         else:
00113             flagfile=[]
00114     
00115     ok = True
00116     for i in range( len( flagfile ) ):
00117         if ( not os.path.exists( flagfile[i] ) ):
00118             casalog.post( 'Unable to locate FLAG file '+ flagfile[i],
00119                           'SEVERE' )
00120             ok = False
00121         
00122     if ( not ok ):
00123         return retValue    
00124 
00125 
00126     # Ok, we've done our preliminary checks, now let's get
00127     # down to business and import our fitsfile.
00128     try:
00129         casalog.post( 'Starting import ...', 'NORMAL' )
00130         importuvfits( fitsfile, vis )
00131     except Exception, instance:
00132         casalog.post( str(instance), 'SEVERE' )
00133         return retValue
00134 
00135 #    mytb, myms, myfg = gentools(['tb', 'ms', 'fg'])
00136     mytb, myms = gentools(['tb', 'ms'])
00137     aflocal = casac.agentflagger()
00138 
00139     # Write history
00140     try:
00141         param_names = importgmrt.func_code.co_varnames[:importgmrt.func_code.co_argcount]
00142         param_vals = [eval(p) for p in param_names]
00143         ok &= write_history(myms, vis, 'importgmrt', param_names,
00144                             param_vals, casalog)
00145     except Exception, instance:
00146         casalog.post("*** Error \'%s\' updating HISTORY" % (instance),
00147                      'WARN')
00148     
00149     # CASA's importuvfits doesn't understand the GMRT antenna
00150     # names very well and/or the FITS files tend to put the
00151     # antenna names in a place that CASA doesn't look for them.
00152     # Luckily the information seems to be in the "station"
00153     # information.  So let's fix up the Name column.
00154     casalog.post( 'Correcting GMRT Antenna names.', 'NORMAL1' )
00155     try:
00156         mytb.open( vis+'/ANTENNA', nomodify=False )
00157         stations = mytb.getcol('STATION')
00158         names    = []
00159         for idx in range( 0, len( stations ) ):
00160             names.append( stations[idx].split(':')[0] )
00161         mytb.putcol( 'NAME', numpy.array( names ) )
00162         mytb.done()
00163     except Exception, instance:
00164         casalog.post( 'Unable to properly name the antennas, but continuing')
00165         casalog.post( str(instance), 'WARN' )
00166 
00167     # If we don't have a flagfile then we are done!
00168     # Yippee awe eh!
00169     if ( len( flagfile ) < 1 ):
00170         return True
00171 
00172     # We've imported let's find out the observation start
00173     # and end times, this will be needed for flagging.
00174     # We use it to find the different days the observation
00175     # was done on.
00176     startObs = ''
00177     endObs   = ''
00178     try: 
00179         myms.open( vis )
00180         trange=myms.range( 'time' )['time']
00181         myms.done()
00182 
00183         # Now have the observation time range in a numpy array.
00184         startObs = qa.time( str( trange[0] )+'s', prec=8, form='ymd' )[0]
00185         endObs   = qa.time( str( trange[1] )+'s', prec=8, form='ymd' )[0]
00186         
00187     except Exception, instance:
00188         casalog.post( 'Unable to find obaservation start/en times', 'SEVERE' )
00189         casalog.post( str(instance), 'SEVERE' )
00190         return retValue
00191     
00192     days=[]
00193     startYY = startObs.split('/')[0]
00194     startMM = startObs.split('/')[1]
00195     startDD = startObs.split('/')[2]
00196 
00197     endYY = endObs.split('/')[0]
00198     endMM =endObs.split('/')[1]
00199     endDD =endObs.split('/')[2]    
00200 
00201     for year in range( int( startYY ), int( endYY ) + 1 ):
00202         for month in range( int( startMM ), int( endMM ) + 1 ):
00203             for day in range( int( startDD ), int( endDD ) + 1 ):
00204                 days.append( str(year) + '/' + str(month) + '/' + str(day) + '/' )
00205         casalog.post( "DAYS: "+str(days), 'DEBUG1')
00206 
00207 
00208     # Read through the flag file(s) and assemble it into
00209     # a data structure.  The hope is to minimize the
00210     # number of calls to "flagdata" to make the flagging
00211     # fast!
00212     #
00213     # The expected file format is what is the flag file
00214     # format used by TVFLAG/UVFLAG in AIPS, which is as follows:
00215     #    1. If a line starts with '!' it's a comment
00216     #    2. 
00217     #
00218     casalog.post( 'Starting the flagging ...', 'NORMAL' )
00219                                           
00220         # First lets save the flag information as we got it from the
00221         # flag file.
00222 #    myfg.open( vis )
00223 #    myfg.saveflagversion( 'none', 'No flagging performed yet' )
00224 #    myfg.done()
00225     aflocal.open( vis )
00226     aflocal.saveflagversion( 'none', 'No flagging performed yet' )
00227     aflocal.done()
00228 
00229     for file in flagfile:
00230         casalog.post( 'Reading flag file '+file, 'NORMAL2' )
00231         FLAG_FILE = open( file, 'r' )
00232         line = FLAG_FILE.readline()
00233         ant_names = []          
00234         
00235         while ( len( line ) > 0 ):
00236             casalog.post( 'Read from flag file: '+str(line), 'DEBUG1' )
00237 
00238             # Default antenna list and time range
00239             antennas  = []
00240             baselines = []
00241             timerange = []
00242 
00243             # Skip comment lines, and the end of file.
00244             if ( len(line) < 1 or ( line[0] == '!' and line[2:6] !=  'PANT' ) ):
00245                                 line = FLAG_FILE.readline()
00246                                 continue
00247                         
00248             # Store the antenna name list
00249             if ( line[0:6] == '! PANT' ):                       
00250                 ant_names= line.split(':')[1].lstrip(' ').split(' ')
00251                 line = FLAG_FILE.readline()
00252                 continue
00253 
00254             # First divide lines up on spaces
00255             casalog.post( 'LINE READ: '+line, 'DEBUG1' )
00256             step1=line.split(' ')   
00257             if ( len( step1 ) < 2 ):
00258                 # We want a line with something on it.
00259                 line = FLAG_FILE.readline()
00260                 continue
00261 
00262             # Parse each bit of the string looking for the
00263             # antenna and timrange information.
00264             for i in range( len( step1 ) ):
00265                 step2=step1[i].split('=')
00266                 
00267                 if ( len( step2 ) != 2 ):
00268                     # We want something with an equals sign in it!
00269                     continue
00270                 casalog.post( 'keyword/value: '+str(step2), 'DEBUG1' )
00271 
00272 
00273                 if ( step2[0].upper().startswith( 'ANT' ) ):
00274                     # Should be a comma separated list which
00275                     # is acceptable to flag data
00276                     #antennas=list(step2[1].split(',')
00277                     antennas=step2[1]
00278                     casalog.post( "antennas: "+antennas, 'DEBUG1')
00279 
00280 
00281                 if ( step2[0].upper().startswith( 'BASEL' ) ):
00282                     # Should be a comma separated list which
00283                     # is acceptable to flag data
00284                     baselines=step2[1]
00285                     casalog.post( "baselines: "+baselines, 'DEBUG1')
00286 
00287                 if ( step2[0].upper() == 'TIMERANGE' ):
00288                     # Time in the FLAG file is in the format
00289                     #   start-day,start-hour,start-min,star-sec,\
00290                     #   end-day,end-hour,end-min,end-sec
00291                     #
00292                     # CASA expects data to be of the form
00293                     #   YYYY/MM/DD/HH:MM:SS~YYYY/MM/DD/HH:MM:SS
00294                     # http://casa.nrao.edu/Memos/msselection/node8.html
00295                     #
00296                     # We use the day # to index into our days list and
00297                     # append the time on the end.  But day # starts at
00298                     # 1 and indexing starts at 0 so we need to adjust.
00299                     times     = step2[1].split(',')
00300                     casalog.post( "time: "+str(times), 'DEBUG1')
00301                     if ( int(times[0]) < 1 and int(times[4]) > len(days) ):
00302                         # All data for this antenna
00303                         timerange = ''
00304                     elif ( int(times[0]) < 1 ):
00305                         # We want all time before the end time
00306                         timerange = '<' + days[int(times[4])-1] \
00307                                     + times[5] + ':' + times[6] + ':' +times[7]
00308                     elif ( int(times[4]) > len(days) ):
00309                         # We want all times after the start time
00310                         timerange = '>' + days[int(times[0])-1] \
00311                                     + times[1]+':'+times[2]+':'+times[3]
00312                     else:
00313                         # We ahve a fully specified time range
00314                         timerange = days[int(times[0])-1] \
00315                                 + times[1]+':'+times[2]+':'+times[3]\
00316                                 +'~'\
00317                                 +days[int(times[4])-1] \
00318                                 +times[5]+':'+times[6]+':'+times[7]
00319                 
00320             # Construct the antenna query string
00321             #
00322             # Verify the syntax of the ANTENNA & BASELINE keywords, 
00323             # if ANTENNA will always have only one.  I suspect we
00324             # with to loop through the baselines only here.
00325                         #
00326                         # OLD CODE FOR creating ANT SELECTION STRING
00327                         #
00328             #antStr=''
00329             #selectAntStr=''
00330             #if ( len( antennas ) > 0 and len( baselines ) > 0 ):
00331             #    # We are selecting baselines.
00332             #    for ant1 in antennas:
00333             #        for ant2 in baselines:
00334             #            if ( len( antStr ) < 1 ):
00335             #                antStr=str(int(ant1)-1)+'&'+str(int(ant2)-1)
00336             #    selectAntStr=antStr
00337             #elif ( len( antennas ) > 0 ):
00338             #    antStr=str(int(antennas)-1)
00339             #    # A hack to make sure we always have some unflagged data.
00340             #    # CASA seg faults if the myfg.setdata() results in no data.
00341             #    selectAntStr+=antStr+','+antennas
00342 
00343             
00344             antStr=''
00345             #if ( len( antennas ) > 0 and len( baselines ) > 0 ):
00346                         #       # We are selecting baselines.
00347                         #       for ant1 in antennas:
00348             #        for ant2 in baselines:
00349             #            if ( len( antStr ) < 1 ):
00350             #                antStr=str(int(ant1)-1)+'&'+str(int(ant2)-1)
00351                         #    selectAntStr=antStr
00352             #elif ( len( antennas ) > 0 ):
00353             if ( len( antennas ) > 0 ):
00354                 antStr = ant_names[ int(antennas)-1 ].strip( ' ' )
00355 
00356             try:
00357 #                casalog.post( "flagdata( "+vis+", mode='manualflag', antenna='"
00358 #                              +antStr+"', timerange='"+timerange+"' )", 'NORMAL')
00359                 casalog.post( "flagdata( "+vis+", mode='manual', antenna='"
00360                               +antStr+"', timerange='"+timerange+"' )", 'NORMAL')
00361                 #flagdata( vis, mode='manualflag', selectdata=True, \
00362                 #            antenna=antStr, timerange=timerange, \
00363                 #            flagbackup=False )
00364 #                myfg.open( vis )
00365 #                #print "myfg.setdata( time='",timerange,"' )"
00366 #                myfg.setdata( time=timerange )
00367 #                #print "myfg.setmanualflags( baseline='",antStr,"', time='",timerange,"' )"
00368 #                myfg.setmanualflags( baseline=antStr, time=timerange, unflag=False, clipexpr='' )
00369 #                myfg.run()
00370 #                myfg.done()
00371                 aflocal.open(vis)
00372                 aflocal.selectdata(timerange=timerange)
00373                 aflocal.parsemanualparameters(antenna=antStr, timerange=timerange)
00374                 aflocal.init()
00375                 aflocal.run(True, True)
00376                 aflocal.done()
00377             except Exception, instance:
00378                 casalog.post( 'Unable to flag data from flag file '+file\
00379                               +'.\nAntennas='+antennas+' and timerange='\
00380                               +timerange, 'WARN' )
00381                 casalog.post( str(instance), 'SEVERE' )
00382                 return retValue
00383 
00384             line = FLAG_FILE.readline()
00385             
00386         FLAG_FILE.close()
00387 
00388         # Save a Flag version so we can revert back to it if we wish
00389 #    myfg.open( vis )
00390 #    myfg.saveflagversion( 'import', 'Flagged the data from the GMRT flag file' )
00391 #    myfg.done()
00392     aflocal.open( vis )
00393     aflocal.saveflagversion( 'import', 'Flagged the data from the GMRT flag file' )
00394     aflocal.done()
00395 
00396     retValue = True
00397     return retValue
00398 
00399 
00400     
00401     
00402