Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_simalma.py: 1%
977 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 19:10 +0000
1import os
2import shutil
3import re
4#import pdb
6from casatools import ctsys
7from casatasks import concat, imregrid, immath, sdimaging, impbcor, simobserve, simanalyze, feather, casalog
8from .simutil import *
9from . import sdbeamutil
11def simalma(
12 project=None,
13 dryrun=None,
14 skymodel=None, inbright=None, indirection=None, incell=None,
15 incenter=None, inwidth=None,
16 complist=None, compwidth=None,
17 ########
18 setpointings=None,
19 ptgfile=None,
20 integration=None, direction=None, mapsize=None,
21 antennalist=None,
22 hourangle=None,
23 totaltime=None,
24 ###
25 tpnant = None,
26 tptime = None,
27 ###
28 pwv=None,
29 image=None,
30 imsize=None, imdirection=None,cell=None,
31 niter=None, threshold=None,
32 graphics=None,
33 verbose=None,
34 overwrite=None
35 ):
37 # Collect a list of parameter values to save inputs
38 in_params = locals()
40 #-------------------------
41 # Create the utility object
42 myutil = simutil(direction)
43 if verbose: myutil.verbose = True
44 msg = myutil.msg
46 try:
48 ###########################
49 # preliminaries
50 ###########################
52 # Predefined parameters
53 pbcoeff = 1.13 ## PB defined as pbcoeff*lambda/d
54 nyquist = 0.5/1.13 ## Nyquist spacing = PB*nyquist
55 maptype_int = 'ALMA'
56 maptype_tp = 'square'
58 # time ratios for 12extended, 12compact,7m, and TP:
59 default_timeratio = [1,0.5,2,4]
61 # pbgridratio_tp = 0.25 # -> would be gridratio_tp = 1/3.4
62 # the grid spacing is defined in terms of lambda/d times these factors
63 # ALMA OT uses lambda/d/sqrt(3)
64 gridratio_int = 1./pl.sqrt(3) # 0.5 is nyquist
65 gridratio_tp = 1./3
67 # number of 12m primary beams to pad the total power image during
68 # the gridding stage (i.e. even larger pad than the padding
69 # added for the observation).
70 tppad_npb = 2.
72 # weight of 7m data relative to 12m data
73 weightratio_7_12 = 0.34
76 # the scale factor to correct expected Gauss PB size to empirical simPB
77 simpb_factor = 0.96
78 caldirection = ""
79 calflux = "0Jy"
80 tpantid = 0
81 t_ground = 270.
82 if pwv > 0:
83 thermalnoise = "tsys-atm"
84 else:
85 thermalnoise = ""
86 leakage = 0.
87 weighting = "briggs"
89 antlist_tp_default="aca.tp.cfg"
91 #----------------------------------------
92 # Save outputs in a directory called "project"
93 fileroot = project
94 # simalma is not supposed to run multiple times.
95 if os.path.exists(fileroot):
96 infomsg = "Project directory, '%s', already exists." % fileroot
97 if overwrite:
98 casalog.post(infomsg)
99 casalog.post("Removing old project directory '%s'" % fileroot)
100 shutil.rmtree(fileroot)
101 else:
102 raise RuntimeError(infomsg)
104 if not os.path.exists(fileroot):
105 os.mkdir(fileroot)
107 concatname=fileroot+"/"+project+".concat"
109 #-------------------------
110 # set up reporting to file, terminal, logger
111 casalog.origin('simalma')
112 if verbose:
113 casalog.filter(level="DEBUG2")
114 v_priority = "INFO"
115 else:
116 v_priority = None
118 simobserr = "simalma caught an exception in task simobserve"
119 simanaerr = "simalma caught an exception in task simanalyze"
121 # open report file
122 myutil.reportfile=fileroot+"/"+project+".simalma.report.txt"
123 myutil.openreport()
127 #----------------------------------------
128 # Get globals to call saveinputs()
129 casalog.post("saveinputs not available in casatasks, skipping saving simalma inputs", priority='WARN')
131 # filename parsing of cfg file here so that the project filenames
132 # can contain the cfg
133 repodir = ctsys.resolve("alma/simmos")
135 #--------------------------
136 # format mapsize
137 if not is_array_type(mapsize):
138 mapsize = [mapsize, mapsize]
139 elif len(mapsize) < 2:
140 mapsize = [mapsize[0], mapsize[0]]
142 #---------------------------
143 # Operation flags
144 addnoise = (thermalnoise != '')
145 # Rectangle setup mode
146 multiptg = (not setpointings) \
147 or (is_array_type(direction) and len(direction) > 1)
148 rectmode = (not multiptg)
150 # Use full model image as a mapsize = ["", ""]
151 fullsize = (len(mapsize[0]) == 0) or (len(mapsize[1]) == 0)
161 ###########################
162 # analyze input sky model
163 ###########################
164 msg("="*60,priority="INFO")
165 msg("Checking the sky model",priority="INFO",origin="simalma")
166 msg(" ",priority="INFO")
168 #----------------------------
169 # Either skymodel or complist should exist
170 if is_array_type(skymodel):
171 if len(skymodel)>1:
172 msg("You have given more than one skymodel - simalma will only use the first one, "+skymodel[0],priority="INFO")
173 skymodel = skymodel[0]
174 skymodel = skymodel.replace('$project',project)
176 if is_array_type(complist):
177 if len(complist)>1:
178 msg("You have given more than one componentlist - simalma will only use the first one, "+componentlist[0],priority="INFO")
179 complist = complist[0]
181 if((not os.path.exists(skymodel)) and (not os.path.exists(complist))):
182 if len(skymodel)>0:
183 msg("Your skymodel '"+skymodel+"' could not be found.",priority="warn")
184 if len(complist)>0:
185 msg("Your complist '"+complist+"' could not be found.",priority="warn")
186 if len(skymodel)==0 and len(complist)==0:
187 msg("At least one of skymodel or complist must be set.",priority="error")
189 else:
190 msg("No sky input found. At least one of skymodel or complist must exist.",priority="error")
193 #-------------------------
194 # Get model_size and model_center
195 # TODO: check if outmodel==inmodel works (just collect info)
196 if os.path.exists(skymodel):
197 components_only=False
198 outmodel = fileroot+"/"+project+"temp.skymodel"
199 model_vals = myutil.modifymodel(skymodel, outmodel, inbright,
200 indirection, incell, incenter,
201 inwidth, -1, False)
202 shutil.rmtree(outmodel)
203 model_refdir = model_vals[0]
204 model_cell = model_vals[1]
205 model_size = model_vals[2]
206 model_nchan = model_vals[3]
207 model_center = model_vals[4]
208 model_width = model_vals[5]
209 del model_vals
210 msg("You will be simulating from sky model image "+skymodel,priority="info")
211 msg(" pixel(cell) size = "+ qa.tos(model_cell[0]),priority="info")
212 msg(" image size = "+ qa.tos(model_size[0]),priority="info")
213 msg(" reference direction = "+model_refdir,priority="info")
214 if model_nchan>1:
215 msg(" # channels = "+ qa.tos(model_nchan),priority="info")
216 msg(" channel width = "+ qa.tos(model_width),priority="info")
217 else:
218 msg(" single channel / continuum image, with bandwidth = "+qa.tos(model_width),priority="info")
220 if os.path.exists(complist):
221 msg(" ",priority="info")
222 msg("You will also be simulating the components in "+complist,priority="info")
223 msg(" These will get added to a copy of the skymodel image.",priority="info")
224 else:
225 # XXX TODO make sure components AND image work here
226 components_only=True
227 compdirs = []
228 cl.open(complist)
229 for i in range(cl.length()):
230 compdirs.append(myutil.dir_m2s(cl.getrefdir(i)))
232 model_refdir, coffs = myutil.average_direction(compdirs)
233 model_center = cl.getspectrum(0)['frequency']['m0']
234 cmax = 0.0014 # ~5 arcsec
235 for i in range(coffs.shape[1]):
236 xc = pl.absolute(coffs[0,i]) # offsets in deg
237 yc = pl.absolute(coffs[1,i])
238 cmax = max(cmax, xc)
239 cmax = max(cmax, yc)
240 model_size = ["%fdeg" % (2*cmax), "%fdeg" % (2*cmax)]
241 cl.done()
242 model_cell = None
243 model_nchan = 1
244 del compdirs, coffs, xc, yc, cmax
245 msg("You will be simulating only from component list "+complist,priority="info")
246 msg("Based on the spatial distribution of components, a sky model image will be generated with these parameters:",priority="info")
247 msg(" image size = "+ qa.tos(model_size[0]),priority="info")
248 msg(" reference direction = "+model_refdir,priority="info")
249 msg("and each call to simobserve will chose a skymodel cell size of 1/20 the expected PSF FWHM (to accurately locate components).",priority="info")
250 msg("Simulation from components-only produces a single channel / continuum observation",priority="info")
251 msg(" with bandwidth = "+compwidth,priority="info")
254 #-----------------------------------
255 # determine mapsize - copied code from simobserve
256 # if the user has not input a map size, then use model_size
257 if len(mapsize) == 0:
258 mapsize = model_size
259 msg("setting map size to "+str(model_size))
260 else:
261 if type(mapsize) == type([]):
262 if len(mapsize[0]) == 0:
263 mapsize = model_size
264 msg("setting map size to "+str(model_size))
266 if components_only:
267 # if map is greater tham model (defined from components)
268 # then use mapsize as modelsize
269 if type(mapsize) == type([]):
270 map_asec = qa.convert(mapsize[0],"arcsec")['value']
271 else:
272 map_asec = qa.convert(mapsize,"arcsec")['value']
273 if type(model_size) == type([]):
274 mod_asec = qa.convert(model_size[0],"arcsec")['value']
275 else:
276 mod_asec = qa.convert(model_size,"arcsec")['value']
277 if map_asec>mod_asec:
278 model_size=["%farcsec" % map_asec,"%farcsec" % map_asec]
281 msg(" ",priority="info")
282 msg("You will be mapping an area of size "+qa.tos(mapsize[0])+','+qa.tos(mapsize[1]))
291 ###########################
292 # Calculate 12-m PB and grid spacing for int and tp
293 ###########################
295 Dant = 12.
296 wave_length = 0.2997924/ qa.convert(qa.quantity(model_center),'GHz')['value']
297 # (wave length)/D_ant in arcsec
298 lambda_D = wave_length / Dant * 3600. * 180 / pl.pi
299 PB12 = qa.quantity(lambda_D*pbcoeff, "arcsec")
300 # Correction factor for PB in simulation
301 # (PS of simulated image is somehow smaller than PB12)
302 PB12sim = qa.mul(PB12, simpb_factor)
304 msg(" primary beam size: %s" % (qa.tos(PB12)),priority="info")
306 # Pointing spacing of observations - define in terms of lambda/d
307 # to avoid ambiguities in primary beam definition.
308 # it would be best for simobserve to accept the shorthand "LD"
309 # but until it does that, we'll divide by pbcoeff and use "PB":
310 # this is a fragile solution since it depends on pbcoeff being
311 # the same in simalma and simobserve.
312 pointingspacing = str(gridratio_int/pbcoeff)+"PB"
313 ptgspacing_tp = qa.tos(qa.mul(PB12sim, (gridratio_tp/pbcoeff)))
324 ###########################
325 # analyze antennalist(s)
326 ###########################
328 if type(antennalist)==type(" "):
329 antennalist=[antennalist]
331 # number of arrays/configurations to run:
332 nconfigs=len(antennalist)
333 msg(" ",priority="info")
334 msg("="*60,priority="info")
335 msg("You are requesting %i configurations:" % nconfigs,origin="simalma",priority="info")
337 configtypes=[] # ALMA, ACA, or ALMASD (latter is special case)
338 cycles=[]
339 resols=[]
340 tp_inconfiglist=False
342 # check filename consistency, get resolution, etc
343 for configfile in antennalist:
345 isalma=0
347 # antennalist should contain ALMA or ACA
348 if configfile.find(";")>=0:
349 telescopename="BYRES"
350 configparms=configfile.split(";")
351 res=configparms[-1]
352 res_arcsec=-1
353 if myutil.isquantity(res):
354 if qa.compare(res,"arcsec"):
355 res_arcsec=qa.convert(res,"arcsec")['value']
357 if res_arcsec>0:
358 resols.append(res_arcsec)
359 else:
360 msg("simalma cannot interpret the antennalist entry '"+configfile+"' as desired array and resolution e.g. ALMA;0.5arcsec",priority="ERROR")
362 else:
363 # we can only verify explicit files right now, not
364 # configs specified as "ALMA;0.5arcsec"
365 #
366 configfile_short = (configfile.split("/"))[-1]
367 # Search order is fileroot/ -> specified path -> repository
368 if len(configfile) > 0:
369 if os.path.exists(fileroot+"/"+configfile):
370 configfile = fileroot + "/" + configfile
371 elif not os.path.exists(configfile) and \
372 os.path.exists(os.path.join(repodir,configfile)):
373 configfile = os.path.join(repodir,configfile)
374 # Now make sure the configfile exists
375 if not os.path.exists(configfile):
376 msg("Couldn't find configfile: %s" % configfile, priority="error")
377 stnx, stny, stnz, stnd, padnames, antnames, telescopename, posobs = myutil.readantenna(configfile)
378 if telescopename=="ALMASD":
379 resols.append(qa.convert(PB12,'arcsec')['value'])
380 else:
381 psfsize = myutil.approxBeam(configfile,qa.convert(qa.quantity(model_center),'GHz')['value'])
382 resols.append(psfsize) # FWHM in arcsec
384 q = re.compile('CYCLE\d?\d')
385 isCycle = q.search(configfile_short.upper())
386 if isCycle:
387 whatCycle = isCycle.group()[-1] # will break if cycle>9
388 cycles.append(whatCycle)
389 else:
390 cycles.append(-1) # -1 is unknown; defaults to full ALMA
392 if configfile_short.upper().find("ALMA") >= 0:
393 if telescopename=="ALMA" or telescopename=="BYRES":
394 configtypes.append("ALMA")
395 isalma=isalma+1
396 else:
397 msg("Inconsistent antennalist entry '"+configfile_short+"' has ALMA in the name but not set as the observatory",priority="error")
398 if configfile_short.upper().find("ACA") >= 0:
399 if telescopename=="ACA" or telescopename=="BYRES":
400 configtypes.append("ACA")
401 isalma=isalma+1
402 elif telescopename=="ALMASD" or telescopename=="BYRES":
403 tp_inconfigfile=True
404 configtypes.append("ALMASD")
405 isalma=isalma+1
406 else:
407 msg("Inconsistent antennalist entry '"+configfile_short+"' has ACA in the name but the observatory is not ACA or ALMASD",priority="error")
408 #if configfile.upper().find("CYCLE0")
411 if isalma==0:
412 s="simalma can't accept antennalist entry '"+configfile_short+"' because it is neither ALMA nor ACA (in the name and the observatory in the file)"
413 msg(s,origin="simalma",priority="error")
414# raise ValueError, s # not ness - msg2 raises the exception
415 if isalma==2:
416 s="simalma doesn't understand your antennalist entry '"+configfile_short+"'"
417 msg(s,origin="simalma",priority="error")
418# raise ValueError,s
420 #-----------------------------
421 # total power parameter:
422 tptime_min=0.
423 if myutil.isquantity(tptime,halt=False):
424 if qa.compare(tptime,'s'):
425 tptime_min=qa.convert(tptime,'min')['value']
426 else:
427 msg("Can't interpret tptime='"
428 +tptime+"' as a time quantity e.g. '3h'",
429 priority="error")
430 else:
431 msg("Can't interpret tptime='"
432 +tptime+"' as a time quantity e.g. '3h'",
433 priority="error")
436 #-----------------------------
437 # is there a mix of cycle0,1,2,etc,full?
438 whatCycle=cycles[0]
439 cyclesInconsistent=False
440 for i in range(nconfigs):
441 if cycles[i]!=whatCycle:
442 cyclesInconsistent=True
445 #-----------------------------
446 # exposure time parameter totaltime
447 totaltime_min=[]
448 if len(totaltime)==1:
449 # scalar input - use defaults
450 totaltime_min0=0.
451 if not myutil.isquantity(totaltime[0],halt=False):
452 raise ValueError("Can't interpret totaltime parameter '"+totaltime[0]+"' as a time quantity - example quantities: '1h', '20min', '3600sec'")
453 if qa.compare(totaltime[0],'s'):
454 totaltime_min0=qa.convert(totaltime[0],'min')['value']
455 else:
456 raise ValueError("Can't convert totaltime parameter '"+totaltime[0]+"' to minutes - example quantities: '1h', '20min','3600sec'")
457 totaltime_min=pl.zeros(nconfigs)
458 # sort by res'l - TP could still be on here
459 resols=pl.array(resols)
460 intorder=resols.argsort()
461 # assume the scalar user input refers to the highest res 12m
462 totaltime_min[intorder[0]]=totaltime_min0
463 for j in intorder[1:]:
464 if configtypes[j]=='ALMA':
465 # lower res 12m
466 totaltime_min[j]=totaltime_min0*default_timeratio[1]
467 elif configtypes[j]=='ACA':
468 # 7m
469 totaltime_min[j]=totaltime_min0*default_timeratio[2]
470 elif configtypes[j]=='ALMASD':
471 # tp
472 totaltime_min[j]=totaltime_min0*default_timeratio[3]
473 else:
474 raise ValueError("configuration types = "+str(configtypes))
475 else:
476 for time in totaltime:
477 if not myutil.isquantity(time,halt=False):
478 raise ValueError("Can't interpret totaltime vector element '"+time+"' as a time quantity - example quantities: '1h', '20min', '3600sec'")
479 if qa.compare(time,'s'):
480 time_min=qa.convert(time,'min')['value']
481 else:
482 raise ValueError("Can't convert totaltime vector element '"+time+"' to minutes - example quantities: '1h', '20min','3600sec'")
483 totaltime_min.append(time_min)
485 if len(totaltime_min)!=len(antennalist):
486 raise ValueError("totaltime must either be the same length vector as antennalist or a scalar")
490 #-----------------------------
491 # print out what's requested.
492 antlist_tp=None
493 for i in range(nconfigs):
494 configfile=antennalist[i]
495 msg(" ",priority="info")
496 msg(" "+configfile+":",priority="info")
498 if configtypes[i]=="ALMA":
499 msg(" 12m ALMA array",priority="info")
500 elif configtypes[i]=="ACA":
501 msg(" 7m ACA array",priority="info")
502 elif configtypes[i]=="ALMASD":
503 msg(" 12m total power observation",priority="info")
504 if tpnant>0:
505 msg(" This antennalist entry will be ignored in favor of the tpnant requesting %d total power antennas" % tpnant,priority="info")
506 antlist_tp=antlist_tp_default
507 if tptime_min>0:
508 msg(" The total map integration time will be %s minutes as per the tptime parameter." % tptime_min,priority="info")
509 else:
510 tptime_min=totaltime_min[i]
511 msg(" The total map integration time will be %s minutes as per the totaltime parameter." % tptime_min,priority="info")
512 else:
513 # Note: assume the user won't put in >1 TP antlist.
514 msg(" Note: it is preferred to specify total power with the tpnant,tptime parameters so that one can also specify the number of antennas.",priority="info")
515 msg(" Specifying the total power array as one of the antenna configurations will only allow you to use a single antenna for the simulation",priority="info")
516 msg(" simalma will proceed with a single total power antenna observing for %d minutes." % totaltime_min[i],priority="info")
517 tpnant=1
518 tptime_min=totaltime_min[i]
519 antlist_tp=configfile
521 # print cycle and warn if mixed
522 if cycles[i]>'-1':
523 msg(" This is a cycle "+cycles[i]+" configuration",priority="info")
524 else:
525 msg(" This is a full ALMA configuration",priority="info")
526 if cyclesInconsistent:
527 msg(" WARNING: Your choices of configurations mix different cycles and/or Full ALMA. Assuming you know what you want.",priority="info")
530 if configtypes[i]=='ALMASD':
531 # imsize check for PB:
532 if not components_only:
533 minsize = min(qa.convert(model_size[0],'arcsec')['value'],\
534 qa.convert(model_size[1],'arcsec')['value'])
535 PB12sec=qa.convert(PB12,"arcsec")['value']
536 if minsize < 2.5*PB12sec:
537 msg(" WARNING: For TP imaging, skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*PB12sec),priority="warn")
538 msg(" Total power imaging when the source is this small is not necessary, and may fail with an error.",priority="warn")
539 del minsize
540 else:
541 # casalog.post resolution for INT, and warn if model_cell too large
542 msg(" approximate synthesized beam FWHM = %f arcsec" % resols[i],priority="info")
543 msg(" (at zenith; the actual beam will depend on declination, hourangle, and uv coverage)",priority="info")
545 if is_array_type(model_cell):
546 cell_asec=qa.convert(model_cell[0],'arcsec')['value']
547 else:
548 cell_asec=qa.convert(model_cell,'arcsec')['value']
549 if cell_asec > 0.2*resols[i]:
550 if cell_asec >= resols[i]:
551 # XXX if not dryrun raise exception here
552 msg(" ERROR: your sky model cell %f arcsec is too large to simulate this beam. Simulation will fail" % cell_asec,priority="info")
553 else:
554 msg(" WARNING: your sky model cell %f arcsec does not sample this beam well. Simulation may be unreliable or clean may fail." % cell_asec,priority="info")
556 msg(" Observation time = %f min" % totaltime_min[i],priority="info")
557 if totaltime_min[i]>360:
558 msg(" WARNING: this is longer than 6hr - simalma won't (yet) split the observation into multiple nights, so your results may be unrealistic",priority="info")
561 if not tp_inconfiglist and tpnant>0:
562 msg("",priority="info")
563 if tptime_min>0:
564 msg("You are also requesting Total Power observations:",priority="info")
565 msg(" %d antennas for %f minutes" % (tpnant,tptime_min),priority="info")
566 # imsize check for PB:
567 if not components_only:
568 minsize = min(qa.convert(model_size[0],'arcsec')['value'],\
569 qa.convert(model_size[1],'arcsec')['value'])
570 PB12sec=qa.convert(PB12,"arcsec")['value']
571 if minsize < 2.5*PB12sec:
572 msg(" WARNING: For TP imaging, skymodel should be larger than 2.5*primary beam. Your skymodel: %.3f arcsec < %.3f arcsec: 2.5*primary beam" % (minsize, 2.5*PB12sec),priority="warn")
573 msg(" Total power imaging when the source is this small is not necessary, and may fail with an error.",priority="warn")
575 else:
576 msg("You have requested %d total power antennas (tpnant), but no finite integration (tptime) -- check your inputs; no Total Power will be simulated." % tpnant,priority="info")
578 ### WORKAROUND for wrong flux in COMP TP simulations
579 if (tpnant > 0) and os.path.exists(complist):
580 idx_min = pl.where(resols == min(resols))[0]
581 idx = idx_min[0] if len(idx_min) > 0 else 0
582 dummy_proj = "gen_skymodel"
583 errmsg = "You requested Single dish simulation with components list.\n"
584 errmsg += "Single dish simulation has flux recovery issue "+\
585 "when using a components list as an input.\n"
586 errmsg += "Please generate compskymodel image first by task "+\
587 "simobserve and use the image as the skymodel input. "
588 errmsg += "Sorry for the inconvenience.\n\n"
589 errmsg += "How to workaround the issue:\n"
590 errmsg += "1. Generate skymodel image by simobserve\n"
591 errmsg += ("\tsimobserve(project='%s', complist='%s', compwidth='%s', "\
592 % (dummy_proj, complist, compwidth))
593 if os.path.exists(skymodel):
594 skysuffix = '.skymodel'
595 errmsg += ( "skymodel='%s', inbright='%s', indirection='%s', " \
596 % (skymodel, inbright, indirection))
597 errmsg += ( "incell='%s', incenter='%s', inwidth='%s', " \
598 % (skymodel, inbright, indirection, incell, incenter, inwidth) )
599 else:
600 skysuffix = '.compskymodel'
601 errmsg += ("setpointings=True, obsmode='', antennalist='%s', thermalnoise='')\n" \
602 % antennalist[idx])
603 errmsg += "2. Use the generated skymodel image in project directory as an input of simalma.\n"
604 errmsg += ("\tsimalma(project='%s', skymodel='%s/%s', complist='', ....)" % \
605 (project, dummy_proj, \
606 get_data_prefix(antennalist[idx], dummy_proj)+skysuffix))
607 msg(errmsg,priority="error")
608 ### End of WORKAROUND
610 # remove tp from configlist
611 antennalist=pl.array(antennalist)
612 configtypes=pl.array(configtypes)
613 totaltime_min=pl.array(totaltime_min)
614 resols=pl.array(resols)
615 z=pl.where(configtypes!='ALMASD')[0]
616 antennalist=antennalist[z]
617 configtypes=configtypes[z]
618 totaltime_min=totaltime_min[z]
619 resols=resols[z]
620 nconfigs=len(antennalist)
621 if nconfigs < 1:
622 msg("No interferometer configuration is requested. At least one interferometer configuration should be selected.", \
623 origin="simalma", priority="error")
626# TODO check model_size against mapsize - separately after this?
641# ###########################
642# # Resolve prefixes of simulation data (as defined in
643# # simobserve and simanalyze)
644#
645# pref_bl = get_data_prefix(antennalist, project)
646# pref_aca = get_data_prefix(antlist_aca, project)
647 # Resolve output names (as defined in simobserve and simanalyze)
648# if addnoise:
649# msname_bl = pref_bl+".noisy.ms"
650# msname_aca = pref_aca+".noisy.ms"
651# msname_tp = pref_tp+".noisy.sd.ms"
652# #imagename_aca = pref_aca+".noisy.image"
653# else:
654# msname_bl = pref_bl+".ms"
655# msname_aca = pref_aca+".ms"
656# msname_tp = pref_tp+".sd.ms"
657# #imagename_aca = pref_aca+".image"
658#
659# imagename_tp = project+".sd.image"
660# imagename_int = project+".concat.image"
661# msname_concat = project+".concat.ms"
662#
663 combimage = project+".feather.image"
664 simana_file = project+".simanalyze.last"
669 ############################################################
670 # run simobserve
671 msg("",priority="info")
672 msg("="*60,priority="info")
673 msg("simalma calls simobserve to simulate each component:",priority="info",origin="simalma")
675 # sort by res'l - save the ptgfile from the max res'l config
676 # for use in TP pointing definition.
677 resols=pl.array(resols)
678 intorder=resols.argsort()
679 pref_bl=get_data_prefix(antennalist[intorder[0]],project)
680 ptgfile_bl=fileroot+"/"+pref_bl+".ptg.txt"
682 step = 0
684 for i in range(nconfigs):
685 if configtypes[i]=="ALMA":
686 s="12m ALMA array"
687 elif configtypes[i]=='ACA':
688 s="7m ACA array"
689 else:
690 s="12m total power map"
692 if configtypes[i]!='ALMASD':
693 step += 1
695 msg("",priority="info")
696 msg("-"*60, origin="simalma", priority="warn")
697 msg(("Step %d: simulating " % step)+s, origin="simalma", priority="warn")
698 msg("-"*60, origin="simalma", priority="warn")
700 # filename prefixes
701 pref = get_data_prefix(antennalist[i], project)
703 obsmode_int = 'int'
704 ptgfile_int = fileroot+"/"+pref+".ptg.txt"
706 task_param = {}
707 task_param['project'] = project
708 task_param['skymodel'] = skymodel
709 task_param['inbright'] = inbright
710 task_param['indirection'] = indirection
711 task_param['incell'] = incell
712 task_param['incenter'] = incenter
713 task_param['inwidth'] = inwidth
714 task_param['complist'] = complist
715 task_param['compwidth'] = compwidth
716 task_param['setpointings'] = setpointings
717 task_param['ptgfile'] = ptgfile
718 task_param['integration'] = integration
719 task_param['direction'] = direction
720 task_param['mapsize'] = mapsize
721 task_param['maptype'] = maptype_int
722 # this is approriate for 7m or 12m since its in terms of PB
723 task_param['pointingspacing'] = pointingspacing
724 task_param['caldirection'] = caldirection
725 task_param['calflux'] = calflux
726 task_param['obsmode'] = obsmode_int
727 task_param['hourangle'] = hourangle
728 task_param['totaltime'] = "%fmin" % totaltime_min[i]
729 task_param['antennalist'] = antennalist[i]
730 task_param['sdantlist'] = ""
731 task_param['sdant'] = 0
732 task_param['thermalnoise'] = thermalnoise
733 task_param['user_pwv'] = pwv
734 task_param['t_ground'] = t_ground
735 task_param['leakage'] = leakage
736 task_param['graphics'] = graphics
737 task_param['verbose'] = verbose
738 task_param['overwrite'] = overwrite
740 msg(get_taskstr('simobserve', task_param), priority="info")
741 if not dryrun:
742 try:
743 simobserve(**task_param)
744 del task_param
745 except:
746 raise RuntimeError(simobserr)
747 finally:
748 casalog.origin('simalma')
750 qimgsize_tp = None
752 if tpnant>0 and tptime_min<=0:
753 raise ValueError("You requested total power (tpnant=%d) but did not specify a valid nonzero tptime" % tpnant)
756 mslist_tp = []
757 if tptime_min>0:
758 ########################################################
759 # ACA-TP simulation - always do this last
760 obsmode_sd = "sd"
761 step += 1
763 msg(" ",priority="info")
764 msg("-"*60, origin="simalma", priority="warn")
765 msg(("Step %d: simulating Total Power" % step), origin="simalma", priority="warn")
766 msg("-"*60, origin="simalma", priority="warn")
768 if antlist_tp is None:
769 antlist_tp = antlist_tp_default
771 pref_tp = get_data_prefix(antlist_tp, project)
772 if addnoise:
773 msname_tp = pref_tp+".noisy.sd.ms"
774 else:
775 msname_tp = pref_tp+".sd.ms"
778 ###########################
779 # Resolve pointing directions of ACA-TP.
780 #
781 # Pointing directions of TP simulation is defined as follows:
782 #
783 # [I] if ALMA-12m maps a rectangle region (rectmode=T),
784 # TP maps slightly larger region than ALMA-12m by adding 1 PB to
785 # mapsize (pointing extent of ALMA-12m).
786 #
787 # [II] if a list of pointing deirections are specified for the
788 # ALMA-12m observation (multiptg=T), TP pointings are defined as
789 # rectangle areas of 2PB x 2PB centered at each pointing direction
790 # of ALMA-12m. However, in some cases, it is more efficient to
791 # just map a rectangle region that covers all ALMA-12m pointings.
792 # In such case, ACA-TP maps a rectangle region whose extent is 2PB
793 # larger than the extent of all ALMA-12m pointings.
794 if rectmode:
795 # Add 1PB to mapsize
796 if fullsize:
797 mapx = qa.add(PB12,model_size[0]) # in the unit same as PB
798 mapy = qa.add(PB12,model_size[1]) # in the unit same as PB
799 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)]
800 msg("The full skymodel is being mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", origin="simalma", priority='warn')
801 else:
802 # mapsize is defined. Add 1 PB to mapsize.
803 mapx = qa.add(qa.quantity(mapsize[0]), PB12)
804 mapy = qa.add(qa.quantity(mapsize[1]), PB12)
805 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)]
806 msg("Only part of the skymodel is being mapped by ALMA 12-m and ACA 7-m arrays. The total power antenna observes 1PB larger extent.", priority='warn')
807 else:
808 # multi-pointing mode
809 npts, pointings, time = myutil.read_pointings(ptgfile_bl)
810 center, offsets = myutil.average_direction(pointings)
811 del time
812 qx = qa.quantity(max(offsets[0])-min(offsets[0]),"deg")
813 qy = qa.quantity(max(offsets[1])-min(offsets[1]),"deg")
814 # map extent to cover all pointings + 2PB
815 mapx = qa.add(qa.mul(PB12,2.),qx) # in the unit same as PB
816 mapy = qa.add(qa.mul(PB12,2.),qy) # in the unit same as PB
817 mapsize_tp = [qa.tos(mapx), qa.tos(mapy)]
818 # number of pointings to map vicinity of each pointings
819 qptgspc_tp = qa.quantity(ptgspacing_tp)
820 dirs_multi_tp = myutil.calc_pointings2(qptgspc_tp,
821 qa.tos(qa.mul(PB12,2.)),
822 "square", pointings[0])
823 npts_multi = npts * len(dirs_multi_tp)
825 msg("Number of pointings to map vicinity of each direction = %d" % npts_multi)
826 del qptgspc_tp, dirs_multi_tp
828 # save imsize for TP image generation
829# qimgsize_tp = [mapx, mapy]
830 # for imaging later, we need to pad even more!
831 mapx_im = qa.add(mapx, qa.mul(PB12,(2*tppad_npb)))
832 mapy_im = qa.add(mapy, qa.mul(PB12,(2*tppad_npb)))
833 qimgsize_tp = [mapx_im, mapy_im]
836 qptgspc_tp = qa.quantity(ptgspacing_tp)
837 pbunit = PB12['unit']
838 # number of pointings to map pointing region
839 # TODO: use calc pointings for consistent calculation
840 npts_rect = int(qa.convert(mapx, pbunit)['value'] \
841 / qa.convert(qptgspc_tp, pbunit)['value']) \
842 * int(qa.convert(mapy, pbunit)['value'] \
843 / qa.convert(qptgspc_tp, pbunit)['value'])
844 msg("Number of pointings to map a rect region = %d" % npts_rect, priority="DEBUG2")
846 if rectmode:
847 dir_tp = direction
848 npts_tp = npts_rect
849 msg("Rectangle mode: The total power antenna observes 1PB larger region compared to ALMA 12-m and ACA 7-m arrays", priority='info')
850 else:
851 if npts_multi < npts_rect:
852 # Map 2PB x 2PB extent centered at each pointing direction
853 # need to get a list of pointings
854 dir_tp = []
855 locsize = qa.mul(2, PB12)
856 for dir in pointings:
857 dir_tp += myutil.calc_pointings2(qa.tos(qptgspc_tp),
858 qa.tos(locsize),
859 "square", dir)
861 mapsize_tp = ["", ""]
862 #npts_tp = npts_multi
863 npts_tp = len(dir_tp)
864 msg("Multi-pointing mode: The total power antenna observes +-1PB of each point", origin="simalma", priority='warn')
865 else:
866 # Map a region that covers all directions
867 dir_tp = center
868 npts_tp = npts_rect
869 msg("Multi-pointing mode: The total power antenna maps a region that covers all pointings", origin="simalma", priority='warn')
870 msg("- Center of poinings: %s" % center, origin="simalma", priority='warn')
871 msg("- Map size: [%s, %s]" % (mapsize_tp[0], mapsize_tp[1]), origin="simalma", priority='warn')
874 # Scale integration time of TP (assure >= 1 visit per direction)
875 tottime_tp = "%dmin" % tptime_min
876 integration_tp = integration
877 ndump = int(qa.convert(tottime_tp, 's')['value']
878 / qa.convert(integration, 's')['value'])
879 msg("Max number of dump in %s (integration %s): %d" % \
880 (tottime_tp, integration, ndump), origin="simalma")
882 if ndump < npts_tp:
883 t_scale = float(ndump)/float(npts_tp)
884 integration_tp = qa.tos(qa.mul(integration, t_scale))
885 msg("Integration time is scaled to cover all pointings in observation time.", origin="simalma", priority='warn')
886 msg("-> Scaled total power integration time: %s" % integration_tp, origin="simalma", priority='warn')
887 ## Sometimes necessary to avoid the effect of round-off error
888 #iunit = qa.quantity(integration_tp)['unit']
889 #intsec = qa.convert(integration_tp,"s")
890 #totsec = intsec['value']*npts_tp#+0.000000001)
891 ##tottime_tp = qa.tos(qa.convert(qa.quantity(totsec, "s"), iunit))
892 #tottime_tp = qa.tos(qa.quantity(totsec, "s"))
894 task_param = {}
895 task_param['project'] = project
896 task_param['skymodel'] = skymodel
897 task_param['inbright'] = inbright
898 task_param['indirection'] = indirection
899 task_param['incell'] = incell
900 task_param['incenter'] = incenter
901 task_param['inwidth'] = inwidth
902 task_param['complist'] = complist
903 task_param['compwidth'] = compwidth
904 task_param['setpointings'] = True
905 task_param['ptgfile'] = '$project.ptg.txt'
906 task_param['integration'] = integration_tp
907 task_param['direction'] = dir_tp
908 task_param['mapsize'] = mapsize_tp
909 task_param['maptype'] = maptype_tp
910 task_param['pointingspacing'] = ptgspacing_tp
911 task_param['caldirection'] = caldirection
912 task_param['calflux'] = calflux
913 task_param['obsmode'] = obsmode_sd
914 task_param['hourangle'] = hourangle
915 task_param['totaltime'] = "%dmin" % tptime_min
916 task_param['antennalist'] = ""
917 task_param['sdantlist'] = antlist_tp
918 task_param['sdant'] = tpantid
919 task_param['thermalnoise'] = thermalnoise
920 task_param['user_pwv'] = pwv
921 task_param['t_ground'] = t_ground
922 task_param['leakage'] = leakage
923 task_param['graphics'] = graphics
924 task_param['verbose'] = verbose
925 task_param['overwrite'] = overwrite
927 msg(" ",priority="info")
928 msg("Simulating %d TP antennas" % tpnant, priority="info")
929 for iant in range(tpnant):
930 task_param['sdant'] = iant
931 task_param['seed'] = int(pl.random()*100000)
932 msg(" ",priority=v_priority)
933 msg("Running TP simulation with sdant = %d" % task_param['sdant'], priority=v_priority)
934 msg(get_taskstr('simobserve', task_param), priority="info")
935 if not dryrun:
936 try:
937 simobserve(**task_param)
938 except:
939 raise RuntimeError(simobserr)
940 finally:
941 casalog.origin('simalma')
942 if tpnant == 1:
943 mslist_tp = [msname_tp]
944 else:
945 # copy MSes
946 orig_tp = fileroot + "/" + msname_tp
947 suffix = (".Ant%d" % iant)
948 msg(" ",priority=v_priority)
949 msg("Renaming '%s' to '%s'" % (orig_tp, orig_tp+suffix), priority=v_priority)
950 if not dryrun:
951 shutil.move(orig_tp, orig_tp+suffix)
953 mslist_tp.append(msname_tp+suffix)
954 # noiseless MS
955 if addnoise:
956 orig_tp = orig_tp.replace(".noisy", "")
957 msg("Renaming '%s' to '%s'" % (orig_tp, orig_tp+suffix), priority=v_priority)
958 if not dryrun:
959 shutil.move(orig_tp, orig_tp+suffix)
961 del task_param
973 ################################################################
974 # Imaging
975 if dryrun: image=True # why not?
977 # for 4.2 print more info
978 v_priority="info"
980 if image:
981 modelimage = ""
982 imagename_tp = project+".sd.image"
983 if tptime_min > 0:
984 ########################################################
985 # Image ACA-TP
986 step += 1
987 msg(" ",priority="info")
988 msg("-"*60, origin="simalma", priority="info")
989 msg("Step %d: generating a total power image. " % step, origin="simalma", priority="info")
990 msg(" WARNING: Optimal gridding parameters are being analyzed by ALMA and may change.",priority="warn",origin="simalma")
991 msg("-"*60, origin="simalma", priority="info")
993 vis_tp = []
994 for msname_tp in mslist_tp:
995 if dryrun:
996 vis_tp.append(fileroot+"/"+msname_tp)
997 elif os.path.exists(fileroot+"/"+msname_tp):
998 vis_tp.append(fileroot+"/"+msname_tp)
999 else:
1000 msg("Total power MS '%s' is not found" \
1001 % msname_tp, origin="simalma", priority="error")
1003 tp_kernel = 'SF'
1004 #tp_kernel = 'GJINC'
1006 # Define imsize to cover TP map region
1007 msg(" ",priority=v_priority)
1008 msg("Defining image size to cover map region of total power simulation", priority=v_priority)
1009 msg("-> The total power map size: [%s, %s]" % \
1010 (qa.tos(qimgsize_tp[0]), qa.tos(qimgsize_tp[1])), \
1011 priority=v_priority)
1012 if tp_kernel.upper() == 'SF':
1013 beamsamp=9
1014 pb_asec = sdbeamutil.primaryBeamArcsec(qa.tos(qa.convert(qa.quantity(model_center),'GHz')),12.0,0.75,10.0)
1015 qcell=qa.quantity(pb_asec/beamsamp, 'arcsec')
1016 cell_tp = [qa.tos(qcell), qa.tos(qcell)]
1017 msg("Using fixed cell size for SF grid: %s" % cell_tp[0], \
1018 priority=v_priority)
1019 elif cell != '':
1020 # user-defined cell size
1021 msg("The user defined cell size: %s" % cell, \
1022 priority=v_priority)
1023 cell_tp = [qa.tos(cell), qa.tos(cell)]
1024 else:
1025 if model_cell == None:
1026 # components only simulation
1027 compmodel = fileroot+"/"+pref_bl+".compskymodel"
1028 msg("getting the cell size of input compskymodel", \
1029 priority=v_priority)
1030 if not os.path.exists(compmodel):
1031 msg("Could not find the skymodel, '%s'" % \
1032 compmodel, priority='error')
1033 # modifymodel just collects info if outmodel==inmodel
1034 model_vals = myutil.modifymodel(compmodel,compmodel,
1035 "","","","","",-1,
1036 flatimage=False)
1037 model_cell = model_vals[1]
1038 model_size = model_vals[2]
1040 cell_tp = [qa.tos(model_cell[0]), qa.tos(model_cell[1])]
1042 #####################################################
1044 imsize_tp = calc_imsize(mapsize=qimgsize_tp, cell=cell_tp)
1046 msg(" ",priority=v_priority)
1047 msg("-> The number of pixels needed to cover the map region: [%d, %d]" % \
1048 (imsize_tp[0], imsize_tp[1]), \
1049 priority=v_priority)
1051 msg("Compare with interferometer image area and adopt the larger one:", \
1052 priority=v_priority)
1053 # Compare with imsize of BL (note: imsize is an intArray)
1054 imsize_bl = []
1055 if is_array_type(imsize) and imsize[0] > 0:
1056 if len(imsize) > 1:
1057 imsize_bl = imsize[0:2]
1058 else:
1059 imsize_bl = [imsize[0], imsize[0]]
1061 if tp_kernel.upper() == 'SF':
1062 if len(imsize_bl) > 0:
1063 if cell != '':
1064 # user defined cell size for INT
1065 imarea = [qa.tos(qa.mul(cell, imsize[0])),
1066 qa.tos(qa.mul(cell, imsize[1]))]
1067 else:
1068 # using model_cell for INT
1069 imarea = [qa.tos(qa.mul(model_cell[0], imsize[0])),
1070 qa.tos(qa.mul(model_cell[1], imsize[1]))]
1071 tmpimsize = calc_imsize(mapsize=imarea, cell=cell_tp)
1072 else:
1073 msg("estimating imsize from input sky model.", \
1074 priority=v_priority)
1075 tmpimsize = calc_imsize(mapsize=model_size, cell=cell_tp)
1076 msg("-> TP imsize to cover interferometrer image area: [%d, %d]" % \
1077 (tmpimsize[0], tmpimsize[1]), \
1078 priority=v_priority)
1079 elif len(imsize_bl) > 0:
1080 # User has defined imsize
1081 tmpimsize = imsize_bl
1082 msg("-> Interferometer imsize (user defined): [%d, %d]" % \
1083 (tmpimsize[0], tmpimsize[1]), \
1084 priority=v_priority)
1085 else:
1086 # the same as input model (calculate from model_size)
1087 msg("estimating imsize of interferometer from input sky model.", \
1088 priority=v_priority)
1089 tmpimsize = calc_imsize(mapsize=model_size, cell=cell_tp)
1090 msg("-> Estimated interferometer imsize (sky model): [%d, %d]" % \
1091 (tmpimsize[0], tmpimsize[1]), \
1092 priority=v_priority)
1095 imsize_tp = [max(imsize_tp[0], tmpimsize[0]), \
1096 max(imsize_tp[1], tmpimsize[1])]
1098 msg("The image pixel size of TP: [%d, %d]" % \
1099 (imsize_tp[0], imsize_tp[1]), \
1100 priority=v_priority)
1102 # Generate TP image
1103 msg(" ",priority=v_priority)
1104 temp_out = fileroot+"/"+imagename_tp + '0'
1105 task_param = dict(infiles=vis_tp, outfile=temp_out, imsize=imsize,
1106 cell=cell_tp, phasecenter=model_refdir,
1107 mode='channel', nchan= model_nchan)
1108 if tp_kernel.upper() == 'SF':
1109 msg("Generating TP image using 'SF' kernel.",\
1110 priority=v_priority)
1111 # Parameters for sdimaging
1112 task_param['gridfunction'] = 'sf'
1113 task_param['convsupport'] = 6
1114 else:
1115 msg("Generating TP image using 'GJinc' kernel.",\
1116 priority=v_priority)
1117 gfac = 2.52 # b in Mangum et al. (2007)
1118 jfac = 1.55 # c in Mangum et al. (2007)
1119 convfac = 1.8 # The conversion factor to get HWHM of kernel roughly equal to qhwhm
1120 kernelfac = 0.7 # ratio of (kernel HWHM)/(TP pointingspacing)
1121 #qfwhm = PB12 # FWHM of GJinc kernel.
1122 #gwidth = qa.tos(qa.mul(gfac/3.,qfwhm))
1123 #jwidth = qa.tos(qa.mul(jfac/3.,qfwhm))
1124 qhwhm = qa.mul(qptgspc_tp, kernelfac) # hwhm of GJinc kernel
1125 gwidth = qa.tos(qa.mul(qhwhm, convfac))
1126 jwidth = qa.tos(qa.mul(jfac/gfac/pl.log(2.),gwidth))
1127 #casalog.post("Kernel parameter: [qhwhm, gwidth, jwidth] = [%s, %s, %s]" % (qa.tos(qhwhm), gwidth, jwidth))
1128 # Parameters for sdimaging
1129 task_param['gridfunction'] = 'gjinc'
1130 task_param['gwidth'] = gwidth
1131 task_param['jwidth'] = jwidth
1133 casalog.post("saveinputs not available in casatasks, skipping saving sdimaging inputs", priority='WARN')
1135 msg("Having set up the gridding parameters, the sdimaging task is called to actually creage the image:",priority=v_priority)
1136 msg(get_taskstr('sdimaging', task_param), priority="info")
1138 if not dryrun:
1139 sdimaging(**task_param)
1140 #del task_param
1141 # TODO: scale TP image
1143 # Scale TP image
1144 if dryrun: #emulate beam calc in sdimaging
1145 bu = sdbeamutil.TheoreticalBeam()
1146 bu.set_antenna("12m", "0.75m")
1147 bu.set_sampling([ptgspacing_tp, ptgspacing_tp], "0.0deg")
1148 bu.set_image_param(task_param['cell'], model_center,task_param['gridfunction'],
1149 task_param['convsupport'] if 'convsupport' in task_param else -1,
1150 -1,
1151 task_param['gwidth'] if 'gwidth' in task_param else -1,
1152 task_param['jwidth'] if 'jwidth' in task_param else -1,
1153 is_alma=True)
1154 #bu.summary()
1155 imbeam = bu.get_beamsize_image()
1156 else:
1157 ia.open(temp_out)
1158 imbeam = ia.restoringbeam()
1159 ia.close()
1160 bmsize = imbeam['major']
1161 beam_area_ratio = qa.getvalue(qa.convert(bmsize, 'arcsec'))**2 \
1162 / qa.getvalue(qa.convert(PB12sim, 'arcsec'))**2
1163 msg(" ",priority=v_priority)
1164 msg("Scaling TP image intensity by beam area before and after gridding: %f" % beam_area_ratio)
1165 task_param = dict(imagename=temp_out, mode='evalexpr',
1166 expr=("IM0*%f" % (beam_area_ratio)),
1167 outfile = fileroot+"/"+imagename_tp)
1168 casalog.post("saveinputs not available in casatasks, skipping saving inmath inputs", priority='WARN')
1170 msg(get_taskstr('immath', task_param), priority="info")
1171 if not dryrun:
1172 immath(**task_param)
1173 del task_param
1175 # Analyze TP image
1176 tpskymodel = fileroot+"/"+pref_tp+".skymodel"
1177 if components_only:
1178 tpskymodel = fileroot+"/"+pref_tp+".compskymodel"
1180 msg(" ",priority="info")
1181 msg("Analyzing TP image", priority="info")
1183 task_param = {}
1184 task_param['project'] = project
1185 task_param['image'] = False
1186 task_param['imagename'] = fileroot+"/"+imagename_tp
1187 task_param['skymodel'] = tpskymodel
1188 task_param['analyze'] = True
1189 task_param['showuv'] = False
1190 task_param['showpsf'] = False
1191 task_param['showconvolved'] = True
1192 task_param['graphics'] = graphics
1193 task_param['verbose'] = verbose
1194 task_param['overwrite'] = overwrite
1195 task_param['dryrun'] = dryrun
1196 task_param['logfile'] = myutil.reportfile
1198 msg(get_taskstr('simanalyze', task_param), priority="info")
1200 try:
1201 myutil.closereport()
1202 simanalyze(**task_param)
1203 del task_param
1204 myutil.openreport()
1205 except:
1206 raise RuntimeError(simanaerr)
1207 finally:
1208 casalog.origin('simalma')
1210 # Back up simanalyze.last file
1211 if os.path.exists(fileroot+"/"+simana_file):
1212 simana_new = imagename_tp.replace(".image",".simanalyze.last")
1213 msg("Back up input parameter set to '%s'" % simana_new, \
1214 priority=v_priority)
1215 shutil.move(fileroot+"/"+simana_file, fileroot+"/"+simana_new)
1217 if not os.path.exists(fileroot+"/"+imagename_tp) and not dryrun:
1218 msg("TP image '%s' is not found" \
1219 % imagename_tp, priority="error")
1220 #modelimage = imagename_aca
1227 ############################################################
1228 # Image each INT array separately
1229 nconfig=len(antennalist)
1231 for i in range(nconfig):
1232 step += 1
1233 msg(" ",priority="info")
1234 msg("-"*60, origin="simalma", priority="info")
1235 msg(("Step %d: imaging and analyzing " % step)+antennalist[i], origin="simalma", priority="info")
1236 msg(" This step is optional, but useful to assess the result from just one configuration.",priority="warn",origin="simalma")
1237 msg(" WARNING: The example clean shown here uses no mask, may diverge, and almost certainly is not optimal.",priority="warn",origin="simalma")
1238 msg(" Users are HIGHLY recommended to use interactive clean masking (in simanalyze or directly in clean)",priority="warn",origin="simalma")
1239 msg(" Auto-masking is under development for use in the ALMA pipeline and will be included here in a future release",priority="warn",origin="simalma")
1240 msg("-"*60, origin="simalma", priority="info")
1242 pref=get_data_prefix(antennalist[i], project)
1243 if addnoise:
1244 msname = pref+".noisy.ms"
1245 imagename_int=pref+".noisy.image"
1246 else:
1247 msname= pref+".ms"
1248 imagename_int=pref+".image"
1250 if dryrun:
1251 vis_int = msname
1252 else:
1253 if os.path.exists(fileroot+"/"+msname):
1254 vis_int = msname
1255 else:
1256 msg("Could not find MS to image, '%s'" \
1257 % msname, origin="simalma", priority="error")
1259# i think simanalye is fixed 20130826
1260# # TMP fix: get correct skymodel file so that simanalyze picks it
1261# if acaratio > 0:
1262# if os.path.exists(tpskymodel):
1263# shutil.move(tpskymodel,tpskymodel+".save")
1264# else:
1265# msg("TP skymodel '%s' is not found" \
1266# % tpskymodel, origin="simalma", priority="error")
1267#
1268# if os.path.exists(acaskymodel+".save"):
1269# shutil.move(acaskymodel+".save",acaskymodel)
1270# else:
1271# msg("ACA skymodel '%s' is not found" \
1272# % acaskymodel+".save", origin="simalma", priority="error")
1274 # dryrun requires feeding cell and imsize from here.
1275 task_param = {}
1276 task_param['project'] = project
1277 task_param['image'] = image
1278 task_param['vis'] = vis_int
1279 task_param['modelimage'] = ""
1280 if cell:
1281 task_param['cell'] = cell
1282 else:
1283 task_param['cell'] = qa.tos(model_cell[0])
1284 if imsize[0]>0:
1285 task_param['imsize'] = imsize
1286 else:
1287 task_param['imsize'] = int(qa.div(model_size[0],model_cell[0])['value'])
1288 if len(imdirection)>0:
1289 task_param['imdirection'] = imdirection
1290 else:
1291 task_param['imdirection'] = model_refdir
1292 task_param['niter'] = niter
1293 task_param['threshold'] = threshold
1294 task_param['weighting'] = weighting
1295 task_param['mask'] = []
1296 task_param['outertaper'] = []
1297 task_param['stokes'] = 'I'
1298 task_param['analyze'] = True
1299 task_param['graphics'] = graphics
1300 task_param['verbose'] = verbose
1301 task_param['overwrite'] = overwrite
1302 task_param['dryrun'] = dryrun
1303 task_param['logfile'] = myutil.reportfile
1305 msg(get_taskstr('simanalyze', task_param), priority="info")
1307 try:
1308 myutil.closereport()
1309 simanalyze(**task_param)
1310 del task_param
1311 myutil.openreport()
1312 except:
1313 raise RuntimeError(simanaerr)
1314 finally:
1315 casalog.origin('simalma')
1319 if nconfig>1:
1320 ############################################################
1321 # concat
1322 step += 1
1323 msg(" ",priority="info")
1324 msg("-"*60, origin="simalma", priority="info")
1325 msg("Step %d: concatenating interferometric visibilities." % step, origin="simalma", priority="info")
1326 msg("-"*60, origin="simalma", priority="info")
1328 weights=pl.zeros(nconfig)+1
1329 z=pl.where(configtypes=='ACA')[0]
1330 if len(z)>0:
1331 weights[z]=weightratio_7_12
1333 mslist=[]
1334 for i in range(nconfig):
1335 pref=get_data_prefix(antennalist[i],project)
1336 if addnoise:
1337 msname = fileroot + "/" + pref+".noisy.ms"
1338 else:
1339 msname = fileroot + "/" + pref+".ms"
1340 mslist.append(msname)
1342 msg("concat(vis="+str(mslist)+",concatvis="+concatname+".ms,visweightscale="+str(weights.tolist()),priority="info")
1343 if not dryrun:
1344 try:
1345 concat(vis=mslist,concatvis=concatname+".ms",visweightscale=weights)
1346 except:
1347 raise RuntimeError(simanaerr)
1348 finally:
1349 casalog.origin('simalma')
1354 ############################################################
1355 # Image ALMA-BL + ACA-7m
1356 step += 1
1357 msg(" ",priority="info")
1358 msg("-"*60, origin="simalma", priority="info")
1359 msg(("Step %d: imaging and analyzing " % step)+concatname+".ms", origin="simalma", priority="info")
1360 msg(" WARNING: The example clean shown here uses no mask, may diverge, and almost certainly is not optimal.",priority="warn",origin="simalma")
1361 msg(" Users are HIGHLY recommended to use interactive clean masking (in simanalyze or directly in clean)",priority="warn",origin="simalma")
1362 msg(" Auto-masking is under development for use in the ALMA pipeline and will be included here in a future release",priority="warn",origin="simalma")
1363 msg("-"*60, origin="simalma", priority="info")
1366 if dryrun:
1367 vis_int = concatname+".ms"
1368 else:
1369 if os.path.exists(fileroot+"/"+concatname+".ms"):
1370 vis_int = fileroot+"/"+concatname+".ms"
1371 elif os.path.exists(concatname+".ms"):
1372 vis_int = concatname+".ms"
1373 else: msg("Could not find MS to image, "+concatname+".ms", origin="simalma", priority="error")
1375 task_param = {}
1376 task_param['project'] = project
1377 task_param['image'] = image
1378 task_param['vis'] = vis_int
1379 task_param['modelimage'] = ""
1380 if cell:
1381 task_param['cell'] = cell
1382 else:
1383 task_param['cell'] = model_cell
1384 if imsize[0]>0:
1385 task_param['imsize'] = imsize
1386 else:
1387 task_param['imsize'] = int(qa.div(model_size[0],model_cell[0])['value'])
1388 if len(imdirection)>0:
1389 task_param['imdirection'] = imdirection
1390 else:
1391 task_param['imdirection'] = model_refdir
1392 task_param['niter'] = niter
1393 task_param['threshold'] = threshold
1394 task_param['weighting'] = weighting
1395 task_param['mask'] = []
1396 task_param['outertaper'] = []
1397 task_param['stokes'] = 'I'
1398 task_param['analyze'] = True
1399 task_param['graphics'] = graphics
1400 task_param['verbose'] = verbose
1401 task_param['overwrite'] = overwrite
1402 task_param['dryrun'] = dryrun
1403 task_param['logfile'] = myutil.reportfile
1405 msg(get_taskstr('simanalyze', task_param), priority="info")
1406 imagename_int=os.path.basename(concatname.rstrip("/"))+".image"
1408 try:
1409 myutil.closereport()
1410 simanalyze(**task_param)
1411 del task_param
1412 myutil.openreport()
1413 except:
1414 raise RuntimeError(simanaerr)
1415 finally:
1416 casalog.origin('simalma')
1425 if tptime_min > 0:
1426 ########################################################
1427 # Combine TP + INT image
1428 step += 1
1429 msg(" ",priority="info")
1430 msg("-"*60, origin="simalma", priority="info")
1431 msg("Step %d: combining a total power and synthesis image. " % step, origin="simalma", priority="info")
1432 msg(" WARNING: feathering the two images is only one way to combine them. ",priority="warn",origin="simalma")
1433 msg(" Using the total power image as a model in cleaning the interferometric visibilities may work better in some circumstances.",priority="warn",origin="simalma")
1434 msg("-"*60, origin="simalma", priority="info")
1436 if os.path.exists(fileroot+"/"+imagename_int) or dryrun:
1437 highimage0 = fileroot+"/"+imagename_int
1438 else:
1439 msg("The synthesized image '%s' is not found" \
1440 % imagename_int, origin="simalma", priority="error")
1441 if (not os.path.exists(fileroot+"/"+imagename_tp)) and (not dryrun):
1442 msg("ACA is requested but total power image '%s' is not found" \
1443 % imagename_tp, origin="simalma", priority="error")
1444 #lowimage = fileroot+"/"+imagename_tp
1446 # Need to manipulate TP image here
1447 outimage0 = fileroot+"/" + combimage+"0"
1448 outimage = fileroot+"/" + combimage
1449 # CAS-13369 -- imclean call in simanalyze is now imtclean
1450 #pbcov = highimage0.rstrip("image") + "flux.pbcoverage"
1451 pbcov = highimage0.rstrip("image") + "pb"
1452 regridimg = fileroot + "/" + imagename_tp + ".regrid"
1453 scaledimg = fileroot + "/" + imagename_tp + ".pbscaled"
1454 lowimage = scaledimg
1456 msg(" ",priority="info")
1457 msg("Regrid total power image to interferometric image grid:",priority=v_priority)
1458 # regrid TP image
1459 msg("inttemplate = imregrid(imagename = '"+highimage0+"', template='get')",priority=v_priority)
1460 if not dryrun:
1461 inttemplate = imregrid(imagename = highimage0, template='get')
1462 msg("imregrid(imagename = '"+fileroot+"/"+imagename_tp+
1463 "',interpolation='cubic',template = inttemplate, output = '"+regridimg+"')",priority="info")
1464 if not dryrun:
1465 imregrid(imagename = fileroot+"/"+imagename_tp,
1466 interpolation="cubic",
1467 template = inttemplate, output = regridimg)
1468 # multiply SD image with INT PB coverage
1469 if not os.path.exists(pbcov):
1470 msg("The flux image '%s' is not found" \
1471 % pbcov, origin="simalma", priority="error")
1473 msg(" ",priority="info")
1474 msg("Multiply total power image by interferometric sensitivity map:",priority=v_priority)
1475# msg("immath(imagename=['"+regridimg+"','"+pbcov+"'],expr='IM1*IM0',outfile='"+scaledimg+"')",priority="info")
1476# if not dryrun:
1477# immath(imagename=[regridimg, pbcov],
1478# expr='IM1*IM0',outfile=scaledimg)
1479#
1480# msg(" ",priority="info")
1481# msg("Set total power image beam and brightness unit:",priority=v_priority)
1482# msg("ia.open('"+fileroot+"/"+imagename_tp+"')",priority="info")
1483# msg("beam_tp = ia.restoringbeam()",priority="info")
1484# msg("bunit_tp = ia.brightnessunit()",priority="info")
1485# msg("ia.close()",priority="info")
1486# msg("ia.open('"+scaledimg+"')",priority="info")
1487# msg("ia.setrestoringbeam(beam=beam_tp)",priority="info")
1488# msg("ia.setbrightnessunit(bunit_tp)",priority="info")
1489# msg("ia.close()",priority="info")
1490#
1491# if not dryrun:
1492# pdb.set_trace()
1493# # restore TP beam and brightness unit
1494# ia.open(fileroot+"/"+imagename_tp)
1495# beam_tp = ia.restoringbeam()
1496# bunit_tp = ia.brightnessunit()
1497# ia.close()
1498# ia.open(scaledimg)
1499# ia.setrestoringbeam(beam=beam_tp)
1500# ia.setbrightnessunit(bunit_tp)
1501# ia.close()
1503 msg("impbcor('"+regridimg+"', '"+pbcov+"', outfile='"+scaledimg+"',mode='multiply')",priority="info")
1504 if not dryrun:
1505 impbcor(regridimg, pbcov, outfile=scaledimg,mode='multiply')
1507 # de-pbcor the INT image
1508 # not needed now that imtclean is used and .image from tclean is flat sky by default, see CAS-13370
1509 #highimage = fileroot+"/"+imagename_int+".pbscaled"
1510 #immath(imagename=[highimage0, pbcov],
1511 # expr='IM1/IM0',outfile=highimage)
1512# msg(" ",priority="info")
1513# msg("Multiply interferometric image by sensitivity map (un-pbcor):",priority="info")
1514# msg("immath(imagename=['"+highimage0+"','"+pbcov+"'],expr='IM1*IM0',outfile='"+highimage+"')",priority="info")
1515# msg("Restore interferometric beam and brightness unit:",priority="info")
1516# msg("ia.open('"+highimage0+"')",priority="info")
1517# msg("beam_int = ia.restoringbeam()",priority="info")
1518# msg("bunit_int = ia.brightnessunit()",priority="info")
1519# msg("ia.close()",priority="info")
1520# msg("ia.open('"+highimage+"')",priority="info")
1521# msg("ia.setrestoringbeam(beam=beam_int)",priority="info")
1522# msg("ia.setbrightnessunit(bunit_int)",priority="info")
1523# msg("ia.close()",priority="info")
1524#
1525# if not dryrun:
1526# immath(imagename=[highimage0, pbcov],
1527# expr='IM1*IM0',outfile=highimage)
1528# # restore INT beam and brightness unit
1529# ia.open(highimage0)
1530# beam_int = ia.restoringbeam()
1531# bunit_int = ia.brightnessunit()
1532# ia.close()
1533# ia.open(highimage)
1534# ia.setrestoringbeam(beam=beam_int)
1535# ia.setbrightnessunit(bunit_int)
1536# ia.close()
1537#
1538# msg("impbcor('"+highimage0+"', '"+pbcov+"', outfile='"+highimage+"',mode='multiply')",priority="info")
1539# if not dryrun:
1540# impbcor(highimage0, pbcov, outfile=highimage,mode='multiply')
1546 # Feathering
1547 task_param = {}
1548 task_param['imagename'] = outimage0
1549 task_param['highres'] = highimage0
1550 task_param['lowres'] = lowimage
1552 msg(" ",priority="info")
1553 msg(get_taskstr('feather', task_param), priority="info")
1554 try:
1555 casalog.post("saveinputs not available in casatasks, skipping saving feather inputs", priority='WARN')
1557 if not dryrun:
1558 feather(**task_param)
1559 del task_param
1561 # This seems not necessary anymore.
1562 ## transfer mask - feather should really do this
1563 #ia.open(outimage0)
1564 #ia.maskhandler('copy',[highimage+":mask0",'mask0'])
1565 #ia.maskhandler('set','mask0')
1566 #ia.done()
1567 except Exception as exc:
1568 raise Exception("simalma caught an exception in task feather: {}".
1569 format(exc))
1570 finally:
1571 if not dryrun: shutil.rmtree(regridimg)
1572 #shutil.rmtree(scaledimg)
1573 casalog.origin('simalma')
1576 # re-pbcor the result
1577 msg(" ",priority="info")
1578 msg("Re-apply the primary beam correction to the feathered result:",priority=v_priority)
1579# msg("immath(imagename=['"+outimage0+"','"+pbcov+"'],expr='IM0/IM1',outfile='"+outimage+"')",priority="info")
1580# if not dryrun:
1581# immath(imagename=[outimage0, pbcov],
1582# expr='IM0/IM1',outfile=outimage)
1584 msg("impbcor('"+outimage0+"', '"+pbcov+"', outfile='"+outimage+"')",priority="info")
1585 if not dryrun:
1586 impbcor(outimage0, pbcov, outfile=outimage)
1592 ########################################################
1593 # Generate Summary Plot
1594 grscreen = False
1595 grfile = False
1596 if not dryrun:
1597 if graphics == "both":
1598 grscreen = True
1599 grfile = True
1600 if graphics == "screen":
1601 grscreen = True
1602 if graphics == "file":
1603 grfile = True
1605 if grscreen or grfile:
1606 if grfile:
1607 file = fileroot + "/" + project + ".combine.png"
1608 else:
1609 file = ""
1611 # check for image pathes
1612 if os.path.exists(skymodel):
1613 flatsky = pref_bl + ".skymodel.flat"
1614 else:
1615 flatsky = pref_bl + ".compskymodel.flat"
1616 if not os.path.exists(fileroot+"/"+flatsky):
1617 raise RuntimeError("Coud not find a skymodel image '%s'" % flatsky)
1619 if not os.path.exists(fileroot+"/"+combimage):
1620 raise RuntimeError("Coud not find the combined image '%s'" % combimage)
1622 if not os.path.exists(fileroot+"/"+imagename_int):
1623 raise RuntimeError("Coud not find the synthesized image '%s'" % imagename_int)
1625 if not os.path.exists(fileroot+"/"+imagename_tp):
1626 raise RuntimeError("Coud not find the total power image '%s'" % (imagename_tp))
1627 # Now the actual plotting
1628 disprange = None
1629 myutil.newfig(multi=[2,2,1],show=grscreen)
1630 # skymodel
1631 #discard = myutil.statim(fileroot+"/"+flatsky,disprange=disprange)
1633 #disprange = []
1634 # generate flat synthesized (7m+12m) image
1635 flatint = fileroot + "/" + imagename_int + ".flat"
1636 myutil.flatimage(fileroot+"/"+imagename_int,verbose=verbose)
1637 if not os.path.exists(flatint):
1638 raise RuntimeError("Failed to generate '%s'" % (flatint))
1640 # generate convolved sky model image
1641 myutil.convimage(fileroot+"/"+flatsky, flatint)
1642 discard = myutil.statim(fileroot+"/"+flatsky+".regrid.conv",disprange=disprange)
1643 shutil.rmtree(fileroot+"/"+flatsky+".regrid")
1644 shutil.rmtree(fileroot+"/"+flatsky+".regrid.conv")
1646 # total power image
1647 flattp = fileroot + "/" + imagename_tp + ".flat"
1648 myutil.flatimage(fileroot+"/"+imagename_tp,verbose=verbose)
1649 #flattp = scaledimg + ".flat"
1650 #myutil.flatimage(scaledimg,verbose=verbose)
1651 if not os.path.exists(flattp):
1652 raise RuntimeError("Failed to generate '%s'" % (flattp))
1653 myutil.nextfig()
1654 discard = myutil.statim(flattp,disprange=disprange)
1655 shutil.rmtree(flattp)
1657 #disprange = []
1658 # display flat synthesized (7m+12m) image
1659 myutil.nextfig()
1660 discard = myutil.statim(flatint,disprange=disprange)
1661 shutil.rmtree(flatint)
1663 # combined image
1664 flatcomb = fileroot + "/" + combimage + ".flat"
1665 myutil.flatimage(fileroot+"/"+combimage,verbose=verbose)
1666 if not os.path.exists(flatcomb):
1667 raise RuntimeError("Failed to generate '%s'" % (flatcomb))
1668 myutil.nextfig()
1669 discard = myutil.statim(flatcomb,disprange=disprange)
1670 myutil.endfig(show=grscreen,filename=file)
1671 shutil.rmtree(flatcomb)
1673 if myutil.isreport():
1674 myutil.closereport()
1675 finalize_tools()
1677 finally:
1678 finalize_tools()
1679 if myutil.isreport():
1680 myutil.closereport()
1683def finalize_tools():
1684 if ia.isopen(): ia.close()
1685 sm.close()
1686 #cl.close() # crashes casa
1688def get_data_prefix(cfgname, project=""):
1689 if str.upper(cfgname[0:4]) == "ALMA":
1690 foo=cfgname.replace(';','_')
1691 else:
1692 foo = cfgname
1693 foo=foo.replace(".cfg","")
1694 sfoo=foo.split('/')
1695 if len(sfoo)>1: foo=sfoo[-1]
1696 return project+"."+foo
1698def calc_imsize(mapsize=None, cell=None):
1699 if mapsize == None:
1700 raise ValueError("mapsize is not defined")
1701 if cell == None:
1702 raise ValueError("cell is not defined")
1703 # get a list of cell size
1704 if is_array_type(cell):
1705 if len(cell) < 2:
1706 cell = [cell[0], cell[0]]
1707 else:
1708 cell = [cell, cell]
1710 for qval in cell:
1711 if not qa.compare(qval, "deg"):
1712 raise TypeError("cell should be an angular size")
1714 qcellx = qa.quantity(cell[0])
1715 qcelly = qa.quantity(cell[1])
1717 # get a list of map size
1718 if is_array_type(mapsize):
1719 if len(mapsize) < 2:
1720 mapsize = [mapsize[0], mapsize[0]]
1721 else:
1722 mapsize = [mapsize, mapsize]
1724 for qval in mapsize:
1725 if not qa.compare(qval, "deg"):
1726 raise TypeError("mapsize should be an angular size")
1728 vsizex = qa.convert(mapsize[0], qcellx['unit'])['value']
1729 vsizey = qa.convert(mapsize[1], qcelly['unit'])['value']
1731 # Calculate the number of pixels to cover the map size
1732 npixx = int(pl.ceil(vsizex/qcellx['value']))
1733 npixy = int(pl.ceil(vsizey/qcelly['value']))
1735 return [npixx, npixy]