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

800 statements  

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

1# sd task for imaging 

2import os 

3import re 

4import shutil 

5import string 

6 

7import numpy 

8 

9from casatasks import casalog 

10from casatools import image, measures, ms, msmetadata, quanta, table 

11 

12from . import mslisthelper, sdbeamutil, sdutil 

13from .task_tsdimaging import conform_mslist 

14 

15 

16@sdutil.sdtask_decorator 

17def sdimaging(infiles, outfile, overwrite, field, spw, antenna, scan, intent, 

18 mode, nchan, start, width, veltype, outframe, 

19 gridfunction, convsupport, truncate, gwidth, jwidth, 

20 imsize, cell, phasecenter, projection, ephemsrcname, 

21 pointingcolumn, restfreq, stokes, minweight, brightnessunit, clipminmax, 

22 # Performances optimization options 

23 enablecache, convertfirst, 

24 interpolation): 

25 with sdimaging_worker(**locals()) as worker: 

26 worker.initialize() 

27 worker.execute() 

28 worker.finalize() 

29 

30 

31def is_string_type(val): 

32 """Return True if the argument is string type.""" 

33 return type(val) in [str, numpy.str_] 

34 

35 

36def smart_remove(path): 

37 if os.path.isdir(path): 

38 shutil.rmtree(path) 

39 else: 

40 os.remove(path) 

41 

42 

43def get_subtable_path(vis, name): 

44 return os.path.join(vis, name) 

45 

46 

47# functions from cleanhelper 

48class cleanhelper_minimal(object): 

49 def __init__(self, imtool='', vis='', sort_index=None, usescratch=False, casalog=casalog): 

50 if((type(vis)==list) and (len(vis)==1)): 

51 vis=vis[0] 

52 #### 

53 if isinstance(vis, str): 

54 vislist = [vis] 

55 else: 

56 vislist = vis 

57 assert not isinstance(imtool, str) 

58 self.im = imtool 

59 self.vislist = vislist 

60 assert sort_index is not None 

61 self.sortedvisindx = sort_index 

62 self.sortedvislist = [self.vislist[i] for i in sort_index] 

63 

64 self.dataspecframe='LSRK' 

65 self.usespecframe='' 

66 self.inframe=False 

67 self._casalog = casalog 

68 

69 def convertvf(self,vf,frame,field,restf,veltype='radio'): 

70 """ 

71 returns doppler(velocity) or frequency in string 

72 currently use first rest frequency 

73 Assume input vf (velocity or fequency in a string) and 

74 output are the same 'frame'. 

75 """ 

76 #pdb.set_trace() 

77 docalcf=False 

78 #if(frame==''): frame='LSRK' 

79 #Use datasepcframe, it is cleanhelper initialized to set 

80 #to LSRK 

81 if(frame==''): frame=self.dataspecframe 

82 qa = quanta() 

83 if(qa.quantity(vf)['unit'].find('m/s') > -1): 

84 docalcf=True 

85 elif(qa.quantity(vf)['unit'].find('Hz') > -1): 

86 docalcf=False 

87 else: 

88 if vf !=0: 

89 raise TypeError("Unrecognized unit for the velocity or frequency parameter") 

90 ##fldinds=ms.msseltoindex(self.vis, field=field)['field'].tolist() 

91 myms = ms() 

92 fldinds=myms.msseltoindex(self.vislist[self.sortedvisindx[0]], field=field)['field'].tolist() 

93 if(len(fldinds) == 0): 

94 fldid0=0 

95 else: 

96 fldid0=fldinds[0] 

97 if restf=='': 

98 #tb.open(self.vis+'/FIELD') 

99 fldtab=get_subtable_path(self.vislist[self.sortedvisindx[0]],'FIELD') 

100 tb = table() 

101 tb.open(fldtab) 

102 nfld = tb.nrows() 

103 if nfld >= fldid0: 

104 srcid=tb.getcell('SOURCE_ID',fldid0) 

105 else: 

106 raise TypeError( "Cannot set REST_FREQUENCY from the data: " + 

107 "no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 ) 

108 tb.close() 

109 # SOUECE_ID in FIELD table = -1 if no SOURCE table 

110 if srcid==-1: 

111 raise TypeError("Rest frequency info is not supplied") 

112 #tb.open(self.vis+'/SOURCE') 

113 sourcetab=get_subtable_path(self.vislist[self.sortedvisindx[0]], 'SOURCE') 

114 tb.open(sourcetab) 

115 tb2=tb.query('SOURCE_ID==%s' % srcid) 

116 tb.close() 

117 nsrc = tb2.nrows() 

118 if nsrc > 0: 

119 rfreq=tb2.getcell('REST_FREQUENCY',0) 

120 else: 

121 raise TypeError( "Cannot set REST_FREQUENCY from the data: "+ 

122 " no SOURCE corresponding field ID=%s, please supply restfreq" % fldid0 ) 

123 tb2.close() 

124 if(rfreq<=0): 

125 raise TypeError("Rest frequency does not seems to be properly set, check the data") 

126 else: 

127 if type(restf)==str: restf=[restf] 

128 if(qa.quantity(restf[0])['unit'].find('Hz') > -1): 

129 rfreq=[qa.convert(qa.quantity(restf[0]),'Hz')['value']] 

130 #print("using user input rest freq=",rfreq) 

131 else: 

132 raise TypeError("Unrecognized unit or type for restfreq") 

133 if(vf==0): 

134 # assume just want to get a restfrequecy from the data 

135 ret=str(rfreq[0])+'Hz' 

136 else: 

137 me = measures() 

138 if(docalcf): 

139 dop=me.doppler(veltype, qa.quantity(vf)) 

140 rvf=me.tofrequency(frame, dop, qa.quantity(rfreq[0],'Hz')) 

141 else: 

142 frq=me.frequency(frame, qa.quantity(vf)) 

143 rvf=me.todoppler(veltype, frq, qa.quantity(rfreq[0],'Hz')) 

144 ret=str(rvf['m0']['value'])+rvf['m0']['unit'] 

145 return ret 

146 

147 

148 def setChannelizeDefault(self,mode,spw,field,nchan,start,width,frame,veltype,phasec, restf,obstime=''): 

149 """ 

150 Determine appropriate values for channelization 

151 parameters when default values are used 

152 for mode='velocity' or 'frequency' or 'channel' 

153 This makes use of ms.cvelfreqs. 

154 """ 

155 ############### 

156 # for debugging 

157 ############### 

158 debug=False 

159 ############### 

160 spectable=get_subtable_path(self.vislist[self.sortedvisindx[0]], "SPECTRAL_WINDOW") 

161 tb = table() 

162 tb.open(spectable) 

163 chanfreqscol=tb.getvarcol('CHAN_FREQ') 

164 chanwidcol=tb.getvarcol('CHAN_WIDTH') 

165 spwframe=tb.getcol('MEAS_FREQ_REF'); 

166 tb.close() 

167 # first parse spw parameter: 

168 # use MSSelect if possible 

169 if len(self.sortedvislist) > 0: 

170 invis = self.sortedvislist[0] 

171 inspw = self.vislist.index(self.sortedvislist[0]) 

172 else: 

173 invis = self.vislist[0] 

174 inspw = 0 

175 myms = ms() 

176 myms.open(invis) 

177 if type(spw)==list: 

178 spw=spw[inspw] 

179 if spw in ('-1', '*', '', ' '): 

180 spw='*' 

181 if field=='': 

182 field='*' 

183 mssel=myms.msseltoindex(vis=invis, spw=spw, field=field) 

184 selspw=mssel['spw'] 

185 selfield=mssel['field'] 

186 chaninds=mssel['channel'].tolist() 

187 chanst0 = chaninds[0][1] 

188 

189 # frame 

190 spw0=selspw[0] 

191 chanfreqs=chanfreqscol['r'+str(spw0+1)].transpose()[0] 

192 chanres = chanwidcol['r'+str(spw0+1)].transpose()[0] 

193 

194 # ascending or desending data frequencies? 

195 # based on selected first spw's first CHANNEL WIDTH 

196 # ==> some MS data may have positive chan width 

197 # so changed to look at first two channels of chanfreq (TT) 

198 #if chanres[0] < 0: 

199 descending = False 

200 if len(chanfreqs) > 1 : 

201 if chanfreqs[1]-chanfreqs[0] < 0: 

202 descending = True 

203 else: 

204 if chanres[0] < 0: 

205 descending = True 

206 

207 # set dataspecframe: 

208 elspecframe=["REST", 

209 "LSRK", 

210 "LSRD", 

211 "BARY", 

212 "GEO", 

213 "TOPO", 

214 "GALACTO", 

215 "LGROUP", 

216 "CMB"] 

217 self.dataspecframe=elspecframe[spwframe[spw0]]; 

218 

219 # set usespecframe: user's frame if set, otherwise data's frame 

220 if(frame != ''): 

221 self.usespecframe=frame 

222 self.inframe=True 

223 else: 

224 self.usespecframe=self.dataspecframe 

225 

226 # some start and width default handling 

227 if mode!='channel': 

228 if width==1: 

229 width='' 

230 if start==0: 

231 start='' 

232 

233 #get restfreq 

234 if restf=='': 

235 fldtab=get_subtable_path(invis,'FIELD') 

236 tb.open(fldtab) 

237 nfld=tb.nrows() 

238 try: 

239 if nfld >= selfield[0]: 

240 srcid=tb.getcell('SOURCE_ID',selfield[0]) 

241 else: 

242 if mode=='velocity': 

243 raise TypeError( "Cannot set REST_FREQUENCY from the data: " + 

244 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] ) 

245 finally: 

246 tb.close() 

247 #SOUECE_ID in FIELD table = -1 if no SOURCE table 

248 if srcid==-1: 

249 if mode=='velocity': 

250 raise TypeError("Rest frequency info is not supplied") 

251 try: 

252 srctab=get_subtable_path(invis, 'SOURCE') 

253 tb.open(srctab) 

254 tb2=tb.query('SOURCE_ID==%s' % srcid) 

255 nsrc = tb2.nrows() 

256 if nsrc > 0 and tb2.iscelldefined('REST_FREQUENCY',0): 

257 rfqs = tb2.getcell('REST_FREQUENCY',0) 

258 if len(rfqs)>0: 

259 restf=str(rfqs[0])+'Hz' 

260 else: 

261 if mode=='velocity': 

262 raise TypeError( "Cannot set REST_FREQUENCY from the data: " + 

263 "REST_FREQUENCY entry for ID %s in SOURCE table is empty, please supply restfreq" % srcid ) 

264 else: 

265 if mode=='velocity': 

266 raise TypeError( "Cannot set REST_FREQUENCY from the data: " + 

267 "no SOURCE corresponding field ID=%s, please supply restfreq" % selfield[0] ) 

268 finally: 

269 tb.close() 

270 tb2.close() 

271 

272 if type(phasec)==list: 

273 inphasec=phasec[0] 

274 else: 

275 inphasec=phasec 

276 if type(inphasec)==str and inphasec.isdigit(): 

277 inphasec=int(inphasec) 

278 #if nchan==1: 

279 # use data chan freqs 

280 # newfreqs=chanfreqs 

281 #else: 

282 # obstime not included here 

283 if debug: print("before ms.cvelfreqs (start,width,nchan)===>",start, width, nchan) 

284 try: 

285 newfreqs=myms.cvelfreqs(spwids=selspw,fieldids=selfield,mode=mode,nchan=nchan, 

286 start=start,width=width,phasec=inphasec, restfreq=restf, 

287 outframe=self.usespecframe,veltype=veltype).tolist() 

288 except: 

289 # ms must be closed here if ms.cvelfreqs failed with an exception 

290 myms.close() 

291 raise 

292 myms.close() 

293 

294 #print(newfreqs) 

295 descendingnewfreqs=False 

296 if len(newfreqs)>1: 

297 if newfreqs[1]-newfreqs[0] < 0: 

298 descendingnewfreqs=True 

299 

300 

301 try: 

302 if((nchan in [-1, "-1", "", " "]) or (start in [-1, "-1", "", " "])): 

303 frange=im.advisechansel(msname=invis, spwselection=spw, fieldid=selfield[0], getfreqrange=True) 

304 startchan=0 

305 endchan=len(newfreqs)-1 

306 if(descendingnewfreqs): 

307 startchan=numpy.min(numpy.where(frange['freqend'] < numpy.array(newfreqs))) 

308 endchan=numpy.min(numpy.where(frange['freqstart'] < numpy.array(newfreqs))) 

309 else: 

310 startchan=numpy.max(numpy.where(frange['freqstart'] > numpy.array(newfreqs))) 

311 endchan=numpy.max(numpy.where(frange['freqend'] > numpy.array(newfreqs))) 

312 if((start not in [-1, "-1", "", " "]) and (mode=="channel")): 

313 startchan=start 

314 if(nchan not in [-1, "-1", "", " "]): 

315 endchan=startchan+nchan-1 

316 newfreqs=(numpy.array(newfreqs)[startchan:endchan]).tolist() 

317 except: 

318 pass 

319 if debug: print("Mode, Start, width after cvelfreqs =",mode, start,width ) 

320 if type(newfreqs)==list and len(newfreqs) ==0: 

321 raise TypeError( "Output frequency grid cannot be calculated: " + 

322 " please check start and width parameters" ) 

323 if debug: 

324 if len(newfreqs)>1: 

325 print("FRAME=",self.usespecframe) 

326 print("newfreqs[0]===>",newfreqs[0]) 

327 print("newfreqs[1]===>",newfreqs[1]) 

328 print("newfreqs[-1]===>",newfreqs[-1]) 

329 print("len(newfreqs)===>",len(newfreqs)) 

330 else: 

331 print("newfreqs=",newfreqs) 

332 

333 # set output number of channels 

334 if nchan ==1: 

335 retnchan=1 

336 else: 

337 if len(newfreqs)>1: 

338 retnchan=len(newfreqs) 

339 else: 

340 retnchan=nchan 

341 newfreqs=chanfreqs.tolist() 

342 

343 # set start parameter 

344 # first analyze data order etc 

345 reverse=False 

346 negativew=False 

347 if descending: 

348 # channel mode case (width always >0) 

349 if width!="" and (type(width)==int or type(width)==float): 

350 if descendingnewfreqs: 

351 reverse=False 

352 else: 

353 reverse=True 

354 elif width=="": #default width 

355 if descendingnewfreqs and mode=="frequency": 

356 reverse=False 

357 else: 

358 reverse=True 

359 

360 elif type(width)==str: 

361 if width.lstrip().find('-')==0: 

362 negativew=True 

363 if descendingnewfreqs: 

364 if negativew: 

365 reverse=False 

366 else: 

367 reverse=True 

368 else: 

369 if negativew: 

370 reverse=True 

371 else: 

372 reverse=False 

373 else: #ascending data 

374 # depends on sign of width only 

375 # with CAS-3117 latest change(rev.15179), velocity start 

376 # means lowest velocity for default width 

377 if width=="" and mode=="velocity": #default width 

378 # ms.cvelfreqs returns correct order so no reversing 

379 reverse=False 

380 elif type(width)==str: 

381 if width.lstrip().find('-')==0: 

382 reverse=True 

383 else: 

384 reverse=False 

385 

386 if reverse: 

387 newfreqs.reverse() 

388 #if (start!="" and mode=='channel') or \ 

389 # (start!="" and type(start)!=int and mode!='channel'): 

390 # for now to avoid inconsistency later in imagecoordinates2 call 

391 # user's start parameter is preserved for channel mode only. 

392 # (i.e. the current code may adjust start parameter for other modes but 

393 # this probably needs to be changed, especially for multiple ms handling.) 

394 if ((start not in [-1, "", " "]) and mode=='channel'): 

395 retstart=start 

396 else: 

397 # default cases 

398 if mode=="frequency": 

399 retstart=str(newfreqs[0])+'Hz' 

400 elif mode=="velocity": 

401 #startfreq=str(newfreqs[-1])+'Hz' 

402 startfreq=(str(max(newfreqs))+'Hz') if(start=="") else (str(newfreqs[-1])+'Hz') 

403 retstart=self.convertvf(startfreq,frame,field,restf,veltype) 

404 elif mode=="channel": 

405 # default start case, use channel selection from spw 

406 retstart=chanst0 

407 

408 # set width parameter 

409 if width!="": 

410 retwidth=width 

411 else: 

412 if nchan==1: 

413 finc = chanres[0] 

414 else: 

415 finc = newfreqs[1]-newfreqs[0] 

416 if debug: print("finc(newfreqs1-newfreqs0)=",finc) 

417 if mode=="frequency": 

418 # It seems that this is no longer necessary... TT 2013-08-12 

419 #if descendingnewfreqs: 

420 # finc = -finc 

421 retwidth=str(finc)+'Hz' 

422 elif mode=="velocity": 

423 # for default width assume it is vel<0 (incresing in freq) 

424 if descendingnewfreqs: 

425 ind1=-2 

426 ind0=-1 

427 else: 

428 ind1=-1 

429 ind0=-2 

430 v1 = self.convertvf(str(newfreqs[ind1])+'Hz',frame,field,restf,veltype=veltype) 

431 v0 = self.convertvf(str(newfreqs[ind0])+'Hz',frame,field,restf,veltype=veltype) 

432 ##v1 = self.convertvf(str(newfreqs[-1])+'Hz',frame,field,restf,veltype=veltype) 

433 ##v0 = self.convertvf(str(newfreqs[-2])+'Hz',frame,field,restf,veltype=veltype) 

434 #v1 = self.convertvf(str(newfreqs[1])+'Hz',frame,field,restf,veltype=veltype) 

435 #v0 = self.convertvf(str(newfreqs[0])+'Hz',frame,field,restf,veltype=veltype) 

436 qa = quanta() 

437 if(qa.lt(v0, v1) and start==""): 

438 ###user used "" as start make sure step is +ve in vel as start is min vel possible for freqs selected 

439 retwidth=qa.tos(qa.sub(v1, v0)) 

440 else: 

441 retwidth = qa.tos(qa.sub(v0, v1)) 

442 else: 

443 retwidth=1 

444 if debug: print("setChan retwidth=",retwidth) 

445 return retnchan, retstart, retwidth 

446 

447 

448class sdimaging_worker(sdutil.sdtask_template_imaging): 

449 def __init__(self, **kwargs): 

450 super(sdimaging_worker, self).__init__(**kwargs) 

451 self.imager_param = {} 

452 self.sorted_idx = [] 

453 self.image_unit = "" 

454 

455 def __exit__(self, exc_type, exc_val, exc_tb): 

456 self.__del__() 

457 if exc_type: 

458 return False 

459 else: 

460 return True 

461 

462 def parameter_check(self): 

463 # outfile check 

464 sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile, 

465 'im', 

466 self.overwrite) 

467 sdutil.assert_outfile_canoverwrite_or_nonexistent(self.outfile + '.weight', 

468 'im', 

469 self.overwrite) 

470 # fix spw 

471 if type(self.spw) == str: 

472 self.spw = self.__format_spw_string(self.spw) 

473 # check unit of start and width 

474 # fix default 

475 if self.mode == 'channel': 

476 if self.start == '': 

477 self.start = 0 

478 if self.width == '': 

479 self.width = 1 

480 else: 

481 if self.start == 0: 

482 self.start = '' 

483 if self.width == 1: 

484 self.width = '' 

485 # fix unit 

486 if self.mode == 'frequency': 

487 myunit = 'Hz' 

488 elif self.mode == 'velocity': 

489 myunit = 'km/s' 

490 else: # channel 

491 myunit = '' 

492 

493 for name in ['start', 'width']: 

494 param = getattr(self, name) 

495 new_param = self.__format_quantum_unit(param, myunit) 

496 if new_param is None: 

497 raise ValueError("Invalid unit for %s in mode %s: %s" % 

498 (name, self.mode, param)) 

499 setattr(self, name, new_param) 

500 

501 casalog.post("mode='%s': start=%s, width=%s, nchan=%d" % 

502 (self.mode, self.start, self.width, self.nchan)) 

503 

504 # check length of selection parameters 

505 if is_string_type(self.infiles): 

506 nfile = 1 

507 self.infiles = [self.infiles] 

508 else: 

509 nfile = len(self.infiles) 

510 

511 for name in ['field', 'spw', 'antenna', 'scanno']: 

512 param = getattr(self, name) 

513 if not self.__check_selection_length(param, nfile): 

514 raise ValueError("Length of %s != infiles." % (name)) 

515 # set convsupport default 

516 if self.convsupport >= 0 and self.gridfunction.upper() != 'SF': 

517 casalog.post( 

518 f"user defined convsupport is ignored for {self.gridfunction} kernel", 

519 priority='WARN') 

520 self.convsupport = -1 

521 

522 def __format_spw_string(self, spw): 

523 """Return formatted spw selection string which is accepted by imager.""" 

524 if type(spw) != str: 

525 raise ValueError("The parameter should be string.") 

526 if spw.strip() == '*': 

527 spw = '' 

528 # WORKAROUND for CAS-6422, i.e., ":X~Y" fails while "*:X~Y" works. 

529 if spw.startswith(":"): 

530 spw = '*' + spw 

531 return spw 

532 

533 def __format_quantum_unit(self, data, unit): 

534 """Format quantity data. 

535 

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

537 input unit. 

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

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

540 """ 

541 my_qa = quanta() 

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

543 return data 

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

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

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

547 return None 

548 

549 def __check_selection_length(self, data, nfile): 

550 """Check length of the data. 

551 

552 Returns true if data is either a string, an array with length 

553 1 or nfile 

554 """ 

555 if not is_string_type(data) and len(data) not in [1, nfile]: 

556 return False 

557 return True 

558 

559 def get_selection_param_for_ms(self, fileid, param): 

560 """Return valid selection string for a certain ms. 

561 

562 Arguments 

563 fileid : file idx in infiles list 

564 param : string (array) selection value 

565 """ 

566 if is_string_type(param): 

567 return param 

568 elif len(param) == 1: 

569 return param[0] 

570 else: 

571 return param[fileid] 

572 

573 def get_selection_idx_for_ms(self, file_idx): 

574 """Return a dictionary of selection indices for i-th MS in infiles. 

575 

576 Argument: file idx in infiles list 

577 """ 

578 if file_idx < len(self.infiles) and file_idx > -1: 

579 vis = self.infiles[file_idx] 

580 field = self.get_selection_param_for_ms(file_idx, self.field) 

581 spw = self.get_selection_param_for_ms(file_idx, self.spw) 

582 spw = self.__format_spw_string(spw) 

583 antenna = self.get_selection_param_for_ms(file_idx, self.antenna) 

584 if antenna == -1: 

585 antenna = '' 

586 scan = self.get_selection_param_for_ms(file_idx, self.scanno) 

587 intent = self.get_selection_param_for_ms(file_idx, self.intent) 

588 my_ms = ms() 

589 sel_ids = my_ms.msseltoindex(vis=vis, spw=spw, field=field, 

590 baseline=antenna, scan=scan) 

591 fieldid = list(sel_ids['field']) if len(sel_ids['field']) > 0 else -1 

592 baseline = self.format_ac_baseline(sel_ids['antenna1']) 

593 scanid = list(sel_ids['scan']) if len(sel_ids['scan']) > 0 else "" 

594 # SPW (need to get a list of valid spws instead of -1) 

595 if len(sel_ids['channel']) > 0: 

596 spwid = [chanarr[0] for chanarr in sel_ids['channel']] 

597 elif spw == "": # No spw selection 

598 my_ms.open(vis) 

599 try: 

600 spwinfo = my_ms.getspectralwindowinfo() 

601 except Exception: 

602 raise 

603 finally: 

604 my_ms.close() 

605 

606 spwid = [int(idx) for idx in spwinfo.keys()] 

607 else: 

608 raise RuntimeError("Invalid spw selction, %s ,for MS %d" % (str(spw), file_idx)) 

609 

610 return {'field': fieldid, 'spw': spwid, 'baseline': baseline, 'scan': scanid, 

611 'intent': intent, 'antenna1': sel_ids['antenna1']} 

612 else: 

613 raise ValueError("Invalid file index, %d" % file_idx) 

614 

615 def format_ac_baseline(self, in_antenna): 

616 """Format auto-correlation baseline string from antenna idx list.""" 

617 # exact match string 

618 if is_string_type(in_antenna): 

619 # return sdutil.convert_antenna_spec_autocorr(in_antenna) 

620 if (len(in_antenna) != 0) and (in_antenna.find('&') == -1) \ 

621 and (in_antenna.find(';') == -1): 

622 in_antenna = + '&&&' 

623 return in_antenna 

624 # single integer -> list of int 

625 if type(in_antenna) == int: 

626 if in_antenna >= 0: 

627 in_antenna = [in_antenna] 

628 else: 

629 return -1 

630 # format auto-corr string from antenna idices. 

631 baseline = '' 

632 for idx in in_antenna: 

633 if len(baseline) > 0: 

634 baseline += ';' 

635 if idx >= 0: 

636 baseline += (str(idx) + '&&&') 

637 return baseline 

638 

639 def compile(self): 

640 # imaging mode 

641 self.imager_param['mode'] = self.mode 

642 

643 # Work on selection of the first table in sorted list 

644 # to get default restfreq and outframe 

645 # chronological sort 

646 sorted_vislist, sorted_timelist = mslisthelper.sort_mslist(self.infiles) 

647 self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist] 

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

649 # conform MS 

650 conform_mslist(sorted_vislist, ignore_columns=[]) 

651 selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0]) 

652 self.__update_subtable_name(self.infiles[self.sorted_idx[0]]) 

653 # field 

654 fieldid = selection_ids['field'][0] if type( 

655 selection_ids['field']) != int else selection_ids['field'] 

656 sourceid = -1 

657 self.open_table(self.field_table) 

658 source_ids = self.table.getcol('SOURCE_ID') 

659 self.close_table() 

660 if self.field == '' or fieldid == -1: 

661 sourceid = source_ids[0] 

662 elif fieldid >= 0 and fieldid < len(source_ids): 

663 sourceid = source_ids[fieldid] 

664 else: 

665 raise ValueError("No valid field in the first MS.") 

666 

667 # restfreq 

668 if self.restfreq == '' and self.source_table != '': 

669 self.open_table(self.source_table) 

670 source_ids = self.table.getcol('SOURCE_ID') 

671 for i in range(self.table.nrows()): 

672 if sourceid == source_ids[i] \ 

673 and self.table.iscelldefined('REST_FREQUENCY', i) \ 

674 and (selection_ids['spw'] == -1 or 

675 self.table.getcell('SPECTRAL_WINDOW_ID', i) in selection_ids['spw']): 

676 rf = self.table.getcell('REST_FREQUENCY', i) 

677 if len(rf) > 0: 

678 self.restfreq = self.table.getcell('REST_FREQUENCY', i)[0] 

679 break 

680 self.close_table() 

681 casalog.post("restfreq set to %s" % self.restfreq, "INFO") 

682 # REST_FREQUENCY column is optional (need retry if not exists) 

683 self.imager_param['restfreq'] = self.restfreq 

684 

685 # 

686 # spw (define representative spw id = spwid_ref) 

687 spwid_ref = selection_ids['spw'][0] \ 

688 if type(selection_ids['spw']) != int else selection_ids['spw'] 

689 # Invalid spw selection should have handled at msselectiontoindex(). 

690 # -1 means all spw are selected. 

691 self.open_table(self.spw_table) 

692 if spwid_ref < 0: 

693 for id in range(self.table.nrows()): 

694 if self.table.getcell('NUM_CHAN', id) > 0: 

695 spwid_ref = id 

696 break 

697 if spwid_ref < 0: 

698 self.close_table() 

699 msg = 'No valid spw id exists in the first table' 

700 raise ValueError(msg) 

701 self.allchannels = self.table.getcell('NUM_CHAN', spwid_ref) 

702 # freq_chan0 = self.table.getcell('CHAN_FREQ', spwid_ref)[0] 

703 # freq_inc0 = self.table.getcell('CHAN_WIDTH', spwid_ref)[0] 

704 # in case rest frequency is not defined yet. 

705 if self.restfreq == '': 

706 self.restfreq = '%fHz' % self.table.getcell('CHAN_FREQ', spwid_ref).mean() 

707 self.imager_param['restfreq'] = self.restfreq 

708 casalog.post("Using mean freq of spw %d as restfreq: %s" % 

709 (spwid_ref, self.restfreq), "INFO") 

710 self.close_table() 

711 self.imager_param['spw'] = -1 # spwid_ref 

712 

713 # outframe (force using the current frame) 

714 self.imager_param['outframe'] = self.outframe 

715 if self.outframe == '': 

716 if len(self.infiles) > 1: 

717 # The default will be 'LSRK' 

718 casalog.post("Multiple MS inputs. The default outframe is set to 'LSRK'") 

719 self.imager_param['outframe'] = 'LSRK' 

720 else: 

721 # get from MS 

722 my_ms = ms() 

723 my_ms.open(self.infiles[0]) 

724 spwinfo = my_ms.getspectralwindowinfo() 

725 my_ms.close() 

726 del my_ms 

727 for key, spwval in spwinfo.items(): 

728 if spwval['SpectralWindowId'] == spwid_ref: 

729 self.imager_param['outframe'] = spwval['Frame'] 

730 casalog.post( 

731 f"Using frequency frame of MS, '{self.imager_param['outframe']}'") 

732 break 

733 if self.imager_param['outframe'] == '': 

734 raise Exception("Internal error of getting frequency frame of spw=%d." % spwid_ref) 

735 else: 

736 casalog.post( 

737 f"Using frequency frame defined by user, '{self.imager_param['outframe']}'") 

738 

739 # brightness unit 

740 if len(self.brightnessunit) > 0: 

741 if self.brightnessunit.lower() == 'k': 

742 self.image_unit = 'K' 

743 elif self.brightnessunit.lower() == 'jy/beam': 

744 self.image_unit = 'Jy/beam' 

745 else: 

746 raise ValueError("Invalid brightness unit, %s" % self.brightnessunit) 

747 

748# # antenna 

749# in_antenna = self.antenna # backup for future use 

750# if type(self.antenna)==int: 

751# if self.antenna >= 0: 

752# self.antenna=str(self.antenna)+'&&&' 

753# else: 

754# if (len(self.antenna) != 0) and (self.antenna.find('&') == -1) \ 

755# and (self.antenna.find(';')==-1): 

756# self.antenna = self.antenna + '&&&' 

757 

758 def _configure_map_property(self): 

759 selection_ids = self.get_selection_idx_for_ms(self.sorted_idx[0]) 

760 

761 # stokes 

762 if self.stokes == '': 

763 self.stokes = 'I' 

764 self.imager_param['stokes'] = self.stokes 

765 

766 # gridfunction 

767 

768 # outfile 

769 if os.path.exists(self.outfile) and self.overwrite: 

770 smart_remove(self.outfile) 

771 if os.path.exists(self.outfile + '.weight') and self.overwrite: 

772 smart_remove(self.outfile + '.weight') 

773 

774 # cell 

775 cell = self.cell 

776 if cell == '' or cell[0] == '': 

777 # Calc PB 

778 grid_factor = 3. 

779 casalog.post( 

780 "The cell size will be calculated using PB size of antennas in the first MS") 

781 qPB = self._calc_PB(selection_ids['antenna1']) 

782 cell = '%f%s' % (qPB['value'] / grid_factor, qPB['unit']) 

783 casalog.post("Using cell size = PB/%4.2F = %s" % (grid_factor, cell)) 

784 

785 (cellx, celly) = sdutil.get_cellx_celly(cell, unit='arcmin') 

786 self.imager_param['cellx'] = cellx 

787 self.imager_param['celly'] = celly 

788 

789 # Calculate Pointing center and extent (if necessary) 

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

791 imsize = sdutil.to_list(self.imsize, int) or \ 

792 sdutil.to_list(self.imsize, numpy.integer) 

793 if imsize is None: 

794 imsize = self.imsize if hasattr(self.imsize, '__iter__') else [self.imsize] 

795 imsize = [int(numpy.ceil(v)) for v in imsize] 

796 casalog.post( 

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

798 priority="WARN") 

799 casalog.post("rounded-up imsize: %s --> %s" % (str(self.imsize), str(imsize))) 

800 

801 phasecenter = self.phasecenter 

802 if self.phasecenter == "" or \ 

803 len(imsize) == 0 or imsize[0] < 1: 

804 map_param = self._get_pointing_extent() 

805 # imsize 

806 if len(imsize) == 0 or imsize[0] < 1: 

807 imsize = self._get_imsize(map_param['width'], map_param['height'], cellx, celly) 

808 if self.phasecenter != "": 

809 casalog.post( 

810 "You defined phasecenter but not imsize. " 

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

812 "but be centered at phasecenter. " 

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

814 "apart from the center of pointings", 

815 priority='WARN') 

816 if imsize[0] > 1024 or imsize[1] > 1024: 

817 casalog.post( 

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

819 "It could take time to generate the image depending on " 

820 "your computer resource. Please wait...", 

821 priority='WARN') 

822 

823 # phasecenter 

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

825 if self.phasecenter == "": 

826 phasecenter = map_param['center'] 

827 

828 # imsize 

829 (nx, ny) = sdutil.get_nx_ny(imsize) 

830 self.imager_param['nx'] = nx 

831 self.imager_param['ny'] = ny 

832 

833 # phasecenter 

834 self.imager_param['phasecenter'] = phasecenter 

835 

836 self.imager_param['movingsource'] = self.ephemsrcname 

837 

838 # channel map 

839 sorted_vislist, _ = mslisthelper.sort_mslist(self.infiles) 

840 self.sorted_idx = [self.infiles.index(vis) for vis in sorted_vislist] 

841 imhelper = cleanhelper_minimal(self.imager, self.infiles, sort_index=self.sorted_idx, casalog=casalog) 

842 spwsel = str(',').join([str(spwid) for spwid in selection_ids['spw']]) 

843 srestf = self.imager_param['restfreq'] \ 

844 if is_string_type(self.imager_param['restfreq']) \ 

845 else "%fHz" % self.imager_param['restfreq'] 

846 (imnchan, imstart, imwidth) = imhelper.setChannelizeDefault( 

847 self.mode, spwsel, self.field, self.nchan, 

848 self.start, self.width, self.imager_param['outframe'], 

849 self.veltype, self.imager_param['phasecenter'], srestf) 

850 del imhelper 

851 

852 # start and width 

853 if self.mode == 'velocity': 

854 startval = [self.imager_param['outframe'], imstart] 

855 widthval = imwidth 

856 elif self.mode == 'frequency': 

857 startval = [self.imager_param['outframe'], imstart] 

858 widthval = imwidth 

859 else: # self.mode==channel 

860 startval = int(self.start) 

861 widthval = int(self.width) 

862 

863 if self.nchan < 0: 

864 self.nchan = self.allchannels 

865 self.imager_param['start'] = startval 

866 self.imager_param['step'] = widthval 

867 self.imager_param['nchan'] = imnchan # self.nchan 

868 

869 self.imager_param['projection'] = self.projection 

870 

871 def execute(self): 

872 # imaging 

873 casalog.post("Start imaging...", "INFO") 

874 if len(self.infiles) == 1: 

875 self.open_imager(self.infiles[0]) 

876 selection_ids = self.get_selection_idx_for_ms(0) 

877 spwsel = self.get_selection_param_for_ms(0, self.spw) 

878 if spwsel.strip() in ['', '*']: 

879 spwsel = selection_ids['spw'] 

880 # TODO: channel selection based on spw 

881 ok = self.imager.selectvis(field=selection_ids['field'], 

882 # spw=selection_ids['spw'], 

883 spw=spwsel, 

884 nchan=-1, start=0, step=1, 

885 baseline=selection_ids['baseline'], 

886 scan=selection_ids['scan'], 

887 intent=selection_ids['intent']) 

888 if not ok: 

889 raise ValueError("Selection is empty: you may want to review this MS selection") 

890 else: 

891 self.close_imager() 

892 # self.sorted_idx.reverse() 

893 for idx in self.sorted_idx.__reversed__(): 

894 name = self.infiles[idx] 

895 selection_ids = self.get_selection_idx_for_ms(idx) 

896 spwsel = self.get_selection_param_for_ms(idx, self.spw) 

897 if spwsel.strip() in ['', '*']: 

898 spwsel = selection_ids['spw'] 

899 # TODO: channel selection based on spw 

900 self.imager.selectvis(vis=name, field=selection_ids['field'], 

901 # spw=selection_ids['spw'], 

902 spw=spwsel, 

903 nchan=-1, start=0, step=1, 

904 baseline=selection_ids['baseline'], 

905 scan=selection_ids['scan'], 

906 intent=selection_ids['intent']) 

907 # need to do this 

908 self.is_imager_opened = True 

909 

910 # it should be called after infiles are registered to imager 

911 self._configure_map_property() 

912 

913 casalog.post(f"Using phasecenter {self.imager_param['phasecenter']}", "INFO") 

914 

915 self.imager.defineimage(**self.imager_param) # self.__get_param()) 

916 self.imager.setoptions(ftmachine='sd', gridfunction=self.gridfunction, freqinterp=self.interpolation) 

917 self.imager.setsdoptions( 

918 pointingcolumntouse=self.pointingcolumn, 

919 convsupport=self.convsupport, 

920 truncate=self.truncate, 

921 gwidth=self.gwidth, 

922 jwidth=self.jwidth, 

923 minweight=0., 

924 clipminmax=self.clipminmax, 

925 enablecache=self.enablecache, 

926 convertfirst=self.convertfirst 

927 ) 

928 

929 # Create images 

930 imgtype_suffix = {'singledish': '', 'coverage': '.weight'} 

931 for img_type, img_suffix in imgtype_suffix.items(): 

932 img_file = self.outfile + img_suffix 

933 msg_fmt = string.Template(f"$state {img_type} image {img_file}") 

934 casalog.post(msg_fmt.substitute(state="Generating"), "INFO") 

935 self.imager.makeimage(type=img_type, image=img_file) 

936 if not os.path.exists(img_file): 

937 raise RuntimeError(msg_fmt.substitute(state="Failed to generate")) 

938 casalog.post(msg_fmt.substitute(state="Generated"), "INFO") 

939 # Close imager 

940 self.close_imager() 

941 

942 # Convert output images to proper output frame and set brightness unit (if necessary) 

943 my_ia = image() 

944 my_ia.open(self.outfile) 

945 csys = my_ia.coordsys() 

946 csys.setconversiontype(spectral=csys.referencecode('spectra')[0]) 

947 my_ia.setcoordsys(csys.torecord()) 

948 if len(self.image_unit) > 0: 

949 casalog.post("Setting brightness unit '%s' to image." % self.image_unit) 

950 my_ia.setbrightnessunit(self.image_unit) 

951 csys.done() 

952 my_ia.close() 

953 

954 # CAS-12984 set brightness unit for weight image to '' 

955 weightfile = self.outfile + imgtype_suffix['coverage'] 

956 my_ia.open(weightfile) 

957 my_ia.setbrightnessunit('') 

958 my_ia.close() 

959 

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

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

962 casalog.post("Start masking the map using minweight = %f" % 

963 self.minweight, "INFO") 

964 my_ia.open(weightfile) 

965 try: 

966 stat = my_ia.statistics(mask="'" + weightfile + "' > 0.0", robust=True) 

967 valid_pixels = stat['npts'] 

968 except RuntimeError as e: 

969 if str(e).find('No valid data found.') >= 0: 

970 valid_pixels = [0] 

971 else: 

972 raise 

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

974 my_ia.close() 

975 casalog.post( 

976 "All pixels weight zero. This indicates no data in MS is in image area. " 

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

978 priority="WARN") 

979 return 

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

981 weight_threshold = median_weight * self.minweight 

982 casalog.post("Median of weight in the map is %f" % median_weight, 

983 "INFO") 

984 casalog.post("Pixels in map with weight <= median(weight)*minweight = %f will be masked." % 

985 (weight_threshold), "INFO") 

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

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

988 # if one wanted to find how many pixel were masked one could easily count the 

989 # number of pixels set to false 

990 # e.g after masking self.outfile below one could just do this 

991 # nmasked_pixels=tb.calc( 

992 # '[select from "'+self.outfile+'"/mask0'+'" giving [nfalse(PagedArray)]]') 

993 my_tb = table() 

994 nmask_pixels = 0 

995 nchan = stat['trc'][3] + 1 

996 casalog.filter('ERROR') # hide the useless message of tb.calc 

997 

998 # doing it by channel to make sure it does not go out of memory 

999 # tab.calc try to load the whole chunk in ram 

1000 for k in range(nchan): 

1001 nmask_pixels += my_tb.calc( 

1002 '[select from "' + weightfile + '" giving [ntrue(map[,,,' + 

1003 str(k) + '] <=' + str(median_weight * self.minweight) + ')]]')['0'][0] 

1004 casalog.filter() # set logging back to normal 

1005 

1006 casalog.filter() # set logging back to normal 

1007 imsize = numpy.prod(my_ia.shape()) 

1008 my_ia.close() 

1009 # Modify default mask 

1010 my_ia.open(self.outfile) 

1011 my_ia.calcmask("'%s'>%f" % (weightfile, weight_threshold), asdefault=True) 

1012 my_ia.close() 

1013 masked_fraction = 100. * (1. - (imsize - nmask_pixels) / float(valid_pixels[0])) 

1014 casalog.post("This amounts to %5.1f %% of the area with nonzero weight." % 

1015 (masked_fraction), "INFO") 

1016 casalog.post( 

1017 f"The weight image '{weightfile}' is returned by this task, " 

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

1019 priority="INFO") 

1020 

1021 # Calculate theoretical beam size 

1022 casalog.post("Calculating image beam size.") 

1023 if self.gridfunction.upper() not in ['SF']: 

1024 casalog.post( 

1025 f"Beam size definition for '{self.gridfunction}' kernel is experimental.", 

1026 priority='WARN') 

1027 casalog.post( 

1028 "You may want to take careful look at the restoring beam in the image.", 

1029 priority='WARN') 

1030 my_msmd = msmetadata() 

1031 # antenna diameter and blockage 

1032 ref_ms_idx = self.sorted_idx[0] 

1033 ref_ms_name = self.infiles[ref_ms_idx] 

1034 selection_ids = self.get_selection_idx_for_ms(ref_ms_idx) 

1035 ant_idx = selection_ids['antenna1'] 

1036 diameter = self._get_average_antenna_diameter(ant_idx) 

1037 my_msmd.open(ref_ms_name) 

1038 ant_name = my_msmd.antennanames(ant_idx) 

1039 my_msmd.close() 

1040 is_alma = False 

1041 for name in ant_name: 

1042 if name[0:2] in ["PM", "DV", "DA", "CM"]: 

1043 is_alma = True 

1044 break 

1045 blockage = "0.75m" if is_alma else "0.0m" # unknown blockage diameter 

1046 # output reference code 

1047 my_ia.open(self.outfile) 

1048 csys = my_ia.coordsys() 

1049 my_ia.close() 

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

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

1052 csys.done() 

1053 # pointing sampling 

1054 ref_ms_spw = self.get_selection_param_for_ms(ref_ms_idx, self.spw) 

1055 ref_ms_field = self.get_selection_param_for_ms(ref_ms_idx, self.field) 

1056 ref_ms_scan = self.get_selection_param_for_ms(ref_ms_idx, self.scanno) 

1057 # xSampling, ySampling, angle = sdutil.get_ms_sampling_arcsec(ref_ms_name, spw=ref_ms_spw, 

1058 # antenna=selection_ids['baseline'], 

1059 # field=ref_ms_field, 

1060 # scan=ref_ms_scan,#timerange='', 

1061 # outref=outref) 

1062 

1063 # obtain sampling interval for beam calculation. 

1064 self.open_imager(ref_ms_name) 

1065 ok = self.imager.selectvis(field=ref_ms_field, 

1066 # spw=selection_ids['spw'], 

1067 spw=ref_ms_spw, 

1068 nchan=-1, start=0, step=1, 

1069 baseline=selection_ids['baseline'], 

1070 scan=ref_ms_scan, 

1071 intent=selection_ids['intent']) 

1072 if len(ant_idx) > 1: 

1073 casalog.post("Using only antenna %s to calculate sampling interval" % ant_name[0]) 

1074 ptg_samp = self.imager.pointingsampling(pattern='raster', ref=outref, 

1075 movingsource=self.ephemsrcname, 

1076 pointingcolumntouse=self.pointingcolumn, 

1077 antenna=('%s&&&' % ant_name[0])) 

1078 self.close_imager() 

1079 my_qa = quanta() 

1080 xSampling, ySampling = my_qa.getvalue(my_qa.convert(ptg_samp['sampling'], 'arcsec')) 

1081 angle = my_qa.getvalue(my_qa.convert(ptg_samp['angle'], "deg"))[0] 

1082 

1083 casalog.post("Detected raster sampling = [{x:f}, {y:f}] arcsec".format( 

1084 x=xSampling, 

1085 y=ySampling 

1086 )) 

1087 

1088 # handling of failed sampling detection 

1089 valid_sampling = True 

1090 sampling = [xSampling, ySampling] 

1091 if abs(xSampling) < 2.2e-3 or not numpy.isfinite(xSampling): 

1092 casalog.post( 

1093 f"Invalid sampling={xSampling} arcsec. " 

1094 f"Using the value of orthogonal direction={ySampling} arcsec", 

1095 priority="WARN") 

1096 sampling = [ySampling] 

1097 angle = 0.0 

1098 valid_sampling = False 

1099 if abs(ySampling) < 1.0e-3 or not numpy.isfinite(ySampling): 

1100 if valid_sampling: 

1101 casalog.post( 

1102 f"Invalid sampling={ySampling} arcsec. " 

1103 f"Using the value of orthogonal direction={xSampling} arcsec", 

1104 priority="WARN") 

1105 sampling = [xSampling] 

1106 angle = 0.0 

1107 valid_sampling = True 

1108 # reduce sampling and cell if it's possible 

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

1110 sampling = [sampling[0]] 

1111 angle = 0.0 

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

1113 cell = [cell[0]] 

1114 if valid_sampling: 

1115 # actual calculation of beam size 

1116 bu = sdbeamutil.TheoreticalBeam() 

1117 bu.set_antenna(diameter, blockage) 

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

1119 bu.set_image_param(cell, self.restfreq, self.gridfunction, 

1120 self.convsupport, self.truncate, self.gwidth, 

1121 self.jwidth, is_alma) 

1122 bu.summary() 

1123 imbeam_dict = bu.get_beamsize_image() 

1124 casalog.post("Setting image beam: major={major}, minor={minor}, pa={pa}".format( 

1125 major=imbeam_dict['major'], 

1126 minor=imbeam_dict['minor'], 

1127 pa=imbeam_dict['pa'] 

1128 )) 

1129 # set beam size to image 

1130 my_ia.open(self.outfile) 

1131 my_ia.setrestoringbeam(**imbeam_dict) 

1132 my_ia.close() 

1133 else: 

1134 # BOTH sampling was invalid 

1135 casalog.post( 

1136 "Could not detect valid raster sampling. " 

1137 "Exitting without setting beam size to image", 

1138 priority='WARN') 

1139 

1140 def _calc_PB(self, antenna): 

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

1142 

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

1144 and rest frequency 

1145 Average antenna diamter and reference frequency are adopted for 

1146 calculation. 

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

1148 """ 

1149 casalog.post("Calculating Primary beam size:") 

1150 # CAS-5410 Use private tools inside task scripts 

1151 my_qa = quanta() 

1152 

1153 pb_factor = 1.175 

1154 # Reference frequency 

1155 ref_freq = self.restfreq 

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

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

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

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

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

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

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

1163 raise Exception(msg) 

1164 # Antenna diameter 

1165 antdiam_ave = self._get_average_antenna_diameter(antenna) 

1166 # Calculate PB 

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

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

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

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

1171 # Summary 

1172 casalog.post("- Antenna diameter: %s m" % D_m) 

1173 casalog.post("- Reference Frequency: %s" % ref_freq) 

1174 casalog.post("PB size = %5.3f * lambda/D = %s" % (pb_factor, my_qa.tos(PB))) 

1175 return PB 

1176 

1177 def _get_imsize(self, width, height, dx, dy): 

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

1179 # CAS-5410 Use private tools inside task scripts 

1180 my_qa = quanta() 

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

1182 my_qa.getvalue(dy))) 

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

1184 my_qa.getvalue(dx))) 

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

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

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

1188 (nx + 1, ny + 1)) 

1189 return (int(nx + 1), int(ny + 1)) 

1190 

1191 def _get_pointing_extent(self): 

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

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

1194 # CAS-5410 Use private tools inside task scripts 

1195 my_qa = quanta() 

1196 ret_dict = {} 

1197 

1198 colname = self.pointingcolumn.upper() 

1199 

1200 # MSs should be registered to imager 

1201 if not self.is_imager_opened: 

1202 raise RuntimeError('Internal error: imager should be opened here.') 

1203 

1204 if self.phasecenter == "": 

1205 # defaut is J2000 

1206 base_mref = 'J2000' 

1207 elif isinstance(self.phasecenter, int) or self.phasecenter.isdigit(): 

1208 # may be field id 

1209 self.open_table(self.field_table) 

1210 base_mref = self.table.getcolkeyword('PHASE_DIR', 'MEASINFO')['Ref'] 

1211 self.close_table() 

1212 else: 

1213 # may be phasecenter is explicitly specified 

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

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

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

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

1218 # DMS string: 9.15.29, -9d15m29s 

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

1220 # composite pattern 

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

1222 items = self.phasecenter.split() 

1223 base_mref = 'J2000' 

1224 for i in items: 

1225 s = i.strip() 

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

1227 base_mref = s 

1228 break 

1229 

1230 mapextent = self.imager.mapextent(ref=base_mref, movingsource=self.ephemsrcname, 

1231 pointingcolumntouse=colname) 

1232 if mapextent['status']: 

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

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

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

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

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

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

1239 

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

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

1242 my_qa.tos(qheight))) 

1243 ret_dict['center'] = scenter 

1244 ret_dict['width'] = qwidth 

1245 ret_dict['height'] = qheight 

1246 else: 

1247 casalog.post( 

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

1249 'due to missing valid data.', 

1250 priority='SEVERE') 

1251 ret_dict['center'] = '' 

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

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

1254 return ret_dict 

1255 

1256 def _get_x_minmax(self, x): 

1257 # assumed the x is in unit of rad. 

1258 pi2 = 2. * numpy.pi 

1259 x = (x % pi2) 

1260 npart = 4 

1261 dlon = pi2 / float(npart) 

1262 pos = [int(v / dlon) for v in x] 

1263 voids = [False for dummy in range(npart)] 

1264 for ipos in range(npart): 

1265 try: 

1266 pos.index(ipos) 

1267 except Exception: 

1268 voids[ipos] = True 

1269 if not any(voids): 

1270 raise Exception( 

1271 "Failed to find global pointing gap. " 

1272 f"The algorithm requires at least 2PI/{npart} of pointing gap") 

1273 rot_pos = [] 

1274 if (not voids[0]) and (not voids[npart - 1]): 

1275 gmax = -1 

1276 for idx in range(npart - 2, 0, -1): 

1277 if voids[idx]: 

1278 gmax = idx 

1279 break 

1280 if gmax < 0: 

1281 raise Exception("Failed to detect gap max") 

1282 rot_pos = range(gmax + 1, npart) 

1283 for idx in range(len(x)): 

1284 x[idx] = (x[idx] - pi2) if pos[idx] in rot_pos else x[idx] 

1285 

1286 return (x.min(), x.max()) 

1287 

1288 def __update_subtable_name(self, msname): 

1289 self.open_table(msname) 

1290 keys = self.table.getkeywords() 

1291 self.close_table() 

1292 self.field_table = sdutil.get_subtable_name(keys['FIELD']) 

1293 self.spw_table = sdutil.get_subtable_name(keys['SPECTRAL_WINDOW']) 

1294 self.source_table = sdutil.get_subtable_name(keys['SOURCE']) 

1295 self.antenna_table = sdutil.get_subtable_name(keys['ANTENNA']) 

1296 self.polarization_table = sdutil.get_subtable_name(keys['POLARIZATION']) 

1297 self.observation_table = sdutil.get_subtable_name(keys['OBSERVATION']) 

1298 self.pointing_table = sdutil.get_subtable_name(keys['POINTING']) 

1299 self.data_desc_table = sdutil.get_subtable_name(keys['DATA_DESCRIPTION']) 

1300 

1301 def _get_average_antenna_diameter(self, antenna): 

1302 my_qa = quanta() 

1303 self.open_table(self.antenna_table) 

1304 try: 

1305 antdiam_unit = self.table.getcolkeyword('DISH_DIAMETER', 'QuantumUnits')[0] 

1306 diams = self.table.getcol('DISH_DIAMETER') 

1307 finally: 

1308 self.close_table() 

1309 

1310 if len(antenna) == 0: 

1311 antdiam_ave = my_qa.quantity(diams.mean(), antdiam_unit) 

1312 else: 

1313 d_ave = sum([diams[idx] for idx in antenna]) / float(len(antenna)) 

1314 antdiam_ave = my_qa.quantity(d_ave, antdiam_unit) 

1315 return antdiam_ave