|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os 00002 import pdb 00003 from taskinit import * 00004 import time 00005 import numpy as np 00006 from cleanhelper import * 00007 from casac import casac 00008 class imagecont(): 00009 def __init__(self, ftmachine='ft', wprojplanes=10, facets=1, pixels=[3600, 3600], cell=['3arcsec', '3arcsec'], spw='', field='', phasecenter='', weight='natural', robust=0.0, stokes='I', npixels=0, uvtaper=False, outertaper=[], timerange='', uvrange='', baselines='', scan='', observation='', gain=0.1, numthreads=-1, pbcorr=False, minpb=0.2, cyclefactor=1.5): 00010 self.im=casac.imager() 00011 self.imperms={} 00012 self.dc=casac.deconvolver() 00013 self.ft=ftmachine 00014 self.origms='' 00015 self.wprojplanes=wprojplanes 00016 self.facets=facets 00017 print 'cell', cell 00018 self.pixels=pixels 00019 self.cell=cell 00020 if(spw==''): 00021 spw='*' 00022 self.spw=spw 00023 self.field=field 00024 self.phCen=phasecenter 00025 self.weight=weight 00026 self.imageparamset=False 00027 self.robust=robust 00028 self.weightnpix=npixels 00029 self.stokes=stokes 00030 self.numthreads=numthreads 00031 self.imagetilevol=1000000 00032 self.visInMem=False 00033 self.painc=360.0 00034 self.pblimit=0.1 00035 self.dopbcorr=False 00036 self.cyclefactor=cyclefactor 00037 self.novaliddata={} 00038 self.applyoffsets=False 00039 self.cfcache='cfcache.dir' 00040 self.epjtablename='' 00041 self.gain=gain 00042 self.uvtaper=uvtaper 00043 self.outertaper=outertaper 00044 self.timerange=timerange 00045 self.uvrange=uvrange 00046 self.baselines=baselines 00047 self.scan=scan 00048 self.observation=observation 00049 self.pbcorr=pbcorr 00050 self.minpb=minpb 00051 00052 def setextraoptions(self, im, cyclefactor=1.5, fluxscaleimage='', scaletype='SAULT'): 00053 im.setmfcontrol(scaletype=scaletype, fluxscale=[fluxscaleimage], 00054 cyclefactor=cyclefactor, minpb=self.minpb) 00055 if(self.pbcorr): 00056 im.setvp(True) 00057 im.setmfcontrol(scaletype='NONE', cyclefactor=cyclefactor, minpb=self.minpb, fluxscale=[fluxscaleimage]) 00058 00059 def setparamcont(self, im, freq, band, singleprec=False): 00060 im.defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], 00061 celly=self.cell[1], phasecenter=self.phCen, mode='frequency', 00062 nchan=1, start=freq, step=band, facets=self.facets, 00063 stokes=self.stokes) 00064 im.weight(type=self.weight, rmode='norm', npixels=self.weightnpix, 00065 robust=self.robust) 00066 im.setoptions(ftmachine=self.ft, wprojplanes=self.wprojplanes, 00067 pastep=self.painc, pblimit=self.pblimit, 00068 cfcachedirname=self.cfcache, dopbgriddingcorrections=self.dopbcorr, 00069 applypointingoffsets=self.applyoffsets, 00070 imagetilevol=self.imagetilevol, 00071 singleprecisiononly=singleprec, numthreads=self.numthreads) 00072 00073 def imagecontmultims(self, msnames=[''], start=0, numchan=1, spw=0, field=0, freq='1.20GHz', band='200MHz', imname='newmodel'): 00074 im=self.im 00075 ###either psf 0 or no channel selected 00076 if(self.novaliddata==True): 00077 return 00078 else: 00079 self.novaliddata=False 00080 if(not self.imageparamset): 00081 try: 00082 conlis=lambda a,n: a if (type(a)==list) else [a]*n 00083 if (type(msnames)==str): 00084 msnames=[msnames] 00085 numms=len(msnames) 00086 start=conlis(start,numms) 00087 numchan=conlis(numchan, numms) 00088 spw=conlis(spw, numms) 00089 field=conlis(field, numms) 00090 numfail=0 00091 for k in range(numms): 00092 try: 00093 im.selectvis(vis=msnames[k], field=field[k], spw=spw[k], nchan=numchan[k], start=start[k], step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00094 except: 00095 numfail+=1 00096 if(numfail==numms): 00097 self.novaliddata=True 00098 except Exception, instance: 00099 ###failed to selectdata 00100 self.novaliddata=True 00101 self.setparamcont(im, freq, band, singleprec=False) 00102 if((len(numchan)==0) or (np.sum(numchan)==0)): 00103 self.novaliddata=True 00104 self.makecontimage(im, self.novaliddata, imname) 00105 self.imageparamset=True 00106 00107 00108 00109 #### uniform and brigg's weight combination functions 00110 def getweightgrid(self, msname='', wgttype='imaging'): 00111 if(self.imperms.has_key(msname) and (not self.novaliddata[msname])): 00112 return self.imperms[msname].getweightgrid(type=wgttype) 00113 def setweightgrid(self, msname='', weight=[], wgttype='imaging'): 00114 if(self.imperms.has_key(msname) and (not self.novaliddata[msname])): 00115 self.imperms[msname].setweightgrid(weight=weight, type=wgttype) 00116 def getftweight(self, msname='', wgtimage=['']): 00117 if(self.imperms.has_key(msname) and (not self.novaliddata[msname])): 00118 self.imperms[msname].getweightgrid(type='ftweight', wgtimages=wgtimage) 00119 #### 00120 def imagecont(self, msname='spw00_4chan351rowTile.ms', start=[0], numchan=[-1], spw='', field=0, freq='1.20GHz', band='200MHz', imname='newmodel', nterms=1, scales=[0]): 00121 #casalog.post('KEYS '+str(self.imperms.keys())) 00122 if(not self.imperms.has_key(msname)): 00123 self.imageparamset=False 00124 im=casac.imager() 00125 self.imperms[msname]=im 00126 self.novaliddata[msname]=False 00127 #casalog.post('MSNAME '+msname) 00128 else: 00129 #casalog.post('reMSNAME '+msname) 00130 im=self.imperms[msname] 00131 self.imageparamset=True 00132 ###either psf 0 or no channel selected 00133 if(self.novaliddata[msname]): 00134 return 00135 #j=start 00136 #end=start+numchan-1 00137 #spwstring=str(spw)+':'+str(start)+'~'+str(end) 00138 #print 'spwstring', spwstring 00139 if(not self.imageparamset): 00140 self.origms=msname 00141 try: 00142 im.selectvis(vis=msname, field=field, spw=spw, nchan=numchan, start=start, step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00143 if(self.uvtaper): 00144 im.filter(type='gaussian', bmaj=self.outertaper[0], 00145 bmin=self.outertaper[1], bpa=self.outertaper[2]) 00146 except Exception, instance: 00147 ###failed to selectdata 00148 self.novaliddata[msname]=True 00149 #### 00150 #imname=imname+'_%02d'%(j) 00151 self.setparamcont(im, freq, band) 00152 if(self.ft=='mosaic'): 00153 self.setextraoptions(im, fluxscaleimage=imname+'.flux', scaletype='SAULT') 00154 if((len(numchan)==0) or (np.sum(numchan)==0)): 00155 self.novaliddata[msname]=True 00156 if(nterms==1): 00157 self.makecontimage(im, self.novaliddata[msname], imname) 00158 else: 00159 self.makemtcontimage(im, imname, nterms, scales, freq) 00160 self.imageparamset=True 00161 00162 def imagecontbychan(self, msname='spw00_4chan351rowTile.ms', start=[0], numchan=[1], spw=[0], field=0, freq='1.20GHz', band='200MHz', imname='newmodel'): 00163 ia=casac.image() 00164 if(type(spw) == int): 00165 spw=[spw] 00166 start=[start] 00167 numchan=[numchan] 00168 totchan=np.sum(numchan) 00169 if(not self.imageparamset): 00170 if(totchan==0): 00171 ###make blanks 00172 self.im.selectvis(vis=msname, field=field, spw=spw, nchan=numchan, start=start, step=1, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation,writeaccess=False) 00173 self.im.defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], celly=self.cell[1], phasecenter=self.phCen, mode='frequency', nchan=1, start=freq, step=band, facets=self.facets) 00174 self.im.make(imname+'.image') 00175 self.im.make(imname+'.residual') 00176 self.im.make(imname+'.model') 00177 self.im.make(imname+'.psf') 00178 self.imageparamset=True 00179 return 00180 del self.im 00181 self.im=[] 00182 del self.novaliddata 00183 self.novaliddata=[] 00184 selkey={} 00185 spwind=0 00186 chancounter=0 00187 for k in range(totchan): 00188 if(chancounter==numchan[spwind]): 00189 chancounter=0 00190 spwind+=1 00191 selkey[k]= {'spw': spw[spwind], 'start':start[spwind]+chancounter} 00192 chancounter += 1 00193 self.im.append(casac.imager()) 00194 self.novaliddata.append(False) 00195 print 'selkey', selkey 00196 else: 00197 if(totchan==0): 00198 return 00199 origname=msname 00200 #pdb.set_trace() 00201 for k in range(totchan): 00202 ###either psf 0 or no channel selected 00203 if(k==0): 00204 image=imname+'.image' 00205 residual=imname+'.residual' 00206 psf=imname+'.psf' 00207 else: 00208 image=imname+'_'+str(k)+'.image' 00209 residual=imname+'_'+str(k)+'.residual' 00210 psf=imname+'_'+str(k)+'.psf' 00211 if(not self.novaliddata[k]): 00212 msname=origname 00213 if(not self.imageparamset): 00214 self.im[k].selectvis(vis=msname, field=field, spw=selkey[k]['spw'], nchan=1, start=selkey[k]['start'], step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00215 self.im[k].defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], celly=self.cell[1], phasecenter=self.phCen, mode='frequency', nchan=1, start=freq, step=band, facets=self.facets) 00216 self.im[k].weight(type=self.weight, rmode='norm', npixels=self.weightnpix, 00217 robust=self.robust) 00218 self.im[k].setoptions(ftmachine=self.ft, wprojplanes=self.wprojplanes, pastep=self.painc, pblimit=self.pblimit, cfcachedirname=self.cfcache, dopbgriddingcorrections=self.dopbcorr, applypointingoffsets=self.applyoffsets, imagetilevol=self.imagetilevol, singleprecisiononly=True, numthreads=self.numthreads) 00219 #im.regionmask(mask='lala.mask', boxes=[[0, 0, 3599, 3599]]) 00220 #im.setmfcontrol(cyclefactor=0.0) 00221 if(not self.imageparamset): 00222 try: 00223 self.im[k].clean(algorithm='clark', gain=self.gain, niter=0, model=imname+'.model', image=image, residual=residual, psfimage=psf) 00224 if(k > 0): 00225 ia.open(imname+'.psf') 00226 ia.calc('"'+imname+'.psf'+'" + "'+psf+'"') 00227 ia.done() 00228 except Exception, instance: 00229 if(string.count(instance.message, 'PSFZero') >0): 00230 self.novaliddata[k]=True 00231 ###make a blank image 00232 if(k==0): 00233 self.im[k].make(image) 00234 self.im[k].make(residual) 00235 else: 00236 raise instance 00237 else: 00238 if(not self.novaliddata[k]): 00239 self.im[k].restore(model=imname+'.model', image=image, residual=residual) 00240 if(k > 0): 00241 ia.open(imname+'.residual') 00242 ia.calc('"'+imname+'.residual'+'" + "'+residual+'"') 00243 ia.done() 00244 ia.open(imname+'.image') 00245 ia.calc('"'+imname+'.image'+'" + "'+residual+'"') 00246 ia.done() 00247 tb.showcache() 00248 00249 ia.open(imname+'.residual') 00250 ia.calc('"'+imname+'.residual'+'"/'+str(totchan)) 00251 ia.done() 00252 ia.open(imname+'.image') 00253 ia.calc('"'+imname+'.image'+'"/'+str(totchan)) 00254 ia.done() 00255 if(not self.imageparamset): 00256 ia.open(imname+'.psf') 00257 ia.calc('"'+imname+'.psf'+'"/'+str(totchan)) 00258 ia.done() 00259 #im.done() 00260 del ia 00261 tb.showcache() 00262 self.imageparamset=False 00263 00264 def cleancont(self, niter=100, alg='clark', thr='0.0mJy', psf='newmodel.psf', dirty='newmodel.dirty', model='newmodel.model', mask='', scales=[0]): 00265 dc=casac.deconvolver() 00266 dc.open(dirty=dirty, psf=psf) 00267 if((alg=='hogbom') or (alg == 'msclean') or (alg == 'multiscale')): 00268 sca=scales if ((alg=='msclean') or (alg=='multiscale')) else [0] 00269 dc.setscales(scalemethod='uservector', uservector=sca) 00270 alg='fullmsclean' 00271 retval={} 00272 if(alg=='clark'): 00273 retval=dc.fullclarkclean(niter=niter, gain=self.gain, threshold=thr, model=model, mask=mask) 00274 else: 00275 retval=dc.clean(algorithm=alg, gain=self.gain, niter=niter, threshold=thr, model=model, mask=mask) 00276 dc.done() 00277 del dc 00278 return retval 00279 00280 def imagechan(self, msname='spw00_4chan351rowTile.ms', start=0, numchan=1, spw=0, field=0, imroot='newmodel', imchan=0, niter=100, alg='clark', thr='0.0mJy', mask='', majcycle=1, scales=[0], fstart='1GHz', width='10MHz', chanchunk=1): 00281 origname=msname 00282 # a=cleanhelper() 00283 imname=imroot+str(imchan) 00284 maskname='' 00285 if(mask != ''): 00286 maskname=imname+'.mask' 00287 if( not self.makemask(inmask=mask, outmask=maskname, imchan=imchan, chanchunk=chanchunk)): 00288 maskname='' 00289 # a.getchanimage(cubeimage=imroot+'.model', outim=imname+'.model', chan=imchan) 00290 im.selectvis(vis=msname, field=field, spw=spw, nchan=numchan, start=start, step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00291 im.defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], celly=self.cell[1], phasecenter=self.phCen, facets=self.facets, mode='frequency', nchan=chanchunk, start=fstart, step=width) 00292 im.weight(type=self.weight, rmode='norm', npixels=self.weightnpix, 00293 robust=self.robust) 00294 00295 im.setoptions(ftmachine=self.ft, wprojplanes=self.wprojplanes, imagetilevol=self.imagetilevol, singleprecisiononly=True, numthreads=self.numthreads) 00296 im.setscales(scalemethod='uservector', uservector=scales) 00297 self.setextraoptions(im, cyclefactor=0.0) 00298 majcycle = majcycle if (niter/majcycle) >0 else 1 00299 for kk in range(majcycle): 00300 im.clean(algorithm=alg, gain=self.gain, niter= (niter/majcycle), threshold=thr, model=imname+'.model', image=imname+'.image', residual=imname+'.residual', mask=maskname, psfimage='') 00301 im.done() 00302 if(maskname != ''): 00303 shutil.rmtree(maskname, True) 00304 # a.putchanimage(cubimage=imroot+'.model', inim=imname+'.model', 00305 # chan=imchan) 00306 # a.putchanimage(cubimage=imroot+'.residual', inim=imname+'.residual', 00307 # chan=imchan) 00308 # a.putchanimage(cubimage=imroot+'.image', inim=imname+'.image', 00309 # chan=imchan) 00310 ###need to clean up here the channel image...will do after testing phase 00311 00312 def imagechan_selfselect(self, msname='spw00_4chan351rowTile.ms', spwids=[0], field=0, imroot='newmodel', imchan=0, niter=100, alg='clark', thr='0.0mJy', mask='', majcycle=1, scales=[0], chanchunk=1, startchan=0): 00313 ###need generate fstart, width from 1 channel width, start, numchan and spw 00314 origname=msname 00315 retval={} 00316 retval['maxresidual']=0.0 00317 retval['iterations']=0 00318 retval['converged']=False 00319 ia.open(imroot+'.model') 00320 csys=ia.coordsys() 00321 ia.done() 00322 im=casac.imager() 00323 fstart=csys.toworld([0,0,0,startchan],'n')['numeric'][3] 00324 fstep=csys.toworld([0,0,0,startchan+1],'n')['numeric'][3]-fstart 00325 fend=fstep*(chanchunk-1)+fstart 00326 print 'fstat bef findchansel ', fstart 00327 spw, start, nchan=self.findchanselLSRK(msname=msname, spw=spwids, 00328 field=field, 00329 numpartition=1, 00330 beginfreq=fstart, endfreq=fend, chanwidth=fstep) 00331 print 'spw, start, nchan', spw, start, nchan 00332 if((len(spw[0])==0) or (len(nchan[0])==0) or (len(start[0]) ==0) ): 00333 return retval 00334 imname=imroot+str(imchan) 00335 maskname='' 00336 if(mask != ''): 00337 maskname=imname+'.mask' 00338 if( not self.makemask(inmask=mask, outmask=maskname, imchan=imchan, chanchunk=chanchunk, startchan=startchan)): 00339 maskname='' 00340 # a.getchanimage(cubeimage=imroot+'.model', outim=imname+'.model', chan=imchan) 00341 im.selectvis(vis=msname, field=field, spw=spw[0], nchan=nchan[0], start=start[0], step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00342 print 'fstart bef def', fstart 00343 if(self.uvtaper): 00344 im.filter(type='gaussian', bmaj=self.outertaper[0], 00345 bmin=self.outertaper[1], bpa=self.outertaper[2]) 00346 im.defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], celly=self.cell[1], phasecenter=self.phCen, facets=self.facets, mode='frequency', nchan=chanchunk, start=str(fstart)+'Hz', step=str(fstep)+'Hz', outframe='LSRK') 00347 im.weight(type=self.weight, rmode='norm', npixels=self.weightnpix, 00348 robust=self.robust) 00349 00350 im.setoptions(ftmachine=self.ft, wprojplanes=self.wprojplanes, imagetilevol=self.imagetilevol, singleprecisiononly=True, numthreads=self.numthreads) 00351 im.setscales(scalemethod='uservector', uservector=scales) 00352 fluxscaleimage=imname+'.flux' if(self.ft=='mosaic') else '' 00353 majcycle = majcycle if (niter/majcycle) >0 else 1 00354 self.setextraoptions(im, fluxscaleimage=fluxscaleimage, cyclefactor=(0.0 if(majcycle >1) else self.cyclefactor)) 00355 #### non simple ft machine...should use mf 00356 if(self.ft != 'ft'): 00357 alg='mf'+alg 00358 try: 00359 00360 for kk in range(majcycle): 00361 retval=im.clean(algorithm=alg, gain=self.gain, niter= (niter/majcycle), threshold=thr, model=imname+'.model', image=imname+'.image', residual=imname+'.residual', mask=maskname, psfimage=imname+'.psf') 00362 im.done() 00363 del im 00364 except Exception as instance: 00365 im.done() 00366 del im 00367 if(string.count(str(instance), 'PSFZero') <1): 00368 raise Exception(instance) 00369 if(maskname != ''): 00370 shutil.rmtree(maskname, True) 00371 return retval 00372 #self.putchanimage(imroot+'.model', imname+'.model', imchan*chanchunk) 00373 #self.putchanimage(imroot+'.residual', imname+'.residual', imchan*chanchunk) 00374 #self.putchanimage(imroot+'.image', imname+'.image', imchan*chanchunk) 00375 00376 00377 00378 def imagechan_new(self, msname='spw00_4chan351rowTile.ms', start=0, numchan=1, spw=0, field=0, cubeim='imagecube', imroot='newmodel', imchan=0, chanchunk=1, niter=100, alg='clark', thr='0.0mJy', mask='', majcycle=1, scales=[0]): 00379 origname=msname 00380 # a=cleanhelper() 00381 imname=imroot+str(imchan) 00382 maskname='' 00383 if(mask != ''): 00384 maskname=imname+'.mask' 00385 if( not self.makemask(inmask=mask, outmask=maskname, imchan=imchan, chanchunk=chanchunk)): 00386 maskname='' 00387 self.getchanimage(inimage=cubeim+'.model', outimage=imname+'.model', chan=imchan*chanchunk, nchan=chanchunk) 00388 im.selectvis(vis=msname, field=field, spw=spw, nchan=numchan, start=start, step=1, datainmemory=self.visInMem, time=self.timerange, uvrange=self.uvrange, baseline=self.baselines, scan=self.scan, observation=self.observation, writeaccess=False) 00389 im.weight(type=self.weight, rmode='norm', npixels=self.weightnpix, 00390 robust=self.robust) 00391 im.defineimage(nx=self.pixels[0], ny=self.pixels[1], cellx=self.cell[0], celly=self.cell[1], phasecenter=self.phCen, facets=self.facets) 00392 im.setoptions(ftmachine=self.ft, wprojplanes=self.wprojplanes, imagetilevol=self.imagetilevol, singleprecisiononly=True, numthreads=self.numthreads) 00393 im.setscales(scalemethod='uservector', uservector=scales) 00394 self.setextraoptions(im, cyclefactor=0.0) 00395 majcycle = majcycle if (niter/majcycle) >0 else 1 00396 for kk in range(majcycle): 00397 im.clean(algorithm=alg, gain=self.gain, niter= (niter/majcycle), threshold=thr, model=imname+'.model', image=imname+'.image', residual=imname+'.residual', 00398 mask=maskname, psfimage='') 00399 im.done() 00400 if(maskname != ''): 00401 shutil.rmtree(maskname, True) 00402 #self.putchanimage(cubeim+'.model', imname+'.model', imchan*chanchunk) 00403 #self.putchanimage(cubeim+'.residual', imname+'.residual', imchan*chanchunk) 00404 #self.putchanimage(cubeim+'.image', imname+'.image', imchan*chanchunk) 00405 00406 00407 @staticmethod 00408 def getallchanmodel(inimage, chanchunk=1): 00409 #tim1=time.time() 00410 ia.open(inimage+'.model') 00411 modshape=ia.shape() 00412 nchan=modshape[3]/chanchunk; 00413 if(nchan >1): 00414 blc=[0,0,modshape[2]-1,0] 00415 trc=[modshape[0]-1,modshape[1]-1,modshape[2]-1,0+chanchunk] 00416 for k in range(nchan): 00417 #timbeg=time.time() 00418 blc[3]=k*chanchunk 00419 trc[3]=(k+1)*chanchunk-1 00420 if((trc[3]) >= modshape[3]): 00421 trc[3]=modshape[3]-1 00422 sbim=ia.subimage(outfile=inimage+str(k)+'.model', region=rg.box(blc,trc), overwrite=True) 00423 sbim.done() 00424 #print 'time taken for ', k, time.time()-timbeg 00425 ia.done() 00426 tb.clearlocks() 00427 00428 @staticmethod 00429 def getchanimage(inimage, outimage, chan, nchan=1): 00430 """ 00431 create a slice of channels image from cubeimage 00432 """ 00433 #pdb.set_trace() 00434 ia.open(inimage) 00435 modshape=ia.shape() 00436 if (modshape[3]==1) or (chan > (modshape[3]-1)) : 00437 return False 00438 if((nchan+chan) < modshape[3]): 00439 endchan= chan+nchan-1 00440 else: 00441 endchan=modshape[3]-1 00442 blc=[0,0,modshape[2]-1,chan] 00443 trc=[modshape[0]-1,modshape[1]-1,modshape[2]-1,endchan] 00444 sbim=ia.subimage(outfile=outimage, region=rg.box(blc,trc), overwrite=True) 00445 sbim.done() 00446 ia.done() 00447 tb.clearlocks() 00448 #casalog.post('getLOCKS: '+ str(inimage)+ ' --- ' + str(tb.listlocks())) 00449 return True 00450 #getchanimage = staticmethod(getchanimage) 00451 def cleanupcubeimages(self, readyputchan, doneputchan, imagename, nchanchunk, chanchunk): 00452 """ 00453 This function will put the True values of readyputchan into the final cubes and set the doneputchan to True 00454 Confused ..read the code 00455 """ 00456 for k in range(nchanchunk): 00457 if(readyputchan[k] and (not doneputchan[k])): 00458 self.putchanimage(imagename+'.model', imagename+str(k)+'.model', k*chanchunk, False) 00459 self.putchanimage(imagename+'.residual', imagename+str(k)+'.residual', k*chanchunk, False) 00460 self.putchanimage(imagename+'.image', imagename+str(k)+'.image', k*chanchunk, False) 00461 doneputchan[k]=True 00462 00463 def cleanupimages(self, readyputchan, imagename, nchanchunk, chanchunk, removefile=True): 00464 self.cleanupmodelimages(readyputchan, imagename, nchanchunk, chanchunk, removefile) 00465 self.cleanupresidualimages(readyputchan, imagename, nchanchunk, chanchunk, removefile) 00466 self.cleanuprestoredimages(readyputchan, imagename, nchanchunk, chanchunk, removefile) 00467 00468 def cleanupmodelimages(self, readyputchan, imagename, nchanchunk, chanchunk, removefile=True): 00469 """ 00470 This function will put model images only 00471 """ 00472 for k in range(nchanchunk): 00473 if(readyputchan[k]): 00474 self.putchanimage(imagename+'.model', imagename+str(k)+'.model', k*chanchunk, removefile) 00475 def cleanupresidualimages(self, readyputchan, imagename, nchanchunk, chanchunk, removefile=True): 00476 """ 00477 This function will put residual images only 00478 """ 00479 for k in range(nchanchunk): 00480 if(readyputchan[k]): 00481 self.putchanimage(imagename+'.residual', imagename+str(k)+'.residual', k*chanchunk, removefile) 00482 00483 def cleanuprestoredimages(self, readyputchan, imagename, nchanchunk, chanchunk, removefile=True): 00484 """ 00485 This function will put residual images only 00486 """ 00487 for k in range(nchanchunk): 00488 if(readyputchan[k]): 00489 self.putchanimage(imagename+'.image', imagename+str(k)+'.image', k*chanchunk, removefile) 00490 00491 @staticmethod 00492 def concatimages(cubeimage, inim, csys=None, removeinfile=True): 00493 if((csys==None) and os.path.exists(cubeimage)): 00494 ia.open(cubeimage) 00495 csys=ia.coordsys() 00496 ia.done() 00497 if(type(inim) != list): 00498 inim=[inim] 00499 ###Temporary bypass of CAS-4423 00500 elbeamo={} 00501 for elim in inim: 00502 ia.open(elim) 00503 beam=ia.restoringbeam() 00504 #if( (len(beam) > 0) and (not beam.has_key('beams'))): 00505 # ia.setrestoringbeam(remove=True) 00506 # nchan=ia.shape()[3] 00507 # for k in range(nchan): 00508 # ia.setrestoringbeam(major=beam['major'], minor=beam['minor'], pa=beam['positionangle'], channel=k, polarization=0) 00509 ia.setrestoringbeam(remove=True) 00510 elbeamo=beam if(not beam.has_key('beams')) else beam['beams']['*0']['*0'] 00511 ia.done() 00512 #### 00513 if(len(inim)==1): 00514 ib=ia.fromimage(ourfile=cubeimage, infile=inim[0], overwrite=True) 00515 else: 00516 ib=ia.imageconcat(outfile=cubeimage, infiles=inim, 00517 axis=3, relax=True, overwrite=True) 00518 ia.done() 00519 ##### CAS-4423 temp 00520 if(len(elbeamo) >0): 00521 ib.setrestoringbeam(beam=elbeamo) 00522 ##### 00523 if(csys != None): 00524 ib.setcoordsys(csys=csys.torecord()) 00525 ib.done() 00526 if(removeinfile): 00527 for k in range(len(inim)): 00528 ia.removefile(inim[k]) 00529 00530 @staticmethod 00531 def putchanimage(cubimage,inim,chans, removeinfile=True): 00532 """ 00533 put channels image back to a pre-exisiting cubeimage 00534 """ 00535 #pdb.set_trace() 00536 if(type(chans) != list): 00537 chans=[chans] 00538 if(type(inim) != list): 00539 inim=[inim] 00540 #if( not os.path.exists(inim[0])): 00541 # return False 00542 #ia.open(inim[0]) 00543 #inimshape=ia.shape() 00544 ############ 00545 #imdata=ia.getchunk() 00546 #immask=ia.getchunk(getmask=True) 00547 ############## 00548 ia.done() 00549 ia.open(cubimage) 00550 cubeshape=ia.shape() 00551 #if inimshape[0:3]!=cubeshape[0:3]: 00552 # return False 00553 k=0 00554 ib,=gentools(['ia']) 00555 for chan in chans: 00556 print 'chan', chan 00557 if(os.path.exists(inim[k])): 00558 ib.open(inim[k]) 00559 inimshape=ib.shape() 00560 #############3 00561 imdata=ib.getchunk() 00562 immask=ib.getchunk(getmask=True) 00563 ############ 00564 ib.done() 00565 blc=[0,0,0,chan] 00566 trc=[inimshape[0]-1,inimshape[1]-1,inimshape[2]-1,chan+inimshape[3]-1] 00567 if( not (cubeshape[3] > (chan+inimshape[3]-1))): 00568 return False 00569 00570 ############ 00571 rg0=ia.setboxregion(blc=blc,trc=trc) 00572 ########### 00573 00574 ######## 00575 ia.putregion(pixels=imdata,pixelmask=immask, region=rg0) 00576 ########### 00577 ###ia.insert(infile=inim[k], locate=blc) 00578 if(removeinfile): 00579 ia.removefile(inim[k]) 00580 k+=1 00581 ia.close() 00582 tb.clearlocks() 00583 #casalog.post('putLOCKS: '+ str(inim)+ ' --- ' + str(tb.listlocks())) 00584 00585 00586 return True 00587 @staticmethod 00588 def putchanimage2(cubimage,inim,chans, notprocesschan, removeinfile=True): 00589 """ 00590 put channels image back to a pre-exisiting cubeimage 00591 """ 00592 if(type(chans) != list): 00593 chans=[chans] 00594 if(type(inim) != list): 00595 inin=[inim] 00596 if(type(notprocesschan) != list): 00597 notprocesschan=[notprocesschan] 00598 #if( not os.path.exists(inim[0])): 00599 # return False 00600 ibig, = gentools(['ia']) 00601 #ia.open(inim[0]) 00602 #inimshape=ia.shape() 00603 ############ 00604 #imdata=ia.getchunk() 00605 #immask=ia.getchunk(getmask=True) 00606 ############## 00607 #ia.done() 00608 00609 ####get the shape of an image 00610 for k in range(len(inim)): 00611 if(os.path.exists(inim[k])): 00612 break 00613 ia.open(inim[k]) 00614 subshape=ia.shape() 00615 ia.done() 00616 ############## 00617 00618 ibig.open(cubimage) 00619 cubeshape=ibig.shape() 00620 chanstart=0 00621 blc=np.array([0,0,0,0]) 00622 #if inimshape[0:3]!=cubeshape[0:3]: 00623 # return False 00624 nchTile=1 00625 if((ibig.summary()['tileshape'][3]%subshape[3] ==0) or (subshape[3] > ibig.summary()['tileshape'][3])): 00626 nchTile=max(ibig.summary()['tileshape'][3], subshape[3]) 00627 else: 00628 nchTile=subshape[3] 00629 trc=np.array([int(cubeshape[0]-1),int(cubeshape[1]-1),int(cubeshape[2]-1),nchTile-1]) 00630 arr=ibig.getchunk(blc=blc.tolist(), trc=trc.tolist()) 00631 nchanchunk=cubeshape[3]/nchTile 00632 nreschan=cubeshape[3]%nchTile 00633 k=0 00634 chanchunk=1 00635 for chan in chans: 00636 if(chan >= chanchunk*nchTile): 00637 ibig.putchunk(arr, blc.tolist()) 00638 blc[3]=chanchunk*nchTile 00639 chanchunk+=1 00640 if(chanchunk< nchanchunk): 00641 trc[3]=chanchunk*nchTile-1 00642 else: 00643 trc[3]=nchanchunk*nchTile+nreschan-1 00644 arr=ibig.getchunk(blc=blc.tolist(), trc=trc.tolist()) 00645 #print 'new slice', blc, trc 00646 if(not notprocesschan[k] and os.path.exists(inim[k])): 00647 ia.open(inim[k]) 00648 inimshape=ia.shape() 00649 arrsub=ia.getchunk() 00650 ia.done() 00651 #blc=[0,0,0,chan] 00652 #trc=[inimshape[0]-1,inimshape[1]-1,inimshape[2]-1,chan+inimshape[3]-1] 00653 #if( not (cubeshape[3] > (chan+inimshape[3]-1))): 00654 # return False 00655 relchan=chan-(chanchunk-1)*nchTile 00656 arr[:,:,:,relchan:relchan+inimshape[3]]=arrsub 00657 #print 'chan', chan 00658 arrsub=[] 00659 ############ 00660 #rg0=ia.setboxregion(blc=blc,trc=trc) 00661 ########### 00662 00663 ######## 00664 #ia.putregion(pixels=imdata,pixelmask=immask, region=rg0) 00665 ########### 00666 #ia.insert(infile=inim[k], locate=blc) 00667 if(removeinfile): 00668 shutil.rmtree(inim[k], True) 00669 k+=1 00670 ibig.putchunk(arr, blc.tolist()) 00671 ibig.unlock() 00672 ibig.close() 00673 tb.clearlocks() 00674 #casalog.post('putLOCKS: '+ str(inim)+ ' --- ' + str(tb.listlocks())) 00675 00676 00677 return True 00678 @staticmethod 00679 def makemask(inmask='', outmask='' , imchan=0, chanchunk=1, startchan=-1): 00680 if(startchan < 0): 00681 startchan=imchan*chanchunk 00682 if(not os.path.exists(inmask)): 00683 return False 00684 ia,rg=gentools(['ia', 'rg']) 00685 ia.open(inmask) 00686 shp=ia.shape() 00687 if(shp[3]==1): 00688 shutil.rmtree(outmask, True) 00689 shutil.copytree(inmask, outmask) 00690 ia.done() 00691 return True 00692 if(shp[3] > (startchan+chanchunk-1)): 00693 mybox=rg.box(blc=[0, 0, 0, startchan], trc=[shp[0]-1, shp[1]-1, shp[2]-1, startchan+chanchunk-1]) 00694 ia.subimage(outfile=outmask, region=mybox, overwrite=True) 00695 ia.done() 00696 return True 00697 return False 00698 00699 @staticmethod 00700 def findchanselLSRK(msname='', field='*', spw='*', numpartition=1, beginfreq=0.0, endfreq=1e12, chanwidth=0.0): 00701 im,ms=gentools(['im', 'ms']) 00702 if(field==''): 00703 field='*' 00704 fieldid=0 00705 fieldids=[] 00706 if(type(field)==str): 00707 fieldids=ms.msseltoindex(vis=msname, field=field)['field'] 00708 elif(type(field)==int): 00709 fieldids=[field] 00710 elif(type(field)==list): 00711 fieldids=field 00712 00713 #im.selectvis(vis=msname, field=field, spw=spw) 00714 partwidth=(endfreq-beginfreq)/float(numpartition) 00715 spwsel=[] 00716 nchansel=[] 00717 startsel=[] 00718 for k in range(numpartition): 00719 time0=time.time() 00720 a={} 00721 a['ms_0']={} 00722 a['ms_0']['spw']=np.array([]) 00723 fid=0 00724 while((a['ms_0']['spw'].tolist()==[]) and (fid < len(fieldids))): 00725 a=im.advisechansel(freqstart=(beginfreq+float(k)*partwidth), freqend=(beginfreq+float(k+1)*partwidth), freqstep=chanwidth, freqframe='LSRK', msname=msname, fieldid=fieldids[fid]) 00726 fid+=1 00727 00728 00729 if(a['ms_0']['spw'].tolist() != []): 00730 spwsel.append(a['ms_0']['spw'].tolist()) 00731 nchansel.append(a['ms_0']['nchan'].tolist()) 00732 startsel.append(a['ms_0']['start'].tolist()) 00733 #print k, 'fstart', (beginfreq+float(k)*partwidth), 'fend' , (beginfreq+float(k+1)*partwidth) 00734 print 'TIme taken for lsrk calc', time.time()-time0 00735 im.done() 00736 del im 00737 return spwsel, startsel, nchansel 00738 00739 #putchanimage=staticmethod(putchanimage) 00740 def makecontimage(self, im, novaliddata, imname): 00741 if(novaliddata==True): 00742 ###make blanks 00743 im.make(imname+'.image') 00744 im.make(imname+'.residual') 00745 im.make(imname+'.model') 00746 im.make(imname+'.psf') 00747 return 00748 if(not self.imageparamset): 00749 try: 00750 im.clean(algorithm='mfclark', gain=self.gain, niter=0, threshold='0.05mJy', 00751 model=imname+'.model', image=imname+'.image', 00752 residual=imname+'.residual', psfimage=imname+'.psf') 00753 except Exception, instance: 00754 if(string.count(instance.message, 'PSFZero') >0): 00755 novaliddata=True 00756 ###make a blank image 00757 im.make(imname+'.image') 00758 else: 00759 raise instance 00760 else: ### else of (if not self.imageparamset) 00761 if(not novaliddata): 00762 #casalog.post('Updating '+msname+' imname '+imname) 00763 im.updateresidual(model=imname+'.model', image=imname+'.image', 00764 residual=imname+'.residual') 00765 #im.clean(algorithm='mfclark', niter=0, threshold='0.05mJy', 00766 # model=imname+'.model', image=imname+'.image', 00767 # residual=imname+'.residual') 00768 #casalog.post('CACHE: '+ str(tb.showcache())) 00769 def makemtcontimage(self, im, imname, nterms, scales, reffreq): 00770 incremental=self.imageparamset 00771 models=[] 00772 psfs=[] 00773 residuals=[] 00774 restoreds=[] 00775 sumwts=[] 00776 npsftaylor = 2 * nterms - 1 00777 for tt in range(nterms): 00778 models.append(imname+'.model.tt'+str(tt)); 00779 residuals.append(imname+'.residual.tt'+str(tt)); 00780 restoreds.append(imname+'.image.tt'+str(tt)); 00781 sumwts.append(imname+'.sumwt.tt'+str(tt)); 00782 if(not os.path.exists(models[tt])): 00783 im.make(models[tt]); 00784 for tt in range(0,npsftaylor): 00785 psfs.append(imname+'.psf.tt'+str(tt)); 00786 #this is a dummy i believe as deconvolution is done else where 00787 im.setscales(scalemethod='uservector',uservector=scales); 00788 im.settaylorterms(ntaylorterms=nterms, 00789 reffreq=(qa.convert(qa.unit(reffreq),"Hz"))['value']) 00790 im.setoptions(ftmachine=self.ft) 00791 im.clean(model=models,image=restoreds,psfimage=psfs[0:nterms], residual=residuals, algorithm='msmfs',niter=-1) 00792 if( incremental == False ) : # first major cycle 00793 im.getweightgrid(type='ftweight', wgtimages=sumwts) 00794 ## rename extra psfs.... 00795 ##im.getextrapsfs(multifieldid=0, psfs = psfs[ntaylor:npsftaylor]) 00796 for p in range( nterms, npsftaylor ): 00797 os.system('mv '+imname+'.TempPsf.'+str(p)+' '+psfs[p]) 00798 #casalog.post('CACHE: '+ str(tb.showcache())) 00799
1.8.0
