Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/sdint_helper.py: 6%
571 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-21 07:43 +0000
1from scipy import fftpack
2import numpy as np
3import shutil
4import os
5import time
6import random
8from casatools import quanta, table, image, componentlist, regionmanager, imager
9from casatasks import casalog, feather
11_ia = image()
12_qa = quanta()
14class SDINT_helper:
16 def __init__(self, tmpFileNumber=0):
17 if type(tmpFileNumber) != int or tmpFileNumber<=0:
18 casalog.post('No valid tmp file tag provided. Will use a random number.', 'DEBUG')
19 tmpFileNumber = random.randint(100000,999999)
21 self.tmpFileTag = str(tmpFileNumber)
22 casalog.post('Initializing SDINT Helper instance with tmp file tag '+self.tmpFileTag)
24 # define tmp file names
25 self.tmpintplane = 'tmp_'+self.tmpFileTag+'_intplane'
26 self.tmpsdplane = 'tmp_'+self.tmpFileTag+'_sdplane'
27 self.tmpjointplane = 'tmp_'+self.tmpFileTag+'_jointplane'
28 self.tmpallplanes = 'tmp_'+self.tmpFileTag+'_*plane'
30 def getTmpFileTag(self):
31 return self.tmpFileTag
33 def getAllTmpFileRegex(self):
34 return self.tmpallplanes
36 def deleteTmpFiles(self):
37 casalog.post('Deleting '+self.tmpallplanes)
38 os.system('rm -rf '+self.tmpallplanes)
41################################################
42 def getFreqAxisIndex(self, myia):
43 try:
44 mysummary = myia.summary(list=False)
45 freqaxis_index = list(mysummary['axisnames']).index('Frequency')
46 except(ValueError):
47 myia.close()
48 casalog.post('The image '+myia.name()+' has no frequency axis. Try adding one with ia.adddegaxis() .', 'SEVERE')
50 return freqaxis_index
52 def getFreqList(self,imname=''):
53 """ Get the list of frequencies for the given image, one for each channel.
55 Returns:
56 list[float] The frequencies for each channel in the image, in Hz.
57 """
59 _ia.open(imname)
60 csys =_ia.coordsys()
61 shp = _ia.shape()
62 freqaxis_index = self.getFreqAxisIndex(_ia)
63 _ia.close()
65 if(csys.axiscoordinatetypes()[freqaxis_index] == 'Spectral'):
66 restfreq = csys.referencevalue()['numeric'][freqaxis_index]
67 freqincrement = csys.increment()['numeric'][freqaxis_index]
68 freqlist = [];
69 for chan in range(0,shp[freqaxis_index]):
70 freqlist.append(restfreq + chan * freqincrement);
71 elif(csys.axiscoordinatetypes()[freqaxis_index] == 'Tabular'):
72 freqlist = (csys.torecord()['tabular2']['worldvalues'])
73 else:
74 casalog.post('Unknown frequency axis. Exiting.','SEVERE');
75 return False;
77 csys.done()
78 return freqlist
80################################################
82 def copy_restoringbeam(self,fromthis='',tothis=''):
84 ia2 = image()
86 freqlist = self.getFreqList(fromthis)
87 casalog.post('Copying restoring beam(s) for '+str(len(freqlist))+' channel(s) from '+fromthis+' to '+tothis)
88 ia2.open(tothis)
89 try:
90 ia2.setrestoringbeam(imagename=fromthis)
91 except Exception as e:
92 casalog.post('Error while copying restoringbeam: '+str(e), 'WARN')
93 finally:
94 ia2.close()
97################################################
99 def feather_residual(self, int_cube, sd_cube, joint_cube, applypb, inparm):
101 if applypb==True:
102 ## Take initial INT_dirty image to flat-sky.
103 self.modify_with_pb(inpcube=int_cube+'.residual',
104 pbcube=int_cube+'.pb',
105 cubewt=int_cube+'.sumwt',
106 chanwt=inparm['chanwt'],
107 action='div',
108 pblimit=inparm['pblimit'],
109 freqdep=True)
111 ## Feather flat-sky INT dirty image with SD image
112 self.feather_int_sd(sdcube=sd_cube+'.residual',
113 intcube=int_cube+'.residual',
114 jointcube=joint_cube+'.residual', ## output
115 sdgain=inparm['sdgain'],
116 dishdia=inparm['dishdia'],
117 usedata=inparm['usedata'],
118 chanwt = inparm['chanwt'])
120 if applypb==True:
121 if inparm['specmode'].count('cube')>0:
122 ## Multiply the new JOINT dirty image by the frequency-dependent PB.
123 fdep_pb = True
124 else:
125 ## Multiply new JOINT dirty image by a common PB to get the effect of conjbeams.
126 fdep_pb = False
127 self.modify_with_pb(inpcube=joint_cube+'.residual',
128 pbcube=int_cube+'.pb',
129 cubewt=int_cube+'.sumwt',
130 chanwt=inparm['chanwt'],
131 action='mult',
132 pblimit=inparm['pblimit'],
133 freqdep=fdep_pb)
136################################################
138 def feather_int_sd(self, sdcube='', intcube='', jointcube='',sdgain=1.0, dishdia=-1, usedata='sdint', chanwt=''):
139 """
140 Run the feather task to combine the SD and INT Cubes.
142 There's a bug in feather for cubes. Hence, do each channel separately.
143 FIX feather and then change this. CAS-5883 is the JIRA ticket that contains a fix for this issue....
145 TODO : Add the effdishdia usage to get freq-indep feathering.
147 """
149 ### Do the feathering.
150 if usedata=='sdint':
151 ## Feather runs in a loop on chans internally, but there are issues with open tablecache images
152 ## Also, no way to set effective dish dia separately for each channel.
153 #feather(imagename = jointcube, highres = intcube, lowres = sdcube, sdfactor = sdgain, effdishdiam=-1)
155 freqlist = self.getFreqList(sdcube)
157 os.system('rm -rf '+jointcube)
158 os.system('cp -r ' + intcube + ' ' + jointcube)
160 _ib = image()
162 _ia.open(jointcube)
163 _ia.set(0.0) ## Initialize this to zero for all planes
165 for i in range(len(freqlist)):
166 if chanwt[i] != 0.0 : ## process only the nonzero channels
167 if(dishdia <=0):
168 casalog.post('Parameter dishdia (SD dish diameter in meters) must be > 0.', 'SEVERE')
170 freqdishdia = dishdia ## * (freqlist[0] / freqlist[i]) # * 0.5
172 self.deleteTmpFiles()
173 self.createplaneimage(imagename=sdcube, outfile=self.tmpsdplane, chanid = str(i));
174 self.createplaneimage(imagename=intcube, outfile=self.tmpintplane, chanid = str(i));
176 #feather(imagename = self.tmpjointplane, highres = self.tmpintplane, lowres = self.tmpsdplane, sdfactor = sdgain, effdishdiam=freqdishdia)
177 # feathering via toolkit
178 try:
179 casalog.post("start Feathering.....")
180 imFea=imager( )
181 imFea.setvp(dovp=True)
182 imFea.setsdoptions(scale=sdgain)
183 imFea.feather(image=self.tmpjointplane, highres=self.tmpintplane, lowres=self.tmpsdplane, effdishdiam=freqdishdia)
184 imFea.done( )
185 del imFea
186 except Exception as instance:
187 casalog.post('*** Error *** %s' % instance, 'ERROR')
188 raise
190 _ib.open(self.tmpjointplane)
191 pixjoint = _ib.getchunk()
192 _ib.close()
193 _ia.putchunk(pixjoint, blc=[0,0,0,i])
195 _ia.close()
197 if usedata=='sd': ## Copy sdcube to joint.
198 os.system('rm -rf '+jointcube)
199 os.system('cp -r ' + sdcube + ' ' + jointcube)
200 elif usedata=='int': ## Copy intcube to joint
201 os.system('rm -rf '+jointcube)
202 os.system('cp -r ' + intcube + ' ' + jointcube)
204################################################
205 def calc_renorm(self, intname='', jointname=''):
206 """
207 Calculate a new .sumwt spectrum containing the peak of the feathered PSF.
208 The PSF and each residual image calculation will be re-normalized by this.
209 This will keep the PSFs in all channels at a peak of 1.
210 """
211 psfname = jointname+'.psf'
212 os.system('cp -r '+intname+'.sumwt ' + jointname + '.sumwt')
213 _ia.open(jointname+'.sumwt')
214 vals = _ia.getchunk()
215 shp = _ia.shape()
216 freqaxis_index = self.getFreqAxisIndex(_ia)
217 _ia.close()
219 if shp[0]>1:
220 casalog.post("WARNING : Cannot use this task with faceting", 'WARN')
222 _ia.open(jointname+'.psf')
223 for i in range(0, shp[freqaxis_index]):
224 onepsf = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
225 vals[0,0,0,i] = np.max(onepsf)
226 _ia.close()
228 _ia.open(jointname+'.sumwt')
229 _ia.putchunk(vals)
230 _ia.close()
232 casalog.post("********************Re-norm with "+str(vals))
235################################################
237 def apply_renorm(self, imname='', sumwtname=''):
238 """
239 Divide each plane of the input image by the sumwt value for that channel
240 """
241 _ia.open(sumwtname)
242 shp = _ia.shape()
243 freqaxis_index = self.getFreqAxisIndex(_ia)
244 vals = _ia.getchunk() ## This is one pixel per channel.
245 _ia.close()
247 casalog.post("******************** Re-norm with "+str(vals))
249 _ia.open(imname)
250 for i in range(0, shp[freqaxis_index]):
251 oneplane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
252 if vals[0,0,0,i]>0.0:
253 normplane = oneplane/vals[0,0,0,i]
254 else:
255 normplane = oneplane.copy()
256 normplane.fill(0.0)
257 _ia.putchunk( normplane , blc=[0,0,0,i] )
258 _ia.close()
262################################################
264 def modify_with_pb(self, inpcube='', pbcube='',cubewt='', chanwt='', action='mult',pblimit=0.2, freqdep=True):
265 """
266 Multiply or divide by the PB
268 Args:
269 inpcube: The cube to be modified. For example: "try.int.cube.model"
270 pbcube: The primary beam to multiply/divide by. For example: "try.int.cube.pb"
271 cubewt: The per-channel weight of the inpcube. For example: "try.int.cube.sumwt"
272 chanwt: List of 0s and 1s, one per channel, to effectively disable the effect of a channel on the resulting images.
273 action: 'mult' or 'div', to multiply by the PB or divide by it.
274 pblimit: For pixels less than this value in the PB, set those same pixels in the inpcube to zero.
275 freqdep: True for channel by channel, False to use a freq-independent PB from the middle of the list before/after deconvolution
276 """
277 casalog.post('Modify with PB : ' + action + ' with frequency dependence ' + str(freqdep))
279 freqlist = self.getFreqList(inpcube)
281 _ia.open(inpcube)
282 shp=_ia.shape()
283 freqaxis_index = self.getFreqAxisIndex(_ia)
284 _ia.close()
286 if(freqaxis_index!=3):
287 casalog.post('The Frequency axis index of '+inpcube+' is '+str(freqaxis_index)+' but modify_with_pb requires index 3.', 'SEVERE')
289 ##############
290 ### Calculate a reference Primary Beam
291 ### Weighted sum of pb cube
293 ##midchan = int(len(freqlist)/2)
294 ##refchan = len(freqlist)-1 ## This assumes ascending frequency ordering in chans.
295 refchan=0
296 _ia.open(pbcube)
297 pbplane = _ia.getchunk(blc=[0,0,0,refchan],trc=[shp[0],shp[1],0,refchan])
298 _ia.close()
299 pbplane.fill(0.0)
301# pbplane = np.zeros( (shp[0],shp[1]), 'float')
303 if freqdep==False:
304 _ia.open(cubewt) # .sumwt
305 cwt = _ia.getchunk()[0,0,0,:]
306 _ia.close()
308 if shp[3] != len(cwt) or len(freqlist) != len(cwt):
309 raise Exception("Modify with PB : Nchan shape mismatch between cube and sumwt.")
311 cwt = cwt * chanwt ## Merge the weights and flags
313 sumchanwt = np.sum(cwt)
315 if sumchanwt==0:
316 raise Exception("Weights are all zero ! ")
318 _ia.open(pbcube)
319 for i in range(len(freqlist)):
320 ## Read the pb per plane
321 pbplane = pbplane + cwt[i] * _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
322 _ia.close()
324 pbplane = pbplane / sumchanwt
326 ##############
329 ## Special-case for setting the PBmask to be same for all freqs
330 if freqdep==False:
331 shutil.copytree(pbcube, pbcube+'_tmpcopy')
333 for i in range(len(freqlist)):
335 ## Read the pb per plane
336 if freqdep==True:
337 _ia.open(pbcube)
338 pbplane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
339 _ia.close()
341 ## Make a tmp pbcube with the same pb in all planes. This is for the mask.
342 if freqdep==False:
343 _ia.open(pbcube+'_tmpcopy')
344 _ia.putchunk(pbplane, blc=[0,0,0,i])
345 _ia.close()
347 _ia.open(inpcube)
348 implane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
350 outplane = pbplane.copy()
351 outplane.fill(0.0)
353 if action=='mult':
354 pbplane[pbplane<pblimit]=0.0
355 outplane = implane * pbplane
356 else:
357 implane[pbplane<pblimit]=0.0
358 pbplane[pbplane<pblimit]=1.0
359 outplane = implane / pbplane
361 _ia.putchunk(outplane, blc=[0,0,0,i])
362 _ia.close()
364# if freqdep==True:
365# ## Set a mask based on frequency-dependent PB
366# self.addmask(inpcube,pbcube,pblimit)
367# else:
368 if freqdep==False:
369 ## Set a mask based on the PB in refchan
370 self.addmask(inpcube,pbcube+'_tmpcopy',pblimit)
371 shutil.rmtree(pbcube+'_tmpcopy')
374################################################
375 def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='replace'):
376 #def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='add'):
377 #def addmask(self, inpimage='',pbimage='',pblimit=0.2, mode='old'):
378 """
379 add pb mask: create a new mask called 'pbmask' and set it as a defualt mask
380 mode: "replace" or "add"
381 relpalce: create a pbmask based on pblimit without account for the exist mask
382 add: create a pbmask based on pblimit and merge with existing default mask
383 """
384 _ia.open(inpimage)
385 defaultmaskname=_ia.maskhandler('default')[0]
386 allmasknames = _ia.maskhandler('get')
387 # casalog.post("defaultmaskname=",defaultmaskname)
388 if mode=='replace':
389 if defaultmaskname!='' and defaultmaskname!='mask0':
390 _ia.calcmask(mask='"'+pbimage+'"'+'>'+str(pblimit), name=defaultmaskname);
392 elif defaultmaskname=='mask0':
393 if 'pbmask' in allmasknames:
394 _ia.maskhandler('delete','pbmask')
395 _ia.calcmask(mask='"'+pbimage+'"'+'>'+str(pblimit), name='pbmask');
397 #elif mode=='add':
398 # After deleting a pixel mask it sometimes leaves it in cache
399 #
400 # _ia.open(inpimage)
401 # if defaultmaskname=='pbmask':
402 # _ia.maskhandler('delete',defaultmaskname)
403 # _ia.close()
404 # _ia.open(inpimage)
406 # _ia.calcmask(mask='mask("'+inpimage+'")||'+'"'+pbimage+'"'+'>'+str(pblimit), name='pbmask');
407 elif mode=='old':
408 # this one create a new mask every time this function is called!
409 #_ia.open(inpimage)
410 _ia.calcmask(mask='mask("'+inpimage+'")||'+'"'+pbimage+'"'+'>'+str(pblimit));
411 else:
412 raise Exception("Unrecongnized value for mode: ",mode)
413 _ia.close()
414 _ia.done()
416################################################
417 def cube_to_taylor_sum(self, cubename='', cubewt='', chanwt='', mtname='',reffreq='1.5GHz',nterms=2,dopsf=False):
418 """
419 Convert Cubes (output of major cycle) to Taylor weighted averages (inputs to the minor cycle)
420 Input : Cube image <cubename>, with channels weighted by image <cubewt>
421 Output : Set of images : <mtname>.tt0, <mtname>.tt1, etc...
422 Algorithm: I_ttN = sum([ I_v * ((f-ref)/ref)**N for f in freqs ])
424 Args:
425 cubename: Name of a cube image to interpret into a set of taylor term .ttN images, eg "try.residual", "joint.cube.psf".
426 cubewt: Name of a .sumwt image that contains the per-channel weighting for the interferometer image.
427 chanwt: List of 0s and 1s, one per channel, to effectively disable the effect of a channel on the resulting images.
428 mtname: The prefix output name, to be concatenated with ".ttN" strings, eg "try_mt.residual", "joint.multiterm.psf"
429 These images should already exist by the time this function is called.
430 It's suggested that this have same suffix as cubename.
431 dopsf: Signals that cubename represents a point source function, should be true if cubename ends with ".psf".
432 If true, then output 2*nterms-1 ttN images.
433 """
435 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value']
437 # casalog.post("&&&&&&&&& REF FREQ : " + str(refnu))
439 pix=[]
441 num_terms=nterms
443 if dopsf==True:
444 num_terms=2*nterms-1
446 for tt in range(0,num_terms):
447 _ia.open(mtname+'.tt'+str(tt))
448 pix.append( _ia.getchunk() )
449 _ia.close()
450 pix[tt].fill(0.0)
452 _ia.open(cubename)
453 shp = _ia.shape()
454 freqaxis_index = self.getFreqAxisIndex(_ia)
455 _ia.close()
457 if(freqaxis_index!=3):
458 casalog.post('The Frequency axis index of '+cubename+' is '+str(freqaxis_index)+' but cube_to_taylor_sum requires index 3.', 'SEVERE')
460 _ia.open(cubewt)
461 cwt = _ia.getchunk()[0,0,0,:]
462 _ia.close()
465 freqlist = self.getFreqList(cubename)
467 if shp[3] != len(cwt) or len(freqlist) != len(cwt):
468 raise Exception("Nchan shape mismatch between cube and sumwt.")
470 cwt = cwt * chanwt ## Merge the weights and flags.
472 sumchanwt = np.sum(cwt) ## This is a weight
474 if sumchanwt==0:
475 raise Exception("Weights are all zero ! ")
476 else:
478 _ia.open(cubename)
479 for i in range(len(freqlist)):
480 wt = (freqlist[i] - refnu)/refnu
481 implane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
482 for tt in range(0,num_terms):
483 pix[tt] = pix[tt] + (wt**tt) * implane * cwt[i]
484 _ia.close()
486 for tt in range(0,num_terms):
487 pix[tt] = pix[tt]/sumchanwt
488# ia.close()
490 for tt in range(0,num_terms):
491 _ia.open(mtname+'.tt'+str(tt))
492 _ia.putchunk(pix[tt])
493 _ia.close()
495################################################
496 def taylor_model_to_cube(self, cubename='', mtname='',reffreq='1.5GHz',nterms=2):
497 """
498 Convert Taylor coefficients (output of minor cycle) to cube (input to major cycle)
499 Input : Set of images with suffix : .tt0, .tt1, etc...
500 Output : Cube
501 """
503 if not os.path.exists(cubename+'.model'):
504 shutil.copytree(cubename+'.psf', cubename+'.model')
505 _ia.open(cubename+'.model')
506 _ia.set(0.0)
507 _ia.setrestoringbeam(remove=True)
508 _ia.setbrightnessunit('Jy/pixel')
509 _ia.close()
512 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value']
514 pix=[]
516 for tt in range(0,nterms):
517 _ia.open(mtname+'.model.tt'+str(tt))
518 pix.append( _ia.getchunk() )
519 _ia.close()
521 _ia.open(cubename+'.model')
522 shp = _ia.shape()
523 _ia.close()
525 implane = pix[0].copy()
527 freqlist = self.getFreqList(cubename+'.psf')
528 _ia.open(cubename+'.model')
529 for i in range(len(freqlist)):
530 wt = (freqlist[i] - refnu)/refnu
531 implane.fill(0.0)
532 for tt in range(0,nterms):
533 implane = implane + (wt**tt) * pix[tt]
534 _ia.putchunk(implane, blc=[0,0,0,i])
535 _ia.close()
537##################################################
540 def calc_sd_residual(self,origcube='', modelcube = '', residualcube = '', psfcube=''):
542 """
543 Residual = Original - ( PSF * Model )
544 """
546 ia_orig = image()
547 ia_model = image()
548 ia_psf = image()
549 ia_res = image()
551 ia_orig.open(origcube)
553 shp = ia_orig.shape()
554 freqaxis_index = self.getFreqAxisIndex(ia_orig)
556 if(freqaxis_index!=3):
557 ia_orig.close()
558 casalog.post('The Frequency axis index of '+origcube+' is '+str(freqaxis_index)+' but calc_sd_residual requires index 3.', 'SEVERE')
560 ia_model.open(modelcube)
561 ia_psf.open(psfcube)
562 ia_res.open(residualcube)
564 for i in range(0,shp[freqaxis_index]): # loop over all channels
566 im_orig = ia_orig.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
568 im_model = ia_model.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
570 im_psf = ia_psf.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
572 smoothedim = self.myconvolve(im_model[:,:,0,0], im_psf[:,:,0,0])
574 if( np.nansum(im_orig)==0.0):
575 smoothedim.fill(0.0)
577 im_residual=im_psf.copy() ## Just to init the shape of this thing
578 im_residual[:,:,0,0] = im_orig[:,:,0,0] - smoothedim
580 ia_res.putchunk(im_residual, blc=[0,0,0,i])
582 ia_orig.close()
583 ia_model.close()
584 ia_psf.close()
585 ia_res.close()
587##################################################
589 def myconvolve(self,im1,im2):
591 t3 = time.time()
593 shp = im1.shape
594 pshp = (shp[0]*2, shp[1]*2)
595 pim1 = np.zeros(pshp,'float')
596 pim2 = np.zeros(pshp,'float')
598 pim1[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2] = im1
599 pim2[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2] = im2
601 fftim1 = fftpack.fftshift( fftpack.fft2( fftpack.ifftshift( pim1 ) ) )
602 fftim2 = fftpack.fftshift( fftpack.fft2( fftpack.ifftshift( pim2 ) ) )
603 fftprod = fftim1*fftim2
604 psmoothedim = np.real(fftpack.fftshift( fftpack.ifft2( fftpack.ifftshift( fftprod ) ) ) )
606 smoothedim = psmoothedim[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2]
608 return smoothedim
610##########################################
612 def regridimage(self, imagename, template, outfile):
613 outia = None
614 _myia = image()
615 _myia.open(template)
616 csys = _myia.coordsys()
617 shape = _myia.shape()
618 _myia.done()
620 _myia.open(imagename)
621 myshape = _myia.shape()
622 if len(myshape)!=len(shape):
623 raise Exception("Image "+imagename+" and template image "+template+" have different number of axes: "
624 +str(len(myshape))+" and "+str(len(shape))+" .")
626 mysummary = _myia.summary(list=False)
627 extractStokes = False
629 for i in [2,3]:
630 if len(myshape)>i:
631 if myshape[i]!=shape[i]:
632 casalog.post("Image "+imagename+" and template image "+template+" have different shape in axis "
633 +str(i)+" ("+mysummary['axisnames'][i]+"): "+str(myshape[i])+" and "+str(shape[i])+" .", 'WARN')
634 if mysummary['axisnames'][i] == 'Stokes' and myshape[i]>1 and shape[i]==1:
635 extractStokes = True
636 try:
637 if extractStokes:
638 casalog.post("Will try to continue by extracting Stokes I from "+imagename+" ...", 'WARN')
639 tmpimage = 'tmp_'+self.tmpFileTag+'_'+imagename+'_StokesIOnly'
640 tmp_csys = _myia.coordsys()
641 myrg = regionmanager()
642 xregion = myrg.frombcs(csys=tmp_csys.torecord(), shape=myshape, stokes='I', stokescontrol="a")
643 tmp_csys.done()
644 _myia.subimage(outfile=tmpimage, region=xregion, dropdeg=False, overwrite=True, stretch=False, keepaxes=[])
645 _myia.close()
646 casalog.post("Success. Created "+tmpimage, 'INFO')
647 _myia.open(tmpimage)
649 outia=_myia.regrid(outfile=outfile,
650 shape=shape,
651 csys=csys.torecord(),
652 axes=[0,1],
653 overwrite=True,
654 asvelocity=False)
656 except Exception as instance:
657 casalog.post("Error while regridding image \'"+imagename+"\':", 'WARN')
658 casalog.post("*** \'%s\'" % (instance), 'WARN')
659 raise
661 finally:
662 csys.done()
663 if outia != None and outia.isopen():
664 outia.done()
665 _myia.done()
668##########################################
670 def createplaneimage(self,imagename, outfile, chanid):
671 """
672 extract a channel plane image
673 """
674 _tmpia=image()
675 _tmprg=regionmanager()
676 outia=None
678 _tmpia.open(imagename)
679 theregion = _tmprg.frombcs(csys=_tmpia.coordsys().torecord(), shape=_tmpia.shape(), chans=chanid)
680 try:
681 outia=_tmpia.subimage(outfile=outfile, region=theregion)
682 except Exception as instance:
683 casalog.post("*** Error \'%s\' in creating subimage" % (instance), 'WARN')
685 _tmpia.close()
686 _tmpia.done()
687 _tmprg.done()
688 if outia != None and outia.isopen():
689 outia.done()
691 def pbcor(self, imagename, pbimage, cutoff, outfile):
692 """
693 pb-correction
694 """
695 outia=None
696 _myia=image()
697 _myia.open(imagename)
699 try:
700 outia = _myia.pbcor(pbimage=pbimage, outfile=outfile, overwrite=True,
701 mode='divide', cutoff=cutoff)
702 except Exception as instance:
703 casalog.post("*** Error \'%s\' in creating pb-corrected image" % (instance), 'WARN')
705 finally:
706 _myia.done()
707 if outia != None and outia.isopen():
708 outia.done()
711 def checkpsf(self, inpsf, refpsf):
712 """
713 check the center of psf if diffent for
714 refpsf center and (shift to refpsf position)
715 in returned psf
716 """
717 tol=0.001
718 allowshift=True
719 _ia.open(inpsf)
720 incsys = _ia.coordsys().torecord()
721 _ia.close()
722 #_ia.done()
723 _tmpia = image()
724 _tmpia.open(refpsf)
725 refcsys = _tmpia.coordsys().torecord()
726 _tmpia.close()
727 #_tmpia.done()
728 # check the field center
729 ramismatch = False
730 decmismatch = False
732 indir = incsys['direction0']
733 refdir = refcsys['direction0']
735 diff_ra = indir['crval'][0] - refdir['crval'][0]
736 diff_dec = indir['crval'][1] - refdir['crval'][1]
738 if diff_ra/refdir['crval'][0] > tol:
739 ramismatch = True
740 if diff_dec/refdir['crval'][1] > tol:
741 decmismatch = True
742 if ramismatch or decmismatch:
743 casalog.post("The position of SD psf is different from the the psf by (diffRA,diffDec)=( %s, %s)." % (diff_ra, diff_dec)
744,'WARN')
745 if allowshift:
746 modsdpsf=inpsf+'_mod'
747 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)
748 shutil.copytree(inpsf, inpsf+'_mod')
749 _ia.open(modsdpsf)
750 thecsys = _ia.coordsys()
751 themodcsysrec = thecsys.torecord()
752 #repalcing ra, dec of the sd psf to those of the int psf
753 themodcsysrec['direction0']['crval'][0] = refdir['crval'][0]
754 themodcsysrec['direction0']['crval'][1] = refdir['crval'][1]
755 thecsys.fromrecord(themodcsysrec)
756 _ia.setcoordsys(thecsys)
757 _ia.close()
758 #_ia.done()
759 else:
760 raise Exception("the center of the psf different from the int psf by (diffRA, diffDec)=(%s,%s)" % (diff_ra, diff_dec))
762 else:
763 casalog.post(" The center of psf coincides with int psf: (diffRA,diffDec)=( %s, %s)" % (diff_ra, diff_dec))
765 #### check for frequency axis
766 _ia.open(inpsf)
767 sdshape = _ia.shape()
768 freqaxis_index1 = self.getFreqAxisIndex(_ia)
769 _ia.close()
770 _ia.open(refpsf)
771 tshape = _ia.shape()
772 freqaxis_index2 = self.getFreqAxisIndex(_ia)
773 _ia.close()
774 if freqaxis_index1 != freqaxis_index2 or sdshape[freqaxis_index1] != tshape[freqaxis_index1]:
775 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.")
777 #return modpsf
779 def create_sd_psf(self, sdimage, sdpsfname ):
780 """
781 If sdpsf="", create an SD_PSF cube using restoringbeam information from the sd image.
782 Start from the regridded SD_IMAGE cube
783 """
784 sdintlib = SDINT_helper()
786 ia2 = image()
787 _cl = componentlist()
788 _rg = regionmanager()
790 ## Get restoringbeam info for all channels
791 ia2.open(sdimage)
792 restbeams = ia2.restoringbeam()
793 shp = ia2.shape()
794 try:
795 mysummary = ia2.summary(list=False)
796 freqaxis_index = list(mysummary['axisnames']).index('Frequency')
797 except(ValueError):
798 ia2.close()
799 casalog.post('The SD image has no frequency axis. Try adding one with ia.adddegaxis() .', 'SEVERE')
800 try:
801 stokesaxis_index = list(mysummary['axisnames']).index('Stokes')
802 nstokes = shp[stokesaxis_index]
803 except(ValueError):
804 casalog.post('The SD image has no Stokes axis.', 'WARN')
805 nstokes = 1
807 csys = ia2.coordsys()
808 ia2.close()
810 # Handle images without per-plane restoring beams
811 if not 'nChannels' in restbeams or restbeams['nChannels'] != shp[freqaxis_index]:
812 casalog.post("The input SD image does not have per-plane-restoring beams. Working around that ...", 'WARN')
814 if shp[freqaxis_index] == 1: # If there is only one channel, just use the one beam we have
815 mybeams = {'beams': {'*0': {'*0': restbeams.copy()}}, 'nChannels':1, 'nStokes': nstokes}
816 restbeams = mybeams
817 elif shp[freqaxis_index] > 1: # create a copy of the SD image with per-plane restoring beams
818 casalog.post("Constructing per-plane restoring beams based on "+sdimage, 'WARN')
819 ia2.open(sdimage)
820 ia2.setrestoringbeam(remove=True)
821 ia2.setrestoringbeam(beam=restbeams, channel=0, polarization=-1)
822 restbeams = ia2.restoringbeam()
823 ia2.close()
824 else:
825 casalog.post('The SD image has a frequency axis of length < 1. Cannot proceed.', 'SEVERE')
827 cdir = csys.torecord()['direction0']
828 compdir = [cdir['system'] , str(cdir['crval'][0])+cdir['units'][0] , str(cdir['crval'][1])+cdir['units'][1] ]
830 ## Make empty SD psf cube from SD image cube
831 os.system('rm -rf '+sdpsfname)
832 os.system('cp -R '+sdimage+' '+sdpsfname)
834 ## Iterate through PSF cube and replace pixels with Gaussians matched to restoringbeam info
836 ia2.open(sdpsfname)
837 for ch in range(0,shp[freqaxis_index]):
838 os.system('rm -rf '+self.tmpsdplane)
839 rbeam = restbeams['beams']['*'+str(ch)]['*0']
841 _cl.close()
842 _cl.addcomponent(flux=1.0, fluxunit='Jy',polarization='Stokes', dir=compdir,
843 shape='Gaussian', majoraxis=rbeam['major'],
844 minoraxis=rbeam['minor'], positionangle=rbeam['positionangle'],
845 spectrumtype='constant') #, freq=str(freqs[ch])+'Hz')
848 implane = ia2.getchunk(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch])
849 implane.fill(0.0)
850 ia2.putchunk(implane, blc=[0,0,0,ch])
851 ia2.modify(model=_cl.torecord(), subtract=False, region=_rg.box(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch]))
852 ## Now, normalize it.
853 implane = ia2.getchunk(blc=[0,0,0,ch],trc=[shp[0],shp[1],0,ch])
854 pmax = np.max(implane)
855 #print(pmax)
856 if pmax>0.0:
857 implane = implane/pmax
858 else:
859 implane.fill(0.0)
860 ia2.putchunk(implane, blc=[0,0,0,ch])
862 ia2.close()
865 def check_coords(self, intres='', intpsf='', intwt = '', sdres='', sdpsf='',sdwt = '', pars=None):
866 """
867 ### (1) Frequency range of cube, data selection range. mtmfs reffreq.
868 ### (2) nchan too small or too large
869 ### (3) sumwt : flagged channels in int cubes
870 ### (4) sd cube empty channels ? Weight image ?
871 """
872 validity=True
874 freqlist = self.getFreqList(intres)
875 chanwt = np.ones( len(freqlist), 'float')
876 if pars['usedata']=='int':
877 pars['chanwt'] = chanwt
878 return validity, pars
880 ## get shapes and gather stats
881 _ia.open(intres)
882 int_shp = _ia.shape()
883 int_stat = _ia.statistics(axes=[0,1])
884 _ia.close()
886 _ia.open(sdres)
887 sd_shp = _ia.shape()
888 sd_stat = _ia.statistics(axes=[0,1])
889 _ia.close()
891 _ia.open(intwt)
892 sumwt = _ia.getchunk()
893 _ia.close()
896 ### For mtmfs minor cycle only
897 if pars['specmode'] in ['mfs','cont'] and pars['deconvolver']=='mtmfs':
898 ##casalog.post('DOING EXTRA CHECKS##############','WARN')
899 #(1) # Check reference frequency w.r.to cube freq range.
900 if pars['reffreq']=='':
901 reffreq =str( ( freqlist[0] + freqlist[ len(freqlist)-1 ] )/2.0 ) + 'Hz'
902 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")
903 ## Modified parameters :
904 pars['reffreq'] = reffreq
905 else:
906 refval = _qa.convert(_qa.quantity( pars['reffreq'] ), 'Hz') ['value']
907 if refval < freqlist[0] or refval > freqlist[ len(freqlist)-1 ] :
908 if len(freqlist)>1:
909 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")
910 validity=False
911 else:
912 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')
915 #(2.1)# Too many channels
916 if len(freqlist) > 50:
917 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")
919 #(2.2)# Too few channels
920 if len(freqlist) < 5 and pars['nterms'] > 1 :
921 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")
922 if len(freqlist) < pars['nterms']:
923 validity=False
925 ### For both cube and mtmfs
926 #(3) ## If there are channels filled with zeros.... create a chanflag to use during 'cube_to_taylor' and 'feathering'
928 ## INT : If some chans are zero, check that sumwt reflects it too.
929 zerochans = np.count_nonzero( int_stat['sum']==0.0 )
930 if zerochans>0:
931 casalog.post('There are '+str(zerochans)+' empty channels in the interferometer cube. These channels will be excluded from the feathering step.')
932 chanwt[ int_stat['sum']==0.0 ] = 0.0
934 ## SD : If some chans are zero.
935 ## Set the wt to zero and use the logical AND of int_wt and sd_wt for feathering?
936 zerochans = np.count_nonzero( sd_stat['sum']==0.0 )
937 if zerochans>0:
938 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. ')
939 chanwt[ sd_stat['sum']==0.0 ] = 0.0
941 ### If there are channels to flag.... list them.
942 if np.count_nonzero( chanwt==0.0 ):
943 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')
944 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')
947 if np.count_nonzero( chanwt != 0.0 ) == 0:
948 casalog.post("There are no channels with data in both the INT and the SD cubes. Cannot proceed","WARN")
949 validity=False
951 pars['chanwt'] = chanwt
952 return validity, pars
955 def setup_cube_params(self,sdcube=''):
956 """
957 Read coordinate system from the SD cube
958 Decide parameters to input into sdintimaging for the INT cube to match.
960 This is a helper method, not currently used in the sdintimaging task.
961 We will add this later (after 6.1), and also remove some parameters from the task level.
962 """
963 _ia.open(sdcube)
964 csys = _ia.coordsys()
965 shp = _ia.shape()
966 ctypes = csys.axiscoordinatetypes()
967 casalog.post("Shape of SD cube : "+str(shp))
968 casalog.post("Coordinate ordering : "+str(ctypes))
969 if len(ctypes) !=4 or ctypes[3] != 'Spectral':
970 casalog.post("The SD cube needs to have 4 axes, in the RA/DEC/Stokes/Spectral order", 'ERROR')
971 _ia.close()
972 return False
973 nchan = shp[3]
974 start = str( csys.referencevalue()['numeric'][3] ) + csys.units()[3]
975 width = str( csys.increment()['numeric'][3]) + csys.units()[3]
976 ## Number of channels
977 casalog.post("nchan = "+str(nchan))
978 ## Start Frequency
979 casalog.post("start = " + start )
980 ## Width
981 casalog.post("width = " + width )
982 ## Test for restoringbeams
983 rbeams = _ia.restoringbeam()
984 #if not rbeams.has_key('nChannels') or rbeams['nChannels']!= shp[3]:
985 if not 'nChannels' in rbeams or rbeams['nChannels']!= shp[3]:
986 casalog.post("The SD Cube needs to have per-plane restoringbeams", 'ERROR')
987 _ia.close()
988 return False
989 else:
990 casalog.post("Found " + str(rbeams['nChannels']) + " per-plane restoring beams")
991 casalog.post("\n(For specmode='mfs' in sdintimaging, please remember to set 'reffreq' to a value within the freq range of the cube)\n")
992 _ia.close()
993 return {'nchan':nchan, 'start':start, 'width':width}
998### Using Old Imager. Does not work for cubes ?
999# def fit_psf_beam(self,msname = '', psfname =''):
1000# _im.open(msname)
1001# _ia.open(psfname)
1002# csys = _ia.coordsys()
1003# rbeam_old = _ia.restoringbeam()
1004# print(rbeam_old)
1005# shp = _ia.shape()
1006# _ia.close()
1007# cellx = csys.increment()['numeric'][0];
1008# celly = csys.increment()['numeric'][1];
1009# _im.defineimage(nx=shp[0],ny=shp[1],cellx=str(cellx)+'rad',celly=str(celly)+'rad',nchan=3)
1010# params =_im.fitpsf(psfname)
1011# print(params)
1012# _im.close()
1013#