|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 import os 00002 import shutil 00003 import numpy 00004 from taskinit import * 00005 from cleanhelper import * 00006 im,cb,ms,tb,fg,me,ia,po,sm,cl,cs,rg,sl,dc,vp,msmd,fi,fn=gentools() 00007 00008 def clean(vis, imagename,outlierfile, field, spw, selectdata, timerange, 00009 uvrange, antenna, scan, observation, mode, gridmode, 00010 wprojplanes, facets, cfcache, painc, aterm, psterm, wbawp, epjtable, 00011 interpolation, 00012 niter, gain, threshold, psfmode, imagermode, ftmachine, mosweight, 00013 scaletype, multiscale, negcomponent, smallscalebias, 00014 interactive, mask, nchan, start, width, outframe, 00015 veltype, imsize, cell, phasecenter, restfreq, stokes, weighting, 00016 robust, uvtaper, outertaper, innertaper, modelimage, restoringbeam, 00017 pbcor, minpb, usescratch, noise, npixels, npercycle, cyclefactor, 00018 cyclespeedup, nterms, reffreq, chaniter, flatnoise, allowchunk): 00019 00020 #Python script 00021 casalog.origin('clean') 00022 casalog.post('nchan='+str(nchan)+' start='+str(start)+' width='+str(width)) 00023 00024 #If using new FT-Machines, do not use the on-the-fly model_data columns. 00025 if gridmode == 'advancedaprojection' and usescratch==False: 00026 casalog.post('Forcing usescratch=True for new FTMs. This is temporary.', 'WARN') 00027 usescratch=True 00028 00029 ####################################################################### 00030 # 00031 # start of the big cube treatment 00032 # 00033 ####################################################################### 00034 #paralist=[vis, imagename,outlierfile, field, spw, selectdata, \ 00035 # timerange, uvrange, antenna, scan, observation, mode, gridmode, \ 00036 # wprojplanes, facets, cfcache, painc, epjtable, \ 00037 # interpolation, niter, gain, threshold, psfmode, \ 00038 # imagermode, ftmachine, mosweight, scaletype, multiscale, \ 00039 # negcomponent, smallscalebias, interactive, mask, nchan, \ 00040 # start, width, outframe, veltype, imsize, cell, \ 00041 # phasecenter, restfreq, stokes, weighting, \ 00042 # robust, uvtaper, outertaper, innertaper, modelimage, \ 00043 # restoringbeam, pbcor, minpb, calready, noise, npixels, \ 00044 # npercycle, cyclefactor, cyclespeedup, nterms, reffreq, \ 00045 # chaniter, flatnoise, allowchunk] 00046 #print paralist 00047 00048 if (spw==''): 00049 spw='*'; 00050 if allowchunk and (mode=='velocity' or mode=='frequency' or mode=='channel') \ 00051 and not interactive: 00052 00053 casalog.post('analysing intended channalization...') 00054 imCln=imtool() 00055 imset=cleanhelper(imCln, vis, usescratch, casalog) 00056 00057 (npage, localstart, localwidth)=imset.setChannelizeNonDefault(mode, 00058 spw, field,nchan,start,width,outframe,veltype,phasecenter, restfreq) 00059 casalog.post('start='+str(localstart)+' width='+str(localwidth)+' nchan='+str(npage)) 00060 del imCln 00061 00062 try: 00063 # estimate the size of the image 00064 import commands 00065 00066 nstokes=len(stokes) 00067 casalog.post('imsize='+str(imsize)+' nstokes='+str(nstokes)) 00068 00069 if len(imsize)==1: 00070 npixel=imsize[0]*imsize[0] 00071 else: 00072 npixel=imsize[0]*imsize[1] 00073 00074 volumn=4*nstokes*npage*npixel 00075 00076 00077 av=cu.hostinfo()['memory']['available'] 00078 casalog.post('mem available: '+str(int(av/1024))+'M') 00079 00080 freemem=commands.getoutput("free") 00081 for line in freemem.splitlines(): 00082 if line.startswith('Mem'): 00083 av=float(line.split()[3])/1024 00084 break 00085 casalog.post('mem free: '+str(int(av))+'M') 00086 00087 nd=volumn/1024./1024.*9 00088 casalog.post('mem needed for single chunck clean: '+str(int(nd))+'M') 00089 00090 chunk=1 00091 tchan=npage 00092 fa=2.5 00093 if not av/nd>fa: 00094 tchan=npage*av/nd/fa 00095 chunk=math.ceil(npage/tchan) 00096 tchan=int(math.floor(npage/chunk)) 00097 chunk=int(chunk) 00098 #print 'tchan', tchan, ' chunk', chunk 00099 00100 if chunk>1 and tchan<npage: 00101 casalog.post('will clean the cube in '+str(chunk)+' chunks') 00102 bigimg=imagename 00103 ta=qa.quantity(str(localstart)) 00104 wd=qa.quantity(str(localwidth)) 00105 subimg=[] 00106 for k in range(chunk): 00107 st=0 00108 bg=k*tchan 00109 ed=bg+tchan-1 00110 imname=bigimg+'_'+str(bg)+'-'+str(ed) 00111 if mode=='channel': 00112 st=localstart+k*tchan 00113 ed=st+tchan-1 00114 if mode=='frequency': 00115 st=qa.convert(ta, 'Hz')['value'] 00116 ed=qa.convert(wd, 'Hz')['value'] 00117 st=str(st+k*tchan*ed)+'Hz' 00118 if mode=='velocity': 00119 st=qa.convert(ta, 'm/s')['value'] 00120 ed=qa.convert(wd, 'm/s')['value'] 00121 st=str(st+k*tchan*ed)+'m/s' 00122 00123 #print imname, tchan, st, localwidth 00124 00125 os.system('rm -rf '+imname+'*') 00126 clean(vis=vis,imagename=imname,outlierfile=outlierfile,field=field, 00127 spw=spw,selectdata=selectdata,timerange=timerange,uvrange=uvrange, 00128 antenna=antenna,scan=scan, observation=str(observation), 00129 mode=mode,gridmode=gridmode, 00130 wprojplanes=wprojplanes,facets=facets,cfcache=cfcache,painc=painc, 00131 psterm=psterm,aterm=aterm,wbawp=wbaw, 00132 epjtable=epjtable,interpolation=interpolation,niter=niter, 00133 gain=gain, 00134 threshold=threshold,psfmode=psfmode,imagermode=imagermode, 00135 ftmachine=ftmachine,mosweight=mosweight,scaletype=scaletype, 00136 multiscale=multiscale,negcomponent=negcomponent, 00137 smallscalebias=smallscalebias,interactive=interactive, 00138 mask=mask,nchan=tchan,start=st,width=localwidth,outframe=outframe, 00139 veltype=veltype,imsize=imsize,cell=cell,phasecenter=phasecenter, 00140 restfreq=restfreq,stokes=stokes,weighting=weighting, 00141 robust=robust,uvtaper=uvtaper,outertaper=outertaper, 00142 innertaper=innertaper,modelimage=modelimage, 00143 restoringbeam=restoringbeam,pbcor=pbcor,minpb=minpb, 00144 usescratch=usescratch,noise=noise,npixels=npixels,npercycle=npercycle, 00145 cyclefactor=cyclefactor,cyclespeedup=cyclespeedup,nterms=nterms, 00146 reffreq=reffreq,chaniter=chaniter,flatnoise=flatnoise, 00147 allowchunk=False) 00148 subimg.append(imname) 00149 00150 for i in ['.image', '.flux', '.model', '.psf', '.residual', '.mask']: 00151 casalog.post('concate '+bigimg+'_*'+i+' to '+bigimg+i) 00152 os.system('rm -rf '+bigimg+i) 00153 inf='' 00154 for j in range(len(subimg)): 00155 inf+=' '+subimg[j]+i 00156 bigim=ia.imageconcat(outfile=bigimg+i, infiles=inf, relax=True) 00157 bigim.done() 00158 ia.close() 00159 os.system('rm -rf '+inf) 00160 00161 return 00162 else: 00163 casalog.post('will clean the cube in a single chunk') 00164 #nchan=npage 00165 #start=localstart 00166 #width=localwidth 00167 #let default channelization handles it its own way 00168 except: 00169 raise 00170 00171 else: 00172 casalog.post('Use default channelization for clean') 00173 ####################################################################### 00174 # 00175 # end of big cube treatment 00176 # 00177 ####################################################################### 00178 00179 00180 casalog.post('clean image: '+str(imagename)) 00181 00182 00183 applyoffsets=False; 00184 pbgridcorrect=True; 00185 padding=1.2; 00186 # 00187 # While the following condition is irrelavent due to the change in 00188 # the previous line, I (SB) am leaving it here till after some of 00189 # us discuss this change more carefully. (Feb. 25, 2010). 00190 # 00191 if (facets > 1): 00192 padding=1.2; 00193 00194 # Handle selectdata explicitly 00195 # (avoid use of hidden globals) 00196 if (not selectdata): 00197 timerange='' 00198 uvrange='' 00199 antenna='' 00200 scan='' 00201 observation = '' 00202 00203 try: 00204 # Create a new imager tool 00205 imCln=imtool(); 00206 00207 ###if usescratch open ms with scratch column 00208 ###if mosweight use scratch columns as there in no 00209 ###mosweight available for no scratch column /readonly ms yet 00210 imset=cleanhelper(imCln, vis, usescratch, casalog) 00211 00212 # multims input only (do sorting of vis list based on spw) 00213 if type(vis)==list: imset.sortvislist(spw,mode,width) 00214 00215 # Check imagename 00216 if((len(imagename) == 0) or 00217 ((type(imagename) == str) and imagename.isspace())): 00218 raise Exception, 'Cannot proceed with blank imagename' 00219 00220 multifield=False 00221 if (type(imagename)==list) & (len(imagename) > 1): 00222 multifield=True 00223 elif (type(phasecenter) == list) and (len(phasecenter) >1): 00224 raise TypeError, 'Number of phasecenters has to be equal to number of images' 00225 00226 # Figure out which FTMachine to use. 00227 localFTMachine = getFTMachine(gridmode, imagermode, mode, wprojplanes, 00228 ftmachine); 00229 00230 casalog.post("FTMachine used is %s "%localFTMachine) 00231 00232 00233 # handle mode='mfs' explicitly 00234 if (mode=='mfs'): 00235 start=0; 00236 outframe='' 00237 00238 #some default value handling for channelization 00239 if (mode=='velocity' or mode=='frequency' or mode=='channel'): 00240 # new version: uses ms.cvelfreqs 00241 (localnchan, localstart, localwidth)=imset.setChannelizeDefault(mode,spw,field,nchan,start,width,outframe,veltype,phasecenter, restfreq) 00242 00243 else: 00244 imset.setspecframe(spw) 00245 localnchan=nchan 00246 localstart=start 00247 localwidth=width 00248 00249 #setup for 'per channel' clean 00250 dochaniter=False 00251 #if interactive and chaniter: 00252 if chaniter: 00253 # if veltype=="optical": 00254 # raise Exception, 'The chaniter=True interactive clean for optical velocity mode is not implemented yet.' 00255 if localnchan > 1: 00256 dochaniter=True 00257 00258 # make a template cube for interactive chanter=T 00259 if dochaniter: 00260 imset.makeTemplateCubes(imagename,outlierfile, field, spw, selectdata, 00261 timerange, uvrange, antenna, scan, str(observation), 00262 mode, facets, cfcache, 00263 interpolation, imagermode, localFTMachine, mosweight, 00264 localnchan, localstart, localwidth, outframe, veltype, 00265 imsize, cell, phasecenter, restfreq, stokes, weighting, 00266 robust, uvtaper, outertaper, innertaper, modelimage, 00267 restoringbeam, usescratch, noise, npixels, padding) 00268 00269 nchaniter=localnchan 00270 # check nchan in templatecube 00271 if type(imagename)==list: 00272 imgname=imagename[0]+'.image' 00273 finalimagename=imagename 00274 else: 00275 try: 00276 if int(imset.imageids.values()[0])==0: 00277 # must be using outlier file with old syntax 00278 imgname=imagename+'_'+imset.imageids.values()[0]+'.image' 00279 finalimagename=[imagename+'_'+img for img in imset.imageids.values()] 00280 except: 00281 imgname=imagename+'.image' 00282 finalimagename=imset.imageids.values() 00283 ia.open(imgname) 00284 if localnchan > ia.shape()[3]: 00285 nchaniter = ia.shape()[3] 00286 ia.close() 00287 #finalimagename=imagename 00288 if type(finalimagename)==str: 00289 finalimagename=[finalimagename] 00290 imset.finalimages=finalimagename 00291 # move the following to a helper func. 00292 # create a temporary directory to put channel images 00293 (freqs,finc,newmode,tmppath)=imset.initChaniter(localnchan,spw,localstart,localwidth,finalimagename,mode) 00294 mode=newmode 00295 else: 00296 nchaniter=1 00297 finalimagename='' 00298 tmppath='' 00299 00300 # loop over channels for per-channel clean 00301 for j in xrange(nchaniter): 00302 if dochaniter: 00303 00304 print "Processing channel %s of %s" % (j+1, nchaniter) 00305 casalog.post("Processing channel %s of %s"% (j+1, nchaniter)) 00306 chaniterParms=imset.setChaniterParms(finalimagename,spw,j,localstart,localwidth,freqs,finc,tmppath) 00307 #chaniterParms=imset.setChaniterParms(finalimagename,spw,j,localstart,width,freqs,finc,tmppath) 00308 imagename=chaniterParms['imagename'] 00309 imnchan=chaniterParms['imnchan'] 00310 chanslice=chaniterParms['chanslice'] 00311 localwidth=chaniterParms['width'] 00312 imstart=chaniterParms['imstart'] 00313 visnchan=chaniterParms['visnchan'] 00314 visstart=chaniterParms['visstart'] 00315 00316 # change to handle multifield masks 00317 maskimage='' 00318 if mask in ([], ''): 00319 mask=[''] 00320 if interactive: 00321 if mask==['']: 00322 if type(imagename) == list: 00323 maskimage = [] 00324 for namoo in imagename: 00325 maskimage.append(namoo + '.mask') 00326 imset.maskimages[namoo] = namoo + '.mask' 00327 elif (type(imagename) == str) and (len(outlierfile) 00328 != 0): 00329 maskimage = imagename + '.mask' 00330 imset.maskimages[imagename] = imagename + '.mask' 00331 00332 #read outlierfile 00333 # - initialize local variables to handle multiple fields 00334 imageids=[] 00335 imsizes=[] 00336 phasecenters=[] 00337 parms={} 00338 rootname='' 00339 # 00340 newformat=False 00341 # make a copy of the input, mask, modelimage (will be modified) 00342 if type(modelimage)==str: 00343 modelimage=[modelimage]; 00344 loc_modelimage= modelimage 00345 if type(mask)==str: 00346 mask=[mask]; 00347 loc_mask = mask 00348 00349 # new handling: 00350 # need to combine task parameter inputs and outlier file input 00351 if len(outlierfile) != 0: 00352 #imsizes,phasecenters,imageids=imset.readoutlier(outlierfile) 00353 # if newfomat = False, outlier file is in old format 00354 f_imageids,f_imsizes,f_phasecenters,f_masks,f_modelimages,parms,newformat=imset.newreadoutlier(outlierfile) 00355 #print "from outlierfile: f_imsizes=",f_imsizes," f_phasecenters=",f_phasecenters,\ 00356 # " f_imageids=",f_imageids," f_masks=",f_masks, " parms=",parms 00357 00358 if type(imagename) == list or newformat or dochaniter: 00359 rootname = '' 00360 else: 00361 rootname = imagename 00362 00363 # combine with the task parameter input 00364 if dochaniter: 00365 # imagename is already combined 00366 if type(imagename)!=list: 00367 raise Exception, "imagename=%s expected to be a list." % imagename 00368 else: 00369 imageids=imagename 00370 if newformat: 00371 if type(imsize[0])==list: 00372 imsizes=imsize 00373 else: 00374 imsizes.append(imsize) 00375 if type(phasecenter)==list: 00376 phasecenters=phaseceneter 00377 else: 00378 phasecenters.append(phasecenter) 00379 else: 00380 if type(imagename) == str: 00381 if newformat: 00382 imageids.append(imagename) 00383 imsizes.append(imsize) 00384 phasecenters.append(phasecenter) 00385 else: 00386 imageids=imagename 00387 imsizes=imsize 00388 phasecenters=phasecenter 00389 # for mask, modelimage task input 00390 # turn them into list or list of list 00391 if(loc_mask == []): ## UU 00392 loc_mask = ['']; ## UU 00393 if type(loc_mask) != list: 00394 loc_mask=[loc_mask] 00395 elif type(loc_mask[0]) != list: 00396 loc_mask=[loc_mask] 00397 00398 if(loc_modelimage == []): ## UU 00399 loc_modelimage = ['']; ## UU 00400 if type(loc_modelimage) != list: 00401 loc_modelimage=[loc_modelimage] 00402 elif type(loc_modelimage[0]) != list: ## UUU and type(imagename) != str: 00403 loc_modelimage=[loc_modelimage] 00404 # add extra bracket to correct matching to image 00405 #elif type(loc_modelimage[0]) != list: 00406 if type(loc_modelimage[0]) != list: 00407 # if (type(imagename)==list and len(imagename) < len(loc_modelimage)) or \ 00408 if (type(imagename)==list and len(imagename) != len(loc_modelimage)) or \ 00409 (type(imagename)==str and len(loc_modelimage)>1): 00410 loc_modelimage=[loc_modelimage] 00411 00412 00413 # now append readoutlier content 00414 for indx, name in enumerate(f_imageids): 00415 if not dochaniter: 00416 imageids.append(name) 00417 imsizes.append(f_imsizes[indx]) 00418 phasecenters.append(f_phasecenters[indx]) 00419 if newformat: 00420 loc_mask.append(f_masks[indx]) 00421 #modelimage.append(f_modelimages[indx]) 00422 loc_modelimage.append(f_modelimages[indx]) 00423 else: 00424 # append empty string list to match the size of modelimage list 00425 if indx!=0: 00426 #modelimage.append(f_modelimages[indx]) 00427 loc_modelimage.append(f_modelimages[indx]) 00428 00429 nfield=len(imageids) 00430 if nfield > 1: 00431 multifield=True 00432 #check if number of elements for each list matches 00433 if len(imsizes) != nfield: 00434 raise Exception, "Mismatch in number of imsizes for %s image fields" % nfield 00435 if len(phasecenters) != nfield: 00436 raise Exception, "Mismatch in number of phasecenters for %s image fields" % nfield 00437 # check of mask and modelimage need to be done later... 00438 00439 # check for dulplicated entry 00440 for imname in imageids: 00441 if (imageids.count(imname)!=1): 00442 raise Exception, "Duplicate entry for imagename=%s" % imname 00443 else: 00444 00445 imsizes=[imsize[0], imsize[0]] if ((len(imsize)==1) and numpy.isscalar(imsize[0])) else imsize 00446 phasecenters=phasecenter 00447 imageids=imagename 00448 casalog.post("imsizes="+str(imsizes)+" imageids="+str(imageids), 'DEBUG1') 00449 00450 00451 ###test image sizes 00452 optsize=[0,0] 00453 tmpsize=imsizes if type(imsizes[0])==list else [imsizes] 00454 for ksize in range(len(tmpsize)): 00455 nksize=len(tmpsize[ksize]) 00456 optsize[0]=cleanhelper.getOptimumSize(tmpsize[ksize][0]) 00457 if nksize==1: # imsize can be a single element 00458 optsize[1]=optsize[0] 00459 else: 00460 optsize[1]=cleanhelper.getOptimumSize(tmpsize[ksize][1]) 00461 if((optsize[0] != tmpsize[ksize][0]) or (nksize!=1 and optsize[1] != tmpsize[ksize][1])): 00462 raise ValueError, str(tmpsize[ksize])+' is not an acceptable imagesize, try '+str(optsize) 00463 # 00464 # Moved getAlgorithm() to here so that multifield is set using outlier file. 00465 # 00466 localAlgorithm = getAlgorithm(psfmode, imagermode, gridmode, mode, 00467 multiscale, multifield, facets, nterms, 00468 'clark'); 00469 00470 00471 ###PBCOR or not 00472 sclt='SAULT' 00473 makepbim=False 00474 if (scaletype=='PBCOR') or (scaletype=='pbcor'): 00475 sclt='NONE' 00476 imCln.setvp(dovp=True) 00477 else: 00478 if imagermode != 'mosaic': 00479 makepbim=True 00480 # scaletype is 'SAULT' so use default sclt 00481 # regardless of pbcor is T/F 00482 #elif pbcor: # mosaic and pbcor=true 00483 # sclt='NONE' # do the division in c++ 00484 ###always setvp for mosaic mode 00485 if(imagermode=='mosaic'): 00486 imCln.setvp(dovp=True) 00487 00488 if not dochaniter: 00489 imnchan=localnchan 00490 chanslice=-1 00491 imstart=localstart 00492 visnchan=-1 00493 visstart=0 00494 00495 # data selection only 00496 imset.datsel(field=field, spw=spw, 00497 timerange=timerange, uvrange=uvrange, 00498 antenna=antenna, 00499 scan=scan, observation=str(observation), 00500 usescratch=usescratch, nchan=visnchan, 00501 start=visstart, width=1) 00502 00503 imset.definemultiimages(rootname=rootname, imsizes=imsizes, 00504 cell=cell, stokes=stokes, mode=mode, 00505 spw=spw, nchan=imnchan, start=imstart, 00506 width=localwidth, restfreq=restfreq, 00507 field=field, phasecenters=phasecenters, 00508 names=imageids, facets=facets, 00509 outframe=outframe, veltype=veltype, 00510 makepbim=makepbim, checkpsf=dochaniter) 00511 00512 # weighting and tapering 00513 imset.datweightfilter(field=field, spw=spw, 00514 timerange=timerange, uvrange=uvrange, 00515 antenna=antenna, scan=scan, 00516 wgttype=weighting, robust=robust, 00517 noise=noise, npixels=npixels, 00518 mosweight=mosweight, 00519 uvtaper=uvtaper, 00520 innertaper=innertaper, 00521 outertaper=outertaper, 00522 usescratch=usescratch, nchan=visnchan, 00523 start=visstart, width=1) 00524 00525 # Do data selection and wieghting,papering all at once. 00526 # This does not work for multims input 00527 # imset.datselweightfilter(field=field, spw=spw, 00528 # timerange=timerange, uvrange=uvrange, 00529 # antenna=antenna, scan=scan, 00530 # wgttype=weighting, robust=robust, 00531 # noise=noise, npixels=npixels, 00532 # mosweight=mosweight, 00533 # innertaper=innertaper, 00534 # outertaper=outertaper, 00535 # calready=calready, nchan=visnchan, 00536 # start=visstart, width=1) 00537 # 00538 # make maskimage 00539 if maskimage=='': 00540 maskimage=imset.imagelist[0]+'.mask' 00541 00542 if not multifield: 00543 imset.makemaskimage(outputmask=maskimage, imagename=imagename, 00544 maskobject=mask, slice=chanslice) 00545 00546 else: 00547 #print "mask=",loc_mask, " chanslice=", chanslice 00548 #imset.makemultifieldmask2(mask,chanslice) 00549 imset.makemultifieldmask2(loc_mask,chanslice,newformat, interactive) 00550 maskimage=[] 00551 #for img in sorted(imset.maskimages): 00552 for img in imset.maskimages.keys(): 00553 maskimage.append(imset.maskimages[img]) 00554 casalog.post('Used mask(s) : ' + str(loc_mask) + ' to create mask image(s) : ' + str(maskimage),'INFO'); 00555 00556 if dochaniter: 00557 imset.checkpsf(chanslice) 00558 if imset.skipclean: # for chaniter=T, and if the channel is flagged. 00559 imset.makeEmptyimages() 00560 casalog.post('No valid data, skip CLEANing this channel','WARN') 00561 00562 imCln.setoptions(ftmachine=localFTMachine, 00563 wprojplanes=wprojplanes, 00564 freqinterp=interpolation, padding=padding, 00565 cfcachedirname=cfcache, pastep=painc, 00566 pblimit=minpb, 00567 epjtablename=epjtable, 00568 applypointingoffsets=applyoffsets, 00569 dopbgriddingcorrections=pbgridcorrect, 00570 psterm=psterm,aterm=aterm,wbawp=wbawp); 00571 00572 00573 ##Set the restoring beam 00574 imset.setrestoringbeam(restoringbeam) 00575 00576 #(1) print "Making image names for ", nterms, " terms and " ,len(imset.imagelist) , " fields"; 00577 modelimages=[] 00578 restoredimage=[] 00579 residualimage=[] 00580 psfimage=[] 00581 fluximage=[] 00582 00583 if(nterms==1): 00584 for k in range(len(imset.imagelist)): 00585 modelimages.append(imset.imagelist[k]+'.model') 00586 restoredimage.append(imset.imagelist[k]+'.image') 00587 residualimage.append(imset.imagelist[k]+'.residual') 00588 psfimage.append(imset.imagelist[k]+'.psf') 00589 if(imagermode=='mosaic'): 00590 fluximage.append(imset.imagelist[k]+'.flux') 00591 else: 00592 for tt in range(0, nterms): 00593 for k in range(len(imset.imagelist)): 00594 # make names of all images 00595 modelimages.append( imset.imagelist[k]+'.model.tt'+str(tt) ) 00596 restoredimage.append(imset.imagelist[k]+'.image.tt'+str(tt)) 00597 residualimage.append(imset.imagelist[k]+'.residual.tt'+str(tt)) 00598 psfimage.append(imset.imagelist[k]+'.psf.tt'+str(tt)) 00599 if(imagermode=='mosaic'): 00600 fluximage.append(imset.imagelist[k]+'.flux.tt' + str(tt)) 00601 00602 # (2) print 'Make sure the model images exist on disk' 00603 for tt in range(0, nterms): 00604 for k in range(len(imset.imagelist)): 00605 if nterms==1: 00606 modname = imset.imagelist[k]+'.model'; 00607 else: 00608 modname = imset.imagelist[k]+'.model.tt'+str(tt) ; 00609 if(not os.path.exists( modname ) ): 00610 if( not os.path.exists( imset.imagelist[k] ) ): 00611 raise Exception, "Internal task error. Model image " + imset.imagelist[k] + " does not exist"; 00612 shutil.copytree( imset.imagelist[k] , modname ); 00613 casalog.post("No model found. Making empty initial model : "+modname); 00614 else: 00615 casalog.post("Found and starting from existing model on disk : "+modname); 00616 00617 00618 # (3) print 'Add user-specified model images to the default ones on disk'; 00619 imset.convertAllModelImages(loc_modelimage, mode, nterms, dochaniter, j, tmppath); 00620 00621 # (4) print "Delete the template images made per field."; 00622 for k in range(len(imset.imagelist)): 00623 ia.open( imset.imagelist[k] ); 00624 ia.remove(verbose=False); 00625 ia.close(); 00626 00627 # Multi-term parameters. 00628 if (mode == "mfs") and (nterms > 1): 00629 if multiscale == []: 00630 multiscale = [0] 00631 00632 # Check that reference-frequency is a valid string 00633 reffreqVal=0.0; 00634 qat=qatool(); 00635 try: 00636 rff=qat.canonical(reffreq); 00637 except Exception, instance: 00638 print '*** Error *** In conversion of reffreq=\'',reffreq,'\' to a numerical value'; 00639 raise Exception, instance 00640 reffreqVal=rff['value']; # This is the frequency in Hz 00641 00642 # Set the number of terms and reference-frequency 00643 imCln.settaylorterms(ntaylorterms=nterms, 00644 reffreq=reffreqVal); 00645 00646 # forbid pbcorrection with msmfs for now 00647 if(pbcor): 00648 raise Exception, 'Primary-beam correction is currently not supported with nterms>1' 00649 00650 casalog.post('Running MS-MFS with '+str(nterms)+' Taylor-terms on dataset : ' + str(vis)); 00651 ########################################################### 00652 00653 if len(multiscale) > 0: 00654 imCln.setscales(scalemethod='uservector', 00655 uservector=multiscale) 00656 imCln.setsmallscalebias(smallscalebias) 00657 00658 if imagermode == 'mosaic': 00659 if localFTMachine == 'ft': 00660 imCln.setmfcontrol(stoplargenegatives=negcomponent, 00661 scaletype=sclt,minpb=minpb,constpb=1.0, 00662 cyclefactor=cyclefactor, 00663 cyclespeedup=cyclespeedup, 00664 # fluxscale=[imagename+'.flux']) 00665 fluxscale=fluximage) 00666 else: 00667 imCln.setmfcontrol(stoplargenegatives=negcomponent, 00668 scaletype=sclt,minpb=minpb, 00669 cyclefactor=cyclefactor, 00670 cyclespeedup=cyclespeedup, 00671 # fluxscale=[imagename+'.flux']) 00672 fluxscale=fluximage, flatnoise=flatnoise) 00673 else: 00674 imCln.setmfcontrol(stoplargenegatives=negcomponent, 00675 cyclefactor=cyclefactor, 00676 cyclespeedup=cyclespeedup, minpb=minpb) 00677 00678 ####after all the mask shenanigans...make sure to use the 00679 ####last mask 00680 if not multifield: 00681 maskimage = imset.outputmask 00682 ##if no effective mask was passed down but interactive 00683 if maskimage == '' and interactive: 00684 if imset.maskimages.has_key(imset.imagelist[0]): 00685 maskimage=imset.maskimages[imset.imagelist[0]] 00686 else: 00687 maskimage = imset.imagelist[0] + '.mask' 00688 imset.maskimages[imset.imagelist[0]] = maskimage 00689 maskimg = mask 00690 if mask == True: 00691 maskimg = minpb 00692 if ((len(maskimage) == 0 or maskimage[0] == '') and 00693 isinstance(maskimg, float) and maskimg > 0.0 and maskimg < 1.0 00694 # and imagermode == 'mosaic' 00695 and not interactive): 00696 casalog.post('Making a mask at primary beam level ' + str(maskimg), 00697 'INFO') 00698 casalog.post('Running clean with niter=0 to get the primary beam coverage', 00699 'INFO') 00700 # Run clean with niter = 0 to get the pbcoverage. 00701 imCln.clean(algorithm=localAlgorithm, niter=0, gain=gain, 00702 threshold=qa.quantity(threshold,'mJy'), 00703 model=modelimages, residual=residualimage, 00704 image=restoredimage, psfimage=psfimage, 00705 mask=maskimage, interactive=False, 00706 npercycle=npercycle) 00707 pbcov_image = imset.imagelist[0] + '.flux' 00708 if localFTMachine == 'mosaic': 00709 pbcov_image += '.pbcoverage' 00710 maskimage = imset.make_mask_from_threshhold(pbcov_image, maskimg) 00711 if not imset.skipclean: 00712 #print "imager.clean() starts" 00713 imCln.clean(algorithm=localAlgorithm, niter=niter, gain=gain, 00714 threshold=qa.quantity(threshold,'mJy'), 00715 model=modelimages, residual=residualimage, 00716 image=restoredimage, psfimage=psfimage, 00717 mask=maskimage, interactive=interactive, 00718 npercycle=npercycle) 00719 00720 #In the interactive mode, deconvlution can be skipped and in that case 00721 #psf is not generated. So check if all psfs are there if not, generate 00722 # if interactive: 00723 # if nterms==1: 00724 # for psfim in psfimage: 00725 # if not os.path.isdir(psfim): 00726 # imCln.approximatepsf(psf=psfim) 00727 # else: 00728 # casalog.post('Multi-term PSFs not made','WARN') 00729 00730 if interactive: 00731 for psfim in psfimage: 00732 if not os.path.isdir(psfim): 00733 if nterms==1: 00734 imCln.approximatepsf(psf=psfim) 00735 else: 00736 casalog.post('Multi-term PSF '+psfim+' not made','WARN') 00737 00738 if dochaniter and not imset.skipclean : 00739 imset.storeCubeImages(finalimagename,imset.imagelist,j,imagermode) 00740 00741 imCln.close() 00742 #################################################################### 00743 if dochaniter: 00744 imset.cleanupTempFiles(tmppath) 00745 imset.imagelist=finalimagename 00746 presdir=os.path.realpath('.') 00747 for k in xrange(len(imset.imagelist)): 00748 newimage=imset.imagelist[k] 00749 if(imset.imagelist[k].count('/') > 0): 00750 newimage=os.path.basename(imset.imagelist[k]) 00751 os.chdir(os.path.dirname(imset.imagelist[k])) 00752 result = '\'' + newimage + '.image' + '\''; 00753 fluxscale_image = '\'' + newimage + '.flux' + '\''; 00754 pbcov_image=fluxscale_image 00755 # convert .flux image, which may not be correct outframe if 00756 # it was made from im tool, to correct outframe 00757 imset.convertImageFreqFrame([fluxscale_image.replace('\'','')]) 00758 if(localFTMachine=='mosaic'): 00759 pbcov_image = '\'' + newimage + '.flux.pbcoverage' + '\''; 00760 residim = '\'' + newimage + '.residual' + '\''; 00761 if (pbcor): 00762 if(sclt != 'NONE'): 00763 ##otherwise its already divided 00764 ia.open(newimage+'.image') 00765 00766 ##this is not needed as mask is set in C++ 00767 #pixmask = fluxscale_image+'>'+str(minpb); 00768 #ia.calcmask(pixmask,asdefault=True); 00769 corrfac=minpb*minpb 00770 if(localFTMachine=='ft'): 00771 corrfac=minpb 00772 pixels='iif('+ pbcov_image+'>'+str(corrfac)+','+ result+'/'+fluxscale_image+', 0)' 00773 00774 ia.calc(pixels=pixels) 00775 ia.close() 00776 ia.open(newimage+'.residual') 00777 pixels='iif('+ fluxscale_image+'>'+str(corrfac)+','+ residim+'/'+fluxscale_image+', 0)' 00778 ia.calc(pixels=pixels) 00779 ia.close() 00780 elif sclt == 'NONE': 00781 ia.open(newimage+'.image') 00782 pixels=result+'*'+fluxscale_image 00783 ia.calc(pixels=pixels) 00784 ia.close() 00785 ia.open(newimage+'.residual') 00786 pixels=residim+'*'+fluxscale_image 00787 ia.calc(pixels=pixels) 00788 ia.close() 00789 00790 os.chdir(presdir) 00791 00792 del imCln 00793 00794 except Exception, instance: 00795 print '*** Error *** ',instance 00796 raise Exception, instance 00797
1.8.0
