|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from taskinit import * 00002 from simutil import * 00003 import os 00004 import re 00005 import pylab as pl 00006 import pdb 00007 00008 def simanalyze( 00009 project=None, 00010 image=None, 00011 imagename=None, 00012 vis=None, modelimage=None, cell=None, imsize=None, imdirection=None, 00013 niter=None, threshold=None, 00014 weighting=None, mask=None, outertaper=None, stokes=None, 00015 analyze=None, 00016 showuv=None, showpsf=None, showmodel=None, 00017 showconvolved=None, showclean=None, showresidual=None, showdifference=None, 00018 showfidelity=None, 00019 graphics=None, 00020 verbose=None, 00021 overwrite=None, 00022 async=False): 00023 00024 # Collect a list of parameter values to save inputs 00025 in_params = locals() 00026 00027 import re 00028 import glob 00029 00030 casalog.origin('simanalyze') 00031 if verbose: casalog.filter(level="DEBUG2") 00032 00033 a = inspect.stack() 00034 stacklevel = 0 00035 for k in range(len(a)): 00036 if (string.find(a[k][1], 'ipython console') > 0): 00037 stacklevel = k 00038 myf = sys._getframe(stacklevel).f_globals 00039 00040 # create the utility object: 00041 util = simutil() 00042 if verbose: util.verbose = True 00043 msg = util.msg 00044 00045 # put output in directory called "project" 00046 fileroot = project 00047 if not os.path.exists(fileroot): 00048 msg(fileroot+" directory doesn't exist - the task expects to find results from creating the datasets there, like the skymodel.",priority="error") 00049 return False 00050 00051 00052 saveinputs = myf['saveinputs'] 00053 saveinputs('simanalyze',fileroot+"/"+project+".simanalyze.last", 00054 myparams=in_params) 00055 00056 if (not image) and (not analyze): 00057 casalog.post("No operation to be done. Exiting from task.", "WARN") 00058 return True 00059 00060 grscreen = False 00061 grfile = False 00062 if graphics == "both": 00063 grscreen = True 00064 grfile = True 00065 if graphics == "screen": 00066 grscreen = True 00067 if graphics == "file": 00068 grfile = True 00069 00070 try: 00071 # if True: 00072 00073 # handle '$project' in modelimage 00074 modelimage = modelimage.replace('$project',project) 00075 00076 # things we need: model_cell, model_direction if user doesn't specify - 00077 # so find those first, and get information using util.modifymodel 00078 # with skymodel=newmodel 00079 00080 components_only=False 00081 00082 # first look for skymodel, if not then compskymodel 00083 skymodels=glob.glob(fileroot+"/"+project+"*.skymodel")+glob.glob(fileroot+"/"+project+"*.newmodel") 00084 nmodels=len(skymodels) 00085 if nmodels>1: 00086 msg("Found %i sky model images:" % nmodels) 00087 for ff in skymodels: 00088 msg(" "+ff) 00089 msg("Using "+skymodels[0]) 00090 if nmodels>=1: 00091 skymodel=skymodels[0] 00092 else: 00093 skymodel="" 00094 00095 if os.path.exists(skymodel): 00096 msg("Sky model image "+skymodel+" found.") 00097 else: 00098 skymodels=glob.glob(fileroot+"/"+project+"*.compskymodel") 00099 nmodels=len(skymodels) 00100 if nmodels>1: 00101 msg("Found %i sky model images:" % nmodels) 00102 for ff in skymodels: 00103 msg(" "+ff) 00104 msg("Using "+skymodels[0]) 00105 if nmodels>=1: 00106 skymodel=skymodels[0] 00107 else: 00108 skymodel="" 00109 00110 if os.path.exists(skymodel): 00111 msg("Sky model image "+skymodel+" found.") 00112 components_only=True 00113 else: 00114 msg("Can't find a model image in your project directory, named skymodel or compskymodel - output image will be been created, but comparison with the input model is not possible.",priority="error") 00115 return False 00116 00117 modelflat = skymodel+".flat" 00118 if not os.path.exists(modelflat): 00119 util.flatimage(skymodel,verbose=verbose) 00120 00121 # modifymodel just collects info if skymodel==newmodel 00122 returnpars = util.modifymodel(skymodel,skymodel, 00123 "","","","","",-1, 00124 flatimage=False) 00125 if not returnpars: 00126 return False 00127 00128 (model_refdir,model_cell,model_size, 00129 model_nchan,model_center,model_width, 00130 model_stokes) = returnpars 00131 00132 00133 cell_asec=qa.convert(model_cell[0],'arcsec')['value'] 00134 00135 00136 ##################################################################### 00137 # clean if desired, use noisy image for further calculation if present 00138 # todo suggest a cell size from psf? 00139 ##################################################################### 00140 if (not image): 00141 user_imagename=imagename 00142 if user_imagename=="default" or len(user_imagename)<=0: 00143 images= glob.glob(fileroot+"/*image") 00144 if len(images)<1: 00145 msg("can't find any image in project directory",priority="error") 00146 return False 00147 if len(images)>1: 00148 msg("found multiple images in project directory") 00149 msg("using "+images[0]) 00150 imagename=images[0] 00151 # trim .image suffix: 00152 imagename= imagename.replace(".image","") 00153 00154 00155 beam_current = False 00156 if image: 00157 00158 # make sure cell is defined 00159 if type(cell) == type([]): 00160 if len(cell) > 0: 00161 cell0 = cell[0] 00162 else: 00163 cell0 = "" 00164 else: 00165 cell0 = cell 00166 if len(cell0) <= 0: 00167 cell = model_cell 00168 if type(cell) == type([]): 00169 if len(cell) == 1: 00170 cell = [cell[0],cell[0]] 00171 else: 00172 cell = [cell,cell] 00173 00174 # cells are positive by convention 00175 cell = [qa.abs(cell[0]),qa.abs(cell[1])] 00176 00177 # and imsize 00178 if type(imsize) == type([]): 00179 if len(imsize) > 0: 00180 imsize0 = imsize[0] 00181 if len(imsize) > 1: 00182 imsize1 = imsize[1] 00183 else: 00184 imsize1 = imsize0 00185 else: 00186 imsize0 = -1 00187 else: 00188 imsize0 = imsize 00189 if imsize0 <= 0: 00190 imsize = [int(pl.ceil(qa.convert(qa.div(model_size[0],cell[0]),"")['value'])), 00191 int(pl.ceil(qa.convert(qa.div(model_size[1],cell[1]),"")['value']))] 00192 else: 00193 imsize=[imsize0,imsize1] 00194 00195 00196 # check for default measurement sets: 00197 default_mslist = glob.glob(fileroot+"/*ms") 00198 n_default=len(default_mslist) 00199 # is the user requesting this ms? 00200 default_requested=[] 00201 for i in range(n_default): default_requested.append(False) 00202 00203 00204 # parse ms parameter and check for existance; 00205 # initial ms list 00206 if vis=="default" or len(vis)==0: 00207 mslist0=default_mslist 00208 else: 00209 mslist0 = vis.split(',') 00210 # verified found ms list 00211 mslist = [] 00212 mstype = [] 00213 00214 mstoimage=[] 00215 tpmstoimage=None 00216 00217 for ms0 in mslist0: 00218 if not len(ms0): continue 00219 ms1 = ms0.replace('$project',project) 00220 00221 # MSes in fileroot/ have priority 00222 if os.path.exists(fileroot+"/"+ms1): 00223 ms1 = fileroot + "/" + ms1 00224 00225 if os.path.exists(ms1): 00226 mslist.append(ms1) 00227 00228 # mark as requested 00229 if default_mslist.count(ms1): 00230 i=default_mslist.index(ms1) 00231 default_requested[i]=True 00232 00233 # check if noisy in name 00234 if re.search('noisy.',ms1): 00235 ms1_raw=str.join("",re.split('noisy.',ms1)) 00236 if default_mslist.count(ms1_raw): 00237 i=default_mslist.index(ms1_raw) 00238 default_requested[i]=True 00239 else: # not noisy 00240 if ms1.endswith(".sd.ms"): 00241 ms1_noisy=re.split('.sd.ms',ms1)[0]+'.noisy.sd.ms' 00242 else: 00243 ms1_noisy=re.split('.ms',ms1)[0]+'.noisy.ms' 00244 if default_mslist.count(ms1_noisy): 00245 i=default_mslist.index(ms1_noisy) 00246 default_requested[i]=True 00247 if vis == "default": continue 00248 msg("You requested "+ms1+" but there is a corrupted (noisy) version of the ms in your project directory - if your intent is to model noisy data you may want to check inputs",priority="warn") 00249 00250 # check if the ms is tp data or not. 00251 if util.ismstp(ms1,halt=False): 00252 mstype.append('TP') 00253 tpmstoimage = ms1 00254 # XXX TODO more than one TP ms will not be handled 00255 # correctly 00256 msg("Found a total power measurement set, %s." % ms1) 00257 else: 00258 mstype.append('INT') 00259 mstoimage.append(ms1) 00260 msg("Found a synthesis measurement set, %s." % ms1) 00261 else: 00262 if verbose: 00263 msg("measurement set "+ms1+" not found -- removing from imaging list",priority="warn") 00264 00265 else: 00266 msg("measurement set "+ms1+" not found -- removing from imaging list") 00267 00268 # check default mslist for unrequested ms: 00269 if verbose: 00270 priority="warn" 00271 else: 00272 priority="info" 00273 for i in range(n_default): 00274 if not default_requested[i]: 00275 msg("Project directory contains "+default_mslist[i]+" but you have not requested to include it in your simulated image.",priority=priority) 00276 00277 00278 00279 # more than one to image? 00280 if len(mstoimage) > 1: 00281 msg("Multiple interferometric ms found:",priority="warn") 00282 for i in range(len(mstoimage)): 00283 msg(" "+mstoimage[i],priority="warn") 00284 msg(" will be concated and simultaneously deconvolved; if something else is desired, please specify vis, or image manually and use image=F",priority="warn") 00285 concatms=project+"/"+project+".concat.ms" 00286 from concat import concat 00287 concat(mstoimage,concatms) 00288 mstoimage=[concatms] 00289 00290 if len(mstoimage) == 0: 00291 if tpmstoimage: 00292 sd_only = True 00293 else: 00294 msg("no measurement sets found to image",priority="warn") 00295 return False 00296 else: 00297 sd_only = False 00298 # get some quantities from the interferometric ms 00299 maxbase=0. 00300 # TODO make work better for multiple MS 00301 for msfile in mstoimage: 00302 tb.open(msfile) 00303 rawdata = tb.getcol("UVW") 00304 tb.done() 00305 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m 00306 psfsize = 0.3/qa.convert(qa.quantity(model_center),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec 00307 minimsize = 8* int(psfsize/cell_asec) 00308 00309 if imsize[0] < minimsize: 00310 msg("The number of image pixel in x-axis, %d, is small to cover 8 x PSF. Setting x pixel number, %d." % (imsize[0], minimsize), priority='warn') 00311 imsize[0] = minimsize 00312 if imsize[1] < minimsize: 00313 msg("The number of image pixel in y-axis, %d, is small to cover 8 x PSF. Setting y pixel number, %d" % (imsize[1], minimsize), priority='warn') 00314 imsize[1] = minimsize 00315 00316 00317 # Do single dish imaging first if tpmstoimage exists. 00318 if tpmstoimage and os.path.exists(tpmstoimage): 00319 msg('creating image from ms: '+tpmstoimage) 00320 if len(mstoimage): 00321 tpimage = project + '.sd.image' 00322 else: 00323 tpimage = project + '.image' 00324 tpimage = fileroot + "/" + tpimage 00325 00326 if len(mstoimage): 00327 if len(modelimage) and tpimage != modelimage and \ 00328 tpimage != fileroot+"/"+modelimage: 00329 msg("modelimage parameter set to "+modelimage+" but also creating a new total power image "+tpimage,priority="warn") 00330 msg("assuming you know what you want, and using modelimage="+modelimage+" in deconvolution",priority="warn") 00331 else: 00332 # This forces to use TP image as a model for clean 00333 if len(modelimage) <= 0: 00334 msg("you are generating total power image "+tpimage+". this is used as a model image for clean",priority="warn") 00335 modelimage = tpimage 00336 00337 # format image size properly 00338 sdimsize = imsize 00339 if not isinstance(imsize,list): 00340 sdimsize = [imsize,imsize] 00341 elif len(imsize) == 1: 00342 sdimsize = [imsize[0],imsize[0]] 00343 00344 im.open(tpmstoimage) 00345 im.selectvis(nchan=model_nchan,start=0,step=1,spw=0) 00346 im.defineimage(mode='channel',nx=sdimsize[0],ny=sdimsize[1],cellx=cell[0],celly=cell[1],phasecenter=model_refdir,nchan=model_nchan,start=0,step=1,spw=0) 00347 #im.setoptions(ftmachine='sd',gridfunction='pb') 00348 im.setoptions(ftmachine='sd',gridfunction='pb') 00349 im.makeimage(type='singledish',image=tpimage) 00350 im.close() 00351 del sdimsize 00352 00353 # For single dish: manually set the primary beam 00354 ia.open(tpimage) 00355 beam = ia.restoringbeam() 00356 if len(beam) == 0: 00357 msg('setting primary beam information to image.') 00358 # !! aveant will only be set if modifymodel or setpointings and in 00359 # any case it will the the aveant of the INTERFM array - we want the SD 00360 tb.open(tpmstoimage+"/ANTENNA") 00361 diams = tb.getcol("DISH_DIAMETER") 00362 tb.done() 00363 aveant = pl.mean(diams) 00364 # model_center should be set even if we didn't predict this execution 00365 pb = 1.2*0.3/qa.convert(qa.quantity(model_center),'GHz')['value']/aveant*3600.*180/pl.pi 00366 beam['major'] = beam['minor'] = qa.quantity(pb,'arcsec') 00367 beam['positionangle'] = qa.quantity(0.0,'deg') 00368 msg('Primary beam: '+str(beam['major'])) 00369 ia.setrestoringbeam(beam=beam) 00370 ia.close() 00371 if sd_only: 00372 beam_current = True 00373 bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2 00374 bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix 00375 else: del beam 00376 #del beam 00377 00378 msg('generation of total power image '+tpimage+' complete.') 00379 # update TP ms name the for following steps 00380 sdmsfile = tpmstoimage 00381 sd_any = True 00382 00383 imagename = re.split('.image$',tpimage)[0] 00384 # End of single dish imaging part 00385 00386 outflat_current = False 00387 convsky_current = False 00388 00389 if image and len(mstoimage) > 0: 00390 00391 # for reruns 00392 foo=mstoimage[0] 00393 foo=foo.replace(".ms","") 00394 foo=foo.replace(project,"") 00395 foo=foo.replace("/","") 00396 project=project+foo 00397 00398 imagename = fileroot + "/" + project 00399 00400 # get nfld, sourcefieldlist, from (interfm) ms if it was not just created 00401 # TODO make work better for multiple mstoimage (figures below) 00402 tb.open(mstoimage[0]+"/SOURCE") 00403 code = tb.getcol("CODE") 00404 sourcefieldlist = pl.where(code=='OBJ')[0] 00405 nfld = len(sourcefieldlist) 00406 tb.done() 00407 msfile = mstoimage[0] 00408 00409 # set cleanmode automatically (for interfm) 00410 if nfld == 1: 00411 cleanmode = "csclean" 00412 else: 00413 cleanmode = "mosaic" 00414 00415 00416 00417 00418 # clean insists on using an existing model if its present 00419 if os.path.exists(imagename+".image"): shutil.rmtree(imagename+".image") 00420 if os.path.exists(imagename+".model"): shutil.rmtree(imagename+".model") 00421 00422 # An image in fileroot/ has priority 00423 if len(modelimage) > 0 and os.path.exists(fileroot+"/"+modelimage): 00424 modelimage = fileroot + "/" + modelimage 00425 msg("Found modelimage, %s." % modelimage) 00426 00427 # in simdata we use imdirection instead of model_refdir 00428 if not util.isdirection(imdirection,halt=False): 00429 imdirection=model_refdir 00430 util.imclean(mstoimage,imagename, 00431 cleanmode,cell,imsize,imdirection, 00432 niter,threshold,weighting, 00433 outertaper,stokes, #sourcefieldlist=sourcefieldlist, 00434 modelimage=modelimage,mask=mask) 00435 00436 # create imagename.flat and imagename.residual.flat: 00437 util.flatimage(imagename+".image",verbose=verbose) 00438 util.flatimage(imagename+".residual",verbose=verbose) 00439 outflat_current = True 00440 00441 msg("done inverting and cleaning") 00442 if not type(cell) == type([]): 00443 cell = [cell,cell] 00444 if len(cell) <= 1: 00445 cell = [qa.quantity(cell[0]),qa.quantity(cell[0])] 00446 else: 00447 cell = [qa.quantity(cell[0]),qa.quantity(cell[1])] 00448 cell = [qa.abs(cell[0]),qa.abs(cell[0])] 00449 00450 # get beam from output clean image 00451 if verbose: msg("getting beam from "+imagename+".image",origin="analysis") 00452 ia.open(imagename+".image") 00453 beam = ia.restoringbeam() 00454 beam_current = True 00455 ia.close() 00456 # model has units of Jy/pix - calculate beam area from clean image 00457 # (even if we are not plotting graphics) 00458 bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2 00459 bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix 00460 msg("synthesized beam area in output pixels = %f" % bmarea) 00461 00462 00463 00464 00465 00466 00467 00468 00469 00470 00471 if image: 00472 # show model, convolved model, clean image, and residual 00473 if grfile: 00474 file = fileroot + "/" + project + ".image.png" 00475 else: 00476 file = "" 00477 else: 00478 mslist=[] 00479 00480 if image and len(mstoimage) > 0: 00481 if grscreen or grfile: 00482 util.newfig(multi=[2,2,1],show=grscreen) 00483 00484 # create regridded and convolved sky model image 00485 util.convimage(modelflat,imagename+".image.flat") 00486 convsky_current = True # don't remake this for analysis in this run 00487 00488 disprange = [] # passing empty list causes return of disprange 00489 00490 # original sky regridded to output pixels but not convolved with beam 00491 discard = util.statim(modelflat+".regrid",disprange=disprange,showstats=False) 00492 util.nextfig() 00493 00494 # convolved sky model - units of Jy/bm 00495 disprange = [] 00496 discard = util.statim(modelflat+".regrid.conv",disprange=disprange) 00497 util.nextfig() 00498 00499 # clean image - also in Jy/beam 00500 # although because of DC offset, better to reset disprange 00501 disprange = [] 00502 discard = util.statim(imagename+".image.flat",disprange=disprange) 00503 util.nextfig() 00504 00505 if len(mstoimage) > 0: 00506 util.nextfig() 00507 00508 # clean residual image - Jy/bm 00509 discard = util.statim(imagename+".residual.flat",disprange=disprange) 00510 util.endfig(show=grscreen,filename=file) 00511 00512 00513 00514 00515 ##################################################################### 00516 # analysis 00517 00518 if analyze: 00519 00520 # set above 00521 # if not image: 00522 # imagename = fileroot + "/" + project 00523 00524 if not os.path.exists(imagename+".image"): 00525 msg("Can't find a simulated image - expecting "+imagename,priority="error") 00526 return False 00527 00528 # we should have skymodel.flat created above 00529 00530 if not image: 00531 if not os.path.exists(imagename+".image"): 00532 msg("you must image before analyzing.",priority="error") 00533 return False 00534 00535 # get beam from output clean image 00536 if verbose: msg("getting beam from "+imagename+".image",origin="analysis") 00537 ia.open(imagename+".image") 00538 beam = ia.restoringbeam() 00539 beam_current = True 00540 ia.close() 00541 # model has units of Jy/pix - calculate beam area from clean image 00542 cell = util.cellsize(imagename+".image") 00543 cell= [ qa.convert(cell[0],'arcsec'), 00544 qa.convert(cell[1],'arcsec') ] 00545 # (even if we are not plotting graphics) 00546 bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2 00547 bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix 00548 msg("synthesized beam area in output pixels = %f" % bmarea) 00549 00550 00551 # flat output:? if the user manually cleaned, this may not exist 00552 outflat = imagename + ".image.flat" 00553 if (not outflat_current) or (not os.path.exists(outflat)): 00554 # create imagename.flat and imagename.residual.flat 00555 util.flatimage(imagename+".image",verbose=verbose) 00556 if os.path.exists(imagename+".residual"): 00557 util.flatimage(imagename+".residual",verbose=verbose) 00558 else: 00559 if showresidual: 00560 msg(imagename+".residual not found -- residual will not be plotted",priority="warn") 00561 showresidual = False 00562 outflat_current = True 00563 00564 # regridded and convolved input:? 00565 if not convsky_current: 00566 util.convimage(modelflat,imagename+".image.flat") 00567 convsky_current = True 00568 00569 # now should have all the flat, convolved etc even if didn't run "image" 00570 00571 # make difference image. 00572 # immath does Jy/bm if image but only if ia.setbrightnessunit("Jy/beam") in convimage() 00573 convolved = modelflat + ".regrid.conv" 00574 difference = imagename + '.diff' 00575 diff_ia = ia.imagecalc(difference, "'%s' - '%s'" % (convolved, outflat), overwrite=True) 00576 diff_ia.setbrightnessunit("Jy/beam") 00577 00578 # get rms of difference image for fidelity calculation 00579 #ia.open(difference) 00580 diffstats = diff_ia.statistics(robust=True, verbose=False,list=False) 00581 diff_ia.close() 00582 del diff_ia 00583 maxdiff = diffstats['medabsdevmed'] 00584 if maxdiff != maxdiff: maxdiff = 0. 00585 if type(maxdiff) != type(0.): 00586 if maxdiff.__len__() > 0: 00587 maxdiff = maxdiff[0] 00588 else: 00589 maxdiff = 0. 00590 # Make fidelity image. 00591 absdiff = imagename + '.absdiff' 00592 calc_ia = ia.imagecalc(absdiff, "max(abs('%s'), %f)" % (difference, 00593 maxdiff/pl.sqrt(2.0)), overwrite=True) 00594 calc_ia.close() 00595 fidelityim = imagename + '.fidelity' 00596 calc_ia = ia.imagecalc(fidelityim, "abs('%s') / '%s'" % (convolved, absdiff), overwrite=True) 00597 calc_ia.close() 00598 msg("fidelity image calculated",origin="analysis") 00599 00600 # scalar fidelity 00601 absconv = imagename + '.absconv' 00602 calc_ia = ia.imagecalc(absconv, "abs('%s')" % convolved, overwrite=True) 00603 if ia.isopen(): ia.close() #probably not necessary 00604 calc_ia.close() 00605 del calc_ia 00606 00607 ia.open(absconv) 00608 modelstats = ia.statistics(robust=True, verbose=False,list=False) 00609 maxmodel = modelstats['max'] 00610 if maxmodel != maxmodel: maxmodel = 0. 00611 if type(maxmodel) != type(0.): 00612 if maxmodel.__len__() > 0: 00613 maxmodel = maxmodel[0] 00614 else: 00615 maxmodel = 0. 00616 ia.close() 00617 scalarfidel = maxmodel/maxdiff 00618 msg("fidelity range (max model / rms difference) = "+str(scalarfidel),origin="analysis") 00619 00620 00621 # now, what does the user want to actually display? 00622 00623 # need MS for showuv and showpsf 00624 if not image: 00625 msfile = fileroot + "/" + project + ".ms" 00626 elif sd_only: 00627 # imaged and single dish only 00628 msfile = tpmstoimage 00629 # if sd_only and os.path.exists(sdmsfile): 00630 # # use TP ms for UV plot if only SD sim, i.e., 00631 # # image=sd_only=T or (image=F=predict_uv and predict_sd=T) 00632 # msfile = sdmsfile 00633 # psf is not available for SD only sim 00634 if os.path.exists(msfile) and util.ismstp(msfile,halt=False): 00635 if showpsf: msg("single dish simulation -- psf will not be plotted",priority='warn') 00636 showpsf = False 00637 if (not image) and (not os.path.exists(msfile)): 00638 if showpsf or showuv: 00639 msg("No image is generated in this run. Default MS, '%s', does not exists -- uv and psf will not be plotted" % msfile,priority='warn') 00640 showpsf = False 00641 showuv = False 00642 00643 00644 # if the order in the task input changes, change it here too 00645 figs = [showuv,showpsf,showmodel,showconvolved,showclean,showresidual,showdifference,showfidelity] 00646 nfig = figs.count(True) 00647 if nfig > 6: 00648 msg("only displaying first 6 selected panels in graphic output",priority="warn") 00649 if nfig <= 0: 00650 return True 00651 if nfig < 4: 00652 multi = [1,nfig,1] 00653 else: 00654 if nfig == 4: 00655 multi = [2,2,1] 00656 else: 00657 multi = [2,3,1] 00658 00659 if grfile: 00660 file = fileroot + "/" + project + ".analysis.png" 00661 else: 00662 file = "" 00663 00664 if grscreen or grfile: 00665 util.newfig(multi=multi,show=grscreen) 00666 00667 # if order in task parameters changes, change here too 00668 if showuv: 00669 # TODO loop over all ms - show all UV including zero 00670 if len(mslist)>1: 00671 msg("Using only "+msfile+" for uv plot",priority="warn") 00672 tb.open(msfile) 00673 rawdata = tb.getcol("UVW") 00674 tb.done() 00675 pl.box() 00676 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m 00677 klam_m = 300/qa.convert(model_center,'GHz')['value'] 00678 pl.plot(rawdata[0,]/klam_m,rawdata[1,]/klam_m,'b,') 00679 pl.plot(-rawdata[0,]/klam_m,-rawdata[1,]/klam_m,'b,') 00680 ax = pl.gca() 00681 ax.yaxis.LABELPAD = -4 00682 pl.xlabel('u[klambda]',fontsize='x-small') 00683 pl.ylabel('v[klambda]',fontsize='x-small') 00684 pl.axis('equal') 00685 # Add zero-spacing (single dish) if not yet plotted 00686 # TODO make this a check over all ms 00687 # if predict_sd and not util.ismstp(msfile,halt=False): 00688 # pl.plot([0.],[0.],'r,') 00689 util.nextfig() 00690 00691 if showpsf: 00692 if image: 00693 psfim = imagename + ".psf" 00694 else: 00695 psfim = project + ".quick.psf" 00696 if not os.path.exists(psfim): 00697 if len(mslist)>1: 00698 msg("Using only "+msfile+" for psf generation",priority="warn") 00699 im.open(msfile) 00700 # TODO spectral parms 00701 im.defineimage(cellx=qa.tos(model_cell[0]),nx=max([minimsize,128])) 00702 if os.path.exists(psfim): 00703 shutil.rmtree(psfim) 00704 im.approximatepsf(psf=psfim) 00705 # beam is set above (even in "analyze" only) 00706 # note that if image, beam has fields 'major' whereas if not, it 00707 # has fields like 'bmaj'. 00708 # beam=im.fitpsf(psf=psfim) 00709 im.done() 00710 ia.open(psfim) 00711 beamcs = ia.coordsys() 00712 beam_array = ia.getchunk(axes=[beamcs.findcoordinate("spectral")[1][0],beamcs.findcoordinate("stokes")[1][0]],dropdeg=True) 00713 nn = beam_array.shape 00714 xextent = nn[0]*cell_asec*0.5 00715 xextent = [xextent,-xextent] 00716 yextent = nn[1]*cell_asec*0.5 00717 yextent = [-yextent,yextent] 00718 flipped_array = beam_array.transpose() 00719 ttrans_array = flipped_array.tolist() 00720 ttrans_array.reverse() 00721 pl.imshow(ttrans_array,interpolation='bilinear',cmap=pl.cm.jet,extent=xextent+yextent,origin="bottom") 00722 psfim.replace(project+"/","") 00723 pl.title(psfim,fontsize="x-small") 00724 b = qa.convert(beam['major'],'arcsec')['value'] 00725 pl.xlim([-3*b,3*b]) 00726 pl.ylim([-3*b,3*b]) 00727 ax = pl.gca() 00728 pl.text(0.05,0.95,"bmaj=%7.1e\nbmin=%7.1e" % (beam['major']['value'],beam['minor']['value']),transform = ax.transAxes,bbox=dict(facecolor='white', alpha=0.7),size="x-small",verticalalignment="top") 00729 ia.close() 00730 util.nextfig() 00731 00732 disprange = [] # first plot will define range 00733 if showmodel: 00734 discard = util.statim(modelflat+".regrid",incell=cell,disprange=disprange,showstats=False) 00735 util.nextfig() 00736 disprange = [] 00737 00738 if showconvolved: 00739 discard = util.statim(modelflat+".regrid.conv") 00740 # if disprange gets set here, it'll be Jy/bm 00741 util.nextfig() 00742 00743 if showclean: 00744 # own scaling because of DC/zero spacing offset 00745 discard = util.statim(imagename+".image.flat") 00746 util.nextfig() 00747 00748 if showresidual: 00749 # it gets its own scaling 00750 discard = util.statim(imagename+".residual.flat") 00751 util.nextfig() 00752 00753 if showdifference: 00754 # it gets its own scaling. 00755 discard = util.statim(imagename+".diff") 00756 util.nextfig() 00757 00758 if showfidelity: 00759 # it gets its own scaling. 00760 discard = util.statim(imagename+".fidelity",showstats=False) 00761 util.nextfig() 00762 00763 util.endfig(show=grscreen,filename=file) 00764 00765 sim_min,sim_max,sim_rms,sim_units = util.statim(imagename+".image.flat",plot=False) 00766 # if not displaying still print stats: 00767 # 20100505 ia.stats changed to return Jy/bm: 00768 msg('Simulation rms: '+str(sim_rms/bmarea)+" Jy/pix = "+ 00769 str(sim_rms)+" Jy/bm",origin="analysis") 00770 msg('Simulation max: '+str(sim_max/bmarea)+" Jy/pix = "+ 00771 str(sim_max)+" Jy/bm",origin="analysis") 00772 #msg('Simulation rms: '+str(sim_rms)+" Jy/pix = "+ 00773 # str(sim_rms*bmarea)+" Jy/bm",origin="analysis") 00774 #msg('Simulation max: '+str(sim_max)+" Jy/pix = "+ 00775 # str(sim_max*bmarea)+" Jy/bm",origin="analysis") 00776 msg('Beam bmaj: '+str(beam['major']['value'])+' bmin: '+str(beam['minor']['value'])+' bpa: '+str(beam['positionangle']['value']),origin="analysis") 00777 00778 00779 00780 # cleanup - delete newmodel, newmodel.flat etc 00781 if os.path.exists(imagename+".image.flat"): 00782 shutil.rmtree(imagename+".image.flat") 00783 if os.path.exists(imagename+".residual.flat"): 00784 shutil.rmtree(imagename+".residual.flat") 00785 if os.path.exists(imagename+".flux"): 00786 shutil.rmtree(imagename+".flux") 00787 absdiff = imagename + '.absdiff' 00788 if os.path.exists(absdiff): 00789 shutil.rmtree(absdiff) 00790 absconv = imagename + '.absconv' 00791 if os.path.exists(absconv): 00792 shutil.rmtree(absconv) 00793 # if os.path.exists(imagename+".diff"): 00794 # shutil.rmtree(imagename+".diff") 00795 if os.path.exists(imagename+".quick.psf") and os.path.exists(imagename+".psf"): 00796 shutil.rmtree(imagename+".quick.psf") 00797 00798 00799 except TypeError, e: 00800 finalize_tools() 00801 #msg("simanalyze -- TypeError: %s" % e,priority="error") 00802 casalog.post("simanalyze -- TypeError: %s" % e, priority="ERROR") 00803 raise TypeError, e 00804 return 00805 except ValueError, e: 00806 finalize_tools() 00807 #print "simanalyze -- OptionError: ", e 00808 casalog.post("simanalyze -- OptionError: %s" % e, priority="ERROR") 00809 raise ValueError, e 00810 return 00811 except Exception, instance: 00812 finalize_tools() 00813 #print '***Error***',instance 00814 casalog.post("simanalyze -- Exception: %s" % instance, priority="ERROR") 00815 raise Exception, instance 00816 return 00817 00818 00819 def finalize_tools(): 00820 if ia.isopen(): ia.close() 00821 im.close() 00822 tb.close()
1.8.0
