Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/sdint_helper.py: 6%
522 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 17:39 +0000
2from scipy import fftpack
3import numpy as np
4import shutil
5import os
6import time
8from casatools import quanta, table, image, regionmanager, imager
9from casatasks import casalog, imsubimage, feather
11_ia = image()
12_qa = quanta()
13_rg = regionmanager()
15_mytb = table()
17class SDINT_helper:
19# def __init__(self):
20# casalog.post('Init Helper')
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')
31 return freqaxis_index
33 def getFreqList(self,imname=''):
34 """ Get the list of frequencies for the given image, one for each channel.
36 Returns:
37 list[float] The frequencies for each channel in the image, in Hz.
38 """
40 _ia.open(imname)
41 csys =_ia.coordsys()
42 shp = _ia.shape()
43 freqaxis_index = self.getFreqAxisIndex()
44 _ia.close()
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;
58 csys.done()
59 return freqlist
61################################################
63 def copy_restoringbeam(self,fromthis='',tothis=''):
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()
76################################################
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.
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....
85 TODO : Add the effdishdia usage to get freq-indep feathering.
87 """
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)
95 freqlist = self.getFreqList(sdcube)
97 os.system('rm -rf '+jointcube)
98 os.system('cp -r ' + intcube + ' ' + jointcube)
100 _ib = image()
102 _ia.open(jointcube)
103 _ia.set(0.0) ## Initialize this to zero for all planes
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')
110 freqdishdia = dishdia ## * (freqlist[0] / freqlist[i]) # * 0.5
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));
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
132 _ib.open('tmp_jointplane')
133 pixjoint = _ib.getchunk()
134 _ib.close()
135 _ia.putchunk(pixjoint, blc=[0,0,0,i])
137 _ia.close()
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)
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()
163 if shp[0]>1:
164 casalog.post("WARNING : Cannot use this task with faceting", 'WARN')
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()
172 _ia.open(jointname+'.sumwt')
173 _ia.putchunk(vals)
174 _ia.close()
176 casalog.post("********************Re-norm with "+str(vals))
179################################################
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()
191 casalog.post("********************Re-norm with "+str(vals))
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()
206################################################
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
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))
223 freqlist = self.getFreqList(inpcube)
225 _ia.open(inpcube)
226 shp=_ia.shape()
227 freqaxis_index = self.getFreqAxisIndex()
228 _ia.close()
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')
233 ##############
234 ### Calculate a reference Primary Beam
235 ### Weighted sum of pb cube
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)
245# pbplane = np.zeros( (shp[0],shp[1]), 'float')
247 if freqdep==False:
248 _ia.open(cubewt) # .sumwt
249 cwt = _ia.getchunk()[0,0,0,:]
250 _ia.close()
252 if shp[3] != len(cwt) or len(freqlist) != len(cwt):
253 raise Exception("Modify with PB : Nchan shape mismatch between cube and sumwt.")
255 cwt = cwt * chanwt ## Merge the weights and flags
257 sumchanwt = np.sum(cwt)
259 if sumchanwt==0:
260 raise Exception("Weights are all zero ! ")
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()
268 pbplane = pbplane / sumchanwt
270 ##############
273 ## Special-case for setting the PBmask to be same for all freqs
274 if freqdep==False:
275 shutil.copytree(pbcube, pbcube+'_tmpcopy')
277 for i in range(len(freqlist)):
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()
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()
291 _ia.open(inpcube)
292 implane = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
294 outplane = pbplane.copy()
295 outplane.fill(0.0)
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
305 _ia.putchunk(outplane, blc=[0,0,0,i])
306 _ia.close()
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')
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);
336 elif defaultmaskname=='mask0':
337 if 'pbmask' in allmasknames:
338 _ia.maskhandler('delete','pbmask')
339 _ia.calcmask(mask='"'+pbimage+'"'+'>'+str(pblimit), name='pbmask');
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)
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()
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 ])
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 """
379 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value']
381 # casalog.post("&&&&&&&&& REF FREQ : " + str(refnu))
383 pix=[]
385 num_terms=nterms
387 if dopsf==True:
388 num_terms=2*nterms-1
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)
396 _ia.open(cubename)
397 shp = _ia.shape()
398 freqaxis_index = self.getFreqAxisIndex()
399 _ia.close()
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')
404 _ia.open(cubewt)
405 cwt = _ia.getchunk()[0,0,0,:]
406 _ia.close()
409 freqlist = self.getFreqList(cubename)
411 if shp[3] != len(cwt) or len(freqlist) != len(cwt):
412 raise Exception("Nchan shape mismatch between cube and sumwt.")
414 cwt = cwt * chanwt ## Merge the weights and flags.
416 sumchanwt = np.sum(cwt) ## This is a weight
418 if sumchanwt==0:
419 raise Exception("Weights are all zero ! ")
420 else:
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]
430 for tt in range(0,num_terms):
431 pix[tt] = pix[tt]/sumchanwt
432# ia.close()
434 for tt in range(0,num_terms):
435 _ia.open(mtname+'.tt'+str(tt))
436 _ia.putchunk(pix[tt])
437 _ia.close()
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 """
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()
456 refnu = _qa.convert( _qa.quantity(reffreq) ,'Hz' )['value']
458 pix=[]
460 for tt in range(0,nterms):
461 _ia.open(mtname+'.model.tt'+str(tt))
462 pix.append( _ia.getchunk() )
463 _ia.close()
465 _ia.open(cubename+'.model')
466 shp = _ia.shape()
467 _ia.close()
469 implane = pix[0].copy()
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()
481##################################################
484 def calc_sd_residual(self,origcube='', modelcube = '', residualcube = '', psfcube=''):
485##, pbcube='', applypb=False, pblimit=0.2):
486 """
487 Residual = Original - ( PSF * Model )
488 """
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)
500 freqlist = self.getFreqList(origcube)
502 _ia.open(origcube)
503 shp = _ia.shape()
504 freqaxis_index = self.getFreqAxisIndex()
505 _ia.close()
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')
510 for i in range(0,len(freqlist)):
512 _ia.open(origcube)
513 im_orig = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
514 _ia.close()
516 _ia.open(modelcube)
517 im_model = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
518 _ia.close()
520 _ia.open(psfcube)
521 im_psf = _ia.getchunk(blc=[0,0,0,i],trc=[shp[0],shp[1],0,i])
522 _ia.close()
524 smoothedim = self.myconvolve(im_model[:,:,0,0], im_psf[:,:,0,0])
526 if( np.nansum(im_orig)==0.0):
527 smoothedim.fill(0.0)
529 im_residual=im_psf.copy() ## Just to init the shape of this thing
530 im_residual[:,:,0,0] = im_orig[:,:,0,0] - smoothedim
532 _ia.open(residualcube)
533 _ia.putchunk(im_residual, blc=[0,0,0,i])
534 _ia.close()
536# ia_orig.close()
537# ia_model.close()
538# ia_residual.close()
539# ia_psf.close()
541##################################################
543 def myconvolve(self,im1,im2):
545 t3 = time.time()
547 shp = im1.shape
548 pshp = (shp[0]*2, shp[1]*2)
549 pim1 = np.zeros(pshp,'float')
550 pim2 = np.zeros(pshp,'float')
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
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 ) ) ) )
560 smoothedim = psmoothedim[shp[0]-shp[0]//2 : shp[0]+shp[0]//2, shp[1]-shp[1]//2 : shp[1]+shp[1]//2]
562 return smoothedim
564##########################################
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()
574 _myia.open(imagename)
575 casalog.post("imagename="+imagename)
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
589 finally:
590 csys.done()
591 if outia != None and outia.isopen():
592 outia.done()
593 _myia.done()
596##########################################
598 def createplaneimage(self,imagename, outfile, chanid):
599 """
600 extract a channel plane image
601 """
602 _tmpia=image()
603 _tmprg=regionmanager()
604 outia=None
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')
613 _tmpia.close()
614 _tmpia.done()
615 _tmprg.done()
616 if outia != None and outia.isopen():
617 outia.done()
619 def pbcor(self, imagename, pbimage, cutoff, outfile):
620 """
621 pb-correction
622 """
623 outia=None
624 _myia=image()
625 _myia.open(imagename)
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')
633 finally:
634 _myia.done()
635 if outia != None and outia.isopen():
636 outia.done()
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
660 indir = incsys['direction0']
661 refdir = refcsys['direction0']
663 diff_ra = indir['crval'][0] - refdir['crval'][0]
664 diff_dec = indir['crval'][1] - refdir['crval'][1]
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))
690 else:
691 casalog.post(" The center of psf coincides with int psf: (diffRA,diffDec)=( %s, %s)" % (diff_ra, diff_dec))
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.")
705 #return modpsf
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
715 _ia = image()
716 _cl = componentlist()
717 _rg = regionmanager()
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
736 csys = _ia.coordsys()
737 _ia.close()
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')
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')
756 cdir = csys.torecord()['direction0']
757 compdir = [cdir['system'] , str(cdir['crval'][0])+cdir['units'][0] , str(cdir['crval'][1])+cdir['units'][1] ]
759 ## Make empty SD psf cube from SD image cube
760 os.system('rm -rf '+sdpsfname)
761 os.system('cp -R '+sdimage+' '+sdpsfname)
763 ## Iterate through PSF cube and replace pixels with Gaussians matched to restoringbeam info
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']
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')
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])
791 _ia.close()
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
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
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()
815 _ia.open(sdres)
816 sd_shp = _ia.shape()
817 sd_stat = _ia.statistics(axes=[0,1])
818 _ia.close()
820 _ia.open(intwt)
821 sumwt = _ia.getchunk()
822 _ia.close()
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')
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")
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
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'
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
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
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')
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
880 pars['chanwt'] = chanwt
881 return validity, pars
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.
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}
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#