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

333 statements  

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

1import shutil 

2import os 

3import string 

4 

5from .parallel.parallel_task_helper import ParallelTaskHelper 

6from casatools import ms as mstool 

7from casatools import table as tbtool 

8from casatools import imager as imtool 

9from casatools import measures, quanta 

10from casatasks import casalog 

11from .mstools import write_history 

12_qa = quanta( ) 

13_me = measures( ) 

14 

15def _checkinternalephemtab(vis, field): 

16 """ 

17 This function checks if there is an ephemeris table attached in the MS for 

18 the field selected. It will returned table names found under FIELD table. 

19 

20 """ 

21 import glob 

22 from casatools import table, ms 

23 _tb = table() 

24 _ms = ms() 

25 fids = _ms.msseltoindex(vis=vis,field=field)['field'] 

26 _tb.open(vis+'/FIELD') 

27 ephemnames = [] 

28 if 'EPHEMERIS_ID' in _tb.colnames(): 

29 for i in fids: 

30 ephemid = _tb.getcell('EPHEMERIS_ID',i) 

31 if ephemid != -1: 

32 ephemnames.append(glob.glob(f'{vis}/FIELD/EPHEM{ephemid}*/')[0]) 

33 _tb.close() 

34 return list(set(ephemnames)) 

35 

36def fixplanets(vis, field, fixuvw=False, direction='', refant=0, reftime='first'): 

37 """ 

38 Fix FIELD, SOURCE, and UVW for given fields based on given direction or pointing 

39 table information 

40 

41 This task's main purpose is to correct observations which were performed 

42 with correct pointing and correlation but for which incorrect direction 

43 information was entered in the FIELD and SOURCE table of the MS. 

44 If you actually want to change the phase center of the visibilties in an MS, 

45 you should use task fixvis. 

46 

47 Input Parameters 

48 vis -- Name of the input visibility set 

49  

50 field -- field selection string 

51  

52 fixuvw -- recalc uvw coordinates? (default: False) 

53 

54 direction -- if set, don't use pointing table but set direction to this value. 

55 The direction can either be given explicitly or as the path 

56 to a JPL Horizons ephemeris (for an example of the format, 

57 ephemerides/JPL-Horizons/ in the local data directory or "External Data" section  

58 of the current CASADocs). 

59 Alternatively, the ephemeris table can also be obtained using 

60 getephemtable task. 

61  

62 example: 'J2000 19h30m00 -40d00m00', default= '' (use pointing table) 

63 

64 refant -- if using pointing table information, use it from this antenna 

65 default: 0 (antenna id 0) 

66 examples: 'DV06' (antenna with name DV06) 

67 3 (antenna id 3) 

68 

69 reftime -- if using pointing table information, use it from this timestamp 

70 default: 'first' 

71 examples: 'median' will use the median timestamp for the given field 

72 using only the unflagged maintable rows  

73 '2012/07/11/08:41:32' will use the given timestamp (must be 

74 within the observaton time) 

75 """ 

76 

77 casalog.origin('fixplanets') 

78 

79 mst = mstool() 

80 tbt = tbtool() 

81 imt = None 

82 

83 try: 

84 fields = mst.msseltoindex(vis=vis,field=field)['field'] 

85 numfields = 0 

86 if(len(fields) == 0): 

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

88 return 

89 

90 tbt.open(vis+"/FIELD") 

91 oldrefcol = [] 

92 if('PhaseDir_Ref' in tbt.colnames()): 

93 oldrefcol = tbt.getcol('PhaseDir_Ref') 

94 tbt.close() 

95 

96 fixplanetstemp = 'fixplanetstemp-'+os.path.basename(vis) # temp file name 

97 fixplanetstemp2 = 'fixplanetstemp2-'+os.path.basename(vis) # temp file name 

98 

99 for fld in fields: 

100 thenewra_rad = 0. 

101 thenewdec_rad = 0. 

102 thenewref = -1 

103 thenewrefstr = '' 

104 theephemeris = '' 

105 

106 if(direction==''): # use information from the pointing table 

107 # find median timestamp for this field in the main table 

108 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

109 thetime = 0 

110 if(reftime.lower()=='median'): 

111 tbt.open(vis) 

112 tttb = tbt.query('FIELD_ID=='+str(fld)+' AND FLAG_ROW==False', 

113 name=fixplanetstemp, columns='TIME') 

114 tttb.close() 

115 tttb = None 

116 tbt.close() 

117 tbt.open(fixplanetstemp) 

118 if(tbt.nrows()>0): 

119 thetime = tbt.getcell('TIME',tbt.nrows()//2) 

120 casalog.post( "MEDIAN TIME "+str(thetime), 'NORMAL') 

121 tbt.close() 

122 else: 

123 casalog.post( "No pointing table rows for field "+field, 'NORMAL') 

124 tbt.close() 

125 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

126 return 

127 elif(reftime.lower()=='first'): 

128 tbt.open(vis) 

129 tttb = tbt.query('FIELD_ID=='+str(fld), name=fixplanetstemp, columns='TIME') 

130 tttb.close() 

131 tttb = None 

132 tbt.close() 

133 tbt.open(fixplanetstemp) 

134 if(tbt.nrows()>0): 

135 thetime = tbt.getcell('TIME',0) 

136 casalog.post( "FIRST TIME "+str(thetime), 'NORMAL') 

137 tbt.close() 

138 else: 

139 casalog.post( "No pointing table rows for field "+field, 'NORMAL') 

140 tbt.close() 

141 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

142 return 

143 else: 

144 try: 

145 myqa = _qa.quantity(reftime) 

146 thetime = _qa.convert(myqa,'s')['value'] 

147 except Exception as instance: 

148 raise TypeError("reftime parameter is not a valid date (e.g. YYYY/MM/DD/hh:mm:ss): %s" % reftime) 

149 tbt.open(vis) 

150 tttb = tbt.query('FIELD_ID=='+str(fld), name=fixplanetstemp, columns='TIME') 

151 tttb.close() 

152 tttb = None 

153 tbt.close() 

154 tbt.open(fixplanetstemp) 

155 if(tbt.nrows()>0): 

156 thefirsttime = tbt.getcell('TIME',0) 

157 thelasttime = tbt.getcell('TIME',tbt.nrows()-1) 

158 tbt.close() 

159 else: 

160 casalog.post( "No pointing table rows for field "+field, 'NORMAL') 

161 tbt.close() 

162 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

163 return 

164 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

165 if (thefirsttime<=thetime and thetime<=thelasttime): 

166 casalog.post( "GIVEN TIME "+reftime+" == "+str(thetime), 'NORMAL') 

167 else: 

168 casalog.post( "GIVEN TIME "+reftime+" == "+str(thetime)+" is not within the observation time (" 

169 +str(thefirsttime)+'s'+" - " 

170 +str(thelasttime)+'s'+")", 'SEVERE') 

171 raise TypeError 

172 shutil.rmtree(fixplanetstemp, ignore_errors=True) 

173 

174 # determine reference antenna 

175 antids = mst.msseltoindex(vis=vis,baseline=refant)['antenna1'] 

176 antid = -1 

177 if(len(antids) == 0): 

178 msg = "Antenna selection returned zero results." 

179 raise RuntimeError(msg) 

180 antid = int(antids[0]) 

181 casalog.post('Using antenna id '+str(antid)+' as reference antenna.', 'NORMAL') 

182 tbt.open(vis+'/ANTENNA') 

183 flgcol = tbt.getcol('FLAG_ROW') 

184 tbt.close() 

185 if(flgcol[antid]==True): 

186 msg = 'Antenna id '+str(antid)+' is flagged. Please select a different one.' 

187 raise RuntimeError(msg) 

188 

189 # get direction for the timestamp from pointing table 

190 shutil.rmtree(fixplanetstemp2, ignore_errors=True) 

191 tbt.open(vis+'/POINTING') 

192 ttb = tbt.query('TRACKING==True AND NEARABS(TIME,'+str(thetime)+',INTERVAL/2.) AND ANTENNA_ID==' 

193 +str(antid), 

194 name=fixplanetstemp2) 

195 nr = ttb.nrows() 

196 ttb.close() 

197 ttb = None 

198 if(nr==0): 

199 shutil.rmtree(fixplanetstemp2, ignore_errors=True) 

200 ttb2 = tbt.query('TRACKING==True AND NEARABS(TIME,'+str(thetime)+',3.) AND ANTENNA_ID==' 

201 +str(antid), # search within 3 seconds 

202 name=fixplanetstemp2) 

203 nr = ttb2.nrows() 

204 ttb2.close() 

205 ttb2 = None 

206 if(nr==0): 

207 casalog.post( "Cannot find any POINTING table rows for antenna "+str(antid) 

208 +" with TRACKING==True within 3 seconds of TIME "+str(thetime), 'NORMAL') 

209 casalog.post( "Will try without requiring TRACKING==True ...", 'NORMAL') 

210 shutil.rmtree(fixplanetstemp2, ignore_errors=True) 

211 ttb3 = tbt.query('NEARABS(TIME,'+str(thetime)+',INTERVAL/2.) AND ANTENNA_ID==' 

212 +str(antid), 

213 name=fixplanetstemp2) 

214 nr = ttb3.nrows() 

215 ttb3.close() 

216 ttb3 = None 

217 if(nr==0): 

218 shutil.rmtree(fixplanetstemp2, ignore_errors=True) 

219 ttb4 = tbt.query('NEARABS(TIME,'+str(thetime)+',3.) AND ANTENNA_ID==' 

220 +str(antid), # search within 3 seconds 

221 name=fixplanetstemp2) 

222 nr = ttb4.nrows() 

223 ttb4.close() 

224 ttb4 = None 

225 if(nr==0): 

226 tbt.close() 

227 msg = ("Cannot find any POINTING table rows for antenna "+str(antid) + 

228 " within 3 seconds of TIME "+str(thetime)) 

229 raise RuntimeError(msg) # give up 

230 tbt.close() 

231 tbt.open(fixplanetstemp2) 

232 thedir = tbt.getcell('DIRECTION',0) 

233 tbt.close() 

234 casalog.post( ' field id '+str(fld)+ ' AZ EL '+str(thedir[0])+" "+str(thedir[1]), 'NORMAL') 

235 thedirme = _me.direction(rf='AZELGEO',v0=_qa.quantity(thedir[0][0], 'rad'), 

236 v1=_qa.quantity(thedir[1][0],'rad')) 

237 # convert AZEL to J2000 coordinates 

238 _me.doframe(_me.epoch(rf='UTC', v0=_qa.quantity(thetime,'s'))) 

239 

240 tbt.open(vis+'/ANTENNA') 

241 thepos = tbt.getcell('POSITION',antid) 

242 theposref = tbt.getcolkeyword('POSITION', 'MEASINFO')['Ref'] 

243 theposunits = tbt.getcolkeyword('POSITION', 'QuantumUnits') 

244 tbt.close() 

245 casalog.post( "Ref. antenna position is "+str(thepos)+' ('+theposunits[0]+', '+theposunits[1]+', ' 

246 +theposunits[2]+')('+theposref+')', 'NORMAL') 

247 _me.doframe(_me.position(theposref, 

248 _qa.quantity(thepos[0],theposunits[0]), 

249 _qa.quantity(thepos[1],theposunits[1]), 

250 _qa.quantity(thepos[2],theposunits[2])) 

251 ) 

252 thedirmemod = _me.measure(v=thedirme, rf='J2000') 

253 # casalog.post(thedirmemod) 

254 thenewra_rad = thedirmemod['m0']['value'] 

255 thenewdec_rad = thedirmemod['m1']['value'] 

256 _me.done() 

257 shutil.rmtree(fixplanetstemp2, ignore_errors=True) 

258 

259 else: # direction is not an empty string, use this instead of the pointing table information 

260 if(type(direction)==str): 

261 dirstr = direction.split(' ') 

262 if len(dirstr)==1: # an ephemeris table was given 

263 if(os.path.exists(dirstr[0])): 

264 if os.path.isfile(dirstr[0]): # it is a file, i.e. not a CASA table 

265 msg = "*** Error when interpreting parameter \'direction\':\n A file is given. Use of the JPL-Horizons "+\ 

266 "MIME format file is deprecated." 

267 raise RuntimeError(msg) 

268 else: # not a file, assume it is a CASA table 

269 theephemeris = dirstr[0] 

270 # add a check if it is going to replace the existing table in the MS 

271 existingephemtab = _checkinternalephemtab(vis, field) 

272 if existingephemtab != []: 

273 casalog.post(f'Will replace existing ephemeris table {existingephemtab} in the MS with {direction}. '+\ 

274 'Ephemeris tables attached in the MS are assumed to be the ones used by the correlator and attaching' +\ 

275 ' a different ephemeris table may lead to scientifically wrong results.','WARN') 

276 

277 casalog.post('Will use ephemeris table '+theephemeris+' with offset (0,0)', 'NORMAL') 

278 

279 thenewra_rad = 0. 

280 thenewdec_rad = 0. 

281 

282 else: 

283 msg = "*** Error when interpreting parameter \'direction\':\n string is neither a direction nor a table." 

284 raise RuntimeError(msg) 

285 else: 

286 if len(dirstr)==2: # a direction without ref frame was given 

287 dirstr = ['J2000', dirstr[0], dirstr[1]] 

288 casalog.post('Assuming reference frame J2000 for parameter \'direction\'', 'NORMAL') 

289 elif not len(dirstr)==3: 

290 msg = 'Incorrect format in parameter \'direction\'' 

291 raise RuntimeError(msg) 

292 

293 # process direction and refframe 

294 try: 

295 thedir = _me.direction(dirstr[0], dirstr[1], dirstr[2]) 

296 thenewra_rad = thedir['m0']['value'] 

297 thenewdec_rad = thedir['m1']['value'] 

298 except Exception as instance: 

299 msg = "*** Error \'%s\' when interpreting parameter \'direction\': " % (instance) 

300 raise RuntimeError(msg) 

301 

302 try: 

303 tbt.open(vis+"/FIELD") 

304 thenewrefindex = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefTypes'].tolist().index(dirstr[0]) 

305 thenewref = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefCodes'][thenewrefindex] 

306 thenewrefstr = dirstr[0] 

307 tbt.close() 

308 except Exception as instance: 

309 casalog.post("PHASE_DIR is not a variable reference column. Will leave reference as is.", 'WARN') 

310 thenewref = -1 

311 if(0<thenewref and thenewref<32): 

312 msg = "*** Error when interpreting parameter \'direction\':\n presently only J2000 and solar system objects are supported." 

313 raise RuntimeError(msg) 

314 #endif 

315 

316 # modify FIELD table 

317 tbt.open(vis+'/FIELD', nomodify=False) 

318 numfields = tbt.nrows() 

319 theolddir = tbt.getcell('PHASE_DIR',fld) 

320 planetname = tbt.getcell('NAME',fld) 

321 

322 casalog.post( 'object: '+planetname, 'NORMAL') 

323 casalog.post( 'old RA, DEC (rad) '+str(theolddir[0])+" " +str(theolddir[1]), 'NORMAL') 

324 if(theephemeris==''): 

325 casalog.post( 'new RA, DEC (rad) '+str(thenewra_rad)+" "+ str(thenewdec_rad), 'NORMAL') 

326 else: 

327 casalog.post( 'new RA, DEC to be taken from ephemeris table '+theephemeris 

328 +' with offsets (rad) '+str(thenewra_rad)+" "+ str(thenewdec_rad), 'NORMAL') 

329 

330 pcol = tbt.getcol('PHASE_DIR') 

331 pcol[0][0][fld] = thenewra_rad 

332 pcol[1][0][fld] = thenewdec_rad 

333 tbt.putcol('PHASE_DIR',pcol) 

334 

335 pcol = tbt.getcol('DELAY_DIR') 

336 pcol[0][0][fld] = thenewra_rad 

337 pcol[1][0][fld] = thenewdec_rad 

338 tbt.putcol('DELAY_DIR',pcol) 

339 

340 pcol = tbt.getcol('REFERENCE_DIR') 

341 pcol[0][0][fld] = thenewra_rad 

342 pcol[1][0][fld] = thenewdec_rad 

343 tbt.putcol('REFERENCE_DIR',pcol) 

344 

345 casalog.post("FIELD table PHASE_DIR, DELAY_DIR, and REFERENCE_DIR columns changed for field "+str(fld)+".", 'NORMAL') 

346 

347 if(theephemeris==''): 

348 

349 if ('EPHEMERIS_ID' in tbt.colnames()) and (tbt.getcell('EPHEMERIS_ID',fld)>=0): # originally an ephemeris was used 

350 eidc = tbt.getcol('EPHEMERIS_ID') 

351 eidc[fld] = -1 

352 tbt.putcol('EPHEMERIS_ID', eidc) # remove the reference to it 

353 casalog.post("FIELD table EPHEMERIS_ID column reset to -1 for field "+str(fld)+".", 'NORMAL') 

354 

355 if(thenewref!=-1): 

356 # modify reference of the direction columns permanently 

357 theoldref = -1 

358 theoldref2 = -1 

359 theoldref3 = -1 

360 try: 

361 pcol = tbt.getcol('PhaseDir_Ref') 

362 theoldref = pcol[fld] 

363 pcol[fld] = thenewref 

364 oldrefcol = pcol 

365 

366 pcol2 = tbt.getcol('DelayDir_Ref') 

367 theoldref2 = pcol2[fld] 

368 pcol2[fld] = thenewref 

369 

370 pcol3 = tbt.getcol('RefDir_Ref') 

371 theoldref3 = pcol3[fld] 

372 pcol3[fld] = thenewref 

373 

374 if not (theoldref==theoldref2 and theoldref==theoldref3): 

375 casalog.post("The three FIELD table direction reference frame entries for field "+str(fld) 

376 +" are not identical in the input data: " 

377 +str(theoldref)+", "+str(theoldref2)+", "+str(theoldref3) 

378 +". Will try to continue ...", 'WARN') 

379 theoldref=-1 

380 else: 

381 oldrefstrindex = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefCodes'].tolist().index(theoldref) 

382 oldrefstr = tbt.getcolkeyword('PHASE_DIR', 'MEASINFO')['TabRefTypes'][oldrefstrindex] 

383 casalog.post("Original FIELD table direction reference frame entries for field "+str(fld) 

384 +" are all "+str(theoldref)+" ("+oldrefstr+")", 'NORMAL') 

385 tbt.putcol('PhaseDir_Ref', pcol) 

386 tbt.putcol('DelayDir_Ref', pcol2) 

387 tbt.putcol('RefDir_Ref', pcol3) 

388 

389 except Exception as instance: 

390 msg = "*** Error \'%s\' when writing reference frames in FIELD table: " % (instance) 

391 raise RuntimeError(msg) 

392 

393 if(theoldref != thenewref): 

394 casalog.post("FIELD table direction reference frame entries for field "+str(fld) 

395 +" set to "+str(thenewref)+" ("+thenewrefstr+")", 'NORMAL') 

396 else: 

397 casalog.post("FIELD table direction reference frame entries for field "+str(fld) 

398 +" unchanged.", 'NORMAL') 

399 else: # we are working with an ephemeris 

400 try: 

401 mst.open(vis, nomodify=False) 

402 # sanitize the internal ephemeris table name 

403 valid_chars = "-+_.%s%s" % (string.ascii_letters, string.digits) 

404 planetname_for_table = ''.join(c if c in valid_chars else '_' for c in planetname) 

405 # make sure theemepheris is given as absolute path to avoid MeasIERS warnings 

406 mst.addephemeris(-1, os.path.abspath(theephemeris), planetname_for_table, fld) # -1 = take the next free ID 

407 except Exception as instance: 

408 raise RuntimeError("*** Error \'%s\' when attaching ephemeris: " % 

409 (instance)) 

410 finally: 

411 mst.close() 

412 

413 

414 if(fixuvw and (len(oldrefcol)!=0) and (thenewref>0)): 

415 # modify reference of phase dir for fixuvw 

416 pcol = tbt.getcol('PhaseDir_Ref') 

417 pcol[fld] = 0 # J2000 

418 tbt.putcol('PhaseDir_Ref', pcol) # J2000 

419 

420 tbt.close() 

421 

422 #modify SOURCE table 

423 newsra_rad = thenewra_rad 

424 newsdec_rad = thenewdec_rad 

425 if(theephemeris!=''): # get the nominal position from the ephemeris 

426 mst.open(vis) 

427 trec = mst.getfielddirmeas('PHASE_DIR',fld) 

428 newsra_rad = trec['m0']['value'] 

429 newsdec_rad = trec['m1']['value'] 

430 mst.close() 

431 tbt.open(vis+'/SOURCE', nomodify=False) 

432 sdir = tbt.getcol('DIRECTION') 

433 newsdir = sdir 

434 sname = tbt.getcol('NAME') 

435 

436 for i in range(0,tbt.nrows()): 

437 if(sname[i]==planetname): 

438 #casalog.post('i old dir ' + i + " " + sdir[0][i] + sdir[1][i]) 

439 newsdir[0][i] = newsra_rad 

440 newsdir[1][i] = newsdec_rad 

441 #casalog.post(' new dir ' + newsdir[0][i] + newsdir[1][i]) 

442 tbt.putcol('DIRECTION', newsdir) 

443 tbt.close() 

444 casalog.post("SOURCE table DIRECTION column changed.", 'NORMAL') 

445 

446 # end for 

447 

448 if(fixuvw): 

449 casalog.post("Fixing the UVW coordinates ...", 'NORMAL') 

450 

451 # similar to fixvis 

452 

453 fldids = [] 

454 for i in range(numfields): 

455 if (i in fields): 

456 fldids.append(i) 

457 

458 imt = imtool() 

459 imt.open(vis, usescratch=False) 

460 imt.calcuvw(fldids, refcode='J2000', reuse=False) 

461 imt.close() 

462 imt = None 

463 

464 if((len(oldrefcol)!=0) and (thenewref>0)): 

465 tbt.open(vis+'/FIELD', nomodify=False) 

466 tbt.putcol('PhaseDir_Ref', oldrefcol) 

467 tbt.close() 

468 

469 else: 

470 casalog.post("UVW coordinates not changed.", 'NORMAL') 

471 

472 if (ParallelTaskHelper.isParallelMS(vis)): 

473 casalog.post("Tidying up the MMS subtables ...", 'NORMAL') 

474 ParallelTaskHelper.restoreSubtableAgreement(vis) 

475 

476 finally: 

477 mst = None 

478 tbt = None 

479 imt = None 

480 

481 # Write history to MS 

482 try: 

483 param_names = fixplanets.__code__.co_varnames[:fixplanets.__code__.co_argcount] 

484 vars = locals() 

485 param_vals = [vars[p] for p in param_names] 

486 write_history(mstool(), vis, 'fixplanets', param_names, 

487 param_vals, casalog) 

488 except Exception as instance: 

489 casalog.post("*** Error \'%s\' updating HISTORY" % (instance), 

490 'WARN') 

491 

492 mst = None 

493 tbt = None