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

637 statements  

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

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('Internal Error: invalid vislist \'{0}\''.format(vislist)) 

657 

658 if isinstance(spw, str): 

659 spwsel = spw 

660 elif hasattr(spw, '__iter__'): 

661 spwsel = spw[0] 

662 else: 

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

664 

665 if isinstance(field, str): 

666 fieldsel = field 

667 elif hasattr(field, '__iter__'): 

668 fieldsel = field[0] 

669 else: 

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

671 

672 with open_ms(vis) as ms: 

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

674 ndx = ms.msselectedindices() 

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

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

677 else: 

678 spwid = None 

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

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

681 else: 

682 fieldid = None 

683 sourceid = None 

684 if fieldid is not None: 

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

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

687 if sourceid < 0: 

688 sourceid = None 

689 if rf is None: 

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

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

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

693 tsel = None 

694 taql = '' 

695 if spwid is not None: 

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

697 if sourceid is not None: 

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

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

700 if len(taql) > 0: 

701 tsel = tb.query(taql) 

702 t = tsel 

703 else: 

704 t = tb 

705 try: 

706 nrow = t.nrows() 

707 if nrow > 0: 

708 for irow in range(nrow): 

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

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

711 if len(rfs) > 0: 

712 rf = rfs[0] 

713 break 

714 finally: 

715 if tsel is not None: 

716 tsel.close() 

717 

718 if rf is None: 

719 if spwid is None: 

720 spwid = 0 

721 # otherwise, return mean frequency of given spectral window 

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

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

724 rf = cf.mean() 

725 

726 assert rf is not None 

727 

728 return rf 

729 

730 

731def set_beam_size(vis, imagename, 

732 field, spw, baseline, scan, intent, timerange, 

733 ephemsrcname, pointingcolumntouse, antenna_name, antenna_diameter, 

734 restfreq, gridfunction, convsupport, truncate, gwidth, jwidth): 

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

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

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

738 log_origin = 'set_beam_size' 

739 

740 with open_ia(imagename) as ia: 

741 csys = ia.coordsys() 

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

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

744 csys.done() 

745 

746 old_tool = OldImagerBasedTools() 

747 sampling_params = old_tool.get_pointing_sampling_params(vis, field, spw, baseline, 

748 scan, intent, timerange, 

749 outref=outref, 

750 movingsource=ephemsrcname, 

751 pointingcolumntouse=pointingcolumntouse, 

752 antenna_name=antenna_name) 

753 qa = quanta() 

754 casalog.post( 

755 f'sampling_params={sampling_params}', 

756 origin=log_origin 

757 ) 

758 xsampling, ysampling = qa.getvalue(qa.convert(sampling_params['sampling'], 'arcsec')) 

759 angle = qa.getvalue(qa.convert(sampling_params['angle'], 'deg'))[0] 

760 

761 casalog.post( 

762 f'Detected raster sampling = [{xsampling:f}, {ysampling:f}] arcsec', 

763 origin=log_origin 

764 ) 

765 

766 # handling of failed sampling detection 

767 valid_sampling = True 

768 # TODO: copy from sdimaging implementation 

769 sampling = [xsampling, ysampling] 

770 if abs(xsampling) < 2.2e-3 or not numpy.isfinite(xsampling): 

771 casalog.post( 

772 f"Invalid sampling={xsampling} arcsec. " 

773 f"Using the value of orthogonal direction={ysampling} arcsec", 

774 priority="WARN", 

775 origin=log_origin 

776 ) 

777 sampling = [ysampling] 

778 angle = 0.0 

779 valid_sampling = False 

780 if abs(ysampling) < 1.0e-3 or not numpy.isfinite(ysampling): 

781 if valid_sampling: 

782 casalog.post( 

783 f"Invalid sampling={ysampling} arcsec. " 

784 f"Using the value of orthogonal direction={xsampling} arcsec", 

785 priority="WARN", 

786 origin=log_origin 

787 ) 

788 sampling = [xsampling] 

789 angle = 0.0 

790 valid_sampling = True 

791 # reduce sampling and cell if it's possible 

792 if len(sampling) > 1 and abs(sampling[0] - sampling[1]) <= 0.01 * abs(sampling[0]): 

793 sampling = [sampling[0]] 

794 angle = 0.0 

795 if cell[0] == cell[1]: 

796 cell = [cell[0]] 

797 if valid_sampling: 

798 # actual calculation of beam size 

799 bu = sdbeamutil.TheoreticalBeam() 

800 bu.set_antenna(antenna_diameter, blockage) 

801 bu.set_sampling(sampling, "%fdeg" % angle) 

802 bu.set_image_param(cell, restfreq, gridfunction, 

803 convsupport, truncate, gwidth, 

804 jwidth, is_alma) 

805 bu.summary() 

806 imbeam_dict = bu.get_beamsize_image() 

807 casalog.post( 

808 f"Setting image beam: " 

809 f"major={imbeam_dict['major']}, " 

810 f"minor={imbeam_dict['minor']}, " 

811 f"pa={imbeam_dict['pa']}", 

812 origin=log_origin 

813 ) 

814 # set beam size to image 

815 with open_ia(imagename) as ia: 

816 ia.setrestoringbeam(**imbeam_dict) 

817 else: 

818 # BOTH sampling was invalid 

819 casalog.post( 

820 "Could not detect valid raster sampling. " 

821 "Exiting without setting beam size to image", 

822 priority='WARN', 

823 origin=log_origin 

824 ) 

825 

826 

827def do_weight_mask(imagename, weightimage, minweight): 

828 # Mask image pixels whose weight are smaller than minweight. 

829 # Weight image should have 0 weight for pixels below < minweight 

830 logger = sdutil.Casalog(origin="do_weight_mask") 

831 logger.post(f"Start masking the map using minweight = {minweight:f}", 

832 priority="INFO") 

833 with open_ia(weightimage) as ia: 

834 try: 

835 stat = ia.statistics(mask="'" + weightimage + "' > 0.0", robust=True) 

836 valid_pixels = stat['npts'] 

837 except RuntimeError as e: 

838 if 'No valid data found.' in str(e): 

839 valid_pixels = [0] 

840 else: 

841 raise e 

842 

843 if len(valid_pixels) == 0 or valid_pixels[0] == 0: 

844 logger.post( 

845 "All pixels weight zero. " 

846 "This indicates no data in MS is in image area. " 

847 "Mask will not be set. Please check your image parameters.", 

848 priority="WARN") 

849 return 

850 median_weight = stat['median'][0] 

851 weight_threshold = median_weight * minweight 

852 logger.post(f"Median of weight in the map is {median_weight:f}", 

853 priority="INFO") 

854 logger.post(f"Pixels in map with weight <= median(weight)*minweight = " 

855 "{weight_threshold:f} will be masked.", 

856 priority="INFO") 

857 # Leaving the original logic to calculate the number of masked pixels via 

858 # product of median of and min_weight (which i don't understand the logic) 

859 

860 # Modify default mask 

861 with open_ia(imagename) as ia: 

862 ia.calcmask("'%s'>%f" % (weightimage, weight_threshold), asdefault=True) 

863 

864 ndim = len(ia.shape()) 

865 _axes = numpy.arange(start=0 if ndim <= 2 else 2, stop=ndim) 

866 

867 try: 

868 collapsed = ia.collapse('npts', axes=_axes) 

869 valid_pixels_after = collapsed.getchunk().sum() 

870 collapsed.close() 

871 except RuntimeError as e: 

872 if 'All selected pixels are masked' in str(e): 

873 valid_pixels_after = 0 

874 else: 

875 raise 

876 

877 masked_fraction = 100. * (1. - valid_pixels_after / float(valid_pixels[0])) 

878 

879 logger.post(f"This amounts to {masked_fraction:5.1f} % " 

880 "of the area with nonzero weight.", 

881 priority="INFO") 

882 logger.post( 

883 f"The weight image '{weightimage}' is returned by this task, " 

884 "if the user wishes to assess the results in detail.", 

885 priority="INFO") 

886 

887 

888def get_ms_column_unit(tb, colname): 

889 col_unit = '' 

890 if colname in tb.colnames(): 

891 cdkw = tb.getcoldesc(colname)['keywords'] 

892 if 'QuantumUnits' in cdkw: 

893 u = cdkw['QuantumUnits'] 

894 if isinstance(u, str): 

895 col_unit = u.strip() 

896 elif isinstance(u, list): 

897 col_unit = u[0].strip() 

898 return col_unit 

899 

900 

901def get_brightness_unit_from_ms(msname): 

902 image_unit = '' 

903 with sdutil.table_manager(msname) as tb: 

904 image_unit = get_ms_column_unit(tb, 'DATA') 

905 if image_unit == '': 

906 image_unit = get_ms_column_unit(tb, 'FLOAT_DATA') 

907 if image_unit.upper() == 'K': 

908 image_unit = 'K' 

909 else: 

910 image_unit = 'Jy/beam' 

911 

912 return image_unit 

913 

914 

915@sdutil.sdtask_decorator 

916def tsdimaging( 

917 # Input data: list of MeasurementSets 

918 infiles, 

919 # Output data: CASA images: their path prefix, overwrite control 

920 outfile, overwrite, 

921 # Select data from input MeasurementSets, by 

922 field, spw, antenna, scan, intent, timerange, 

923 # Output images definition: frequency axis 

924 outframe, # velocity frame 

925 mode, nchan, start, width, veltype, # gridding type 

926 specmode, # Doppler handling 

927 interpolation, # interpolation mode 

928 # Output images definition: spatial axes 

929 pointingcolumn, convertfirst, 

930 projection, 

931 imsize, cell, phasecenter, 

932 # Output images definition: stokes axis 

933 stokes, 

934 # Gridder parameters 

935 gridfunction, convsupport, truncate, gwidth, jwidth, 

936 clipminmax, 

937 # Single-dish image: mask control 

938 minweight, 

939 # Single-dish image: metadata 

940 brightnessunit, 

941 restfreq # rest frequency to assign to image 

942 ): 

943 

944 origin = 'tsdimaging' 

945 imager = None 

946 

947 try: # Create the Single-Dish Image 

948 # Validate brightnessunit parameter CAS-11503 

949 image_unit = brightnessunit.lower().capitalize() 

950 if not image_unit in ['', 'K', 'Jy/beam']: 

951 raise ValueError(f"Invalid brightness unit: {brightnessunit}") 

952 

953 # Handle outfile and overwrite parameters 

954 output_path_prefix = outfile.rstrip('/') 

955 singledish_image_path = output_path_prefix + image_suffix 

956 if os.path.exists(singledish_image_path): 

957 if overwrite == False: 

958 raise RuntimeError( 

959 f"Output file exists: '{singledish_image_path}'" 

960 ) 

961 else: 

962 # delete existing images 

963 if os.path.exists(singledish_image_path): 

964 casalog.post(f"Removing '{singledish_image_path}'") 

965 _remove_image(singledish_image_path) 

966 assert not os.path.exists(singledish_image_path) 

967 for _suffix in associate_suffixes: 

968 path_to_remove = output_path_prefix + _suffix 

969 if os.path.exists(path_to_remove): 

970 casalog.post(f"Removing '{path_to_remove}'") 

971 _remove_image(path_to_remove) 

972 assert not os.path.exists(path_to_remove) 

973 

974 # Tweak spw parameter into _spw 

975 if isinstance(spw, str): 

976 _spw = '*' + spw if spw.startswith(':') else spw 

977 else: 

978 _spw = ['*' + v if v.startswith(':') else v for v in spw] 

979 

980 # Handle image spectral axis parameters 

981 imnchan, imstart, imwidth = _configure_spectral_axis( 

982 mode, nchan, start, width, restfreq 

983 ) 

984 

985 # Handle image restfreq parameter's default value 

986 _restfreq = _get_restfreq_if_empty( 

987 infiles, _spw, field, restfreq 

988 ) 

989 

990 # Handle gridder parameters 

991 # convert type of task's default values 

992 # to ones supported by the Synthesis Imager framework 

993 gtruncate = _handle_grid_defaults(truncate) 

994 ggwidth = _handle_grid_defaults(gwidth) 

995 gjwidth = _handle_grid_defaults(jwidth) 

996 

997 # handle infiles parameter 

998 # sort input data to get results consistent with old sdimaging task 

999 _sorted = sort_vis( 

1000 infiles, _spw, mode, imwidth, field, 

1001 antenna, scan, intent, timerange 

1002 ) 

1003 sorted_vis = _sorted[0] 

1004 sorted_field = _sorted[1] 

1005 sorted_spw = _sorted[2] 

1006 sorted_antenna = _sorted[3] 

1007 sorted_scan = _sorted[4] 

1008 sorted_intent = _sorted[5] 

1009 sorted_timerange = _sorted[6] 

1010 

1011 def antenna_to_baseline(s): 

1012 if len(s) == 0: 

1013 return s 

1014 elif len(s) > 3 and s.endswith('&&&'): 

1015 return s 

1016 else: 

1017 return '{0}&&&'.format(s) 

1018 

1019 if isinstance(sorted_antenna, str): 

1020 sorted_baseline = antenna_to_baseline(sorted_antenna) 

1021 else: 

1022 sorted_baseline = [antenna_to_baseline(a) for a in sorted_antenna] 

1023 

1024 # Handle image geometric parameters 

1025 _ephemsrcname = '' 

1026 ephem_sources = ['MERCURY', 'VENUS', 'MARS', 'JUPITER', 'SATURN', 

1027 'URANUS', 'NEPTUNE', 'PLUTO', 'SUN', 'MOON', 'TRACKFIELD'] 

1028 if isinstance(phasecenter, str) and phasecenter.strip().upper() in ephem_sources: 

1029 _ephemsrcname = phasecenter 

1030 _imsize, _cell, _phasecenter = _handle_image_params( 

1031 imsize, cell, phasecenter, sorted_vis, 

1032 sorted_field, sorted_spw, sorted_antenna, 

1033 sorted_scan, sorted_intent, sorted_timerange, 

1034 _restfreq, pointingcolumn, _ephemsrcname 

1035 ) 

1036 

1037 # Set up PySynthesisImager input parameters 

1038 # - List all parameters that you need here 

1039 # - Defaults will be assumed for unspecified parameters 

1040 # - Nearly all parameters are identical to that in the task. Please look at the 

1041 # list of parameters under __init__ using " help ImagerParameters " ) 

1042 casalog.post('*** Creating paramList ***', origin=origin) 

1043 paramList = ImagerParameters( 

1044 # input file name 

1045 msname=sorted_vis, # 'sdimaging.ms', 

1046 # data selection 

1047 field=sorted_field, # '', 

1048 spw=sorted_spw, # '0', 

1049 timestr=sorted_timerange, 

1050 antenna=sorted_baseline, 

1051 scan=sorted_scan, 

1052 state=sorted_intent, 

1053 # image parameters 

1054 imagename=output_path_prefix, # 'try2', 

1055 nchan=imnchan, # 1024, 

1056 start=imstart, # '0', 

1057 width=imwidth, # '1', 

1058 outframe=outframe, 

1059 veltype=veltype, 

1060 restfreq=_restfreq, 

1061 phasecenter=_phasecenter, # 'J2000 17:18:29 +59.31.23', 

1062 imsize=_imsize, # [75,75], 

1063 cell=_cell, # ['3arcmin', '3arcmin'], 

1064 projection=projection, 

1065 stokes=stokes, 

1066 specmode=specmode, 

1067 gridder='singledish', 

1068 # single dish specific parameters 

1069 gridfunction=gridfunction, 

1070 convsupport=convsupport, 

1071 truncate=gtruncate, 

1072 gwidth=ggwidth, 

1073 jwidth=gjwidth, 

1074 pointingcolumntouse=pointingcolumn, 

1075 convertfirst=convertfirst, 

1076 minweight=minweight, 

1077 clipminmax=clipminmax, 

1078 # normalizer 

1079 normtype='flatsky', 

1080 pblimit=1e-16, # TODO: explain why 1e-16 ? 

1081 interpolation=interpolation 

1082 ) 

1083 

1084 # Construct the PySynthesisImager object, with all input parameters 

1085 casalog.post('*** Creating imager object ***', origin=origin) 

1086 imager = PySynthesisImager(params=paramList) 

1087 

1088 # Initialize PySynthesisImager "modules" required for single-dish imaging 

1089 # - Pick only the modules you will need later on. For example, to only make 

1090 # the PSF, there is no need for the deconvolver or iteration control modules. 

1091 casalog.post('*** Initializing imagers ***', origin=origin) 

1092 # This is where the underlying C++ synthesis imager is created 

1093 imager.initializeImagers() 

1094 casalog.post('*** Initializing normalizers ***', origin=origin) 

1095 imager.initializeNormalizers() 

1096 

1097 # Create Single-Dish images 

1098 casalog.post('*** Creating single-dish images ***', origin=origin) 

1099 imager.makeSdImage() 

1100 casalog.post('*** Created single-dish images ***', origin=origin) 

1101 

1102 finally: # Close tools and rename Synthesis Imager's residual image 

1103 casalog.post('*** Cleaning up tools ***', origin=origin) 

1104 if imager is not None: 

1105 imager.deleteTools() 

1106 # Change image suffix from .residual to .image 

1107 residual_image_path = output_path_prefix + residual_suffix 

1108 if os.path.exists(residual_image_path): 

1109 os.rename(residual_image_path, singledish_image_path) 

1110 

1111 # Set single-dish image's beam size 

1112 # TODO: re-define related functions in the new tool framework (sdms?) 

1113 # ms_index = 0 

1114 rep_ms = _get_param(0, infiles) 

1115 rep_field = _get_param(0, field) 

1116 rep_spw = _get_param(0, _spw) 

1117 rep_antenna = _get_param(0, antenna) 

1118 rep_scan = _get_param(0, scan) 

1119 rep_intent = _get_param(0, intent) 

1120 rep_timerange = _get_param(0, timerange) 

1121 if len(rep_antenna) > 0: 

1122 baseline = '{0}&&&'.format(rep_antenna) 

1123 else: 

1124 baseline = '*&&&' 

1125 with open_ms(rep_ms) as ms: 

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

1127 ndx = ms.msselectedindices() 

1128 antenna_index = ndx['antenna1'][0] 

1129 with sdutil.table_manager(os.path.join(rep_ms, 'ANTENNA')) as tb: 

1130 antenna_name = tb.getcell('NAME', antenna_index) 

1131 antenna_diameter = tb.getcell('DISH_DIAMETER', antenna_index) 

1132 casalog.post(f"Setting single-dish image's beam") 

1133 set_beam_size( 

1134 rep_ms, singledish_image_path, 

1135 rep_field, rep_spw, baseline, rep_scan, rep_intent, rep_timerange, 

1136 _ephemsrcname, pointingcolumn, antenna_name, antenna_diameter, 

1137 _restfreq, gridfunction, convsupport, truncate, gwidth, jwidth 

1138 ) 

1139 

1140 # Set single-dish image's brightness unit (CAS-11503) 

1141 if len(image_unit) == 0: 

1142 image_unit = get_brightness_unit_from_ms(rep_ms) 

1143 if len(image_unit) > 0: 

1144 with open_ia(singledish_image_path) as ia: 

1145 casalog.post(f"Setting single-dish image's brightness unit to '{image_unit}'") 

1146 ia.setbrightnessunit(image_unit) 

1147 

1148 # Update single-dish image's mask: mask low weight pixels 

1149 weight_image_path = output_path_prefix + weight_suffix 

1150 casalog.post(f"Creating weight image mask") 

1151 do_weight_mask(singledish_image_path, weight_image_path, minweight) 

1152 

1153 # Delete images systematically generated by the Synthesis Imager 

1154 # which are either not required or currently useless 

1155 # in the context of single-dish imaging 

1156 # CAS-10891 

1157 _remove_image(output_path_prefix + '.sumwt') 

1158 # CAS-10893 

1159 # TODO: remove the following line once the 'correct' SD 

1160 # PSF image based on primary beam can be generated 

1161 _remove_image(output_path_prefix + '.psf')