|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 import shutil 00003 00004 def fixvis(vis, outputvis='',field='', refcode='', reuse=True, phasecenter='', datacolumn='all'): 00005 """ 00006 Input Parameters 00007 vis -- Name of the input visibility set 00008 00009 outputvis -- Name of the output visibility set, default: same as vis 00010 00011 field -- field selection string 00012 00013 refcode -- Reference frame to convert to, 00014 default: the refcode of PHASE_DIR in the FIELD table 00015 example: 'B1950' 00016 00017 reuse -- base recalculation on existing UVW coordinates? default=True 00018 ignored if parameter 'phasecenter' is set 00019 00020 phasecenter -- if set to a valid direction: change the phase center for the given field 00021 to this value 00022 example: 'J2000 9h25m00s 05d12m00s' 00023 If given without the equinox, e.g. '0h01m00s 00d12m00s', the parameter 00024 is interpreted as a pair of offsets in RA and DEC to the present phasecenter. 00025 00026 datacolumn -- when applying a phase center shift, modify visibilities only in this/these column(s) 00027 default: 'all' (DATA, CORRECTED, and MODEL) 00028 example: 'DATA,CORRECTED' (will not modify MODEL) 00029 00030 Examples: 00031 00032 fixvis('NGC3256.ms','NGC3256-fixed.ms') 00033 will recalculate the UVW coordinates for all fields based on the existing 00034 phase center information in the FIELD table. 00035 00036 fixvis('0925+0512.ms','0925+0512-fixed.ms','0925+0512', '', 'J2000 9h25m00s 05d12m00s') 00037 will set the phase center for field '0925+0512' to the given direction and recalculate 00038 the UVW coordinates. 00039 """ 00040 casalog.origin('fixvis') 00041 retval = True 00042 try: 00043 if(vis==outputvis or outputvis==''): 00044 casalog.post('Will overwrite original MS ...', 'NORMAL') 00045 outputvis = vis 00046 else: 00047 casalog.post('Copying original MS to outputvis ...', 'NORMAL') 00048 shutil.rmtree(outputvis, ignore_errors=True) 00049 shutil.copytree(vis, outputvis) 00050 00051 # me is also used, but not in a state-altering way. 00052 tbt, myms, myim = gentools(['tb', 'ms', 'im']) 00053 00054 if(field==''): 00055 field='*' 00056 00057 fields = myms.msseltoindex(vis=outputvis,field=field)['field'] 00058 00059 if(len(fields) == 0): 00060 casalog.post( "Field selection returned zero results.", 'WARN') 00061 return False 00062 00063 #determine therefcode, the reference frame to be used for the output UVWs 00064 tbt.open(outputvis+"/FIELD") 00065 numfields = tbt.nrows() 00066 therefcode = 'J2000' 00067 ckwdict = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO') 00068 tbt.close() 00069 if(refcode==''): 00070 if(ckwdict.has_key('TabRefTypes')): # we have a variable reference column 00071 therefcode = 'J2000' # always use "J2000" 00072 else: # not a variable reference column 00073 therefcode = ckwdict['Ref'] 00074 else: # a refcode was given, enforce validity 00075 if not (type(refcode)==str): 00076 casalog.post('Invalid refcode '+str(refcode), 'SEVERE') 00077 return False 00078 if(ckwdict.has_key('TabRefTypes')): # variable ref column 00079 refcodelist = ckwdict['TabRefTypes'].tolist() 00080 ref = 0 00081 if not (refcode in refcodelist): 00082 casalog.post('Invalid refcode '+refcode, 'SEVERE') 00083 return False 00084 else: # not a variable reference column 00085 if not (refcode in get_validcodes()): 00086 casalog.post('Invalid refcode '+refcode, 'SEVERE') 00087 return False 00088 #endif 00089 therefcode = refcode 00090 #end if 00091 00092 if(phasecenter==''): # we are only modifying the UVW coordinates 00093 casalog.post('Will leave phase centers unchanged.', 'NORMAL') 00094 casalog.post("Recalculating the UVW coordinates ...", 'NORMAL') 00095 00096 fldids = [] 00097 for i in xrange(numfields): 00098 if (i in fields): 00099 fldids.append(i) 00100 00101 # 00102 myim.open(outputvis, usescratch=False) 00103 myim.calcuvw(fields=fldids, refcode=therefcode, reuse=reuse) 00104 myim.close() 00105 else: # we are modifying UVWs and visibilities 00106 ## if observation: 00107 ## casalog.post('Modifying the phase tracking center(s) is imcompatible', 'SEVERE') 00108 ## casalog.post('with operating on only a subset of the observation IDs', 'SEVERE') 00109 ## return False 00110 00111 if(type(phasecenter) != str): 00112 casalog.post("Invalid phase center.", 'SEVERE') 00113 return False 00114 00115 try: 00116 theoldref, theoldrefstr, ckwdict, isvarref, flddict = get_oldref(outputvis, tbt) 00117 except ValueError: 00118 return False 00119 00120 # for the case of a non-variable reference column and several selected fields 00121 commonoldrefstr = '' 00122 00123 # handle the datacolumn parameter 00124 if (not type(datacolumn)==str): 00125 casalog.post("Invalid parameter datacolumn", 'SEVERE') 00126 elif datacolumn=='' or datacolumn.lower()=='all': 00127 datacolumn='all' 00128 casalog.post("Will modify the visibilities in all columns", 'NORMAL') 00129 else: 00130 # need to check datacolumn before any part of the MS gets modified 00131 wantedcolumns = datacolumn.split(',') 00132 tbt.open(outputvis) 00133 thecolumns = tbt.colnames() 00134 tbt.close() 00135 for col in wantedcolumns: 00136 if not (col.lower() in ['data','observed','corrected', 'corrected_data','model','model_data']): 00137 casalog.post("Invalid datacolumn: \""+col+"\"", 'SEVERE') 00138 return False 00139 if (col.lower()=='observed'): 00140 col = 'DATA' 00141 if (col.lower()=='corrected'): 00142 col = 'CORRECTED_DATA' 00143 if (col.lower()=='model'): 00144 col = 'MODEL_DATA' 00145 if not col.upper() in thecolumns: 00146 casalog.post("Datacolumn "+col+" not present", 'SEVERE') 00147 return False 00148 casalog.post("Will only modify the visibilities in the columns "+datacolumn.upper(), 'NORMAL') 00149 00150 for fld in fields: 00151 commonoldrefstr = modify_fld_vis(fld, outputvis, tbt, myim, 00152 commonoldrefstr, phasecenter, 00153 therefcode, reuse, numfields, 00154 ckwdict, theoldref, theoldrefstr, 00155 isvarref, flddict, datacolumn) 00156 if commonoldrefstr == False: 00157 return False 00158 #endif change phasecenter 00159 00160 # Write history to output MS 00161 try: 00162 param_names = fixvis.func_code.co_varnames[:fixvis.func_code.co_argcount] 00163 param_vals = [eval(p) for p in param_names] 00164 retval &= write_history(myms, outputvis, 'fixvis', param_names, param_vals, 00165 casalog) 00166 except Exception, instance: 00167 casalog.post("*** Error \'%s\' updating HISTORY" % (instance), 'WARN') 00168 except Exception, instance: 00169 casalog.post("*** Error \'%s\' " % (instance), 'SEVERE') 00170 retval = False 00171 00172 return retval 00173 00174 def get_validcodes(code=None): 00175 """Returns a list of valid refcodes.""" 00176 if not get_validcodes._hc: # Is it not already cached? 00177 # Making it once per session is often enough. 00178 validcodes = me.listcodes(me.direction('J2000', '0','0')) 00179 get_validcodes._normal = validcodes['normal'].tolist() 00180 get_validcodes._extra = validcodes['extra'].tolist() 00181 if code == 'extra': 00182 return get_validcodes._extra 00183 elif code == 'normal': 00184 return get_validcodes._normal 00185 else: 00186 return get_validcodes._normal + get_validcodes._extra 00187 get_validcodes._hc = False 00188 00189 def is_valid_refcode(refcode): 00190 """Is refcode usable?""" 00191 return refcode in ['J2000', 'B1950', 'B1950_VLA', 'HADEC', 'ICRS'] \ 00192 + get_validcodes('extra') 00193 00194 def log_phasecenter(oldnew, refstr, ra, dec): 00195 """Post a phase center to the logger, along with whether it is old or new.""" 00196 casalog.post(oldnew + ' phasecenter RA, DEC ' + refstr + ' ' 00197 + qa.time(qa.quantity(ra, 'rad'), 10)[0] # 10 digits precision 00198 + " " + qa.angle(qa.quantity(dec, 'rad'), 10)[0], 'NORMAL') 00199 casalog.post(' RA, DEC (rad) ' + refstr + ' ' 00200 + str(ra) + " " + str(dec), 'NORMAL') 00201 00202 def get_oldref(outputvis, tbt): 00203 """Returns the original reference code, string, ckwdict, and isvarref.""" 00204 theoldref = -1 00205 theoldrefstr = '' 00206 tbt.open(outputvis + "/FIELD") 00207 ckwdict = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO') 00208 flddict = {} 00209 colstoget = ['PHASE_DIR', 'NAME'] 00210 if ckwdict.has_key('TabRefTypes') and ckwdict.has_key('TabRefCodes'): 00211 colstoget.append('PhaseDir_Ref') 00212 for c in colstoget: 00213 flddict[c] = tbt.getcol(c) 00214 if flddict['PHASE_DIR'].shape[1] > 1: 00215 casalog.post('The input PHASE_DIR column is of order ' 00216 + str(flddict['PHASE_DIR'].shape[1] - 1), 'WARN') 00217 casalog.post('Orders > 0 are poorly tested.', 'WARN') 00218 flddict['PHASE_DIR'] = flddict['PHASE_DIR'].transpose((2, 0, 1)) 00219 tbt.close() 00220 if(ckwdict.has_key('TabRefTypes') and ckwdict.has_key('TabRefCodes')): 00221 isvarref = True 00222 else: 00223 isvarref = False 00224 theoldrefstr = ckwdict['Ref'] 00225 return theoldref, theoldrefstr, ckwdict, isvarref, flddict 00226 00227 def modify_fld_vis(fld, outputvis, tbt, myim, commonoldrefstr, phasecenter, 00228 therefcode, reuse, numfields, ckwdict, theoldref, 00229 theoldrefstr, isvarref, flddict, datacol): 00230 """Modify the UVW and visibilities of field fld.""" 00231 viaoffset = False 00232 thenewra_rad = 0. 00233 thenewdec_rad = 0. 00234 thenewref = -1 00235 theolddir = flddict['PHASE_DIR'][fld] 00236 fieldname = flddict['NAME'][fld] 00237 if(ckwdict.has_key('TabRefTypes') and ckwdict.has_key('TabRefCodes')): 00238 # determine string name of the phase dir reference frame 00239 theoldref = flddict['PhaseDir_Ref'][fld] 00240 refcodestrlist = ckwdict['TabRefTypes'].tolist() 00241 refcodelist = ckwdict['TabRefCodes'].tolist() 00242 if not (theoldref in refcodelist): 00243 casalog.post('Invalid refcode in FIELD column PhaseDir_Ref: ' + 00244 str(theoldref), 'SEVERE') 00245 return False 00246 theoldrefstr = refcodestrlist[refcodelist.index(theoldref)] 00247 00248 if not isvarref: 00249 if not (commonoldrefstr == ''): 00250 theoldrefstr = commonoldrefstr 00251 else: 00252 commonoldrefstr = theoldrefstr 00253 00254 00255 theoldphasecenter = theoldrefstr + ' ' + \ 00256 qa.time(qa.quantity(theolddir[0], 'rad'), 14)[0] + ' ' + \ 00257 qa.angle(qa.quantity(theolddir[1],'rad'), 14)[0] 00258 00259 if not is_valid_refcode(theoldrefstr): 00260 casalog.post('Refcode for FIELD column PHASE_DIR is valid but not yet supported: ' 00261 + theoldrefstr, 'WARN') 00262 casalog.post('Output MS may not be valid.', 'WARN') 00263 00264 casalog.post('field: ' + fieldname, 'NORMAL') 00265 log_phasecenter('old', theoldrefstr, theolddir[0], theolddir[1]) 00266 00267 if(therefcode != 'J2000'): 00268 casalog.post( 00269 "When changing phase center, can only write new UVW coordinates in J2000.", 00270 'WARN') 00271 therefcode = 'J2000' 00272 if reuse: 00273 casalog.post("When changing phase center, UVW coordinates will be recalculated.", 00274 'NORMAL') 00275 reuse = False 00276 00277 dirstr = parse_phasecenter(phasecenter, isvarref, theoldref, theoldrefstr, theolddir) 00278 if not dirstr: 00279 return False 00280 00281 if(isvarref): 00282 thenewrefindex = ckwdict['TabRefTypes'].tolist().index(dirstr[0]) 00283 thenewref = ckwdict['TabRefCodes'][thenewrefindex] 00284 thenewrefstr = dirstr[0] 00285 else: # not a variable ref col 00286 validcodes = get_validcodes() 00287 if dirstr[0] in validcodes: 00288 thenewref = validcodes.index(dirstr[0]) 00289 thenewrefstr = dirstr[0] 00290 else: 00291 casalog.post('Invalid refcode ' + dirstr[0], 'SEVERE') 00292 return False 00293 if(dirstr[0] != ckwdict['Ref']): 00294 allselected = True 00295 for i in xrange(numfields): 00296 if not (i in fields): 00297 allselected = False 00298 break 00299 if numfields > 1 and not allselected: 00300 casalog.post("You have not selected all " + str(numfields) 00301 + " fields and PHASE_DIR is not a variable reference column.\n" 00302 + " Please use split or provide phase dir in " + ckwdict['Ref'] 00303 + ".", 'SEVERE') 00304 return False 00305 else: 00306 casalog.post( 00307 "The direction column reference frame in the FIELD table will be changed from " 00308 + ckwdict['Ref'] + " to " + dirstr[0], 'NORMAL') 00309 #endif isvarref 00310 00311 try: 00312 thedir = me.direction(thenewrefstr, dirstr[1], dirstr[2]) 00313 thenewra_rad = thedir['m0']['value'] 00314 thenewdec_rad = thedir['m1']['value'] 00315 except Exception, instance: 00316 casalog.post("*** Error \'%s\' when interpreting parameter \'phasecenter\': " 00317 % (instance), 'SEVERE') 00318 return False 00319 00320 if not is_valid_refcode(thenewrefstr): 00321 casalog.post('Refcode for the new phase center is valid but not yet supported: ' 00322 + thenewrefstr, 'WARN') 00323 casalog.post('Output MS may not be valid.', 'WARN') 00324 00325 if(theolddir[0] >= qa.constants('pi')['value']): # old RA uses range 0 ... 2 pi, not -pi ... pi 00326 while (thenewra_rad < 0.): # keep RA positive in order not to confuse the user 00327 thenewra_rad += 2. * qa.constants('pi')['value'] 00328 00329 log_phasecenter('new', thenewrefstr, thenewra_rad, thenewdec_rad) 00330 00331 # modify FIELD table 00332 tbt.open(outputvis + '/FIELD', nomodify=False) 00333 pcol = tbt.getcol('PHASE_DIR') 00334 pcol[0][0][fld] = thenewra_rad 00335 pcol[1][0][fld] = thenewdec_rad 00336 tbt.putcol('PHASE_DIR', pcol) 00337 00338 casalog.post("FIELD table PHASE_DIR column changed for field " + str(fld) + ".", 00339 'NORMAL') 00340 00341 if(thenewref != -1): 00342 # modify reference of the phase dir column; check also the 00343 # other direction columns 00344 theoldref2 = -1 00345 theoldref3 = -1 00346 if(isvarref): 00347 pcol = tbt.getcol('PhaseDir_Ref') 00348 #theoldref was already determined further above 00349 #theoldref = pcol[fld] 00350 pcol[fld] = thenewref 00351 00352 pcol2 = tbt.getcol('DelayDir_Ref') 00353 theoldref2 = pcol2[fld] 00354 00355 pcol3 = tbt.getcol('RefDir_Ref') 00356 theoldref3 = pcol3[fld] 00357 00358 if(theoldref != thenewref): 00359 tbt.putcol('PhaseDir_Ref', pcol) 00360 casalog.post( 00361 "FIELD table phase center direction reference frame for field " 00362 + str(fld) + " set to " + str(thenewref) + " (" 00363 + thenewrefstr + ")", 'NORMAL') 00364 if not (thenewref == theoldref2 and thenewref == theoldref3): 00365 casalog.post( 00366 "*** The three FIELD table direction reference frame entries for field " 00367 + str(fld) 00368 + " will not be identical in the output data: " 00369 + str(thenewref) + ", " + str(theoldref2) + ", " 00370 + str(theoldref3), 'NORMAL') 00371 if not (theoldref == theoldref2 and theoldref == theoldref3): 00372 casalog.post( 00373 "*** The three FIELD table direction reference frame entries for field " 00374 + str(fld) 00375 + " were not identical in the input data either: " 00376 + str(theoldref) + ", " + str(theoldref2) 00377 + ", " + str(theoldref3), 'NORMAL') 00378 else: 00379 casalog.post( 00380 "FIELD table direction reference frame entries for field " + str(fld) 00381 + " unchanged.", 'NORMAL') 00382 00383 else: # not a variable reference column 00384 tmprec = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO') 00385 if(theoldrefstr != thenewrefstr): 00386 tmprec['Ref'] = thenewrefstr 00387 tbt.putcolkeyword('PHASE_DIR', 'MEASINFO', tmprec) 00388 casalog.post( 00389 "FIELD table phase center direction reference frame changed from " 00390 + theoldrefstr + " to " + thenewrefstr, 'NORMAL') 00391 tbt.close() 00392 00393 fldids = [] 00394 phdirs = [] 00395 for i in xrange(numfields): 00396 if (i==fld): 00397 fldids.append(i) 00398 phdirs.append(theoldphasecenter) 00399 00400 # 00401 myim.open(outputvis, usescratch=False) 00402 myim.fixvis(fields=fldids, phasedirs=phdirs, refcode=therefcode, datacolumn=datacol) 00403 myim.close() 00404 return commonoldrefstr 00405 00406 def parse_phasecenter(phasecenter, isvarref, ref, refstr, theolddir): 00407 dirstr = phasecenter.split(' ') 00408 if len(dirstr) == 2: # interpret phasecenter as an offset 00409 casalog.post("No equinox given in parameter \'phasecenter\': " 00410 + phasecenter, 'NORMAL') 00411 casalog.post("Interpreting it as pair of offsets in (RA,DEC) ...", 00412 'NORMAL') 00413 00414 if(isvarref and ref > 31): 00415 casalog.post('*** Refcode in FIELD column PhaseDir_Ref is a solar system object: ' 00416 + refstr, 'NORMAL') 00417 casalog.post( 00418 '*** Will use the nominal entry in the PHASE_DIR column to calculate new phase center', 00419 'NORMAL') 00420 00421 qra = qa.quantity(theolddir[0], 'rad') 00422 qdec = qa.quantity(theolddir[1], 'rad') 00423 qraoffset = qa.quantity(dirstr[0]) 00424 qdecoffset = qa.quantity(dirstr[1]) 00425 if not qa.isangle(qdecoffset): 00426 casalog.post("Invalid phasecenter parameter. DEC offset must be an angle.", 00427 'SEVERE') 00428 return False 00429 qnewdec = qa.add(qdec,qdecoffset) 00430 qnewra = qra 00431 if qa.canonical(qraoffset)['unit'] == 'rad': 00432 casalog.post( 00433 "RA offset is an angle (not a time). Will divide by cos(DEC) to compute time offset.") 00434 if(qa.cos(qnewdec)['value'] == 0.): 00435 casalog.post( 00436 "Resulting DEC is at celestial pole. Will ignore RA offset.", 'WARN') 00437 else: 00438 qraoffset = qa.div(qraoffset, qa.cos(qnewdec)) 00439 qnewra = qa.add(qnewra, qraoffset) 00440 else: 00441 if not (qa.canonical(qraoffset)['unit'] == 's'): 00442 casalog.post( 00443 "Invalid phasecenter parameter. RA offset must be an angle or a time.", 00444 'SEVERE') 00445 return False 00446 qraoffset = qa.convert(qraoffset, 'deg') 00447 qnewra = qa.add(qnewra, qraoffset) 00448 dirstr = [refstr, qa.time(qnewra,12)[0], qa.angle(qnewdec,12)[0]] 00449 elif not len(dirstr) == 3: 00450 casalog.post('Incorrectly formatted parameter \'phasecenter\': ' 00451 + phasecenter, 'SEVERE') 00452 return False 00453 return dirstr
1.8.0
