Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_simanalyze.py: 55%

1import os 

2import re 

3import glob 

4import pylab as pl 

5 

6from casatools import table, image, imager 

7from casatasks import casalog, sdimaging, imregrid, immath, concat, feather 

8from . import sdbeamutil 

9from .simutil import * 

10 

11tb = table( ) 

12ia = image( ) 

13im = imager( ) 

14 

15 # the global tb, ia, and im tools are used 

16 

17def simanalyze( 

18 project=None, 

19 image=None, 

20 # if image==False: 

21 imagename=None, 

22 skymodel=None, 

23 # else: 

24 vis=None, modelimage=None, imsize=None, imdirection=None, cell=None, 

25 interactive=None, niter=None, threshold=None, 

26 weighting=None, mask=None, outertaper=None, pbcor=None, stokes=None, 

27 featherimage=None, 

28 # endif 

29 analyze=None, 

30 showuv=None, showpsf=None, showmodel=None, 

31 showconvolved=None, showclean=None, showresidual=None, showdifference=None, 

32 showfidelity=None, 

33 graphics=None, 

34 verbose=None, 

35 overwrite=None, 

36 dryrun=False, 

37 logfile=None 

38 ): 

39 

40#def simanalyze(project='sim', image=True, imagename='default', skymodel='', vis='default', modelimage='', imsize=[128, 128], imdirection='', cell='', interactive=False, niter=0, threshold='0.1mJy', weighting='natural', mask=[], outertaper=[''], stokes='I', featherimage='', analyze=False, showuv=True, showpsf=True, showmodel=True, showconvolved=False, showclean=True, showresidual=False, showdifference=True, showfidelity=True, graphics='both', verbose=False, overwrite=True, dryrun=False): 

41 

42 

43 # Collect a list of parameter values to save inputs 

44 in_params = locals() 

45 

46 casalog.origin('simanalyze') 

47 if verbose: casalog.filter(level="DEBUG2") 

48 

49 # create the utility object:  

50 myutil = simutil() 

51 if logfile: 

52 myutil.reportfile=logfile 

53 myutil.openreport() 

54 if verbose: myutil.verbose = True 

55 msg = myutil.msg 

56 

57 # put output in directory called "project" 

58 fileroot = project 

59 if not os.path.exists(fileroot): 

60 msg(fileroot+" directory doesn't exist - the task expects to find results from creating the datasets there, like the skymodel.",priority="error") 

61 # msg should raise an exception for priority=error 

62 

63 

64 casalog.post("saveinputs not available in casatasks, skipping saving simanalyze inputs", priority='INFO') 

65 

66 if (not image) and (not analyze): 

67 casalog.post("No operation to be done. Exiting from task.", "WARN") 

68 return 

69 

70 grscreen = False 

71 grfile = False 

72 if graphics == "both": 

73 grscreen = True 

74 grfile = True 

75 if graphics == "screen": 

76 grscreen = True 

77 if graphics == "file": 

78 grfile = True 

79 

80 try: 

81 

82 # Predefined parameters  

83 pbcoeff = 1.13 ## PB defined as pbcoeff*lambda/d 

84 

85 

86 # handle '$project' in modelimage 

87 modelimage = modelimage.replace('$project',project) 

88 featherimage = featherimage.replace('$project',project) 

89 

90 #========================= 

91 # things we need: model_cell, model_direction if user doesn't specify -  

92 # so find those first, and get information using util.modifymodel 

93 # with skymodel=newmodel 

94 

95 

96 # we need to parse either the mslist or imagename (if image=False)  

97 # first, so that we can pick the appropriate skymodel,  

98 # if there are several. 

99 skymodel_searchstring="NOT SPECIFIED" 

100 

101 if not (image or dryrun): 

102 user_imagename=imagename 

103 if user_imagename=="default" or len(user_imagename)<=0: 

104 images= glob.glob(fileroot+"/*image") 

105 if len(images)<1: 

106 emsg = "can't find any image in project directory" 

107 raise RuntimeError(emsg) 

108 if len(images)>1: 

109 msg("found multiple images in project directory",priority="warn") 

110 msg("using "+images[0],priority="warn") 

111 imagename=images[0] 

112 # trim .image suffix: 

113 imagename= imagename.replace(".image","") 

114 

115 # if the user hasn't specified a sky model image, we can try to  

116 # see if their imagename contains something like the project and  

117 # configuration, as it would if simobserve created it. 

118 user_skymodel=skymodel 

119 if not os.path.exists(user_skymodel): 

120 if os.path.exists(fileroot+"/"+user_skymodel): 

121 user_skymodel=fileroot+"/"+user_skymodel 

122 elif len(user_skymodel)>0: 

123 raise RuntimeError("Can't find your specified skymodel "+user_skymodel) 

124 # try to strip a searchable identifier 

125 tmpstring=user_skymodel.split("/")[-1] 

126 skymodel_searchstring=tmpstring.replace(".image","") 

127 

128 

129 if image: 

130 # check for default measurement sets: 

131 default_mslist = glob.glob(fileroot+"/*ms") 

132 n_default=len(default_mslist) 

133 # is the user requesting this ms? 

134 default_requested=[] 

135 for i in range(n_default): default_requested.append(False) 

136 

137 # parse ms parameter and check for existance; 

138 # initial ms list 

139 if vis=="default" or len(vis)==0: 

140 mslist0=default_mslist 

141 else: 

142 mslist0 = vis.split(',') 

143 # verified found ms list 

144 mslist = [] 

145 mstype = [] 

146 

147 mstoimage=[] 

148 tpmstoimage=None 

149 

150 for ms0 in mslist0: 

151 if not len(ms0): continue 

152 ms1 = ms0.replace('$project',project) 

153 

154 # MSes in fileroot/ have priority 

155 if os.path.exists(fileroot+"/"+ms1): 

156 ms1 = fileroot + "/" + ms1 

157 

158 if os.path.exists(ms1)or dryrun: 

159 mslist.append(ms1) 

160 

161 # mark as requested 

162 if default_mslist.count(ms1): 

163 i=default_mslist.index(ms1) 

164 default_requested[i]=True 

165 

166 # check if noisy in name 

167 if re.search('noisy.',ms1): 

168 ms1_raw=str.join("",re.split('noisy.',ms1)) 

169 if default_mslist.count(ms1_raw): 

170 i=default_mslist.index(ms1_raw) 

171 default_requested[i]=True 

172 else: # not noisy 

173 if ms1.endswith(".sd.ms"): 

174 ms1_noisy=re.split('.sd.ms',ms1)[0]+'.noisy.sd.ms' 

175 else: 

176 ms1_noisy=re.split('.ms',ms1)[0]+'.noisy.ms' 

177 if default_mslist.count(ms1_noisy): 

178 i=default_mslist.index(ms1_noisy) 

179 default_requested[i]=True 

180 if vis == "default": continue 

181 msg("You requested "+ms1+" but there is a 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") 

182 

183 # check if the ms is tp data or not. 

184 if dryrun: 

185 # HACK 

186 mstype.append('INT') 

187 mstoimage.append(ms1) 

188 

189 elif myutil.ismstp(ms1,halt=False): 

190 mstype.append('TP') 

191 tpmstoimage = ms1 

192 # XXX TODO more than one TP ms will not be handled 

193 # correctly 

194 msg("Found a total power measurement set, %s." % ms1,origin='simanalyze') 

195 else: 

196 mstype.append('INT') 

197 mstoimage.append(ms1) 

198 msg("Found a synthesis measurement set, %s." % ms1,origin='simanalyze') 

199 else: 

200 msg("measurement set "+ms1+" not found -- removing from imaging list") 

201 

202 # check default mslist for unrequested ms: 

203 for i in range(n_default): 

204 if not default_requested[i]: 

205 msg("Project directory contains "+default_mslist[i]+" but you have not requested to include it in your simulated image.") 

206 

207 

208 if not mstoimage and len(tpmstoimage) == 0: 

209 raise RuntimeError("No MS found to image") 

210 

211 # now try to parse the mslist for an identifier string that  

212 # we can use to find the right skymodel if there are several 

213 if len(mstoimage) == 0 and len(tpmstoimage) > 0: 

214 tmpstring = tpmstoimage.split("/")[-1] 

215 else: 

216 tmpstring=(mstoimage[0]).split("/")[-1] 

217 skymodel_searchstring=tmpstring.replace(".ms","") 

218 

219 

220 

221 # more than one to image? 

222 if len(mstoimage) > 1: 

223 msg("Multiple interferometric ms found:",priority="info",origin='simanalyze') 

224 for i in range(len(mstoimage)): 

225 msg(" "+mstoimage[i],priority="info",origin='simanalyze') 

226 msg(" will be concated and simultaneously deconvolved; if something else is desired, please specify vis, or image manually and use image=F",priority="info",origin='simanalyze') 

227 concatms=project+"/"+project+".concat.ms" 

228 weights = get_concatweights(mstoimage) 

229 msg(" concat("+str(mstoimage)+",concatvis='"+concatms+"',visweightscale="+str(weights)+")",origin='simanalyze') 

230 if not dryrun: 

231 concat(mstoimage,concatvis=concatms,visweightscale=weights) 

232 mstoimage=[concatms] 

233 

234 

235 

236 #======================================================== 

237 # now we can search for skymodel, and if there are several,  

238 # pick the one that is closest to either the imagename,  

239 # or the first MS if there are several MS to image. 

240 

241 components_only=False 

242 

243 # first look for skymodel, if not then compskymodel 

244 skymodels=glob.glob(fileroot+"/"+project+"*.skymodel")+glob.glob(fileroot+"/"+project+"*.newmodel") 

245 nmodels=len(skymodels) 

246 skymodel_index=0 

247 if nmodels>1: 

248 msg("Found %i sky model images:" % nmodels,origin='simanalyze') 

249 # use the skymodel_searchstring to try to pick the right one 

250 # print them out for the user while we're at it. 

251 for i in range(nmodels): 

252 msg(" "+skymodels[i]) 

253 if skymodels[i].count(skymodel_searchstring)>0: 

254 skymodel_index=i 

255 msg("Using skymodel "+skymodels[skymodel_index],origin='simanalyze') 

256 if nmodels>=1: 

257 skymodel=skymodels[skymodel_index] 

258 else: 

259 skymodel="" 

260 

261 if os.path.exists(skymodel): 

262 msg("Sky model image "+skymodel+" found.",origin='simanalyze') 

263 else: 

264 skymodels=glob.glob(fileroot+"/"+project+"*.compskymodel") 

265 nmodels=len(skymodels) 

266 if nmodels>1: 

267 msg("Found %i sky model images:" % nmodels,origin='simanalyze') 

268 for ff in skymodels: 

269 msg(" "+ff) 

270 msg("Using "+skymodels[0],origin='simanalyze') 

271 if nmodels>=1: 

272 skymodel=skymodels[0] 

273 else: 

274 skymodel="" 

275 

276 if os.path.exists(skymodel): 

277 msg("Sky model image "+skymodel+" found.",origin='simanalyze') 

278 components_only=True 

279 elif not dryrun: 

280 msg("Can't find a model image in your project directory, named skymodel or compskymodel - output image will be created, but comparison with the input model is not possible.",priority="warn",origin='simanalyze') 

281 analyze=False 

282 

283 modelflat = skymodel+".flat" 

284 

285 if os.path.exists(skymodel): 

286 if not (os.path.exists(modelflat) or dryrun): 

287 myutil.flatimage(skymodel,verbose=verbose) 

288 

289 # modifymodel just collects info if skymodel==newmodel 

290 (model_refdir,model_cell,model_size, 

291 model_nchan,model_specrefval,model_specrefpix,model_width, 

292 model_stokes) = myutil.modifymodel(skymodel,skymodel, 

293 "","","","","",-1, 

294 flatimage=False) 

295 

296 cell_asec=qa.convert(model_cell[0],'arcsec')['value'] 

297 

298 ##################################################################### 

299 # clean if desired, use noisy image for further calculation if present 

300 # todo suggest a cell size from psf? 

301 ##################################################################### 

302 if image: 

303 

304 # make sure cell is defined 

305 if is_array_type(cell): 

306 if len(cell) > 0: 

307 cell0 = cell[0] 

308 else: 

309 cell0 = "" 

310 else: 

311 cell0 = cell 

312 if len(cell0)<=0: 

313 cell = model_cell 

314 if is_array_type(cell): 

315 if len(cell) == 1: 

316 cell = [cell[0],cell[0]] 

317 else: 

318 cell = [cell,cell] 

319 

320 # cells are positive by convention 

321 cell = [qa.abs(cell[0]),qa.abs(cell[1])] 

322 

323 # and imsize 

324 if is_array_type(imsize): 

325 if len(imsize) > 0: 

326 imsize0 = imsize[0] 

327 if len(imsize) > 1: 

328 imsize1 = imsize[1] 

329 else: 

330 imsize1 = imsize0 

331 else: 

332 imsize0 = -1 

333 else: 

334 imsize0 = imsize 

335 if imsize0 <= 0: 

336 imsize = [int(pl.ceil(qa.convert(qa.div(model_size[0],cell[0]),"")['value'])), 

337 int(pl.ceil(qa.convert(qa.div(model_size[1],cell[1]),"")['value']))] 

338 else: 

339 imsize=[imsize0,imsize1] 

340 

341 

342 if len(mstoimage) == 0: 

343 if tpmstoimage: 

344 sd_only = True 

345 else: 

346 msg("no measurement sets found to image",priority="error",origin='simanalyze') 

347 else: 

348 sd_only = False 

349 # get some quantities from the interferometric ms 

350 # TODO use something like aU.baselineStats for this, and the 90% baseline 

351 maxbase=0. 

352 if len(mstoimage)>1 and dryrun: 

353 msg("imaging multiple ms not possible in dryrun mode",priority="warn",origin="simanalyze") 

354 # TODO make work better for multiple MS 

355 for msfile in mstoimage: 

356 if os.path.exists(msfile): 

357 tb.open(msfile) 

358 rawdata = tb.getcol("UVW") 

359 tb.done() 

360 maxbase = max([max(rawdata[0,]),max(rawdata[1,])]) # in m 

361 psfsize = 0.3/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/maxbase*3600.*180/pl.pi # lambda/b converted to arcsec 

362 minimsize = 8* int(psfsize/cell_asec) 

363 elif dryrun: 

364 minimsize = min(imsize) 

365 psfsize = qa.mul(cell[0],3) # HACK 

366 else: 

367 raise RuntimeError(mstoimage+" not found.") 

368 

369 if imsize[0] < minimsize: 

370 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',origin='simanalyze') 

371 imsize[0] = minimsize 

372 if imsize[1] < minimsize: 

373 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',origin='simanalyze') 

374 imsize[1] = minimsize 

375 

376 tpimage=None 

377 # Do single dish imaging first if tpmstoimage exists. 

378 if tpmstoimage and os.path.exists(tpmstoimage): 

379 msg('creating image from ms: '+tpmstoimage,origin='simanalyze') 

380 #if len(mstoimage): 

381 # tpimage = project + '.sd.image' 

382 #else: 

383 # tpimage = project + '.image' 

384 tpimage = project + '.sd.image' 

385 tpimage = fileroot + "/" + tpimage 

386 

387 if len(mstoimage): 

388 if len(modelimage) and tpimage != modelimage and \ 

389 tpimage != fileroot+"/"+modelimage: 

390 msg("modelimage parameter set to "+modelimage+" but also creating a new total power image "+tpimage,priority="warn",origin='simanalyze') 

391 msg("assuming you know what you want, and using modelimage="+modelimage+" in deconvolution",priority="warn",origin='simanalyze') 

392 elif len(featherimage) and tpimage != featherimage and \ 

393 tpimage != fileroot+"/"+featherimage: 

394 msg("featherimage parameter set to "+featherimage+" but also creating a new total power image "+tpimage,priority="warn",origin='simanalyze') 

395 msg("assuming you know what you want, and using featherimage="+featherimage+" in feather",priority="warn",origin='simanalyze') 

396 

397 # Get PB size of TP Antenna 

398 # !! aveant will only be set if modifymodel or setpointings and in  

399 # any case it will the the aveant of the INTERFM array - we want the SD 

400 if os.path.exists(tpmstoimage): 

401 # antenna diameter 

402 tb.open(tpmstoimage+"/ANTENNA") 

403 diams = tb.getcol("DISH_DIAMETER") 

404 tb.close() 

405 aveant = pl.mean(diams) 

406 # theoretical antenna beam size 

407 pb_asec = sdbeamutil.primaryBeamArcsec(qa.tos(qa.convert(qa.quantity(model_specrefval),'GHz')),aveant,(0.75 if aveant==12.0 else 0.0),10.0) 

408 elif dryrun: 

409 aveant = 12.0 

410 pb_asec = pbcoeff*0.29979/qa.convert(qa.quantity(model_specrefval),'GHz')['value']/aveant*3600.*180/pl.pi 

411 else: 

412 raise RuntimeError(tpmstoimage+" not found.") 

413 

414 # default PSF from PB of antenna 

415 imbeam = {'major': qa.quantity(pb_asec,'arcsec'), 

416 'minor': qa.quantity(pb_asec,'arcsec'), 

417 'positionangle': qa.quantity(0.0,'deg')} 

418 

419 # Common imaging parameters 

420 sdim_param = dict(infiles=[tpmstoimage], overwrite=overwrite, 

421 phasecenter=model_refdir, mode='channel', 

422 nchan=model_nchan, start=0, width=1) 

423 

424 if True: #SF gridding 

425 msg("Generating TP image using 'SF' kernel.",origin='simanalyze') 

426 beamsamp = 9 

427 sfcell_asec = pb_asec/beamsamp 

428 sfcell = qa.tos(qa.quantity(sfcell_asec, "arcsec")) 

429 cell_asec = [qa.convert(cell[0],"arcsec")['value'], 

430 qa.convert(cell[1],"arcsec")['value']] 

431 if cell_asec[0] > sfcell_asec or \ 

432 cell_asec[1] > sfcell_asec: 

433 # imregrid() may not work properly for regrid of 

434 # small to large cell 

435 msg("The requested cell size is too large to invoke SF gridding. Please set cell size <= %f arcsec or grid TP MS '%s' manually" % (sfcell_asec, tpmstoimage),priority="error",origin='simanalyze') 

436 

437 sfsupport = 6 

438 temp_out = tpimage+"0" 

439 temp_cell = [sfcell, sfcell] 

440 # too small - is imsize too small to start with? 

441 # needs to cover all pointings. 

442 temp_imsize = [int(pl.ceil(cell_asec[0]/sfcell_asec*imsize[0])), 

443 int(pl.ceil(cell_asec[1]/sfcell_asec*imsize[1]))] 

444 msg("Using predefined algorithm to define grid parameters.",origin='simanalyze') 

445 msg("SF gridding summary",origin='simanalyze') 

446 msg("- Antenna primary beam: %f arcsec" % pb_asec,origin='simanalyze') 

447 msg("- Image pixels per antenna PB (predefined): %f" % beamsamp,origin='simanalyze') 

448 msg("- Cell size (arcsec): [%s, %s]" % (temp_cell[0], temp_cell[1]),origin='simanalyze') 

449 msg("- Imsize to cover final TP image area: [%d, %d] (type: %s)" % (temp_imsize[0], temp_imsize[1], type(temp_imsize[0])),origin='simanalyze') 

450 msg("- convolution support: %d" % sfsupport,origin='simanalyze') 

451 # kernel specific imaging parameters 

452 sdim_param['gridfunction'] = 'SF' 

453 sdim_param['convsupport'] = sfsupport 

454 sdim_param['outfile'] = temp_out 

455 sdim_param['imsize'] = temp_imsize 

456 sdim_param['cell'] = temp_cell 

457 

458 msg(get_taskstr('sdimaging', sdim_param), priority="info") 

459 if not dryrun: 

460 sdimaging(**sdim_param) 

461 if not os.path.exists(temp_out): 

462 raise RuntimeError("TP imaging failed.") 

463 

464 # Scale image by convolved beam / antenna primary beam 

465 ia.open(temp_out) 

466 imbeam = ia.restoringbeam() 

467 ia.close() 

468 try: 

469 beam_area_ratio = qa.getvalue(qa.convert(imbeam['major'], "arcsec")) \ 

470 * qa.getvalue(qa.convert(imbeam['minor'], "arcsec")) \ 

471 / pb_asec**2 

472 except KeyError: 

473 # this shouldn't hit when simutil.imtclean sets restoringbeam='common' 

474 msg("using center channel values to calculate beam area ratio", 

475 origin='simanalyze', priority='info') 

476 chan_index = int(imbeam['nChannels']/2) 

477 pol_index = 0 

478 beam_area_ratio = qa.getvalue(qa.convert(imbeam['major'], "arcsec")) \ 

479 * qa.getvalue(qa.convert(imbeam['minor'], "arcsec")) \ 

480 / pb_asec**2 

481 msg("Scaling TP image intensity by %f." % (beam_area_ratio),origin='simanalyze') 

482 temp_in = temp_out 

483 temp_out = temp_out + ".scaled" 

484 immath(imagename=temp_in, mode='evalexpr', expr="IM0*%f" % (beam_area_ratio), 

485 outfile=temp_out) 

486 if not os.path.exists(temp_out): 

487 raise RuntimeError("TP image scaling failed.") 

488 

489 # Regrid TP image to final resolution 

490 msg("Regridding TP image to final resolution",origin='simanalyze') 

491 msg("- cell size (arecsec): [%s, %s]" % (cell[0], cell[1]),origin='simanalyze') 

492 msg("- imsize: [%d, %d]" % (imsize[0], imsize[1]),origin='simanalyze') 

493 if not dryrun: 

494 ia.open(temp_out) 

495 newcsys = ia.coordsys() 

496 ia.close() 

497 dir_idx = newcsys.findcoordinate("direction")['world'] 

498 newcsys.setreferencepixel([imsize[0]/2., imsize[1]/2.], 

499 type="direction") 

500 incr = newcsys.increment(type='direction')['numeric'] 

501 newincr = [incr[0]*cell_asec[0]/sfcell_asec, 

502 incr[1]*cell_asec[1]/sfcell_asec,] 

503 newcsys.setincrement(newincr, type="direction") 

504 # 

505 sdtemplate = imregrid(imagename=temp_out, template="get") 

506 sdtemplate['csys'] = newcsys.torecord() 

507 for idx in range(len(dir_idx)): 

508 sdtemplate['shap'][ dir_idx[idx] ] = imsize[idx] 

509 imregrid(imagename=temp_out, interpolation="cubic", 

510 template=sdtemplate, output=tpimage, 

511 overwrite=overwrite) 

512 del newcsys, sdtemplate, incr, newincr, dir_idx 

513 del temp_out, temp_cell, temp_imsize, sfcell_asec, cell_asec 

514 else: #PB grid 

515 msg("Generating TP image using 'PB' kernel.",origin='simanalyze') 

516 # Final TP cell and image size. 

517 # imsize and cell are already int and quantum arrays 

518 sdimsize = imsize 

519 sdcell = [qa.tos(cell[0]), qa.tos(cell[1])] 

520 ### TODO: need to set phasecenter properly based on imdirection 

521 # kernel specific imaging parameters 

522 sdim_param['gridfunction'] = 'PB' 

523 sdim_param['outfile'] = tpimage 

524 sdim_param['imsize'] = sdimsize 

525 sdim_param['cell'] = sdcell 

526 

527 msg(get_taskstr('sdimaging', sdim_param), priority="info") 

528 if not dryrun: 

529 sdimaging(**sdim_param) 

530 del sdimsize, sdcell 

531 # TODO: Define PSF of image here 

532 # for now use default  

533 

534 # get image beam size form TP image 

535 if os.path.exists(tpimage): 

536 ia.open(tpimage) 

537 beam = ia.restoringbeam() 

538 ia.close() 

539 

540 if sd_only: 

541 try: 

542 bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2 

543 except KeyError: 

544 # this shouldn't hit when simutil.imtclean sets restoringbeam='common' 

545 msg("using center channel values to calculate beam area", 

546 origin='simanalyze', priority='info') 

547 chan_index = int(beam['nChannels']/2) 

548 pol_index = 0 

549 bmarea = (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'] * 

550 beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'] * 

551 1.1331) #arcsec2 

552 bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix 

553 else: del beam 

554 #del beam 

555 

556 msg('generation of total power image '+tpimage+' complete.',origin='simanalyze') 

557 # update TP ms name the for following steps 

558 sdmsfile = tpmstoimage 

559 sd_any = True 

560 

561 imagename = re.split('.image$',tpimage)[0] 

562 # End of single dish imaging part 

563 

564 outflat_current = False 

565 convsky_current = False 

566 

567 if image and len(mstoimage) > 0: 

568 

569 # for reruns 

570 foo=mstoimage[0] 

571 foo=foo.replace(".ms","") 

572 foo=foo.replace(project,"") 

573 foo=foo.replace("/","") 

574 project=project+foo 

575 

576 imagename = fileroot + "/" + project 

577 

578 # get nfld, sourcefieldlist, from (interfm) ms if it was not just created 

579 # TODO make work better for multiple mstoimage (figures below) 

580 if os.path.exists(mstoimage[0]): 

581 tb.open(mstoimage[0]+"/SOURCE") 

582 code = tb.getcol("CODE") 

583 sourcefieldlist = pl.where(code=='OBJ')[0] 

584 nfld = len(sourcefieldlist) 

585 tb.done() 

586 elif dryrun: 

587 nfld=1 # HACK 

588 msfile = mstoimage[0] 

589 

590 # set cleanmode automatically (for interfm) 

591 if nfld == 1: 

592 cleanmode = "csclean" 

593 gridder = "standard" 

594 else: 

595 cleanmode = "mosaic" 

596 gridder = "mosaic" 

597 

598 # keep the same default minor cycle algorithm for imtclean 

599 # as previously used in imclean 

600 deconvolver = "clark" 

601 

602 

603 

604 # clean insists on using an existing model if its present 

605 if os.path.exists(imagename+".image"): shutil.rmtree(imagename+".image") 

606 if os.path.exists(imagename+".model"): shutil.rmtree(imagename+".model") 

607 

608 # An image in fileroot/ has priority 

609 if len(modelimage) > 0 and os.path.exists(fileroot+"/"+modelimage): 

610 modelimage = fileroot + "/" + modelimage 

611 msg("Found modelimage, %s." % modelimage,origin='simanalyze') 

612 

613 # in simdata we use imdirection instead of model_refdir 

614 if not myutil.isdirection(imdirection,halt=False): 

615 imdirection=model_refdir 

616 

617 myutil.imtclean(mstoimage,imagename, 

618 gridder,deconvolver, 

619 cell,imsize,imdirection, 

620 interactive,niter,threshold, 

621 weighting,outertaper,pbcor,stokes, 

622 modelimage=modelimage,mask=mask, 

623 dryrun=dryrun) 

624 ## Disabled as part of CAS-12502 

625 #myutil.imclean(mstoimage,imagename, 

626 # cleanmode,cell,imsize,imdirection, 

627 # interactive,niter,threshold,weighting, 

628 # outertaper,pbcor,stokes,  

629 # #sourcefieldlist=sourcefieldlist, 

630 # modelimage=modelimage,mask=mask,dryrun=dryrun) 

631 

632 

633 # create imagename.flat and imagename.residual.flat: 

634 if not dryrun: 

635 myutil.flatimage(imagename+".image",verbose=verbose) 

636 myutil.flatimage(imagename+".residual",verbose=verbose) 

637 outflat_current = True 

638 

639 # feather 

640 if featherimage: 

641 if not os.path.exists(featherimage): 

642 raise RuntimeError("Could not find featherimage "+featherimage) 

643 else: 

644 featherimage="" 

645 if tpimage: 

646 # if you set modelimage, then it won't force tpimage into  

647 # featherimage. this could be hard to explain  

648 # to the user. 

649 if os.path.exists(tpimage) and not os.path.exists(modelimage): 

650 featherimage=tpimage 

651 

652 

653 if os.path.exists(featherimage): 

654 msg("feathering the interfermetric image "+imagename+".image with "+featherimage,origin='simanalyze',priority="info") 

655 # TODO call with params? 

656 msg("feather('"+imagename+".feather.image','"+imagename+".image','"+featherimage+"')",priority="info") 

657 if not dryrun: 

658 feather(imagename+".feather.image",imagename+".image",featherimage) 

659 # copy residual flat image 

660 shutil.copytree(imagename+".residual.flat",imagename+".feather.residual.flat") 

661 imagename=imagename+".feather" 

662 # but replace combined flat image 

663 myutil.flatimage(imagename+".image",verbose=verbose) 

664 

665 

666 

667 if verbose: 

668 msg(" ") 

669 msg("done inverting and cleaning",origin='simanalyze') 

670 if not is_array_type(cell): 

671 cell = [cell,cell] 

672 if len(cell) <= 1: 

673 cell = [qa.quantity(cell[0]),qa.quantity(cell[0])] 

674 else: 

675 cell = [qa.quantity(cell[0]),qa.quantity(cell[1])] 

676 cell = [qa.abs(cell[0]),qa.abs(cell[0])] 

677 

678 # get beam from output clean image 

679 if verbose: msg("getting beam from "+imagename+".image",origin='simanalyze') 

680 if os.path.exists(imagename+".image"): 

681 ia.open(imagename+".image") 

682 beam = ia.restoringbeam() 

683 ia.close() 

684 # model has units of Jy/pix - calculate beam area from clean image 

685 # (even if we are not plotting graphics) 

686 try: 

687 bmarea = beam['major']['value']*beam['minor']['value']*1.1331 #arcsec2 

688 except KeyError: 

689 # this shouldn't hit when simutil.imtclean sets restoringbeam='common' 

690 msg("using center channel values to calculate beam area", 

691 origin='simanalyze', priority='info') 

692 chan_index = int(beam['nChannels']/2) 

693 pol_index = 0 

694 bmarea = (beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['major']['value'] * 

695 beam['beams']['*'+str(chan_index)]['*'+str(pol_index)]['minor']['value'] * 

696 1.1331) #arcsec2 

697 bmarea = bmarea/(cell[0]['value']*cell[1]['value']) # bm area in pix 

698 msg("synthesized beam area in output pixels = %f" % bmarea,origin='simanalyze') 

699 

700 

701 

702 

703 

704 

705 

706 

707 

708 if image: 

709 # show model, convolved model, clean image, and residual 

710 if grfile: 

711 file = fileroot + "/" + project + ".image.png" 

712 else: 

713 file = "" 

714 else: 

715 mslist=[] 

716 

717 if dryrun: 

718 grscreen=False 

719 grfile=False 

720 analyze=False 

721 

722 if image and len(mstoimage) > 0: 

723 if grscreen or grfile: 

724 myutil.newfig(multi=[2,2,1],show=grscreen) 

725 

726 # create regridded and convolved sky model image 

727 myutil.convimage(modelflat,imagename+".image.flat") 

728 convsky_current = True # don't remake this for analysis in this run 

729 

730 disprange = [] # passing empty list causes return of disprange 

731 

732 # original sky regridded to output pixels but not convolved with beam 

733 discard = myutil.statim(modelflat+".regrid",disprange=disprange,showstats=False)