Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_fixvis.py: 5%
304 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 19:10 +0000
1import shutil
3from casatools import quanta, measures, table, ms, imager
4from casatasks import casalog
5from .mstools import write_history
7_myqa = quanta( )
8_myme = measures( )
10def fixvis(vis, outputvis='',field='', refcode='', reuse=True, phasecenter='', distances='', datacolumn='all'):
11 """
12 Input Parameters
13 vis -- Name of the input visibility set
15 outputvis -- Name of the output visibility set, default: same as vis
17 field -- field selection string
19 refcode -- Reference frame to convert to,
20 default: the refcode of PHASE_DIR in the FIELD table
21 example: 'B1950'
23 reuse -- base recalculation on existing UVW coordinates? default=True
24 ignored if parameter 'phasecenter' is set
26 phasecenter -- if set to a valid direction: change the phase center for the given field
27 to this value
28 example: 'J2000 9h25m00s 05d12m00s'
29 If given without the equinox, e.g. '0h01m00s 00d12m00s', the parameter
30 is interpreted as a pair of offsets in RA and DEC to the present phasecenter.
32 distances -- (experimental) List of the distances (as quanta) of the fields selected by field
33 to be used for refocussing.
34 If empty, the distances of all fields are assumed to be infinity.
35 If not a list but just a single value is given, this is applied to
36 all fields.
37 default: []
38 examples: ['2E6km', '3E6km'] '15au'
40 datacolumn -- when applying a phase center shift, modify visibilities only in this/these column(s)
41 default: 'all' (DATA, CORRECTED, and MODEL)
42 example: 'DATA,CORRECTED' (will not modify MODEL)
44 Examples:
46 fixvis('NGC3256.ms','NGC3256-fixed.ms')
47 will recalculate the UVW coordinates for all fields based on the existing
48 phase center information in the FIELD table.
50 fixvis('0925+0512.ms','0925+0512-fixed.ms','0925+0512', '', 'J2000 9h25m00s 05d12m00s')
51 will set the phase center for field '0925+0512' to the given direction and recalculate
52 the UVW coordinates.
53 """
55 # Note: this is duplicated in task_cvel, and really needing CASA-wide harmonization
56 # (CAS-12871)
57 def copy_ms(src, dest):
58 """ This is a MMS-safe copy of an MMS tree directory.
59 :param src: path to the source MS
60 :param dest: path to the destination MS
61 """
62 shutil.copytree(src, dest, symlinks=True)
64 casalog.origin('fixvis')
65 casalog.post(
66 '**ALERT:** Deprecated in CASA 5.9/6.3. Please use task **phaseshift** instead.', 'WARN'
67 )
69 if vis == outputvis or outputvis == '':
70 casalog.post('Will overwrite original MS ...', 'NORMAL')
71 outputvis = vis
72 else:
73 casalog.post('Copying original MS to outputvis ...', 'NORMAL')
75 shutil.rmtree(outputvis, ignore_errors=True)
76 copy_ms(vis, outputvis)
78 tbt = table( )
79 myms = ms( )
80 myim = imager( )
82 if field == '' or isinstance(field,list) and len(field) == 0:
83 field='*'
85 fields = myms.msseltoindex(vis=outputvis,field=field)['field']
87 if len(fields) == 0:
88 casalog.post( "Field selection returned zero results.", 'WARN')
89 return
91 thedistances = []
93 if distances == "":
94 distances = []
95 elif distances != []:
96 if type(distances) == str and _myqa.isquantity(distances):
97 thedist = _myqa.canonical(distances)
98 if thedist['unit'] == 'm': # a length
99 for f in fields: # put nfields copies into the list
100 thedistances.append(thedist['value'])
101 else:
102 msg = "Parameter distances needs to contain quanta with units of length."
103 raise ValueError(msg)
105 elif type(distances) == list:
106 if len(fields) != len(distances):
107 msg = "You selected "+str(len(fields))+" fields but gave "+str(len(distances))+" distances,"
108 raise ValueError(msg)
109 else:
110 for d in distances:
111 if _myqa.isquantity(d):
112 thedist = _myqa.canonical(d)
113 if thedist['unit'] == 'm': # a length
114 thedistances.append(thedist['value'])
115 else:
116 msg = "Parameter distances needs to contain quanta with units of length."
117 raise ValueError(msg)
118 else:
119 msg = "Invalid parameter distances."
120 raise ValueError(msg)
123 if thedistances != []:
124 casalog.post('Will refocus to the given distances: '+str(distances), 'NORMAL')
126 #determine therefcode, the reference frame to be used for the output UVWs
127 tbt.open(outputvis+"/FIELD")
128 numfields = tbt.nrows()
129 therefcode = 'J2000'
130 ckwdict = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')
131 tbt.close()
132 if refcode == '':
133 if 'TabRefTypes' in ckwdict: # we have a variable reference column
134 therefcode = 'J2000' # always use "J2000"
135 else: # not a variable reference column
136 therefcode = ckwdict['Ref']
137 else: # a refcode was given, enforce validity
138 if not (type(refcode)==str):
139 msg = 'Invalid refcode '+str(refcode)
140 raise RuntimeError(msg)
141 if 'TabRefTypes' in ckwdict: # variable ref column
142 refcodelist = ckwdict['TabRefTypes'].tolist()
143 ref = 0
144 if not (refcode in refcodelist):
145 msg = 'Invalid refcode '+refcode
146 raise RuntimeError(msg)
147 else: # not a variable reference column
148 if not (refcode in get_validcodes()):
149 msg = 'Invalid refcode '+refcode
150 raise RuntimeError(msg)
151 #endif
152 therefcode = refcode
153 #end if
155 if phasecenter == '': # we are only modifying the UVW coordinates
156 casalog.post('Will leave phase centers unchanged.', 'NORMAL')
157 casalog.post("Recalculating the UVW coordinates ...", 'NORMAL')
159 fldids = []
160 for i in range(numfields):
161 if (i in fields):
162 fldids.append(i)
164 #
165 myim.open(outputvis, usescratch=False)
166 myim.calcuvw(fields=fldids, refcode=therefcode, reuse=reuse)
167 myim.close()
168 else: # we are modifying UVWs and visibilities
169 ## if observation:
170 ## casalog.post('Modifying the phase tracking center(s) is imcompatible', 'SEVERE')
171 ## casalog.post('with operating on only a subset of the observation IDs', 'SEVERE')
172 ## return False
174 if type(phasecenter) != str:
175 msg = "Invalid phase center."
176 raise ValueError(msg)
178 theoldref, theoldrefstr, ckwdict, isvarref, flddict = get_oldref(outputvis, tbt)
180 # for the case of a non-variable reference column and several selected fields
181 commonoldrefstr = ''
183 # handle the datacolumn parameter
184 if (not type(datacolumn)==str):
185 casalog.post("Invalid parameter datacolumn", 'SEVERE')
186 elif datacolumn=='' or datacolumn.lower()=='all':
187 datacolumn='all'
188 casalog.post("Will modify the visibilities in all columns", 'NORMAL')
189 else:
190 # need to check datacolumn before any part of the MS gets modified
191 wantedcolumns = datacolumn.split(',')
192 tbt.open(outputvis)
193 thecolumns = tbt.colnames()
194 tbt.close()
195 for col in wantedcolumns:
196 if not (col.lower() in ['data','observed','corrected', 'corrected_data','model','model_data']):
197 msg = "Invalid datacolumn: \""+col+"\""
198 raise RuntimeError(msg)
199 if (col.lower()=='observed'):
200 col = 'DATA'
201 if (col.lower()=='corrected'):
202 col = 'CORRECTED_DATA'
203 if (col.lower()=='model'):
204 col = 'MODEL_DATA'
205 if not col.upper() in thecolumns:
206 msg = "Datacolumn "+col+" not present"
207 raise RuntimeError(msg)
209 casalog.post("Will only modify the visibilities in the columns "+datacolumn.upper(), 'NORMAL')
211 for fld in fields:
212 allselected = True
213 for i in range(numfields):
214 if not (i in fields):
215 allselected = False
216 break
218 commonoldrefstr = modify_fld_vis(fld, outputvis, tbt, myim,
219 commonoldrefstr, phasecenter,
220 therefcode, reuse, numfields,
221 ckwdict, theoldref, theoldrefstr,
222 isvarref, flddict, datacolumn,
223 allselected, thedistances)
224 if commonoldrefstr == False:
225 raise RuntimeError('Failure in modify_fld_vis)')
226 #endif change phasecenter
228 # Write history to output MS
229 try:
230 param_names = fixvis.__code__.co_varnames[:fixvis.__code__.co_argcount]
231 local_vars = locals( )
232 param_vals = [local_vars[p] for p in param_names]
234 write_history(myms, outputvis, 'fixvis', param_names, param_vals,
235 casalog)
236 except Exception as instance:
237 casalog.post("*** Error \'%s\' updating HISTORY" % (instance), 'WARN')
239def get_validcodes(code=None):
240 """Returns a list of valid refcodes."""
241 if not get_validcodes._hc: # Is it not already cached?
242 # Making it once per session is often enough.
243 validcodes = _myme.listcodes(_myme.direction('J2000', '0','0'))
244 get_validcodes._normal = validcodes['normal'].tolist()
245 get_validcodes._extra = validcodes['extra'].tolist()
246 if code == 'extra':
247 return get_validcodes._extra
248 elif code == 'normal':
249 return get_validcodes._normal
250 else:
251 return get_validcodes._normal + get_validcodes._extra
253get_validcodes._hc = False
255def is_valid_refcode(refcode):
256 """Is refcode usable?"""
257 return refcode in ['J2000', 'B1950', 'B1950_VLA', 'HADEC', 'ICRS'] \
258 + get_validcodes('extra')
260def log_phasecenter(oldnew, refstr, ra, dec):
261 """Post a phase center to the logger, along with whether it is old or new."""
262 casalog.post(oldnew + ' phasecenter RA, DEC ' + refstr + ' '
263 + _myqa.time(_myqa.quantity(ra, 'rad'), 10)[0] # 10 digits precision
264 + " " + _myqa.angle(_myqa.quantity(dec, 'rad'), 10)[0], 'NORMAL')
265 casalog.post(' RA, DEC (rad) ' + refstr + ' '
266 + str(ra) + " " + str(dec), 'NORMAL')
268def get_oldref(outputvis, tbt):
269 """Returns the original reference code, string, ckwdict, and isvarref."""
270 theoldref = -1
271 theoldrefstr = ''
272 tbt.open(outputvis + "/FIELD")
273 ckwdict = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')
274 flddict = {}
275 colstoget = ['PHASE_DIR', 'NAME']
276 if 'TabRefTypes' in ckwdict and 'TabRefCodes' in ckwdict:
277 colstoget.append('PhaseDir_Ref')
278 for c in colstoget:
279 flddict[c] = tbt.getcol(c)
280 if flddict['PHASE_DIR'].shape[1] > 1:
281 casalog.post('The input PHASE_DIR column is of order '
282 + str(flddict['PHASE_DIR'].shape[1] - 1), 'WARN')
283 casalog.post('Orders > 0 are poorly tested.', 'WARN')
284 flddict['PHASE_DIR'] = flddict['PHASE_DIR'].transpose((2, 0, 1))
285 tbt.close()
286 if 'TabRefTypes' in ckwdict and 'TabRefCodes' in ckwdict:
287 isvarref = True
288 else:
289 isvarref = False
290 theoldrefstr = ckwdict['Ref']
291 return theoldref, theoldrefstr, ckwdict, isvarref, flddict
293def modify_fld_vis(fld, outputvis, tbt, myim, commonoldrefstr, phasecenter,
294 therefcode, reuse, numfields, ckwdict, theoldref,
295 theoldrefstr, isvarref, flddict, datacol, allselected,
296 thedistances):
297 """Modify the UVW and visibilities of field fld."""
298 viaoffset = False
299 thenewra_rad = 0.
300 thenewdec_rad = 0.
301 thenewref = -1
302 theolddir = flddict['PHASE_DIR'][fld]
303 fieldname = flddict['NAME'][fld]
304 if 'TabRefTypes' in ckwdict and 'TabRefCodes' in ckwdict:
305 # determine string name of the phase dir reference frame
306 theoldref = flddict['PhaseDir_Ref'][fld]
307 refcodestrlist = ckwdict['TabRefTypes'].tolist()
308 refcodelist = ckwdict['TabRefCodes'].tolist()
309 if not (theoldref in refcodelist):
310 msg = 'Invalid refcode in FIELD column PhaseDir_Ref: ' + str(theoldref)
311 raise RuntimeError(msg)
312 theoldrefstr = refcodestrlist[refcodelist.index(theoldref)]
314 if not isvarref:
315 if not (commonoldrefstr == ''):
316 theoldrefstr = commonoldrefstr
317 else:
318 commonoldrefstr = theoldrefstr
321 theoldphasecenter = theoldrefstr + ' ' + \
322 _myqa.time(_myqa.quantity(theolddir[0], 'rad'), 14)[0] + ' ' + \
323 _myqa.angle(_myqa.quantity(theolddir[1],'rad'), 14)[0]
325 if not is_valid_refcode(theoldrefstr):
326 casalog.post('Refcode for FIELD column PHASE_DIR is valid but not yet supported: '
327 + theoldrefstr, 'WARN')
328 casalog.post('Output MS may not be valid.', 'WARN')
330 casalog.post('field: ' + fieldname, 'NORMAL')
331 log_phasecenter('old', theoldrefstr, theolddir[0], theolddir[1])
333 if therefcode != 'J2000':
334 casalog.post(
335 "When changing phase center, can only write new UVW coordinates in J2000.",
336 'WARN')
337 therefcode = 'J2000'
338 if reuse:
339 casalog.post("When changing phase center, UVW coordinates will be recalculated.",
340 'NORMAL')
341 reuse = False
343 dirstr = parse_phasecenter(phasecenter, isvarref, theoldref, theoldrefstr, theolddir)
344 if not dirstr:
345 raise RuntimeError('Failed to parse phasecenter')
347 if isvarref:
348 thenewrefindex = ckwdict['TabRefTypes'].tolist().index(dirstr[0])
349 thenewref = ckwdict['TabRefCodes'][thenewrefindex]
350 thenewrefstr = dirstr[0]
351 else: # not a variable ref col
352 validcodes = get_validcodes()
353 if dirstr[0] in validcodes:
354 thenewref = validcodes.index(dirstr[0])
355 thenewrefstr = dirstr[0]
356 else:
357 msg = 'Invalid refcode ' + dirstr[0]
358 raise RuntimeError(msg)
359 if dirstr[0] != ckwdict['Ref']:
360 if numfields > 1 and not allselected:
361 msg = ("You have not selected all " + str(numfields) +
362 " fields and PHASE_DIR is not a variable reference column.\n"
363 " Please use split or provide phase dir in " + ckwdict['Ref']
364 + ".")
365 raise RuntimeError(msg)
366 else:
367 casalog.post(
368 "The direction column reference frame in the FIELD table will be changed from "
369 + ckwdict['Ref'] + " to " + dirstr[0], 'NORMAL')
370 #endif isvarref
372 try:
373 thedir = _myme.direction(thenewrefstr, dirstr[1], dirstr[2])
374 thenewra_rad = thedir['m0']['value']
375 thenewdec_rad = thedir['m1']['value']
376 except Exception as instance:
377 msg = "*** Error \'%s\' when interpreting parameter \'phasecenter\': " % (instance)
378 raise RuntimeError(msg)
380 if not is_valid_refcode(thenewrefstr):
381 casalog.post('Refcode for the new phase center is valid but not yet supported: '
382 + thenewrefstr, 'WARN')
383 casalog.post('Output MS may not be valid.', 'WARN')
385 if theolddir[0] >= _myqa.constants('pi')['value']: # old RA uses range 0 ... 2 pi, not -pi ... pi
386 while (thenewra_rad < 0.): # keep RA positive in order not to confuse the user
387 thenewra_rad += 2. * _myqa.constants('pi')['value']
389 log_phasecenter('new', thenewrefstr, thenewra_rad, thenewdec_rad)
391 # modify FIELD table
392 tbt.open(outputvis + '/FIELD', nomodify=False)
393 pcol = tbt.getcol('PHASE_DIR')
394 pcol[0][0][fld] = thenewra_rad
395 pcol[1][0][fld] = thenewdec_rad
396 tbt.putcol('PHASE_DIR', pcol)
398 casalog.post("FIELD table PHASE_DIR column changed for field " + str(fld) + ".",
399 'NORMAL')
401 if thenewref != -1:
402 # modify reference of the phase dir column; check also the
403 # other direction columns
404 theoldref2 = -1
405 theoldref3 = -1
406 if isvarref:
407 pcol = tbt.getcol('PhaseDir_Ref')
408 #theoldref was already determined further above
409 #theoldref = pcol[fld]
410 pcol[fld] = thenewref
412 pcol2 = tbt.getcol('DelayDir_Ref')
413 theoldref2 = pcol2[fld]
415 pcol3 = tbt.getcol('RefDir_Ref')
416 theoldref3 = pcol3[fld]
418 if theoldref != thenewref:
419 tbt.putcol('PhaseDir_Ref', pcol)
420 casalog.post(
421 "FIELD table phase center direction reference frame for field "
422 + str(fld) + " set to " + str(thenewref) + " ("
423 + thenewrefstr + ")", 'NORMAL')
424 if not (thenewref == theoldref2 and thenewref == theoldref3):
425 casalog.post(
426 "*** The three FIELD table direction reference frame entries for field "
427 + str(fld)
428 + " will not be identical in the output data: "
429 + str(thenewref) + ", " + str(theoldref2) + ", "
430 + str(theoldref3), 'NORMAL')
431 if not (theoldref == theoldref2 and theoldref == theoldref3):
432 casalog.post(
433 "*** The three FIELD table direction reference frame entries for field "
434 + str(fld)
435 + " were not identical in the input data either: "
436 + str(theoldref) + ", " + str(theoldref2)
437 + ", " + str(theoldref3), 'NORMAL')
438 else:
439 casalog.post(
440 "FIELD table direction reference frame entries for field " + str(fld)
441 + " unchanged.", 'NORMAL')
443 else: # not a variable reference column
444 tmprec = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')
445 if theoldrefstr != thenewrefstr:
446 tmprec['Ref'] = thenewrefstr
447 tbt.putcolkeyword('PHASE_DIR', 'MEASINFO', tmprec)
448 casalog.post(
449 "FIELD table phase center direction reference frame changed from "
450 + theoldrefstr + " to " + thenewrefstr, 'NORMAL')
451 tbt.close()
453 fldids = []
454 phdirs = []
455 for i in range(numfields):
456 if (i==fld):
457 fldids.append(i)
458 phdirs.append(theoldphasecenter)
460 if thedistances==[]:
461 thedistances = 0. # the default value
463 #
464 myim.open(outputvis, usescratch=False)
465 myim.fixvis(fields=fldids, phasedirs=phdirs, refcode=therefcode, datacolumn=datacol, distances=thedistances)
466 myim.close()
467 return commonoldrefstr
469def parse_phasecenter(phasecenter, isvarref, ref, refstr, theolddir):
470 dirstr = phasecenter.split(' ')
471 if len(dirstr) == 2: # interpret phasecenter as an offset
472 casalog.post("No equinox given in parameter \'phasecenter\': "
473 + phasecenter, 'NORMAL')
474 casalog.post("Interpreting it as pair of offsets in (RA,DEC) ...",
475 'NORMAL')
477 if isvarref and ref > 31:
478 casalog.post('*** Refcode in FIELD column PhaseDir_Ref is a solar system object: '
479 + refstr, 'NORMAL')
480 casalog.post(
481 '*** Will use the nominal entry in the PHASE_DIR column to calculate new phase center',
482 'NORMAL')
484 qra = _myqa.quantity(theolddir[0], 'rad')
485 qdec = _myqa.quantity(theolddir[1], 'rad')
486 qraoffset = _myqa.quantity(dirstr[0])
487 qdecoffset = _myqa.quantity(dirstr[1])
488 if not _myqa.isangle(qdecoffset):
489 msg = "Invalid phasecenter parameter. DEC offset must be an angle."
490 raise RuntimeError(msg)
491 qnewdec = _myqa.add(qdec,qdecoffset)
492 qnewra = qra
493 ishms = (_myqa.canonical(qraoffset)['unit'] == 'rad') and (('h' in dirstr[0] and 'm' in dirstr[0] and 's' in dirstr[0]) or (dirstr[0].count(':')==2))
494 if (_myqa.canonical(qraoffset)['unit'] == 'rad') and not ishms:
495 casalog.post(
496 "RA offset is an angle (not a time). Will divide by cos(DEC) to compute time offset.")
497 if _myqa.cos(qnewdec)['value'] == 0.0:
498 casalog.post(
499 "Resulting DEC is at celestial pole. Will ignore RA offset.", 'WARN')
500 else:
501 qraoffset = _myqa.div(qraoffset, _myqa.cos(qnewdec))
502 qnewra = _myqa.add(qnewra, qraoffset)
503 else:
504 if not ((_myqa.canonical(qraoffset)['unit'] == 's') or ishms):
505 msg = "Invalid phasecenter parameter. RA offset must be an angle or a time."
506 raise RuntimeError(msg)
507 # RA offset was given as a time; apply as is
508 qraoffset = _myqa.convert(qraoffset, 'deg')
509 qnewra = _myqa.add(qnewra, qraoffset)
511 dirstr = [refstr, _myqa.time(qnewra,12)[0], _myqa.angle(qnewdec,12)[0]]
513 elif not len(dirstr) == 3:
514 msg = 'Incorrectly formatted parameter \'phasecenter\': ' + phasecenter
515 raise RuntimeError(msg)
516 return dirstr