Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_fixvis.py: 70%

304 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1import shutil 

2 

3from casatools import quanta, measures, table, ms, imager 

4from casatasks import casalog 

5from .mstools import write_history 

6 

7_myqa = quanta( ) 

8_myme = measures( ) 

9 

10def fixvis(vis, outputvis='',field='', refcode='', reuse=True, phasecenter='', distances='', datacolumn='all'): 

11 """ 

12 Input Parameters 

13 vis -- Name of the input visibility set 

14  

15 outputvis -- Name of the output visibility set, default: same as vis 

16 

17 field -- field selection string 

18 

19 refcode -- Reference frame to convert to, 

20 default: the refcode of PHASE_DIR in the FIELD table 

21 example: 'B1950' 

22  

23 reuse -- base recalculation on existing UVW coordinates? default=True 

24 ignored if parameter 'phasecenter' is set 

25 

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. 

31 

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' 

39 

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) 

43 

44 Examples: 

45 

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. 

49 

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 """ 

54 

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) 

63 

64 casalog.origin('fixvis') 

65 casalog.post( 

66 '**ALERT:** Deprecated in CASA 5.9/6.3. Please use task **phaseshift** instead.', 'WARN' 

67 ) 

68 

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') 

74 

75 shutil.rmtree(outputvis, ignore_errors=True) 

76 copy_ms(vis, outputvis) 

77 

78 tbt = table( ) 

79 myms = ms( ) 

80 myim = imager( ) 

81 

82 if field == '' or isinstance(field,list) and len(field) == 0: 

83 field='*' 

84 

85 fields = myms.msseltoindex(vis=outputvis,field=field)['field'] 

86 

87 if len(fields) == 0: 

88 casalog.post( "Field selection returned zero results.", 'WARN') 

89 return 

90 

91 thedistances = [] 

92 

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) 

104 

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) 

121 

122 

123 if thedistances != []: 

124 casalog.post('Will refocus to the given distances: '+str(distances), 'NORMAL') 

125 

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 

154 

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') 

158 

159 fldids = [] 

160 for i in range(numfields): 

161 if (i in fields): 

162 fldids.append(i) 

163 

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 

173 

174 if type(phasecenter) != str: 

175 msg = "Invalid phase center." 

176 raise ValueError(msg) 

177 

178 theoldref, theoldrefstr, ckwdict, isvarref, flddict = get_oldref(outputvis, tbt) 

179 

180 # for the case of a non-variable reference column and several selected fields  

181 commonoldrefstr = '' 

182 

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) 

208 

209 casalog.post("Will only modify the visibilities in the columns "+datacolumn.upper(), 'NORMAL') 

210 

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 

217 

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 

227 

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] 

233 

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') 

238 

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 

252 

253get_validcodes._hc = False 

254 

255def is_valid_refcode(refcode): 

256 """Is refcode usable?""" 

257 return refcode in ['J2000', 'B1950', 'B1950_VLA', 'HADEC', 'ICRS'] \ 

258 + get_validcodes('extra') 

259 

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') 

267 

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 

292 

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)] 

313 

314 if not isvarref: 

315 if not (commonoldrefstr == ''): 

316 theoldrefstr = commonoldrefstr 

317 else: 

318 commonoldrefstr = theoldrefstr 

319 

320 

321 theoldphasecenter = theoldrefstr + ' ' + \ 

322 _myqa.time(_myqa.quantity(theolddir[0], 'rad'), 14)[0] + ' ' + \ 

323 _myqa.angle(_myqa.quantity(theolddir[1],'rad'), 14)[0] 

324 

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') 

329 

330 casalog.post('field: ' + fieldname, 'NORMAL') 

331 log_phasecenter('old', theoldrefstr, theolddir[0], theolddir[1]) 

332 

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 

342 

343 dirstr = parse_phasecenter(phasecenter, isvarref, theoldref, theoldrefstr, theolddir) 

344 if not dirstr: 

345 raise RuntimeError('Failed to parse phasecenter') 

346 

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 

371 

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) 

379 

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') 

384 

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'] 

388 

389 log_phasecenter('new', thenewrefstr, thenewra_rad, thenewdec_rad) 

390 

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) 

397 

398 casalog.post("FIELD table PHASE_DIR column changed for field " + str(fld) + ".", 

399 'NORMAL') 

400 

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 

411 

412 pcol2 = tbt.getcol('DelayDir_Ref') 

413 theoldref2 = pcol2[fld] 

414 

415 pcol3 = tbt.getcol('RefDir_Ref') 

416 theoldref3 = pcol3[fld] 

417 

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') 

442 

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() 

452 

453 fldids = [] 

454 phdirs = [] 

455 for i in range(numfields): 

456 if (i==fld): 

457 fldids.append(i) 

458 phdirs.append(theoldphasecenter) 

459 

460 if thedistances==[]: 

461 thedistances = 0. # the default value 

462 

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 

468 

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') 

476 

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') 

483 

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) 

510 

511 dirstr = [refstr, _myqa.time(qnewra,12)[0], _myqa.angle(qnewdec,12)[0]] 

512 

513 elif not len(dirstr) == 3: 

514 msg = 'Incorrectly formatted parameter \'phasecenter\': ' + phasecenter 

515 raise RuntimeError(msg) 

516 return dirstr