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

522 statements  

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

1 

2from scipy import fftpack 

3import numpy as np 

4import shutil 

5import os 

6import time 

7 

8from casatools import quanta, table, image, regionmanager, imager 

9from casatasks import casalog, imsubimage, feather 

10 

11_ia = image() 

12_qa = quanta() 

13_rg = regionmanager() 

14 

15_mytb = table() 

16 

17class SDINT_helper: 

18 

19# def __init__(self): 

20# casalog.post('Init Helper') 

21 

22################################################ 

23 def getFreqAxisIndex(self): 

24 try: 

25 mysummary = _ia.summary(list=False) 

26 freqaxis_index = list(mysummary['axisnames']).index('Frequency') 

27 except(ValueError): 

28 _ia.close() 

29 casalog.post('The image '+_ia.name()+' has no frequency axis. Try adding one with ia.adddegaxis() .', 'SEVERE') 

30 

31 return freqaxis_index 

32 

33 def getFreqList(self,imname=''): 

34 """ Get the list of frequencies for the given image, one for each channel. 

35 

36 Returns: 

37 list[float] The frequencies for each channel in the image, in Hz. 

38 """ 

39 

40 _ia.open(imname) 

41 csys =_ia.coordsys() 

42 shp = _ia.shape() 

43 freqaxis_index = self.getFreqAxisIndex() 

44 _ia.close() 

45 

46 if(csys.axiscoordinatetypes()[freqaxis_index] == 'Spectral'): 

47 restfreq = csys.referencevalue()['numeric'][freqaxis_index]#/1.0e+09; # convert more generally.. 

48 freqincrement = csys.increment()['numeric'][freqaxis_index]# /1.0e+09; 

49 freqlist = []; 

50 for chan in range(0,shp[freqaxis_index]): 

51 freqlist.append(restfreq + chan * freqincrement); 

52 elif(csys.axiscoordinatetypes()[freqaxis_index] == 'Tabular'): 

53 freqlist = (csys.torecord()['tabular2']['worldvalues']) # /1.0e+09; 

54 else: 

55 casalog.post('Unknown frequency axis. Exiting.','SEVERE'); 

56 return False; 

57 

58 csys.done() 

59 return freqlist 

60 

61################################################ 

62 

63 def copy_restoringbeam(self,fromthis='',tothis=''): 

64 

65 freqlist = self.getFreqList(fromthis) 

66 # casalog.post(freqlist) 

67 for i in range(len(freqlist)): 

68 _ia.open(fromthis); 

69 beam = _ia.restoringbeam(channel = i); 

70 _ia.close() 

71 # casalog.post(beam) 

72 _ia.open(tothis) 

73 _ia.setrestoringbeam(beam = beam, channel = i, polarization = 0); 

74 _ia.close() 

75 

76################################################ 

77 

78 def feather_int_sd(self, sdcube='', intcube='', jointcube='',sdgain=1.0, dishdia=-1, usedata='sdint', chanwt=''): 

79 """ 

80 Run the feather task to combine the SD and INT Cubes.  

81  

82 There's a bug in feather for cubes. Hence, do each channel separately. 

83 FIX feather and then change this. CAS-5883 is the JIRA ticket that contains a fix for this issue....  

84 

85 TODO : Add the effdishdia usage to get freq-indep feathering. 

86 

87 """ 

88 

89 ### Do the feathering. 

90 if usedata=='sdint': 

91 ## Feather runs in a loop on chans internally, but there are issues with open tablecache images 

92 ## Also, no way to set effective dish dia separately for each channel. 

93 #feather(imagename = jointcube, highres = intcube, lowres = sdcube, sdfactor = sdgain, effdishdiam=-1) 

94 

95 freqlist = self.getFreqList(sdcube) 

96 

97 os.system('rm -rf '+jointcube) 

98 os.system('cp -r ' + intcube + ' ' + jointcube) 

99 

100 _ib = image() 

101 

102 _ia.open(jointcube) 

103 _ia.set(0.0) ## Initialize this to zero for all planes 

104 

105 for i in range(len(freqlist)): 

106 if chanwt[i] != 0.0 : ## process only the nonzero channels 

107 if(dishdia <=0): 

108 casalog.post('Parameter dishdia (SD dish diameter in meters) must be > 0.', 'SEVERE') 

109 

110 freqdishdia = dishdia ## * (freqlist[0] / freqlist[i]) # * 0.5 

111 

112 os.system('rm -rf tmp_*') 

113 #imsubimage(imagename = sdcube, outfile = 'tmp_sdplane', chans = str(i)); 

114 #imsubimage(imagename = intcube, outfile = 'tmp_intplane', chans = str(i)); 

115 self.createplaneimage(imagename=sdcube, outfile='tmp_sdplane', chanid = str(i)); 

116 self.createplaneimage(imagename=intcube, outfile='tmp_intplane', chanid = str(i)); 

117 

118 #feather(imagename = 'tmp_jointplane', highres = 'tmp_intplane', lowres = 'tmp_sdplane', sdfactor = sdgain, effdishdiam=freqdishdia) 

119 # feathering via toolkit 

120 try: 

121 casalog.post("start Feathering.....") 

122 imFea=imager( ) 

123 imFea.setvp(dovp=True) 

124 imFea.setsdoptions(scale=sdgain) 

125 imFea.feather(image='tmp_jointplane',highres='tmp_intplane',lowres='tmp_sdplane', effdishdiam=freqdishdia) 

126 imFea.done( ) 

127 del imFea 

128 except Exception as instance: 

129 casalog.post('*** Error *** %s' % instance, 'ERROR') 

130 raise 

131 

132 _ib.open('tmp_jointplane') 

133 pixjoint = _ib.getchunk() 

134 _ib.close() 

135 _ia.putchunk(pixjoint, blc=[0,0,0,i]) 

136 

137 _ia.close() 

138 

139 if usedata=='sd': 

140 ## Copy sdcube to joint. 

141 os.system('rm -rf '+jointcube) 

142 os.system('cp -r ' + sdcube + ' ' + jointcube) 

143 if usedata=='int': 

144 ## Copy intcube to joint 

145 os.system('rm -rf '+jointcube) 

146 os.system('cp -r ' + intcube + ' ' + jointcube) 

147 

148################################################ 

149 def calc_renorm(self, intname='', jointname=''): 

150 """ 

151 Calculate a new .sumwt spectrum containing the peak of the feathered PSF. 

152 The PSF and each residual image calculation will be re-normalized by this. 

153 This will keep the PSFs in all channels at a peak of 1. 

154 """ 

155 psfname = jointname+'.psf' 

156 os.system('cp -r '+intname+'.sumwt ' + jointname + '.sumwt') 

157 _ia.open(jointname+'.sumwt') 

158 vals = _ia.getchunk() 

159 shp = _ia.shape() 

160 freqaxis_index = self.getFreqAxisIndex() 

161 _ia.close() 

162 

163 if shp[0]>1: 

164 casalog.post("WARNING : Cannot use this task with faceting", 'WARN') 

165 

166 _ia.open(jointname+'.psf') 

167 for i in range(0, shp[freqaxis_index]): 

168 onepsf = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

169 vals[0,0,0,i] = np.max(onepsf) 

170 _ia.close() 

171 

172 _ia.open(jointname+'.sumwt') 

173 _ia.putchunk(vals) 

174 _ia.close() 

175 

176 casalog.post("********************Re-norm with "+str(vals)) 

177 

178 

179################################################ 

180 

181 def apply_renorm(self, imname='', sumwtname=''): 

182 """ 

183 Divide each plane of the input image by the sumwt value for that channel 

184 """ 

185 _ia.open(sumwtname) 

186 shp = _ia.shape() 

187 freqaxis_index = self.getFreqAxisIndex() 

188 vals = _ia.getchunk() ## This is one pixel per channel. 

189 _ia.close() 

190 

191 casalog.post("********************Re-norm with "+str(vals)) 

192 

193 _ia.open(imname) 

194 for i in range(0, shp[freqaxis_index]): 

195 oneplane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

196 if vals[0,0,0,i]>0.0: 

197 normplane = oneplane/vals[0,0,0,i] 

198 else: 

199 normplane = oneplane.copy() 

200 normplane.fill(0.0) 

201 _ia.putchunk( normplane , blc=[0,0,0,i] ) 

202 _ia.close() 

203 

204 

205 

206################################################ 

207 

208 def modify_with_pb(self, inpcube='', pbcube='',cubewt='', chanwt='', action='mult',pblimit=0.2, freqdep=True): 

209 """ 

210 Multiply or divide by the PB 

211 

212 Args: 

213 inpcube: The cube to be modified. For example: "try.int.cube.model" 

214 pbcube: The primary beam to multiply/divide by. For example: "try.int.cube.pb" 

215 cubewt: The per-channel weight of the inpcube. For example: "try.int.cube.sumwt" 

216 chanwt: List of 0s and 1s, one per channel, to effectively disable the effect of a channel on the resulting images. 

217 action: 'mult' or 'div', to multiply by the PB or divide by it. 

218 pblimit: For pixels less than this value in the PB, set those same pixels in the inpcube to zero. 

219 freqdep: True for channel by channel, False to use a freq-independent PB from the middle of the list before/after deconvolution 

220 """ 

221 casalog.post('Modify with PB : ' + action + ' with frequency dependence ' + str(freqdep)) 

222 

223 freqlist = self.getFreqList(inpcube) 

224 

225 _ia.open(inpcube) 

226 shp=_ia.shape() 

227 freqaxis_index = self.getFreqAxisIndex() 

228 _ia.close() 

229 

230 if(freqaxis_index!=3): 

231 casalog.post('The Frequency axis index of '+inpcube+' is '+str(freqaxis_index)+' but modify_with_pb requires index 3.', 'SEVERE') 

232 

233 ############## 

234 ### Calculate a reference Primary Beam 

235 ### Weighted sum of pb cube 

236 

237 ##midchan = int(len(freqlist)/2) 

238 ##refchan = len(freqlist)-1 ## This assumes ascending frequency ordering in chans. 

239 refchan=0 

240 _ia.open(pbcube) 

241 pbplane = _ia.getchunk(blc=[0,0,0,refchan],trc=[shp[0],shp[1],0,refchan]) 

242 _ia.close() 

243 pbplane.fill(0.0) 

244 

245# pbplane = np.zeros( (shp[0],shp[1]), 'float') 

246 

247 if freqdep==False: 

248 _ia.open(cubewt) # .sumwt 

249 cwt = _ia.getchunk()[0,0,0,:] 

250 _ia.close() 

251 

252 if shp[3] != len(cwt) or len(freqlist) != len(cwt): 

253 raise Exception("Modify with PB : Nchan shape mismatch between cube and sumwt.") 

254 

255 cwt = cwt * chanwt ## Merge the weights and flags 

256 

257 sumchanwt = np.sum(cwt) 

258 

259 if sumchanwt==0: 

260 raise Exception("Weights are all zero ! ") 

261 

262 for i in range(len(freqlist)): 

263 ## Read the pb per plane 

264 _ia.open(pbcube) 

265 pbplane = pbplane + cwt[i] * _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

266 _ia.close() 

267 

268 pbplane = pbplane / sumchanwt 

269 

270 ############## 

271 

272 

273 ## Special-case for setting the PBmask to be same for all freqs 

274 if freqdep==False: 

275 shutil.copytree(pbcube, pbcube+'_tmpcopy') 

276 

277 for i in range(len(freqlist)): 

278 

279 ## Read the pb per plane 

280 if freqdep==True: 

281 _ia.open(pbcube) 

282 pbplane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

283 _ia.close() 

284 

285 ## Make a tmp pbcube with the same pb in all planes. This is for the mask. 

286 if freqdep==False: 

287 _ia.open(pbcube+'_tmpcopy') 

288 _ia.putchunk(pbplane, blc=[0,0,0,i]) 

289 _ia.close() 

290 

291 _ia.open(inpcube) 

292 implane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

293 

294 outplane = pbplane.copy() 

295 outplane.fill(0.0) 

296 

297 if action=='mult': 

298 pbplane[pbplane<pblimit]=0.0 

299 outplane = implane * pbplane 

300 else: 

301 implane[pbplane<pblimit]=0.0 

302 pbplane[pbplane<pblimit]=1.0 

303 outplane = implane / pbplane 

304 

305 _ia.putchunk(outplane, blc=[0,0,0,i]) 

306 _ia.close() 

307 

308# if freqdep==True: 

309# ## Set a mask based on frequency-dependent PB 

310# self.addmask(inpcube,pbcube,pblimit) 

311# else: 

312 if freqdep==False: 

313 ## Set a mask based on the PB in refchan 

314 self.addmask(inpcube,pbcube+'_tmpcopy',pblimit) 

315 shutil.rmtree(pbcube+'_tmpcopy') 

316 

317 

318################################################ 

319 def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='replace'): 

320 #def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='add'): 

321 #def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='old'): 

322 """ 

323 add pb mask: create a new mask called 'pbmask' and set it as a defualt mask  

324 mode: "replace" or "add" 

325 relpalce: create a pbmask based on pblimit without account for the exist mask 

326 add: create a pbmask based on pblimit and merge with existing default mask  

327 """ 

328 _ia.open(inpimage) 

329 defaultmaskname=_ia.maskhandler('default')[0] 

330 allmasknames = _ia.maskhandler('get') 

331 # casalog.post("defaultmaskname=",defaultmaskname) 

332 if mode=='replace': 

333 if defaultmaskname!='' and defaultmaskname!='mask0': 

334 _ia.calcmask(mask='"'+pbimage+'"'+'>'+str(pblimit), name=defaultmaskname); 

335 

336 elif defaultmaskname=='mask0': 

337 if 'pbmask' in allmasknames: 

338 _ia.maskhandler('delete','pbmask') 

339 _ia.calcmask(mask='"'+pbimage+'"'+'>'+str(pblimit), name='pbmask'); 

340 

341 #elif mode=='add': 

342 # After deleting a pixel mask it sometimes leaves it in cache  

343 #  

344 # _ia.open(inpimage) 

345 # if defaultmaskname=='pbmask': 

346 # _ia.maskhandler('delete',defaultmaskname) 

347 # _ia.close() 

348 # _ia.open(inpimage) 

349 

350 # _ia.calcmask(mask='mask("'+inpimage+'")||'+'"'+pbimage+'"'+'>'+str(pblimit), name='pbmask'); 

351 elif mode=='old': 

352 # this one create a new mask every time this function is called! 

353 #_ia.open(inpimage) 

354 _ia.calcmask(mask='mask("'+inpimage+'")||'+'"'+pbimage+'"'+'>'+str(pblimit)); 

355 else: 

356 raise Exception("Unrecongnized value for mode: ",mode) 

357 _ia.close() 

358 _ia.done() 

359 

360################################################ 

361 def cube_to_taylor_sum(self, cubename='', cubewt='', chanwt='', mtname='',reffreq='1.5GHz',nterms=2,dopsf=False): 

362 """ 

363 Convert Cubes (output of major cycle) to Taylor weighted averages (inputs to the minor cycle) 

364 Input : Cube image <cubename>, with channels weighted by image <cubewt> 

365 Output : Set of images : <mtname>.tt0, <mtname>.tt1, etc... 

366 Algorithm: I_ttN = sum([ I_v * ((f-ref)/ref)**N for f in freqs ]) 

367 

368 Args: 

369 cubename: Name of a cube image to interpret into a set of taylor term .ttN images, eg "try.residual", "joint.cube.psf". 

370 cubewt: Name of a .sumwt image that contains the per-channel weighting for the interferometer image. 

371 chanwt: List of 0s and 1s, one per channel, to effectively disable the effect of a channel on the resulting images. 

372 mtname: The prefix output name, to be concatenated with ".ttN" strings, eg "try_mt.residual", "joint.multiterm.psf" 

373 These images should already exist by the time this function is called. 

374 It's suggested that this have same suffix as cubename. 

375 dopsf: Signals that cubename represents a point source function, should be true if cubename ends with ".psf". 

376 If true, then output 2*nterms-1 ttN images. 

377 """ 

378 

379 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value'] 

380 

381 # casalog.post("&&&&&&&&& REF FREQ : " + str(refnu)) 

382 

383 pix=[] 

384 

385 num_terms=nterms 

386 

387 if dopsf==True: 

388 num_terms=2*nterms-1 

389 

390 for tt in range(0,num_terms): 

391 _ia.open(mtname+'.tt'+str(tt)) 

392 pix.append( _ia.getchunk() ) 

393 _ia.close() 

394 pix[tt].fill(0.0) 

395 

396 _ia.open(cubename) 

397 shp = _ia.shape() 

398 freqaxis_index = self.getFreqAxisIndex() 

399 _ia.close() 

400 

401 if(freqaxis_index!=3): 

402 casalog.post('The Frequency axis index of '+cubename+' is '+str(freqaxis_index)+' but cube_to_taylor_sum requires index 3.', 'SEVERE') 

403 

404 _ia.open(cubewt) 

405 cwt = _ia.getchunk()[0,0,0,:] 

406 _ia.close() 

407 

408 

409 freqlist = self.getFreqList(cubename) 

410 

411 if shp[3] != len(cwt) or len(freqlist) != len(cwt): 

412 raise Exception("Nchan shape mismatch between cube and sumwt.") 

413 

414 cwt = cwt * chanwt ## Merge the weights and flags.  

415 

416 sumchanwt = np.sum(cwt) ## This is a weight 

417 

418 if sumchanwt==0: 

419 raise Exception("Weights are all zero ! ") 

420 else: 

421 

422 for i in range(len(freqlist)): 

423 wt = (freqlist[i] - refnu)/refnu 

424 _ia.open(cubename) 

425 implane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

426 _ia.close() 

427 for tt in range(0,num_terms): 

428 pix[tt] = pix[tt] + (wt**tt) * implane * cwt[i] 

429 

430 for tt in range(0,num_terms): 

431 pix[tt] = pix[tt]/sumchanwt 

432# ia.close() 

433 

434 for tt in range(0,num_terms): 

435 _ia.open(mtname+'.tt'+str(tt)) 

436 _ia.putchunk(pix[tt]) 

437 _ia.close() 

438 

439################################################ 

440 def taylor_model_to_cube(self, cubename='', mtname='',reffreq='1.5GHz',nterms=2): 

441 """ 

442 Convert Taylor coefficients (output of minor cycle) to cube (input to major cycle) 

443 Input : Set of images with suffix : .tt0, .tt1, etc... 

444 Output : Cube 

445 """ 

446 

447 if not os.path.exists(cubename+'.model'): 

448 shutil.copytree(cubename+'.psf', cubename+'.model') 

449 _ia.open(cubename+'.model') 

450 _ia.set(0.0) 

451 _ia.setrestoringbeam(remove=True) 

452 _ia.setbrightnessunit('Jy/pixel') 

453 _ia.close() 

454 

455 

456 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value'] 

457 

458 pix=[] 

459 

460 for tt in range(0,nterms): 

461 _ia.open(mtname+'.model.tt'+str(tt)) 

462 pix.append( _ia.getchunk() ) 

463 _ia.close() 

464 

465 _ia.open(cubename+'.model') 

466 shp = _ia.shape() 

467 _ia.close() 

468 

469 implane = pix[0].copy() 

470 

471 freqlist = self.getFreqList(cubename+'.psf') 

472 for i in range(len(freqlist)): 

473 wt = (freqlist[i] - refnu)/refnu 

474 implane.fill(0.0) 

475 for tt in range(0,nterms): 

476 implane = implane + (wt**tt) * pix[tt] 

477 _ia.open(cubename+'.model') 

478 _ia.putchunk(implane, blc=[0,0,0,i]) 

479 _ia.close() 

480 

481################################################## 

482 

483 

484 def calc_sd_residual(self,origcube='', modelcube = '', residualcube = '', psfcube=''): 

485##, pbcube='', applypb=False, pblimit=0.2): 

486 """ 

487 Residual = Original - ( PSF * Model ) 

488 """ 

489 

490# ia_orig = iatool() 

491# ia_model = iatool() 

492# ia_residual = iatool() 

493# ia_psf = iatool() 

494# 

495# ia_orig.open(origcube) 

496# ia_model.open(modelcube) 

497# ia_residual.open(residualcube) 

498# ia_psf.open(psfcube) 

499 

500 freqlist = self.getFreqList(origcube) 

501 

502 _ia.open(origcube) 

503 shp = _ia.shape() 

504 freqaxis_index = self.getFreqAxisIndex() 

505 _ia.close() 

506 

507 if(freqaxis_index!=3): 

508 casalog.post('The Frequency axis index of '+origcube+' is '+str(freqaxis_index)+' but calc_sd_residual requires index 3.', 'SEVERE') 

509 

510 for i in range(0,len(freqlist)): 

511 

512 _ia.open(origcube) 

513 im_orig = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

514 _ia.close() 

515 

516 _ia.open(modelcube) 

517 im_model = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

518 _ia.close() 

519 

520 _ia.open(psfcube) 

521 im_psf = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i]) 

522 _ia.close() 

523 

524 smoothedim = self.myconvolve(im_model[:,:,0,0], im_psf[:,:,0,0]) 

525 

526 if( np.nansum(im_orig)==0.0): 

527 smoothedim.fill(0.0) 

528 

529 im_residual=im_psf.copy() ## Just to init the shape of this thing 

530 im_residual[:,:,0,0] = im_orig[:,:,0,0] - smoothedim 

531 

532 _ia.open(residualcube) 

533 _ia.putchunk(im_residual, blc=[0,0,0,i]) 

534 _ia.close() 

535 

536# ia_orig.close() 

537# ia_model.close() 

538# ia_residual.close() 

539# ia_psf.close() 

540 

541################################################## 

542 

543 def myconvolve(self,im1,im2): 

544 

545 t3 = time.time() 

546 

547 shp = im1.shape 

548 pshp = (shp[0]*2, shp[1]*2) 

549 pim1 = np.zeros(pshp,'float') 

550 pim2 = np.zeros(pshp,'float') 

551 

552 pim1[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2] = im1 

553 pim2[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2] = im2 

554 

555 fftim1 = fftpack.fftshift( fftpack.fft2( fftpack.ifftshift( pim1 ) ) ) 

556 fftim2 = fftpack.fftshift( fftpack.fft2( fftpack.ifftshift( pim2 ) ) ) 

557 fftprod = fftim1*fftim2 

558 psmoothedim = np.real(fftpack.fftshift( fftpack.ifft2( fftpack.ifftshift( fftprod ) ) ) ) 

559 

560 smoothedim = psmoothedim[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2] 

561 

562 return smoothedim 

563 

564########################################## 

565 

566 def regridimage(self, imagename, template, outfile): 

567 outia = None 

568 _myia = image() 

569 _myia.open(template) 

570 csys = _myia.coordsys() 

571 shape = _myia.shape() 

572 _myia.done() 

573 

574 _myia.open(imagename) 

575 casalog.post("imagename="+imagename) 

576 

577 

578 try: 

579 outia=_myia.regrid(outfile=outfile, 

580 shape=shape, 

581 csys=csys.torecord(), 

582 axes=[0,1], 

583 overwrite=True, 

584 asvelocity=False) 

585 except Exception as instance: 

586 casalog.post("*** Error \'%s\' in regridding image" % (instance), 'WARN') 

587 raise 

588 

589 finally: 

590 csys.done() 

591 if outia != None and outia.isopen(): 

592 outia.done() 

593 _myia.done() 

594 

595 

596########################################## 

597 

598 def createplaneimage(self,imagename, outfile, chanid): 

599 """ 

600 extract a channel plane image  

601 """ 

602 _tmpia=image() 

603 _tmprg=regionmanager() 

604 outia=None 

605 

606 _tmpia.open(imagename) 

607 theregion = _tmprg.frombcs(csys=_tmpia.coordsys().torecord(), shape=_tmpia.shape(), chans=chanid) 

608 try: 

609 outia=_tmpia.subimage(outfile=outfile, region=theregion) 

610 except Exception as instance: 

611 casalog.post("*** Error \'%s\' in creating subimage" % (instance), 'WARN') 

612 

613 _tmpia.close() 

614 _tmpia.done() 

615 _tmprg.done() 

616 if outia != None and outia.isopen(): 

617 outia.done() 

618 

619 def pbcor(self, imagename, pbimage, cutoff, outfile): 

620 """ 

621 pb-correction  

622 """ 

623 outia=None 

624 _myia=image() 

625 _myia.open(imagename) 

626 

627 try: 

628 outia = _myia.pbcor(pbimage=pbimage, outfile=outfile, overwrite=True, 

629 mode='divide', cutoff=cutoff) 

630 except Exception as instance: 

631 casalog.post("*** Error \'%s\' in creating pb-corrected image" % (instance), 'WARN') 

632 

633 finally: 

634 _myia.done() 

635 if outia != None and outia.isopen(): 

636 outia.done() 

637 

638 

639 def checkpsf(self, inpsf, refpsf): 

640 """ 

641 check the center of psf if diffent for  

642 refpsf center and (shift to refpsf position) 

643 in returned psf 

644 """ 

645 tol=0.001 

646 allowshift=True 

647 _ia.open(inpsf) 

648 incsys = _ia.coordsys().torecord() 

649 _ia.close() 

650 #_ia.done() 

651 _tmpia = image() 

652 _tmpia.open(refpsf) 

653 refcsys = _tmpia.coordsys().torecord() 

654 _tmpia.close() 

655 #_tmpia.done()  

656 # check the field center 

657 ramismatch = False 

658 decmismatch = False 

659 

660 indir = incsys['direction0'] 

661 refdir = refcsys['direction0'] 

662 

663 diff_ra = indir['crval'][0] - refdir['crval'][0] 

664 diff_dec = indir['crval'][1] - refdir['crval'][1] 

665 

666 if diff_ra/refdir['crval'][0] > tol: 

667 ramismatch = True 

668 if diff_dec/refdir['crval'][1] > tol: 

669 decmismatch = True 

670 if ramismatch or decmismatch: 

671 casalog.post("The position of SD psf is different from the the psf by (diffRA,diffDec)=( %s, %s)." % (diff_ra, diff_dec) 

672,'WARN') 

673 if allowshift: 

674 modsdpsf=inpsf+'_mod' 

675 casalog.post("Modifying the input SD psf, "+inpsf+" by shifting the field center of sd psf to that of int psf. Modified SD psf image:"+modsdpsf) 

676 shutil.copytree(inpsf, inpsf+'_mod') 

677 _ia.open(modsdpsf) 

678 thecsys = _ia.coordsys() 

679 themodcsysrec = thecsys.torecord() 

680 #repalcing ra, dec of the sd psf to those of the int psf 

681 themodcsysrec['direction0']['crval'][0] = refdir['crval'][0] 

682 themodcsysrec['direction0']['crval'][1] = refdir['crval'][1] 

683 thecsys.fromrecord(themodcsysrec) 

684 _ia.setcoordsys(thecsys) 

685 _ia.close() 

686 #_ia.done() 

687 else: 

688 raise Exception("the center of the psf different from the int psf by (diffRA, diffDec)=(%s,%s)" % (diff_ra, diff_dec)) 

689 

690 else: 

691 casalog.post(" The center of psf coincides with int psf: (diffRA,diffDec)=( %s, %s)" % (diff_ra, diff_dec)) 

692 

693 #### check for frequency axis 

694 _ia.open(inpsf) 

695 sdshape = _ia.shape() 

696 freqaxis_index1 = self.getFreqAxisIndex() 

697 _ia.close() 

698 _ia.open(refpsf) 

699 tshape = _ia.shape() 

700 freqaxis_index2 = self.getFreqAxisIndex() 

701 _ia.close() 

702 if freqaxis_index1 != freqaxis_index2 or sdshape[freqaxis_index1] != tshape[freqaxis_index1]: 

703 raise Exception("The frequency axis of the input SD image and the interferometer template do not match and cannot be regridded. This is because when there are per-plane restoring beams, a regrid along the frequency axis cannot be defined at optimal accuracy. Please re-evaluate the SD image and psf onto a frequency grid that matches the interferometer frequency grid, and then retry.") 

704 

705 #return modpsf  

706 

707 def create_sd_psf(self, sdimage, sdpsfname ): 

708 """ 

709 If sdpsf="", create an SD_PSF cube using restoringbeam information from the sd image.  

710 Start from the regridded SD_IMAGE cube 

711 """ 

712 sdintlib = SDINT_helper() 

713 from casatools import image, componentlist, regionmanager 

714 

715 _ia = image() 

716 _cl = componentlist() 

717 _rg = regionmanager() 

718 

719 ## Get restoringbeam info for all channels 

720 _ia.open(sdimage) 

721 restbeams = _ia.restoringbeam() 

722 shp = _ia.shape() 

723 try: 

724 mysummary = _ia.summary(list=False) 

725 freqaxis_index = list(mysummary['axisnames']).index('Frequency') 

726 except(ValueError): 

727 _ia.close() 

728 casalog.post('The SD image has no frequency axis. Try adding one with ia.adddegaxis() .', 'SEVERE') 

729 try: 

730 stokesaxis_index = list(mysummary['axisnames']).index('Stokes') 

731 nstokes = shp[stokesaxis_index] 

732 except(ValueError): 

733 casalog.post('The SD image has no Stokes axis.', 'WARN') 

734 nstokes = 1 

735 

736 csys = _ia.coordsys() 

737 _ia.close() 

738 

739 # Handle images without per-plane restoring beams 

740 if not 'nChannels' in restbeams or restbeams['nChannels'] != shp[freqaxis_index]: 

741 casalog.post("The input SD image does not have per-plane-restoring beams. Working around that ...", 'WARN') 

742 

743 if shp[freqaxis_index] == 1: # If there is only one channel, just use the one beam we have 

744 mybeams = {'beams': {'*0': {'*0': restbeams.copy()}}, 'nChannels':1, 'nStokes': nstokes} 

745 restbeams = mybeams 

746 elif shp[freqaxis_index] > 1: # create a copy of the SD image with per-plane restoring beams 

747 casalog.post("Constructing per-plane restoring beams based on "+sdimage, 'WARN') 

748 _ia.open(sdimage) 

749 _ia.setrestoringbeam(remove=True) 

750 _ia.setrestoringbeam(beam=restbeams, channel=0, polarization=-1) 

751 restbeams = _ia.restoringbeam() 

752 _ia.close() 

753 else: 

754 casalog.post('The SD image has a frequency axis of length < 1. Cannot proceed.', 'SEVERE') 

755 

756 cdir = csys.torecord()['direction0'] 

757 compdir = [cdir['system'] , str(cdir['crval'][0])+cdir['units'][0] , str(cdir['crval'][1])+cdir['units'][1] ] 

758 

759 ## Make empty SD psf cube from SD image cube 

760 os.system('rm -rf '+sdpsfname) 

761 os.system('cp -R '+sdimage+' '+sdpsfname) 

762 

763 ## Iterate through PSF cube and replace pixels with Gaussians matched to restoringbeam info 

764 

765 _ia.open(sdpsfname) 

766 for ch in range(0,shp[freqaxis_index]): 

767 os.system('rm -rf tmp_sdplane') 

768 rbeam = restbeams['beams']['*'+str(ch)]['*0'] 

769 

770 _cl.close() 

771 _cl.addcomponent(flux=1.0, fluxunit='Jy',polarization='Stokes', dir=compdir, 

772 shape='Gaussian', majoraxis=rbeam['major'], 

773 minoraxis=rbeam['minor'], positionangle=rbeam['positionangle'], 

774 spectrumtype='constant') #, freq=str(freqs[ch])+'Hz') 

775 

776 

777 implane = _ia.getchunk(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch]) 

778 implane.fill(0.0) 

779 _ia.putchunk(implane, blc=[0,0,0,ch]) 

780 _ia.modify(model=_cl.torecord(), subtract=False, region=_rg.box(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch])) 

781 ## Now, normalize it. 

782 implane = _ia.getchunk(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch]) 

783 pmax = np.max(implane) 

784 #print(pmax) 

785 if pmax>0.0: 

786 implane = implane/pmax 

787 else: 

788 implane.fill(0.0) 

789 _ia.putchunk(implane, blc=[0,0,0,ch]) 

790 

791 _ia.close() 

792 

793 

794 def check_coords(self, intres='', intpsf='', intwt = '', sdres='', sdpsf='',sdwt = '', pars=None): 

795 """ 

796 ### (1) Frequency range of cube, data selection range. mtmfs reffreq. 

797 ### (2) nchan too small or too large 

798 ### (3) sumwt : flagged channels in int cubes 

799 ### (4) sd cube empty channels ? Weight image ?  

800 """ 

801 validity=True 

802 

803 freqlist = self.getFreqList(intres) 

804 chanwt = np.ones( len(freqlist), 'float') 

805 if pars['usedata']=='int': 

806 pars['chanwt'] = chanwt 

807 return validity, pars 

808 

809 ## get shapes and gather stats 

810 _ia.open(intres) 

811 int_shp = _ia.shape() 

812 int_stat = _ia.statistics(axes=[0,1]) 

813 _ia.close() 

814 

815 _ia.open(sdres) 

816 sd_shp = _ia.shape() 

817 sd_stat = _ia.statistics(axes=[0,1]) 

818 _ia.close() 

819 

820 _ia.open(intwt) 

821 sumwt = _ia.getchunk() 

822 _ia.close() 

823 

824 

825 ### For mtmfs minor cycle only 

826 if pars['specmode'] in ['mfs','cont'] and pars['deconvolver']=='mtmfs': 

827 ##casalog.post('DOING EXTRA CHECKS##############','WARN') 

828 #(1) # Check reference frequency w.r.to cube freq range.  

829 if pars['reffreq']=='': 

830 reffreq =str( ( freqlist[0] + freqlist[ len(freqlist)-1 ] )/2.0 ) + 'Hz' 

831 casalog.post('The reference frequency for MFS is calculated as the middle of the cube frequency range, irrespective of flagged channels : '+reffreq,'INFO', "task_sdintimaging") 

832 ## Modified parameters :  

833 pars['reffreq'] = reffreq 

834 else: 

835 refval = _qa.convert(_qa.quantity( pars['reffreq'] ), 'Hz') ['value'] 

836 if refval < freqlist[0] or refval > freqlist[ len(freqlist)-1 ] : 

837 if len(freqlist)>1: 

838 casalog.post('The specified reffreq for MFS imaging ('+str(refval)+' Hz) is outside the frequency range of the specified Cube image for the major cycle ('+str(freqlist[0])+' Hz - '+str(freqlist[ len(freqlist)-1 ])+' Hz).\n Please specify a reffreq within the cube frequency range or leave it as an empty string to auto-calculate the middle of the range.','WARN', "task_sdintimaging") 

839 validity=False 

840 else: 

841 casalog.post('The specified reffreq for MFS imaging ('+str(refval)+' Hz) is not exactly the same as the frequency of the selected interferometric data for the major cycle ('+str(freqlist[0])+' Hz).\n We will ignore this for now.', 'WARN') 

842 

843 

844 #(2.1)# Too many channels 

845 if len(freqlist) > 50: 

846 casalog.post('The cube for major cycles has '+str(len(freqlist))+' channels. For wideband continuum imaging, it may be possible to reduce the number of channels to (say) one per spectral window to preserve frequency dependent intensity and weight information but also minimize the number of channels in the image cubes. MFS imaging will be performed within each channel. This will reduce the sizes of the image cubes as well as the compute time used for feathering each plane separately. Note that a minimum of nterms=' + str(pars['nterms']) + ' channels is required for an accurate polynomial fit, but where possible at least 5 to 10 channels that span the frequency range are prefered in order to properly encode frequency dependent intensity and weights.', "WARN", "task_sdintimaging") 

847 

848 #(2.2)# Too few channels  

849 if len(freqlist) < 5 and pars['nterms'] > 1 : 

850 casalog.post('The cube for the major cycle has only '+str(len(freqlist))+' channels. A minimum of nterms = ' + str(pars['nterms']) + ' channels is required for an accurate polynomial fit, but where possible at least 5 to 10 channels that span the frequency range are prefered in order to properly encode frequency dependent intensity and weights.','WARN', "task_sdintimaging") 

851 if len(freqlist) < pars['nterms']: 

852 validity=False 

853 

854 ### For both cube and mtmfs  

855 #(3) ## If there are channels filled with zeros.... create a chanflag to use during 'cube_to_taylor' and 'feathering'  

856 

857 ## INT : If some chans are zero, check that sumwt reflects it too.  

858 zerochans = np.count_nonzero( int_stat['sum']==0.0 ) 

859 if zerochans>0: 

860 casalog.post('There are '+str(zerochans)+' empty channels in the interferometer cube. These channels will be excluded from the feathering step.') 

861 chanwt[ int_stat['sum']==0.0 ] = 0.0 

862 

863 ## SD : If some chans are zero.  

864 ## Set the wt to zero and use the logical AND of int_wt and sd_wt for feathering? 

865 zerochans = np.count_nonzero( sd_stat['sum']==0.0 ) 

866 if zerochans>0: 

867 casalog.post('There are '+str(zerochans)+' empty channels in the single dish image cube. These channels will be excluded from the feathering step. NOTE : We do not yet use SD weights. ') 

868 chanwt[ sd_stat['sum']==0.0 ] = 0.0 

869 

870 ### If there are channels to flag.... list them.  

871 if np.count_nonzero( chanwt==0.0 ): 

872 casalog.post('The following channel weights/flags will be used in the feather step and minor cycle. Zero indicates channels that are empty in either the INT or SD input cubes and which will be excluded from the joint reconstruction. : ' + str(chanwt), 'INFO') 

873 casalog.post('There are channels that are filled with zeros either in the INT cube or the SD cube or both, and they will be ignored from the joint reconstruction. Please search the log file for the string "channel weights/flags" to find a listing of channels that are being used','WARN') 

874 

875 

876 if np.count_nonzero( chanwt != 0.0 ) == 0: 

877 casalog.post("There are no channels with data in both the INT and the SD cubes. Cannot proceed","WARN") 

878 validity=False 

879 

880 pars['chanwt'] = chanwt 

881 return validity, pars 

882 

883 

884 def setup_cube_params(self,sdcube=''): 

885 """ 

886 Read coordinate system from the SD cube 

887 Decide parameters to input into sdintimaging for the INT cube to match.  

888 

889 This is a helper method, not currently used in the sdintimaging task. 

890 We will add this later (after 6.1), and also remove some parameters from the task level. 

891 """ 

892 _ia.open(sdcube) 

893 csys = _ia.coordsys() 

894 shp = _ia.shape() 

895 ctypes = csys.axiscoordinatetypes() 

896 casalog.post("Shape of SD cube : "+str(shp)) 

897 casalog.post("Coordinate ordering : "+str(ctypes)) 

898 if len(ctypes) !=4 or ctypes[3] != 'Spectral': 

899 casalog.post("The SD cube needs to have 4 axes, in the RA/DEC/Stokes/Spectral order", 'ERROR') 

900 _ia.close() 

901 return False 

902 nchan = shp[3] 

903 start = str( csys.referencevalue()['numeric'][3] ) + csys.units()[3] 

904 width = str( csys.increment()['numeric'][3]) + csys.units()[3] 

905 ## Number of channels 

906 casalog.post("nchan = "+str(nchan)) 

907 ## Start Frequency 

908 casalog.post("start = " + start ) 

909 ## Width 

910 casalog.post("width = " + width ) 

911 ## Test for restoringbeams 

912 rbeams = _ia.restoringbeam() 

913 #if not rbeams.has_key('nChannels') or rbeams['nChannels']!= shp[3]: 

914 if not 'nChannels' in rbeams or rbeams['nChannels']!= shp[3]: 

915 casalog.post("The SD Cube needs to have per-plane restoringbeams", 'ERROR') 

916 _ia.close() 

917 return False 

918 else: 

919 casalog.post("Found " + str(rbeams['nChannels']) + " per-plane restoring beams") 

920 casalog.post("\n(For specmode='mfs' in sdintimaging, please remember to set 'reffreq' to a value within the freq range of the cube)\n") 

921 _ia.close() 

922 return {'nchan':nchan, 'start':start, 'width':width} 

923 

924 

925 

926 

927### Using Old Imager. Does not work for cubes ?  

928# def fit_psf_beam(self,msname = '', psfname =''): 

929# _im.open(msname)  

930# _ia.open(psfname) 

931# csys = _ia.coordsys() 

932# rbeam_old = _ia.restoringbeam() 

933# print(rbeam_old) 

934# shp = _ia.shape() 

935# _ia.close() 

936# cellx = csys.increment()['numeric'][0]; 

937# celly = csys.increment()['numeric'][1]; 

938# _im.defineimage(nx=shp[0],ny=shp[1],cellx=str(cellx)+'rad',celly=str(celly)+'rad',nchan=3)  

939# params =_im.fitpsf(psfname) 

940# print(params) 

941# _im.close()  

942#