Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/task_tsdimaging.py: 9%

1# sd task for imaging 

2 

3import collections 

4import contextlib 

5import os 

6import re 

7import shutil 

8import time 

9 

10import numpy 

11 

12from functools import partial 

13 

14from casatasks import casalog 

15from casatools import image, imager 

16from casatools import ms as mstool 

17from casatools import quanta 

18 

19from . import mslisthelper, sdbeamutil, sdutil 

20# (1) Import the python application layer 

21from .imagerhelpers.imager_base import PySynthesisImager 

22from .imagerhelpers.input_parameters import ImagerParameters 

23 

24image_suffix = '.image' 

25residual_suffix = '.residual' 

26weight_suffix = '.weight' 

27associate_suffixes = ['.psf', '.sumwt', weight_suffix, residual_suffix] 

28 

29 

30@contextlib.contextmanager 

31def open_ia(imagename): 

32 ia = image() 

33 ia.open(imagename) 

34 try: 

35 yield ia 

36 finally: 

37 ia.close() 

38 

39 

40@contextlib.contextmanager 

41def open_ms(vis): 

42 ms = mstool() 

43 ms.open(vis) 

44 try: 

45 yield ms 

46 finally: 

47 ms.close() 

48 

49 

50class SelectionHandler(object): 

51 def __init__(self, sel): 

52 self.sel = sel 

53 if isinstance(self.sel, str): 

54 self.selector = self._select0 

55 elif len(self.sel) == 1: 

56 self.selector = self._select1 

57 else: 

58 self.selector = self._select2 

59 

60 def __call__(self, i): 

61 return self.selector(i) 

62 

63 def _select0(self, i): 

64 return self.sel 

65 

66 def _select1(self, i): 

67 return self.sel[0] 

68 

69 def _select2(self, i): 

70 return self.sel[i] 

71 

72 

73class OldImagerBasedTools(object): 

74 def __init__(self): 

75 self.imager = imager() 

76 

77 @contextlib.contextmanager 

78 def open_old_imager(self, vis): 

79 try: 

80 self.imager.open(vis) 

81 yield self.imager 

82 finally: 

83 self.imager.close() 

84 

85 @contextlib.contextmanager 

86 def open_and_select_old_imager(self, vislist, field, spw, antenna, scan, intent, timerange): 

87 if isinstance(vislist, str): 

88 with self.open_old_imager(vislist) as im: 

89 im.selectvis(field=field, 

90 spw=spw, 

91 nchan=-1, 

92 start=0, 

93 step=1, 

94 baseline=antenna, 

95 scan=scan, 

96 intent=intent, 

97 time=timerange) 

98 yield im 

99 else: 

100 fieldsel = SelectionHandler(field) 

101 spwsel = SelectionHandler(spw) 

102 antennasel = SelectionHandler(antenna) 

103 scansel = SelectionHandler(scan) 

104 intentsel = SelectionHandler(intent) 

105 timerangesel = SelectionHandler(timerange) 

106 try: 

107 for i in range(len(vislist)): 

108 vis = vislist[i] 

109 _field = fieldsel(i) 

110 _spw = spwsel(i) 

111 _antenna = antennasel(i) 

112 _scan = scansel(i) 

113 _intent = intentsel(i) 

114 _timerangesel = timerangesel(i) 

115 if len(_antenna) == 0: 

116 _baseline = _antenna 

117 elif len(_antenna) < 4 or _antenna[:-3] != '&&&': 

118 _baseline = _antenna + '&&&' 

119 else: 

120 _baseline = _antenna 

121 self.imager.selectvis( 

122 vis, field=_field, spw=_spw, nchan=-1, start=0, step=1, 

123 baseline=_baseline, scan=_scan, intent=_intent, time=_timerangesel) 

124 yield self.imager 

125 finally: 

126 self.imager.close() 

127 

128 def test(self, vis): 

129 with self.open_old_imager(vis): 

130 casalog.post('test') 

131 raise RuntimeError('ERROR!') 

132 

133 def get_pointing_sampling_params(self, vis, field, spw, baseline, scan, intent, timerange, 

134 outref, movingsource, pointingcolumntouse, antenna_name): 

135 with self.open_old_imager(vis) as im: 

136 im.selectvis(field=field, 

137 spw=spw, 

138 nchan=-1, 

139 start=0, 

140 step=1, 

141 baseline=baseline, 

142 scan=scan, 

143 intent=intent, 

144 time=timerange) 

145 sampling_params = im.pointingsampling(pattern='raster', 

146 ref=outref, 

147 movingsource=movingsource, 

148 pointingcolumntouse=pointingcolumntouse, 

149 antenna='{0}&&&'.format(antenna_name)) 

150 return sampling_params 

151 

152 def get_map_extent(self, vislist, field, spw, antenna, scan, intent, timerange, 

153 ref, movingsource, pointingcolumntouse): 

154 

155 with self.open_and_select_old_imager(vislist=vislist, field=field, 

156 spw=spw, antenna=antenna, scan=scan, 

157 intent=intent, timerange=timerange) as im: 

158 map_param = im.mapextent(ref=ref, movingsource=movingsource, 

159 pointingcolumntouse=pointingcolumntouse) 

160 return map_param 

161 

162 

163def check_conformance(mslist, check_result): 

164 """Check conformance of input MS list 

165 

166 Check conformance of input MS list, particularlly existence of 

167 WEIGHT_SPECTRUM column. 

168 

169 Args: 

170 mslist (list): list of names for input MS 

171 check_result (dict): result of conformance check. 

172 see mslisthelper.check_mslist for detail about 

173 the structure of check_result. 

174 

175 Returns: 

176 dict: Per-column set of names for MS that needs to be 

177 edited to resolve the conformance. Top level dict 

178 has two keys, "remove" and "add", which indicate 

179 the operation to be applied to the columns. 

180 """ 

181 process_dict = { 

182 'remove': collections.defaultdict(set), 

183 'add': collections.defaultdict(set) 

184 } 

185 for name, summary in check_result.items(): 

186 if 'Main' in summary: 

187 missingcol_list = [summary['Main']['missingcol_{}'.format(x)] for x in ['a', 'b']] 

188 ms_a = mslist[0] 

189 ms_b = name 

190 

191 # "remove" operation: 

192 # - MS list is opposite order to missingcol_list 

193 # - if WEIGHT_SPECTRUM column is missing in one MS, 

194 # the column should be removed from another MS. 

195 for c, t in zip(missingcol_list, [ms_b, ms_a]): 

196 if 'WEIGHT_SPECTRUM' in c: 

197 process_dict['remove']['WEIGHT_SPECTRUM'].add(t) 

198 # "add" operation: 

199 # - MS list is same order as missingcol_list 

200 # - if CORRECTED_DATA column is missing in one MS, 

201 # the column should be added to that MS. 

202 for c, t in zip(missingcol_list, [ms_a, ms_b]): 

203 if 'CORRECTED_DATA' in c: 

204 process_dict['add']['CORRECTED_DATA'].add(t) 

205 return process_dict 

206 

207 

208def report_conformance(mslist, column_name, process_set): 

209 """Report conformance of input MS 

210 

211 Report conformance of input MS, particularlly on the existence 

212 of the column given by column_name. 

213 

214 Args: 

215 mslist (list): list of names for input MS 

216 column_name (str): name of the column 

217 process_set (set): set of names of MS that need to be edited 

218 """ 

219 if len(process_set) > 0: 

220 casalog.post('Detected non-conformance of {} column in input list of MSes.'.format(column_name), priority='WARN') 

221 cols = ['exists?', 'MS name'] 

222 header = ' '.join(cols) 

223 casalog.post('', priority='WARN') 

224 casalog.post('Summary of existence of {}:'.format(column_name), priority='WARN') 

225 casalog.post(header, priority='WARN') 

226 casalog.post('-' * len(header), priority='WARN') 

227 for name in mslist: 

228 basename = os.path.basename(name.rstrip('/')) 

229 exists = 'YES' if name in process_set else 'NO' 

230 row = '{:^7s} {:<s}'.format(exists, basename) 

231 casalog.post(row, priority='WARN') 

232 

233 

234def fix_conformance(process_dict): 

235 """Resolve non-conformance by removing WEIGHT_SPECTRUM 

236 

237 Two non-conformances are fixed, WEIGHT_SPECTRUM and 

238 CORRECTED_DATA. Remove WEIGHT_SPECTRUM column from, or 

239 add CORRECTED_DATA to the MS provided by process_dict. 

240 Backup is created with the name: 

241 

242 <original_name>.sdimaging.backup-<timestamp> 

243 

244 Args: 

245 process_dict (dict): per-operation ("remove" and "add") 

246 key-value pair of column name and 

247 list of names for MS to be edited 

248 

249 Returns: 

250 dict: mapping of original MS name and the name of backup 

251 """ 

252 backup_list = {} 

253 process_list = set() 

254 for v in process_dict.values(): 

255 for w in v.values(): 

256 process_list = process_list.union(w) 

257 for name in process_list: 

258 basename = os.path.basename(name.rstrip('/')) 

259 timestamp = time.strftime('%Y%m%dT%H%M%S', time.gmtime()) 

260 backup_name = basename + '.sdimaging.backup-{}'.format(timestamp) 

261 with sdutil.table_manager(name) as tb: 

262 tb.copy(backup_name, deep=True, returnobject=True).close() 

263 backup_list[name] = backup_name 

264 casalog.post('Copy of "{}" has been saved to "{}"'.format(name, backup_name), priority='WARN') 

265 

266 for colname, msnames in process_dict['remove'].items(): 

267 for name in msnames: 

268 casalog.post('{} will be removed from "{}"'.format(colname, name), priority='WARN') 

269 with sdutil.table_manager(name, nomodify=False) as tb: 

270 if colname in tb.colnames(): 

271 tb.removecols(colname) 

272 

273 for colname, msnames in process_dict['add'].items(): 

274 for name in msnames: 

275 casalog.post('{} will be added to "{}"'.format(colname, name), priority='WARN') 

276 with sdutil.calibrater_manager(name, addmodel=False, addcorr=True): 

277 pass 

278 return backup_list 

279 

280 

281def conform_mslist(mslist, ignore_columns=['CORRECTED_DATA']): 

282 """Make given set of MS data conform 

283 

284 Here, only conformance on the existence of WEIGHT_SPECTRUM 

285 is checked and resolved because non-conformance of WEIGHT_SPECTRUM, 

286 i.e. some MS have the column while others don't, could cause 

287 the task to crash. If non-conformance is detected, all existing 

288 WEIGHT_SPECTRUM columns are removed. This opration modifies input 

289 MS so data will be backed up with the name: 

290 

291 <original_name>.sdimaging.backup-<timestamp> 

292 

293 Args: 

294 mslist (list): list of names for input MS 

295 """ 

296 check_result = mslisthelper.check_mslist(mslist, testcontent=False) 

297 process_dict = check_conformance(mslist, check_result) 

298 fix_dict = dict() 

299 for op, op_dict in process_dict.items(): 

300 fix_dict[op] = dict() 

301 for col, process_set in op_dict.items(): 

302 if col not in ignore_columns: 

303 report_conformance(mslist, col, process_set) 

304 fix_dict[op][col] = process_set 

305 fix_conformance(fix_dict) 

306 

307 

308def sort_vis(vislist, spw, mode, width, field, antenna, scan, intent, timerange): 

309 """Sort the given MeasurementSet path(s) by their earliest data-taking time, in increasing order. 

310 

311 Return a 7-tuple where 

312 * the first entry holds the re-ordered paths 

313 * the remaining entries hold their corresponding data selection parameters 

314 * FIXME: input parameters mode and width are not used 

315 """ 

316 if isinstance(vislist, str) or len(vislist) == 1: 

317 return vislist, field, spw, antenna, scan, intent, timerange 

318 # chronological sort 

319 sorted_vislist, sorted_timelist = mslisthelper.sort_mslist(vislist) 

320 _vislist = list(vislist) 

321 sorted_idx = [_vislist.index(vis) for vis in sorted_vislist] 

322 mslisthelper.report_sort_result(sorted_vislist, sorted_timelist, sorted_idx, mycasalog=casalog) 

323 # conform MS 

324 conform_mslist(sorted_vislist) 

325 fieldsel = SelectionHandler(field) 

326 sorted_field = [fieldsel(i) for i in sorted_idx] 

327 spwsel = SelectionHandler(spw) 

328 sorted_spw = [spwsel(i) for i in sorted_idx] 

329 antennasel = SelectionHandler(antenna) 

330 sorted_antenna = [antennasel(i) for i in sorted_idx] 

331 scansel = SelectionHandler(scan) 

332 sorted_scan = [scansel(i) for i in sorted_idx] 

333 intentsel = SelectionHandler(intent) 

334 sorted_intent = [intentsel(i) for i in sorted_idx] 

335 timerangesel = SelectionHandler(timerange) 

336 sorted_timerange = [timerangesel(i) for i in sorted_idx] 

337 return sorted_vislist, sorted_field, sorted_spw, sorted_antenna, sorted_scan, sorted_intent, sorted_timerange 

338 

339 

340def _configure_spectral_axis(mode, nchan, start, width, restfreq): 

341 # fix default 

342 if mode == 'channel': 

343 if start == '': 

344 start = 0 

345 if width == '': 

346 width = 1 

347 else: 

348 if start == 0: 

349 start = '' 

350 if width == 1: 

351 width = '' 

352 # fix unit 

353 if mode == 'frequency': 

354 myunit = 'Hz' 

355 elif mode == 'velocity': 

356 myunit = 'km/s' 

357 else: # channel 

358 myunit = '' 

359 

360 tmp_start = _format_quantum_unit(start, myunit) 

361 if tmp_start is None: 

362 raise ValueError("Invalid unit for %s in mode %s: %s" % ('start', mode, start)) 

363 start = tmp_start 

364 if mode == 'channel': 

365 start = int(start) 

366 tmp_width = _format_quantum_unit(width, myunit) 

367 if tmp_width is None: 

368 raise ValueError("Invalid unit for %s in mode %s: %s" % ('width', mode, width)) 

369 width = tmp_width 

370 if mode == 'channel': 

371 width = int(width) 

372 

373 # TODO: work for nchan 

374 imnchan = nchan 

375 imstart = start 

376 imwidth = width 

377 return imnchan, imstart, imwidth 

378 

379 

380def _format_quantum_unit(data, unit): 

381 """Format quantity data. 

382 

383 Returns False if data has an unit which in not a variation of 

384 input unit. 

385 Otherwise, returns input data as a quantum string. The input 

386 unit is added to the return value if no unit is in data. 

387 """ 

388 my_qa = quanta() 

389 if data == '' or my_qa.compare(data, unit): 

390 return data 

391 if my_qa.getunit(data) == '': 

392 casalog.post("No unit specified. Using '%s'" % unit) 

393 return '%f%s' % (data, unit) 

394 return None 

395 

396 

397def _handle_grid_defaults(value): 

398 ret = '' 

399 if isinstance(value, int) or isinstance(value, float): 

400 ret = str(value) 

401 elif isinstance(value, str): 

402 ret = value 

403 return ret 

404 

405 

406def _calc_PB(vis, antenna_id, restfreq): 

407 """Calculate the primary beam size of antenna. 

408 

409 Calculate the primary beam size of antenna, using dish diamenter 

410 and rest frequency 

411 Average antenna diamter and reference frequency are adopted for 

412 calculation. 

413 The input argument should be a list of antenna IDs. 

414 """ 

415 logger = sdutil.Casalog(origin="_calc_PB") 

416 logger.post("Calculating Primary beam size:") 

417 # CAS-5410 Use private tools inside task scripts 

418 my_qa = quanta() 

419 

420 pb_factor = 1.175 

421 # Reference frequency 

422 ref_freq = restfreq 

423 if type(ref_freq) in [float, numpy.float64]: 

424 ref_freq = my_qa.tos(my_qa.quantity(ref_freq, 'Hz')) 

425 if not my_qa.compare(ref_freq, 'Hz'): 

426 msg = "Could not get the reference frequency. " + \ 

427 "Your data does not seem to have valid one in selected field.\n" + \ 

428 "PB is not calculated.\n" + \ 

429 "Please set restreq or cell manually to generate an image." 

430 raise RuntimeError(msg) 

431 # Antenna diameter 

432 with sdutil.table_manager(os.path.join(vis, 'ANTENNA')) as tb: 

433 antdiam_ave = tb.getcell('DISH_DIAMETER', antenna_id) 

434 # antdiam_ave = self._get_average_antenna_diameter(antenna) 

435 # Calculate PB 

436 wave_length = 0.2997924 / my_qa.convert(my_qa.quantity(ref_freq), 'GHz')['value'] 

437 D_m = my_qa.convert(antdiam_ave, 'm')['value'] 

438 lambda_D = wave_length / D_m * 3600. * 180 / numpy.pi 

439 PB = my_qa.quantity(pb_factor * lambda_D, 'arcsec') 

440 # Summary 

441 logger.post(f"- Antenna diameter: {D_m} m") 

442 logger.post(f"- Reference Frequency: {ref_freq}") 

443 logger.post(f"PB size = {pb_factor:5.3f} * lambda/D = {my_qa.tos(PB)}") 

444 return PB 

445 

446 

447def _get_imsize(width, height, dx, dy): 

448 casalog.post("Calculating pixel size.") 

449 # CAS-5410 Use private tools inside task scripts 

450 my_qa = quanta() 

451 ny = numpy.ceil((my_qa.convert(height, my_qa.getunit(dy))['value'] / 

452 my_qa.getvalue(dy))) 

453 nx = numpy.ceil((my_qa.convert(width, my_qa.getunit(dx))['value'] / 

454 my_qa.getvalue(dx))) 

455 casalog.post("- Map extent: [%s, %s]" % (my_qa.tos(width), my_qa.tos(height))) 

456 casalog.post("- Cell size: [%s, %s]" % (my_qa.tos(dx), my_qa.tos(dy))) 

457 casalog.post("Image pixel numbers to cover the extent: [%d, %d] (projected)" % 

458 (nx + 1, ny + 1)) 

459 return [int(nx + 1), int(ny + 1)] 

460 

461 

462def _get_pointing_extent(phasecenter, vislist, field, spw, antenna, scan, intent, timerange, 

463 pointingcolumntouse, ephemsrcname): 

464 # MS selection is ignored. This is not quite right. 

465 casalog.post("Calculating map extent from pointings.") 

466 # CAS-5410 Use private tools inside task scripts 

467 my_qa = quanta() 

468 ret_dict = {} 

469 

470 if isinstance(vislist, str): 

471 vis = vislist 

472 else: 

473 vis = vislist[0] 

474 

475 # colname = pointingcolumntouse.upper() 

476 

477 if phasecenter == "": 

478 # defaut is J2000 

479 base_mref = 'J2000' 

480 elif isinstance(phasecenter, int) or phasecenter.isdigit(): 

481 # may be field id 

482 with sdutil.table_manager(os.path.join(vis, 'FIELD')) as tb: 

483 base_mref = tb.getcolkeyword('PHASE_DIR', 'MEASINFO')['Ref'] 

484 else: 

485 # may be phasecenter is explicitly specified 

486 # numeric value: 3.14, -.3e1, etc. 

487 numeric_pattern = r'[-+]?([0-9]+(.[0-9]*)?|\.[0-9]+)([eE]-?[0-9])?' 

488 # HMS string: 9:15:29, -9h15m29 

489 hms_pattern = r'[-+]?[0-9]+[:h][0-9]+[:m][0-9.]+s?' 

490 # DMS string: 9.15.29, -9d15m29s 

491 dms_pattern = r'[-+]?[0-9]+[.d][0-9]+[.m][0-9.]+s?' 

492 # composite pattern 

493 pattern = fr'^({numeric_pattern}|{hms_pattern}|{dms_pattern})$' 

494 items = phasecenter.split() 

495 base_mref = 'J2000' 

496 for i in items: 

497 s = i.strip() 

498 if re.match(pattern, s) is None: 

499 base_mref = s 

500 break 

501 

502 t = OldImagerBasedTools() 

503 mapextent = t.get_map_extent(vislist, field, spw, antenna, scan, intent, timerange, 

504 ref=base_mref, movingsource=ephemsrcname, 

505 pointingcolumntouse=pointingcolumntouse) 

506 # mapextent = self.imager.mapextent(ref=base_mref, movingsource=ephemsrcname, 

507 # pointingcolumntouse=colname) 

508 if mapextent['status']: 

509 qheight = my_qa.quantity(mapextent['extent'][1], 'rad') 

510 qwidth = my_qa.quantity(mapextent['extent'][0], 'rad') 

511 qcent0 = my_qa.quantity(mapextent['center'][0], 'rad') 

512 qcent1 = my_qa.quantity(mapextent['center'][1], 'rad') 

513 scenter = '%s %s %s' % (base_mref, my_qa.formxxx(qcent0, 'hms'), 

514 my_qa.formxxx(qcent1, 'dms')) 

515 

516 casalog.post("- Pointing center: %s" % scenter) 

517 casalog.post("- Pointing extent: [%s, %s] (projected)" % (my_qa.tos(qwidth), 

518 my_qa.tos(qheight))) 

519 ret_dict['center'] = scenter 

520 ret_dict['width'] = qwidth 

521 ret_dict['height'] = qheight 

522 else: 

523 casalog.post( 

524 'Failed to derive map extent from the MSs registered to the imager probably ' 

525 'due to mising valid data.', 

526 priority='SEVERE') 

527 ret_dict['center'] = '' 

528 ret_dict['width'] = my_qa.quantity(0.0, 'rad') 

529 ret_dict['height'] = my_qa.quantity(0.0, 'rad') 

530 return ret_dict 

531 

532 

533def _handle_image_params(imsize, cell, phasecenter, 

534 vislist, field, spw, antenna, scan, intent, timerange, 

535 restfreq, pointingcolumntouse, ephemsrcname): 

536 logger = sdutil.Casalog(origin="_handle_image_params") 

537 # round-up imsize 

538 _imsize = sdutil.to_list(imsize, int) or sdutil.to_list(imsize, numpy.integer) 

539 if _imsize is None: 

540 _imsize = imsize if hasattr(imsize, '__iter__') else [imsize] 

541 _imsize = [int(numpy.ceil(v)) for v in _imsize] 

542 logger.post( 

543 "imsize is not integers. force converting to integer pixel numbers.", 

544 priority="WARN") 

545 logger.post("rounded-up imsize: %s --> %s" % (str(imsize), str(_imsize))) 

546 

547 # calculate cell based on PB if it is not given 

548 _cell = cell 

549 if _cell == '' or _cell[0] == '': 

550 # calc PB 

551 if isinstance(vislist, str): 

552 vis = vislist 

553 else: 

554 vis = vislist[0] 

555 if isinstance(antenna, str): 

556 antsel = antenna 

557 else: 

558 antsel = antenna[0] 

559 if antsel == '': 

560 antenna_id = 0 

561 else: 

562 if len(antsel) > 3 and antsel[:-3] == '&&&': 

563 baseline = antsel 

564 else: 

565 baseline = antsel + '&&&' 

566 with open_ms(vis) as ms: 

567 ms.msselect({'baseline': baseline}) 

568 ndx = ms.msselectedindices() 

569 antenna_id = ndx['antenna1'][0] 

570 grid_factor = 3. 

571 logger.post("The cell size will be calculated using PB size of antennas in the first MS") 

572 qpb = _calc_PB(vis, antenna_id, restfreq) 

573 _cell = '%f%s' % (qpb['value'] / grid_factor, qpb['unit']) 

574 logger.post("Using cell size = PB/%4.2F = %s" % (grid_factor, _cell)) 

575 

576 # Calculate Pointing center and extent (if necessary) 

577 _phasecenter = phasecenter 

578 if _phasecenter == '' or len(_imsize) == 0 or _imsize[0] < 1: 

579 # return a dictionary with keys 'center', 'width', 'height' 

580 map_param = _get_pointing_extent(_phasecenter, vislist, field, spw, antenna, scan, intent, 

581 timerange, pointingcolumntouse, ephemsrcname) 

582 # imsize 

583 (cellx, celly) = sdutil.get_cellx_celly(_cell, unit='arcmin') 

584 if len(_imsize) == 0 or _imsize[0] < 1: 

585 _imsize = _get_imsize(map_param['width'], map_param['height'], cellx, celly) 

586 if _phasecenter != "": 

587 logger.post( 

588 "You defined phasecenter but not imsize. " 

589 "The image will cover as wide area as pointing in MS extends, " 

590 "but be centered at phasecenter. " 

591 "This could result in a strange image if your phasecenter is " 

592 "apart from the center of pointings", 

593 priority='WARN') 

594 if _imsize[0] > 1024 or _imsize[1] > 1024: 

595 logger.post( 

596 "The calculated image pixel number is larger than 1024. " 

597 "It could take time to generate the image depending on your computer resource. " 

598 "Please wait...", 

599 priority='WARN') 

600 

601 # phasecenter 

602 # if empty, it should be determined here... 

603 if _phasecenter == "": 

604 _phasecenter = map_param['center'] 

605 

606 return _imsize, _cell, _phasecenter 

607 

608 

609def _get_param(ms_index, param): 

610 if isinstance(param, str): 

611 return param 

612 elif hasattr(param, '__iter__'): 

613 if len(param) == 1: 

614 return param[0] 

615 else: 

616 return param[ms_index] 

617 else: 

618 raise RuntimeError('Invalid parameter') 

619 

620 

621def _remove_image(imagename): 

622 if os.path.exists(imagename): 

623 if os.path.isdir(imagename): 

624 shutil.rmtree(imagename) 

625 elif os.path.isfile(imagename): 

626 os.remove(imagename) 

627 else: 

628 # could be a symlink 

629 os.remove(imagename) 

630 

631 

632def _get_restfreq_if_empty(vislist, spw, field, restfreq): 

633 qa = quanta() 

634 rf = None 

635 # if restfreq is nonzero float value, return it 

636 if isinstance(restfreq, float): 

637 if restfreq != 0.0: 

638 rf = restfreq 

639 # if restfreq is valid frequency string, return it 

640 # numeric string is interpreted as a value in the unit of Hz 

641 elif isinstance(restfreq, str): 

642 q = qa.convert(qa.quantity(restfreq), 'Hz') 

643 if q['unit'] == 'Hz' and q['value'] > 0.0: 

644 rf = restfreq 

645 # if restfreq is valid quantity, return it 

646 elif isinstance(restfreq, dict): 

647 q = qa.convert(restfreq, 'Hz') 

648 if q['unit'] == 'Hz' and q['value'] > 0.0: 

649 rf = restfreq 

650 

651 if isinstance(vislist, str): 

652 vis = vislist 

653 elif hasattr(vislist, '__iter__'): 

654 vis = vislist[0] 

655 else: 

656 raise RuntimeError( 

657 'Internal Error: invalid vislist \'{0}\''.format(vislist)) 

658 

659 if isinstance(spw, str): 

660 spwsel = spw 

661 elif hasattr(spw, '__iter__'): 

662 spwsel = spw[0] 

663 else: 

664 raise RuntimeError( 

665 'Internal Error: invalid spw selection \'{0}\''.format(spw)) 

666 

667 if isinstance(field, str): 

668 fieldsel = field 

669 elif hasattr(field, '__iter__'): 

670 fieldsel = field[0] 

671 else: 

672 raise RuntimeError('Internal Error: invalid field selection \'{0}\''.format(field)) 

673 

674 with open_ms(vis) as ms: 

675 ms.msselect({'spw': spwsel, 'field': fieldsel}) 

676 ndx = ms.msselectedindices() 

677 if len(ndx['spw']) > 0: 

678 spwid = ndx['spw'][0] 

679 else: 

680 spwid = None 

681 if len(ndx['field']) > 0: 

682 fieldid = ndx['field'][0] 

683 else: 

684 fieldid = None 

685 sourceid = None 

686 if fieldid is not None: 

687 with sdutil.table_manager(os.path.join(vis, 'FIELD')) as tb: 

688 sourceid = tb.getcell('SOURCE_ID', fieldid) 

689 if sourceid < 0: 

690 sourceid = None 

691 if rf is None: 

692 # if restfrequency is defined in SOURCE table, return it 

693 with sdutil.table_manager(os.path.join(vis, 'SOURCE')) as tb: 

694 if 'REST_FREQUENCY' in tb.colnames(): 

695 tsel = None 

696 taql = '' 

697 if spwid is not None: 

698 taql = 'SPECTRAL_WINDOW_ID == {0}'.format(spwid) 

699 if sourceid is not None: 

700 delimiter = '&&' if len(taql) > 0 else '' 

701 taql += '{0}SOURCE_ID == {1}'.format(delimiter, sourceid) 

702 if len(taql) > 0: 

703 tsel = tb.query(taql) 

704 t = tsel 

705 else: 

706 t = tb 

707 try: 

708 nrow = t.nrows() 

709 if nrow > 0: 

710 for irow in range(nrow): 

711 if t.iscelldefined('REST_FREQUENCY', irow): 

712 rfs = t.getcell('REST_FREQUENCY', irow) 

713 if len(rfs) > 0: 

714 rf = rfs[0] 

715 break 

716 finally: 

717 if tsel is not None: 

718 tsel.close() 

719 

720 if rf is None: 

721 if spwid is None: 

722 spwid = 0 

723 # otherwise, return mean frequency of given spectral window 

724 with sdutil.table_manager(os.path.join(vis, 'SPECTRAL_WINDOW')) as tb: 

725 cf = tb.getcell('CHAN_FREQ', spwid) 

726 rf = cf.mean() 

727 

728 assert rf is not None 

729 

730 return rf 

731 

732 

733def set_beam_size(vis, imagename, 

734 field, spw, baseline, scan, intent, timerange, 

735 ephemsrcname, pointingcolumntouse, 

736 antenna_name, antenna_diameter, 

737 restfreq, 

738 gridfunction, convsupport, truncate, gwidth, jwidth): 

739 """Set estimated beam size to the image.""" 

740 is_alma = antenna_name[0:2] in ['PM', 'DV', 'DA', 'CM'] 

741 blockage = '0.75m' if is_alma else '0.0m' 

742 log_origin = 'set_beam_size' 

743 

744 with open_ia(imagename) as ia: 

745 csys = ia.coordsys() 

746 outref = csys.referencecode('direction')[0] 

747 cell = list(csys.increment(type='direction', format='s')['string']) 

748 csys.done() 

749 

750 old_tool = OldImagerBasedTools() 

751 sampling_params = old_tool.get_pointing_sampling_params( 

752 vis, field, spw, baseline, 

753 scan, intent, timerange, 

754 outref=outref, 

755 movingsource=ephemsrcname, 

756 pointingcolumntouse=pointingcolumntouse, 

757 antenna_name=antenna_name) 

758 qa = quanta() 

759 casalog.post( 

760 f'sampling_params={sampling_params}', 

761 origin=log_origin 

762 ) 

763 xsampling, ysampling = qa.getvalue(qa.convert(sampling_params['sampling'],