|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os 00002 import shutil 00003 import traceback 00004 import pdb 00005 import numpy as np 00006 import sys 00007 from cleanhelper import * 00008 from taskinit import * 00009 im,cb,ms,tb,fg,me,ia,po,sm,cl,cs,rg,sl,dc,vp,msmd,fi,fn=gentools() 00010 00011 def csvclean(vis, imagename,field, spw, advise, mode, nchan, width, imsize, cell, phasecenter, niter, weighting, restoringbeam, interactive): 00012 00013 """ This task does an invert of the visibilities and deconvolve in the 00014 image plane. It does not do a uvdata subtraction (aka Cotton-Schwab 00015 major cycle) of model visibility as in clean. - For ALMA Commissioning 00016 00017 vis -- Name of input visibility file 00018 default: none; example: vis='ngc5921.ms' 00019 00020 imagename -- Name of output CASA image. (only the prefix) 00021 default: none; example: imagename='m2' 00022 output images are: 00023 m2.image; cleaned and restored image 00024 With or without primary beam correction 00025 m2dirty.image; dirty image 00026 m2psf.image; point-spread function (dirty beam) 00027 m2.model; image of clean components 00028 m2.mask; image containing clean regions, when interative=True 00029 00030 00031 field -- Select fields in MS. Use field id(s) or field name(s). 00032 ['go listobs' to obtain the list id's or names] 00033 default: ''= all fields 00034 If field string is a non-negative integer, it is assumed to 00035 be a field index otherwise, it is assumed to be a 00036 field name 00037 field='0~2'; field ids 0,1,2 00038 field='0,4,5~7'; field ids 0,4,5,6,7 00039 field='3C286,3C295'; field named 3C286 and 3C295 00040 field = '3,4C*'; field id 3, all names starting with 4C 00041 00042 spw --Select spectral window/channels 00043 NOTE: This selects the data passed as the INPUT to mode 00044 default: ''=all spectral windows and channels 00045 spw='0~2,4'; spectral windows 0,1,2,4 (all channels) 00046 spw='0:5~61'; spw 0, channels 5 to 61 00047 spw='<2'; spectral windows less than 2 (i.e. 0,1) 00048 spw='0,10,3:3~45'; spw 0,10 all channels, spw 3, 00049 channels 3 to 45. 00050 spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each. 00051 spw='0:0~10;15~60'; spectral window 0 with channels 00052 0-10,15-60 00053 spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10, 00054 spw 1, channels 20-30, and spw 2, channels, 1,2 and 3 00055 00056 imsize -- Image pixel size (x,y). DOES NOT HAVE TO BE A POWER OF 2 00057 default = [256,256]; example: imsize=[350,350] 00058 imsize = 500 is equivalent to [500,500] 00059 Avoid odd-numbered imsize. 00060 00061 cell -- Cell size (x,y) 00062 default= '1.0arcsec'; 00063 example: cell=['0.5arcsec,'0.5arcsec'] or 00064 cell=['1arcmin', '1arcmin'] 00065 cell = '1arcsec' is equivalent to ['1arcsec','1arcsec'] 00066 NOTE:cell = 2.0 => ['2arcsec', '2arcsec'] 00067 00068 phasecenter -- direction measure or fieldid for the mosaic center 00069 default: '' => first field selected ; example: phasecenter=6 00070 or phasecenter='J2000 19h30m00 -40d00m00' 00071 00072 niter -- Maximum number iterations, 00073 if niter=0, then no CLEANing is done ("invert" only) 00074 default: 500; example: niter=5000 00075 00076 weighting -- Weighting to apply to visibilities: 00077 default='natural'; example: weighting='uniform'; 00078 Options: 'natural','uniform','briggs', 00079 'superuniform','briggsabs','radial' 00080 00081 restoringbeam -- Output Gaussian restoring beam for CLEAN image 00082 [bmaj, bmin, bpa] elliptical Gaussian restoring beam 00083 default units are in arc-seconds for bmaj,bmin, degrees 00084 for bpa default: restoringbeam=[]; Use PSF calculated 00085 from dirty beam. 00086 example: restoringbeam=['10arcsec'] or restorinbeam='10arcsec', circular Gaussian. 00087 FWHM 10 arcseconds example: 00088 restoringbeam=['10.0','5.0','45.0deg'] 10"x5" 00089 at 45 degrees 00090 00091 interactive -- Create a mask interactively or not. 00092 default=False; example: interactive=True 00093 The viewer will open with the image displayed. Select the 00094 region for the mask and double click in the middle of it. 00095 00096 """ 00097 00098 00099 #Python script 00100 00101 try: 00102 00103 casalog.origin('csvclean') 00104 ms = casac.ms() 00105 00106 00107 parsummary = 'vis="'+str(vis)+'", imagename="'+str(imagename)+'", ' 00108 parsummary += 'field="'+str(field)+'", spw="'+str(spw)+'", ' 00109 parsummary = 'cell="'+str(cell)+'",' 00110 parsummary = 'phasecenter='+str(phasecenter)+',' 00111 parsummary += 'imsize='+str(imsize)+', niter='+str(niter)+', ' 00112 parsummary += 'weighting="'+str(weighting)+'", ' 00113 parsummary += 'restoringbeam="'+str(restoringbeam)+'", ' 00114 parsummary += 'interactive='+str(interactive)+'' 00115 casalog.post(parsummary,'INFO') 00116 00117 # if (not (type(vis)==str) & (os.path.exists(vis))): 00118 # raise Exception, 'Visibility data set not found - please verify the name' 00119 if ((type(vis)==str) & (os.path.exists(vis))): 00120 ms.open(vis) 00121 else: 00122 raise Exception, 'Visibility data set not found - please verify the name' 00123 if(not advise): 00124 if (imagename == ""): 00125 # ms.close() 00126 raise Exception, "Must provide output image name in parameter imagename." 00127 00128 if os.path.exists(imagename): 00129 # ms.close() 00130 raise Exception, "Output image %s already exists - will not overwrite." % imagename 00131 00132 if (field == ''): 00133 field = '*' 00134 00135 if (spw == ''): 00136 spw = '*' 00137 00138 if ((type(imsize)==int)): 00139 imsize=[imsize,imsize] 00140 00141 if ((len(imsize)==1)): 00142 imsize=[imsize[0],imsize[0]] 00143 00144 nx = imsize[0] 00145 ny = imsize[1] 00146 00147 if ((type(cell)==int) | (type(cell)==float) | (type(cell)==str)): 00148 cell=[cell,cell] 00149 00150 if ((len(cell)==1)): 00151 cell=[cell[0],cell[0]] 00152 00153 cellx=cell[0] 00154 celly=cell[1] 00155 if((type(cell[0])==int) or (type(cell[0])==float)): 00156 cellx=qa.quantity(cell[0], 'arcsec') 00157 celly=qa.quantity(cell[1], 'arcsec') 00158 00159 if(type(phasecenter)==str): 00160 ### blank means take field[0] 00161 if (phasecenter==''): 00162 fieldoo=field 00163 if(fieldoo==''): 00164 fieldoo='0' 00165 phasecenter=int(ms.msseltoindex(vis,field=fieldoo)['field'][0]) 00166 else: 00167 tmppc=phasecenter 00168 try: 00169 if(len(ms.msseltoindex(vis, field=phasecenter)['field']) > 0): 00170 tmppc = int(ms.msseltoindex(vis, 00171 field=phasecenter)['field'][0]) 00172 ##succesful must be string like '0' or 'NGC*' 00173 except Exception, instance: 00174 #failed must be a string type J2000 18h00m00 10d00m00 00175 tmppc = phasecenter 00176 phasecenter = tmppc 00177 00178 00179 if restoringbeam == [''] or len(restoringbeam) == 0: 00180 # calculate from fit below 00181 bmaj = '' 00182 bmin = '' 00183 bpa = '' 00184 else: 00185 if (type(restoringbeam)==str): 00186 restoringbeam=[restoringbeam,restoringbeam,'0deg'] 00187 if (type(restoringbeam)==list and (len(restoringbeam)==1)): 00188 restoringbeam=[restoringbeam[0],restoringbeam[0],'0deg'] 00189 if (type(restoringbeam)==list and (len(restoringbeam)==2)): 00190 restoringbeam=[restoringbeam[0],restoringbeam[1],'0deg'] 00191 if (type(restoringbeam)==list and (len(restoringbeam)==2)): 00192 restoringbeam=[restoringbeam[0],restoringbeam[1],restoringbeam[2]] 00193 00194 if(qa.quantity(restoringbeam[0])['unit'] == ''): 00195 restoringbeam[0]=restoringbeam[0]+'arcsec' 00196 if(qa.quantity(restoringbeam[1])['unit'] == ''): 00197 restoringbeam[1]=restoringbeam[1]+'arcsec' 00198 if(qa.quantity(restoringbeam[2])['unit'] == ''): 00199 restoringbeam[2]=restoringbeam[2]+'deg' 00200 00201 bmaj = restoringbeam[0] 00202 bmin = restoringbeam[1] 00203 bpa = restoringbeam[2] 00204 00205 # Create output names based on imagename parameter 00206 dirtyim = imagename+'dirty.image' 00207 psfim = imagename+'psf.image' 00208 modelname = imagename+'.model' 00209 imname = imagename+'.image' 00210 00211 # Make sure all tables and images are closed 00212 # ms.close() 00213 00214 # Add scratch columns if they don't exist 00215 #tb.open(vis) 00216 #hasit = tb.colnames().count('CORRECTED_DATA')>0 00217 #tb.close() 00218 #if not hasit: 00219 # cb.open(vis) 00220 # cb.close() 00221 00222 # make the dirty image and psf 00223 00224 im.open(vis, usescratch=True) 00225 im.selectvis(spw=spw, field=field) 00226 spwsel=ms.msseltoindex(vis=vis, spw=spw)['spw'] 00227 ch=ms.msseltoindex(vis=vis, spw=spw)['channel'] 00228 if(nchan < 1): 00229 nchan=0 00230 for k in range(len(spwsel)): 00231 nchan += ch[k,2]-ch[k,1]+1 00232 nchan=nchan/width 00233 if(nchan < 1): 00234 nchan=1 00235 if(advise): 00236 tb.open(vis+'/SPECTRAL_WINDOW') 00237 allreffreq=tb.getcol('REF_FREQUENCY') 00238 reffreq=0.0 00239 if(len(allreffreq) > 1): 00240 reffreq=0.0; 00241 for f in allreffreq: 00242 reffreq+=f 00243 reffreq=reffreq/float(len(allreffreq)) 00244 else: 00245 reffreq=allreffreq[0] 00246 tb.done() 00247 tb.open(vis+'/ANTENNA') 00248 diams=tb.getcol('DISH_DIAMETER') 00249 diam=np.min(diams) 00250 tb.done() 00251 fov=qa.quantity(3.0e8/reffreq/diam, 'rad') 00252 adv=im.advise(fieldofview=fov) 00253 cellx=qa.tos(adv[2], prec=4) 00254 celly=cellx 00255 myf = sys._getframe(len(inspect.stack())-1).f_globals 00256 myf['cell']=[cellx,cellx] 00257 myf['imsize']=[adv[1], adv[1]] 00258 nx=ny=adv[1] 00259 myf['advise']=False 00260 return 00261 redopsf=True 00262 redokounter=0 00263 immode='mfs' 00264 start=ch[0,1] 00265 if(mode=='cube'): 00266 immode='channel' 00267 if(width >1): 00268 start=start+ width/2 00269 while(redopsf): 00270 im.defineimage(nx=nx, ny=ny, cellx=cellx, celly=celly, phasecenter=phasecenter, spw=spwsel.tolist(), mode=immode, start=start, step=width, nchan=nchan) 00271 im.weight(weighting) 00272 try: 00273 im.makeimage(type='corrected', image=dirtyim) 00274 im.makeimage(type='psf', image=psfim) 00275 ###make an empty model 00276 im.make(modelname) 00277 if((redokounter==2) and (np.min(nx,ny) > 25)): 00278 #pdb.set_trace() 00279 ia.open(psfim) 00280 csys=ia.coordsys() 00281 rg.setcoordinates(csys=csys.torecord()) 00282 shp=ia.shape() 00283 blc=['10pix', '10pix', '0pix', '0pix'] 00284 trc=[str(shp[0]-10)+'pix',str(shp[1]-10)+'pix', 00285 str(shp[2]-1)+'pix', str(shp[3]-1)+'pix'] 00286 reg=rg.wbox(blc=blc, trc=trc) 00287 ia.set(pixels=0, region=rg.complement(reg)) 00288 ia.done() 00289 #im.done() 00290 00291 # Calculate bpa, bmin, bmaj if not given 00292 if restoringbeam == [''] or len(restoringbeam) == 0: 00293 cx = nx/2 00294 cy = ny/2 00295 box = '' 00296 00297 #if (nx > 100 and ny > 100): 00298 # rrg = [cx-10, cy-10, cx+10, cy+10] 00299 # box = '%s,%s,%s,%s'%(rrg[0],rrg[1],rrg[2],rrg[3]) 00300 #ia.open(psfim) 00301 #shp=ia.shape() 00302 #coords = ia.fitcomponents(box=box) 00303 #ia.close() 00304 #if(coords['converged'] == True): 00305 coords=im.fitpsf(psfim) 00306 if(coords[0]): 00307 bmaj=coords[1] 00308 bmin=coords[2] 00309 bpa=coords[3] 00310 redopsf=False 00311 else: 00312 redopsf=True 00313 nx=nx+1 00314 ny=ny+1 00315 redokounter += 1 00316 if(redokounter==3): 00317 casalog.post('Failed to find a decent psf','SEVERE') 00318 return False 00319 else: 00320 casalog.post('Trying new image with 1 extra pixel','WARN') 00321 00322 except : 00323 redopsf=True 00324 nx=nx+1 00325 ny=ny+1 00326 redokounter += 1 00327 if(redokounter==3): 00328 casalog.post('Failed to find a decent psf','SEVERE') 00329 im.done() 00330 return False 00331 else: 00332 casalog.post('Trying new image with 1 extra pixel','WARN') 00333 im.done() 00334 parsummary = 'restoringbeam values = [\'%s\',\'%s\',\'%s\']'%(qa.tos(bmaj),qa.tos(bmin),qa.tos(bpa)) 00335 casalog.post(parsummary,'INFO') 00336 00337 # Make a mask 00338 maskname='' 00339 if(interactive): 00340 maskname=imagename+'.mask' 00341 if(os.path.exists(maskname)): 00342 ia.open(dirtyim) 00343 csys=ia.coordsys() 00344 shp=ia.shape() 00345 ia.done() 00346 ia.open(maskname) 00347 ia.regrid(outfile='__tmpmask__', shape=shp, csys=csys.torecord(), axes=[0,1]) 00348 ia.remove(True) 00349 ia.done() 00350 shutil.move('__tmpmask__', maskname) 00351 im.drawmask(dirtyim, maskname) 00352 00353 # use deconvolver to do image plane deconvolution 00354 # using a mask image as the mask 00355 dc.open(dirty=dirtyim, psf=psfim) 00356 # NOTE: use the parameter mask which can be an empty 00357 # string if no mask 00358 dc.clean(niter=niter, model=modelname, mask=maskname) 00359 00360 # create the restored image 00361 if restoringbeam == [''] or len(restoringbeam) == 0: 00362 dc.restore(model=modelname, image=imname, bmaj=bmaj, bmin=bmin, bpa=bpa) 00363 else: 00364 dc.restore(model=modelname, image=imname, bmaj=restoringbeam[0], bmin=restoringbeam[1], bpa=restoringbeam[2]) 00365 00366 dc.done() 00367 return True 00368 00369 except Exception, instance: 00370 im.close() 00371 ia.close() 00372 dc.close() 00373 print '*** Error *** ',instance 00374 casalog.post("Error ...", 'SEVERE') 00375 traceback.print_exc() 00376 raise Exception, instance 00377 00378 00379
1.8.0
