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

247 statements  

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

1import os 

2import shutil 

3 

4from casatools import table, quanta 

5from casatools import ms as mstool 

6from casatasks import casalog 

7from .mstools import write_history 

8 

9qa = quanta() 

10 

11def cvel(vis, outputvis, 

12 passall, field, spw, selectdata, antenna, timerange, scan, array, 

13 mode, nchan, start, width, interpolation, 

14 phasecenter, restfreq, outframe, veltype, hanning): 

15 

16 """ regrid an MS to a new spectral window / channel structure or frame 

17  

18 vis -- Name of input visibility file 

19 default: none; example: vis='ngc5921.ms'  

20  

21 outputvis -- Name of output measurement set (required) 

22 default: none; example: vis='ngc5921-regridded.ms'  

23 

24 passall -- if False, data not meeting the selection on field or spw  

25 is omitted; if True, data not meeting the selection  

26 is passed through without modification 

27 default: False; example:  

28 field='NGC5921' 

29 passall=False : only data from NGC5921 is included in output MS,  

30 no data from other fields (e.g. 1331+305) is included 

31 passall=True : data from NGC5921 is transformed by cvel, all other  

32 fields are passed through unchanged  

33  

34 field -- Select fields in MS. Use field id(s) or field name(s). 

35 ['go listobs' to obtain the list id's or names] 

36 default: ''= all fields 

37 If field string is a non-negative integer, it is assumed to 

38 be a field index otherwise, it is assumed to be a  

39 field name 

40 field='0~2'; field ids 0,1,2 

41 field='0,4,5~7'; field ids 0,4,5,6,7 

42 field='3C286,3C295'; field named 3C286 and 3C295 

43 field = '3,4C*'; field id 3, all names starting with 4C 

44  

45 spw --Select spectral window/channels 

46 NOTE: This selects the data passed as the INPUT to mode 

47 default: ''=all spectral windows and channels 

48 spw='0~2,4'; spectral windows 0,1,2,4 (all channels) 

49 spw='0:5~61'; spw 0, channels 5 to 61 

50 spw='<2'; spectral windows less than 2 (i.e. 0,1) 

51 spw='0,10,3:3~45'; spw 0,10 all channels, spw 3,  

52 channels 3 to 45. 

53 spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each. 

54 spw='0:0~10;15~60'; spectral window 0 with channels  

55 0-10,15-60 

56 spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10, 

57 spw 1, channels 20-30, and spw 2, channels, 1,2 and 3 

58  

59 selectdata -- Other data selection parameters 

60 default: True 

61  

62 selectdata=True expandable parameters 

63  

64 antenna -- Select data based on antenna/baseline 

65 default: '' (all) 

66 If antenna string is a non-negative integer, it is  

67 assumed to be an antenna index, otherwise, it is 

68 considered an antenna name. 

69 antenna='5&6'; baseline between antenna index 5 and  

70 index 6. 

71 antenna='VA05&VA06'; baseline between VLA antenna 5  

72 and 6. 

73 antenna='5&6;7&8'; baselines 5-6 and 7-8 

74 antenna='5'; all baselines with antenna index 5 

75 antenna='05'; all baselines with antenna number 05  

76 (VLA old name) 

77 antenna='5,6,9'; all baselines with antennas 5,6,9  

78 index numbers 

79  

80 timerange -- Select data based on time range: 

81 default = '' (all); examples, 

82 timerange = 'YYYY/MM/DD/hh:mm:ss~YYYY/MM/DD/hh:mm:ss' 

83 Note: if YYYY/MM/DD is missing date defaults to first  

84 day in data set 

85 timerange='09:14:0~09:54:0' picks 40 min on first day 

86 timerange= '25:00:00~27:30:00' picks 1 hr to 3 hr  

87 30min on NEXT day 

88 timerange='09:44:00' pick data within one integration  

89 of time 

90 timerange='>10:24:00' data after this time 

91  

92 scan -- Scan number range. 

93 default: '' (all) 

94 example: scan='1~5' 

95 Check 'go listobs' to insure the scan numbers are in  

96 order. 

97  

98 array -- Select data by (sub)array indices 

99 default: '' (all); example: 

100 array='0~2'; arrays 0 to 2  

101  

102 mode -- Frequency Specification: 

103 NOTE: See examples below: 

104 default: 'channel' 

105 mode = 'channel'; Use with nchan, start, width to specify 

106 output spw. See examples below 

107 mode = 'velocity', means channels are specified in  

108 velocity. 

109 mode = 'frequency', means channels are specified in  

110 frequency. 

111  

112 mode expandable parameters  

113 Start, width are given in units of channels, frequency  

114 or velocity as indicated by mode  

115 nchan -- Number of channels in output spw 

116 default: -1 = all channels; example: nchan=3 

117 start -- Start input channel (relative-0) 

118 default=0; example: start=5 

119 width -- Output channel width in units of the input 

120 channel width (>1 indicates channel averaging) 

121 default=1; example: width=4 

122 interpolation -- Interpolation method 

123 default = 'linear' 

124 

125 examples: 

126 spw = '0,1'; mode = 'channel' 

127 will produce a single spw containing all channels in spw  

128 0 and 1 

129 spw='0:5~28^2'; mode = 'channel' 

130 will produce a single spw made with channels  

131 (5,7,9,...,25,27) 

132 spw = '0'; mode = 'channel': nchan=3; start=5; width=4 

133 will produce an spw with 3 output channels 

134 new channel 1 contains data from channels (5+6+7+8) 

135 new channel 2 contains data from channels (9+10+11+12) 

136 new channel 3 contains data from channels (13+14+15+16) 

137 spw = '0:0~63^3'; mode='channel'; nchan=21; start = 0;  

138 width = 1 

139 will produce an spw with 21 channels 

140 new channel 1 contains data from channel 0 

141 new channel 2 contains data from channel 2 

142 new channel 21 contains data from channel 61 

143 spw = '0:0~40^2'; mode = 'channel'; nchan = 3; start =  

144 5; width = 4 

145 will produce an spw with three output channels 

146 new channel 1 contains channels (5,7) 

147 new channel 2 contains channels (13,15) 

148 new channel 3 contains channels (21,23) 

149  

150 phasecenter -- direction measure or fieldid for the mosaic center 

151 default: '' => first field selected ; example: phasecenter=6 

152 or phasecenter='J2000 19h30m00 -40d00m00' 

153  

154 restfreq -- Specify rest frequency to use for output image 

155 default='' Occasionally it is necessary to set this (for 

156 example some VLA spectral line data). For example for 

157 NH_3 (1,1) put restfreq='23.694496GHz' 

158  

159 outframe -- output reference frame (not case-sensitive) 

160 possible values: LSRK, LSRD, BARY, GALACTO, LGROUP, CMB, GEO, TOPO or SOURCE 

161 default='' (keep original reference frame) ; example: outframe='BARY'  

162  

163 veltype -- definition of velocity (in mode) 

164 default = 'radio' 

165  

166 hanning -- if true, Hanning smooth frequency channel data before regridding 

167 to remove Gibbs ringing 

168 default = False 

169  

170 """ 

171 

172 # Note: this is duplicated in task_cvel, and really needing CASA-wide harmonization 

173 # (CAS-12871) 

174 def copy_ms(src, dest): 

175 """ This is a MMS-safe copy of an MMS tree directory. 

176 :param src: path to the source MS 

177 :param dest: path to the destination MS 

178 """ 

179 shutil.copytree(src, dest, symlinks=True) 

180 

181 #Python script 

182 

183 # make ms, table tools local 

184 _ms = mstool() 

185 _tb = table() 

186 

187 try: 

188 casalog.origin('cvel') 

189 

190 if not ((type(vis)==str) & (os.path.exists(vis))): 

191 raise ValueError('Visibility data set not found - please verify the name') 

192 

193 assert outputvis != '', "Must provide output data set name in parameter outputvis." 

194 assert not os.path.exists(outputvis), "Output MS %s already exists - will not overwrite." % outputvis 

195 assert not os.path.exists(outputvis+".flagversions"), \ 

196 "The flagversions \"%s.flagversions\" for the output MS already exist. Please delete." % outputvis 

197 

198 # Handle selectdata explicitly 

199 # (avoid hidden globals) 

200 if (selectdata==False): 

201 timerange='' 

202 array='' 

203 antenna='' 

204 scan='' 

205 

206 if(type(antenna) == list): 

207 antenna = ', '.join([str(ant) for ant in antenna]) 

208 

209 if (field == ''): 

210 field = '*' 

211 

212 if (spw == ''): 

213 spw = '*' 

214 

215 if(passall and spw=='*' and field=='*'): 

216 # all spws and fields selected, nothing to pass through 

217 passall = False 

218 

219 doselection = True 

220 if(field=='*' and spw=='*' and 

221 (not selectdata or (selectdata and antenna=='' and timerange=='' and scan=='' and array=='')) 

222 ): 

223 doselection = False 

224 

225 # open input MS 

226 _ms.open(vis) 

227 

228 if(type(phasecenter)==str): 

229 ### blank means take field 0 

230 if (phasecenter==''): 

231 phasecenter=_ms.msseltoindex(vis,field=field)['field'][0] 

232 else: 

233 tmppc=phasecenter 

234 try: 

235 if(len(_ms.msseltoindex(vis, field=phasecenter)['field']) > 0): 

236 tmppc=_ms.msseltoindex(vis, field=phasecenter)['field'][0] 

237 ##succesful must be string like '0' or 'NGC*' 

238 except Exception: 

239 ##failed must be a string 'J2000 18h00m00 10d00m00' 

240 tmppc=phasecenter 

241 phasecenter=tmppc 

242 

243 newphasecenter = phasecenter 

244 phasecentername = phasecenter 

245 if not (type(phasecenter)==str): 

246 # first check that this field will be in the output MS 

247 if not (phasecenter in _ms.msseltoindex(vis,field=field)['field']): 

248 message = "Field id " + str(phasecenter) 

249 message += " was selected as phasecenter but is not among the fields selected for output: " 

250 message += str(_ms.msseltoindex(vis,field=field)['field']) 

251 _ms.close() 

252 raise RuntimeError(message) 

253 

254 _tb.open(vis+"/FIELD") 

255 try: 

256 # get the name for later display 

257 phasecentername = _tb.getcell('NAME', phasecenter) + " (original field " + str(phasecenter) 

258 _tb.close() 

259 # phase center id will change if we select on fields, 

260 # the name column is only optionally filled and not necessarily unique. 

261 # But _ms.msseltoindex(vis,field=field)['field'] returns the old field ids 

262 # in the order in which they will occur in the new field table. 

263 # => need to get index from there as new phase center ID 

264 newphasecenter = (_ms.msseltoindex(vis,field=field)['field']).tolist().index(phasecenter) 

265 phasecentername += ", new field " + str(newphasecenter) + ")" 

266 except: 

267 _tb.close() 

268 message = "phasecenter field id " + str(phasecenter) + " is not valid" 

269 _ms.close() 

270 raise RuntimeError(message) 

271 

272 if(mode=='frequency'): 

273 ## reset the default values 

274 if(start==0): 

275 start = "" 

276 if(width==1): 

277 width = "" 

278 ##check that start and width have units if they are non-empty 

279 if not(start==""): 

280 if (qa.quantity(start)['unit'].find('Hz') < 0): 

281 _ms.close() 

282 raise TypeError("start parameter is not a valid frequency quantity " %start) 

283 if(not(width=="") and (qa.quantity(width)['unit'].find('Hz') < 0)): 

284 _ms.close() 

285 raise TypeError("width parameter %s is not a valid frequency quantity " %width) 

286 

287 elif(mode=='velocity'): 

288 ## reset the default values 

289 if(start==0): 

290 start = "" 

291 if(width==1): 

292 width = "" 

293 ##check that start and width have units if they are non-empty 

294 if not(start==""): 

295 if (qa.quantity(start)['unit'].find('m/s') < 0): 

296 _ms.close() 

297 raise TypeError("start parameter %s is not a valid velocity quantity " %start) 

298 if(not(width=="") and (qa.quantity(width)['unit'].find('m/s') < 0)): 

299 _ms.close() 

300 raise TypeError("width parameter %s is not a valid velocity quantity " %width) 

301 

302 elif(mode=='channel' or mode=='channel_b'): 

303 if((type(width) != int) or (type(start) != int)): 

304 _ms.close() 

305 raise TypeError("start and width have to be integers with mode = %s" %mode) 

306 

307 

308 ## check if preaveraging is necessary 

309 dopreaverage=False 

310 preavwidth = [1] 

311 

312 thespwsel = _ms.msseltoindex(vis=vis, spw=spw)['spw'] 

313 thefieldsel = _ms.msseltoindex(vis=vis, field=field)['field'] 

314 outgrid = _ms.cvelfreqs(spwids=thespwsel, fieldids=thefieldsel, 

315 mode=mode, nchan=nchan, start=start, width=width, 

316 phasec=newphasecenter, restfreq=restfreq, 

317 outframe=outframe, veltype=veltype, verbose=False) 

318 if(len(outgrid)>1): 

319 tmpavwidth = [] 

320 for thespw in thespwsel: 

321 outgridw1 = _ms.cvelfreqs(spwids=[thespw], fieldids=thefieldsel, 

322 mode='channel', nchan=-1, start=0, width=1, # native width 

323 phasec=newphasecenter, restfreq=restfreq, 

324 outframe=outframe, veltype=veltype, verbose=False) 

325 if(len(outgridw1)>1): 

326 widthratio = abs((outgrid[1]-outgrid[0])/(outgridw1[1]-outgridw1[0])) 

327 if(widthratio>=2.0): # do preaverage 

328 tmpavwidth.append( int(widthratio+0.001) ) 

329 dopreaverage=True 

330 else: 

331 tmpavwidth.append(1) 

332 

333 if dopreaverage: 

334 preavwidth = tmpavwidth 

335 

336 _ms.close() 

337 

338 # if in channel mode and preaveraging,  

339 if(dopreaverage and (mode=='channel' or mode=='channel_b')): 

340 if(max(preavwidth)==1): 

341 dopreaverage = False 

342 else: 

343 # switch to frequency mode 

344 mode = 'frequency' 

345 if(width>0): 

346 start = str(outgrid[0]/1E6)+'MHz' 

347 width = str((outgrid[1]-outgrid[0]))+'Hz' 

348 else: 

349 start = str(outgrid[len(outgrid)-1]/1E6)+'MHz' 

350 width = str(-(outgrid[1]-outgrid[0]))+'Hz' 

351 casalog.post("After pre-averaging, will switch to frequency mode with start="+start+", width = "+width, 'INFO') 

352 

353 if(dopreaverage and hanning and max(preavwidth)>2): 

354 casalog.post("NOTE: You have requested Hanning smoothing and at the same time you have chosen\n" 

355 +"a large width parameter which makes pre-averaging necessary.\n" 

356 +"The Hanning-smoothing may be redundant in this context.\n", 'WARN') 

357 

358 if dopreaverage: 

359 # Past this point we know we are going to 'dopreaverage' 

360 # CAS-9798 

361 raise RuntimeError('ERROR: cvel does not regrid properly for channel ' 

362 'widths > or = 2x the native channel width, please use ' 

363 'mstransform, clean, or tclean for larger regridding. ' 

364 'All earlier versions of CASA also have this issue.') 

365 

366 # determine parameter "datacolumn" 

367 _tb.open(vis) 

368 allcols = _tb.colnames() 

369 _tb.close() 

370 dpresent = ('DATA' in allcols) 

371 mpresent = ('MODEL_DATA' in allcols) 

372 cpresent = ('CORRECTED_DATA' in allcols) 

373 if (dpresent and cpresent): 

374 datacolumn = 'all' 

375 elif (dpresent and not cpresent): 

376 datacolumn = 'data' 

377 elif (cpresent and not dpresent): 

378 datacolumn = 'corrected_data' 

379 elif (mpresent and not cpresent and not dpresent): 

380 datacolumn = 'model_data' 

381 else: 

382 raise RuntimeError("Neither DATA nor CORRECTED_DATA nor MODEL_DATA column present. Cannot proceed.") 

383 

384 if(doselection and not dopreaverage): 

385 casalog.post("Creating selected SubMS ...", 'INFO') 

386 _ms.open(vis) 

387 _ms.split(outputms=outputvis, field=field, 

388 spw=spw, step=[1], 

389 baseline=antenna, subarray=array, 

390 timebin='-1s', time=timerange, 

391 whichcol=datacolumn, 

392 scan=scan, uvrange="") 

393 _ms.close() 

394 

395 elif(dopreaverage and not doselection): 

396 if(hanning): 

397 casalog.post("Creating working copy for Hanning-smoothing ...", 'INFO') 

398 shutil.rmtree(outputvis+'TMP',ignore_errors=True) 

399 copy_ms(vis, outputvis+'TMP') 

400 _ms.open(outputvis+'TMP', nomodify=False) 

401 _ms.hanningsmooth(datacolumn=datacolumn) 

402 _ms.close() 

403 hanning = False 

404 _ms.open(outputvis+'TMP') 

405 else: 

406 _ms.open(vis) 

407 

408 casalog.post("Creating preaveraged SubMS using widths "+str(preavwidth), 'INFO') 

409 _ms.split(outputms=outputvis, 

410 whichcol=datacolumn, step=preavwidth) 

411 _ms.close() 

412 

413 elif(doselection and dopreaverage): 

414 if(hanning): 

415 casalog.post("Creating selected working copy for Hanning-smoothing ...", 'INFO') 

416 shutil.rmtree(outputvis+'TMP',ignore_errors=True) 

417 _ms.open(vis) 

418 _ms.split(outputms=outputvis+'TMP', field=field, 

419 spw=spw, step=[1], 

420 baseline=antenna, subarray=array, 

421 timebin='-1s', time=timerange, 

422 whichcol=datacolumn, 

423 scan=scan, uvrange="") 

424 _ms.close() 

425 _ms.open(outputvis+'TMP', nomodify=False) 

426 _ms.hanningsmooth(datacolumn=datacolumn) 

427 _ms.close() 

428 hanning = False 

429 _ms.open(outputvis+'TMP') 

430 casalog.post("Creating preaveraged SubMS using widths "+str(preavwidth), 'INFO') 

431 _ms.split(outputms=outputvis, 

432 whichcol=datacolumn, step=preavwidth) 

433 _ms.close() 

434 else: 

435 casalog.post("Creating selected, preaveraged SubMS using widths "+str(preavwidth), 'INFO') 

436 _ms.open(vis) 

437 _ms.split(outputms=outputvis, field=field, 

438 spw=spw, step=preavwidth, 

439 baseline=antenna, subarray=array, 

440 timebin='-1s', time=timerange, 

441 whichcol=datacolumn, 

442 scan=scan, uvrange="") 

443 _ms.close() 

444 

445 else: 

446 # no selection or preaveraging necessary, just copy 

447 casalog.post("Creating working copy ...", 'INFO') 

448 shutil.rmtree(outputvis,ignore_errors=True) 

449 copy_ms(vis, outputvis) 

450 

451 # Combine and if necessary regrid it 

452 _ms.open(outputvis, nomodify=False) 

453 

454 message = "Using " + phasecentername + " as common direction for the output reference frame." 

455 casalog.post(message, 'INFO') 

456 

457 if not _ms.cvel(mode=mode, nchan=nchan, start=start, width=width, 

458 interp=interpolation, 

459 phasec=newphasecenter, restfreq=restfreq, 

460 outframe=outframe, veltype=veltype, hanning=hanning): 

461 _ms.close() 

462 raise RuntimeError("Error in regridding step ...") 

463 _ms.close() 

464 

465 # deal with the passall option 

466 temp_suffix = ".deselected" 

467 if (passall): 

468 # determine range of fields 

469 fieldsel = _ms.msseltoindex(vis=vis, field=field)['field'] 

470 _tb.open(vis+"/FIELD") 

471 nfields = _tb.nrows() 

472 _tb.close() 

473 fielddesel = "" 

474 for i in range(0,nfields): 

475 if not (i in fieldsel): 

476 if not (fielddesel == ""): 

477 fielddesel += "," 

478 fielddesel += str(i) 

479 

480 # determine range of SPWs 

481 spwsel = _ms.msseltoindex(vis=vis, spw=spw)['spw'] 

482 _tb.open(vis+"/SPECTRAL_WINDOW") 

483 nspws = _tb.nrows() 

484 _tb.close() 

485 spwdesel = "" 

486 for i in range(0,nspws): 

487 if not (i in spwsel): 

488 if not (spwdesel == ""): 

489 spwdesel += "," 

490 spwdesel += str(i) 

491 

492 if not (fielddesel == "" and spwdesel == ""): 

493 # split out the data not selected by the conditions on field and spw 

494 # from the original MS and join it to the output MS  

495 

496 # need to do this in two steps 

497 # I) field == "*" and deselected spws 

498 if not (spwdesel == ""): 

499 _ms.open(vis) 

500 casalog.post("Passing through data with", 'INFO') 

501 casalog.post(" spws: " + spwdesel, 'INFO') 

502 

503 _ms.split(outputms=outputvis+temp_suffix, field='*', 

504 spw=spwdesel, step=[1], 

505 baseline=antenna, subarray=array, 

506 timebin='-1s', time=timerange, 

507 whichcol=datacolumn, 

508 scan=scan, uvrange="") 

509 _ms.close() 

510 

511 # join with the deselected part 

512 _ms.open(outputvis, nomodify=False) 

513 rval = _ms.concatenate(msfile = outputvis+temp_suffix) 

514 _ms.close() 

515 shutil.rmtree(outputvis+temp_suffix) 

516 if not rval: 

517 raise RuntimeError("Error in attaching passed-through data ...") 

518 

519 # II) deselected fields and selected spws 

520 if not (fielddesel == ""): 

521 _ms.open(vis) 

522 casalog.post("Passing through data with", 'INFO') 

523 casalog.post(" fields: " + fielddesel, 'INFO') 

524 casalog.post(" spws: " + spw, 'INFO') 

525 

526 _ms.split(outputms=outputvis+temp_suffix, field=fielddesel, 

527 spw=spw, step=[1], 

528 baseline=antenna, subarray=array, 

529 timebin='-1s', time=timerange, 

530 whichcol=datacolumn, 

531 scan=scan, uvrange="") 

532 _ms.close() 

533 

534 # join with the deselected part 

535 _ms.open(outputvis, nomodify=False) 

536 rval = _ms.concatenate(msfile = outputvis+temp_suffix) 

537 _ms.close() 

538 shutil.rmtree(outputvis+temp_suffix) 

539 if not rval: 

540 raise RuntimeError("Error in attaching passed-through data ...") 

541 

542 # Write history to output MS 

543 try: 

544 param_names = cvel.__code__.co_varnames[:cvel.__code__.co_argcount] 

545 local_vars = locals() 

546 param_vals = [local_vars[p] for p in param_names] 

547 

548 write_history(mstool(), outputvis, 'cvel', param_names, 

549 param_vals, casalog) 

550 except Exception as instance: 

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

552 'WARN') 

553 

554 finally: 

555 # delete temp output (comment out for debugging) 

556 if os.path.exists(outputvis+".spwCombined"): 

557 casalog.post("Deleting temporary output files ...", 'INFO') 

558 shutil.rmtree(outputvis+".spwCombined")