Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/task_wvrgcal.py: 82%
217 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-01 07:19 +0000
1import os
2import sys
3import numpy as np
5from casatasks import casalog
6from casatools import ctsys, ms, quanta, calibrater, wvr
8def wvrgcal(vis=None, caltable=None, toffset=None, segsource=None,
9 sourceflag=None, tie=None, nsol=None, disperse=None,
10 wvrflag=None, statfield=None, statsource=None, smooth=None,
11 scale=None, spw=None, wvrspw=None,
12 reversespw=None, cont=None, maxdistm=None,
13 minnumants=None, mingoodfrac=None, usefieldtab=None,
14 refant=None, offsetstable=None, rseed=None):
15 """
16 Generate a gain table based on Water Vapour Radiometer data.
17 Returns a dictionary containing the RMS of the path length variation
18 with time towards that antenna (RMS) and the discrepency between the RMS
19 path length calculated separately for different WVR channels (Disc.).
21 vis -- Name of input visibility file
23 default: none; example: vis='ngc5921.ms'
25 caltable -- Name of output gain calibration table
26 default: none; example: caltable='ngc5921.wvr'
28 toffset -- Time offset (sec) between interferometric and WVR data
30 default: 0 (ALMA default for cycle 1, for cycle 0 it was -1)
32 segsource -- Do a new coefficient calculation for each source
34 default: True
36 tie -- Prioritise tieing the phase of these sources as well as possible
37 (requires segsource=True)
38 default: [] example: ['3C273,NGC253', 'IC433,3C279']
40 sourceflag -- Flag the WVR data for these source(s) as bad and do not produce corrections for it
41 (requires segsource=True)
42 default: [] (none) example: ['3C273']
44 nsol -- Number of solutions for phase correction coefficients during this observation.
45 By default only one set of coefficients is generated for the entire observation.
46 If more sets are requested, then they will be evenly distributed in time throughout'
47 the observation. Values > 1 require segsource=False.
48 default: 1
50 disperse -- Apply correction for dispersion
51 default: False
53 wvrflag -- Regard the WVR data for these antenna(s) as bad and use interpolated values instead
54 default: [] (none) example: ['DV03','DA05','PM02']
56 statfield -- Compute the statistics (Phase RMS, Disc) on this field only
57 default: '' (all)
59 statsource -- Compute the statistics (Phase RMS, Disc) on this source only
60 default: '' (all)
62 smooth -- Smooth the calibration solution on the given timescale
63 default: '' (no smoothing), example: '3s' smooth on a timescale of 3 seconds
65 scale -- Scale the entire phase correction by this factor
66 default: 1. (no scaling)
68 spw -- List of the spectral window IDs for which solutions should be saved into the caltable
69 default: [] (all spectral windows), example [17,19,21,23]
71 wvrspw -- List of the spectral window IDs from which the WVR data should be taken
72 default: [] (all WVR spectral windows), example [0]
74 reversespw -- Reverse the sign of the correction for the listed SPWs
75 (only needed for early ALMA data before Cycle 0)
76 default: '' (none), example: reversespw='0~2,4'; spectral windows 0,1,2,4
78 cont -- Estimate the continuum (e.g., due to clouds)
79 default: False
81 maxdistm -- maximum distance (m) an antenna may have to be considered for being part
82 of the <=3 antenna set for interpolation of a solution for a flagged antenna
83 default: 500
85 minnumants -- minimum number of near antennas required for interpolation
86 default: 2
88 mingoodfrac -- If the fraction of unflagged data for an antenna is below this value (0. to 1.),
89 the antenna is flagged.
90 default: 0.8
92 usefieldtab -- derive the antenna AZ/EL values from the FIELD rather than the POINTING table
93 default: False
95 refant -- use the WVR data from this antenna for calculating the dT/dL parameters (can give ranked list)
96 default: '' (use the first good or interpolatable antenna),
97 examples: 'DA45' - use DA45
98 ['DA45','DV51'] - use DA45 and if that is not good, use DV51 instead
100 offsetstable -- subtract the temperature offsets in this table from the WVR measurements before
101 using them to calculate the phase corrections
102 default: '' (do not apply any offsets)
103 examples: 'uid___A002_Xabd867_X2277.cloud_offsets' use the given table
105 rseed -- set random seed (integer) for the wvrgcal fitting routine to this specific value
106 default: 0 - use internal default value
107 example: 54321
109 """
111 # make ms tool local
112 myms = ms()
113 myqa = quanta()
115 ## parameters which are different in format between wvrgcal and wvr.gcal:
116 # reverse: only exists in wvr.gcal
117 # reversespw: string list in wvrgcal, list wvr.gcal
118 # wvrflag: list in wvrgcal, single string in wvr.gcal
119 # statfield: single string in wvrgcal, list in wvr.gcal
120 # statsource: single string in wvrgcal, list in wvr.gcal
121 # refant: list in wvrgcal, list in single string in wvr.gcal
122 # rseed: only exists in wvr.gcal
124 try:
125 casalog.origin('wvrgcal')
127 if not (type(vis)==str) or not (os.path.exists(vis)):
128 raise Exception('Visibility data set not found - please verify the name')
130 if (caltable == ""):
131 raise Exception("Must provide output calibration table name in parameter caltable.")
133 if os.path.exists(caltable):
134 raise Exception("Output caltable "+caltable+" already exists - will not overwrite.")
136 outdir = os.path.dirname(caltable)
137 if outdir == '':
138 outdir = '.'
139 if not os.access(outdir, os.W_OK):
140 raise Exception("Don't have write permission for output directory "+outdir)
142 vispar = vis
144 smoothpar = 1 # this is for the internal smoothing of wvr.gcal(), which we don't use
145 smoothing = -1
146 if (type(smooth)==str and smooth!=''):
147 smoothing = myqa.convert(myqa.quantity(smooth), 's')['value']
148 outputpar = caltable + '_unsmoothed' # as intermediate name before smoothing
149 else:
150 outputpar = caltable
152 toffsetpar = toffset
154 nsolpar = 1
155 if nsol>1:
156 if not segsource:
157 nsolpar = nsol
158 else:
159 raise Exception("In order to use nsol>1, segsource must be set to False.")
161 segsourcepar = segsource
163 sourceflagpar = []
164 if segsource and (len(sourceflag)>0):
165 sourceflagpar = sourceflag
166 for src in sourceflag:
167 if not (type(src)==int or type(src)==str) or src=='':
168 raise Exception("List elements of parameter sourceflag must be int or non-empty string.")
170 tiepar = []
171 if segsource and (len(tie)>0):
172 tiepar = tie
173 for i in range(0,len(tie)):
174 src = tie[i]
175 if not (type(src)==str) or src=='':
176 raise Exception("List elements of parameter tie must be non-emptystrings.")
178 spwpar = spw
179 if (len(spw)>0):
180 for myspw in spw:
181 if not (type(myspw)==int) or myspw<0:
182 raise Exception("List elements of parameter spw must be int and >=0.")
184 wvrspwpar = wvrspw
185 if (len(wvrspw)>0):
186 for myspw in wvrspw:
187 if not (type(myspw)==int) or myspw<0:
188 raise Exception("List elements of parameter wvrspw must be int and >=0.")
190 reversespwpar = []
191 if not (reversespw==''):
192 spws = myms.msseltoindex(vis=vis,spw=reversespw)['spw']
193 for id in spws:
194 reversespwpar.append(id)
196 dispersepar = disperse
197 if disperse:
198 dispdirpath = os.getenv('WVRGCAL_DISPDIR', '')
199 if not os.path.exists(dispdirpath+'/libair-ddefault.csv'):
200 path1 = dispdirpath
201 dispdirpath = ctsys.resolve("alma/wvrgcal")
202 if not os.path.exists(dispdirpath+'/libair-ddefault.csv'):
203 raise Exception("Dispersion table libair-ddefault.csv not found in path "\
204 +"given by WVRGCAL_DISPDIR nor in \""+dispdirpath+"\"")
206 os.putenv('WVRGCAL_DISPDIR', dispdirpath)
208 casalog.post('Using dispersion table '+dispdirpath+'/libair-ddefault.csv')
210 contpar = cont
211 if cont and segsource:
212 raise Exception("cont and segsource are not permitted to be True at the same time.")
214 usefieldtabpar = usefieldtab
216 offsetspar = offsetstable
218 wvrflagpar = ""
219 if (len(wvrflag)>0):
220 for ant in wvrflag:
221 if not (type(ant)==int or type(ant)==str):
222 raise Exception("List elements of parameter wvrflag must be int or string.")
223 if (ant != ''):
224 if len(wvrflagpar)>0:
225 wvrflagpar += ","
226 wvrflagpar += str(ant)
228 refantpar = ""
229 if (type(refant)!=list):
230 refant = [refant]
231 if (len(refant)>0):
232 for ant in refant:
233 if not (type(ant)==int or type(ant)==str):
234 raise Exception("Parameter refant must be int or string or a list of them.")
235 if (ant != ''):
236 if len(refantpar)>0:
237 refantpar += ","
238 refantpar += str(ant)
240 statfieldpar = []
241 if not (statfield==None or statfield=="") and type(statfield)==str:
242 statfieldpar = [statfield]
244 statsourcepar = []
245 if not (statsource==None or statsource=="") and type(statsource)==str:
246 statsourcepar = [statsource]
248 scalepar = scale
250 maxdistmpar = maxdistm
252 minnumantspar = minnumants
254 mingoodfracpar = mingoodfrac
256 rseedpar = 0
257 if type(rseed)==int and rseed>=0:
258 rseedpar = rseed
259 elif not (rseed==None or rseed==""):
260 raise Exception("Parameter rseed must be an integer >= 0 (the value 0 will use the internal default seed).")
262 casalog.post('Running wvr.gcal ...')
265 templogfile = 'wvrgcal_tmp_'+str(np.random.randint(1E6,1E8))
266 if not os.access(".", os.W_OK):
267 import tempfile
268 templogfile = tempfile.gettempdir()+"/"+templogfile
270 os.system('rm -rf '+templogfile)
272 mywvr = wvr()
274 rval = mywvr.gcal(vis=vispar,
275 output=outputpar,
276 toffset=toffsetpar,
277 nsol=nsolpar,
278 segsource=segsourcepar,
279 reverse=False, # only reverse those SPWs in reversespwpar
280 reversespw=reversespwpar,
281 disperse=dispersepar,
282 cont=contpar,
283 wvrflag=wvrflagpar,
284 sourceflag=sourceflagpar,
285 statfield=statfieldpar,
286 statsource=statsourcepar,
287 tie=tiepar,
288 smooth=smoothpar,
289 scale=scalepar,
290 maxdistm=maxdistmpar,
291 minnumants=minnumantspar,
292 mingoodfrac=mingoodfracpar,
293 usefieldtab=usefieldtabpar,
294 spw=spwpar,
295 wvrspw=wvrspwpar,
296 refant=refantpar,
297 offsets=offsetspar,
298 rseed=rseedpar,
299 logfile=templogfile)
301 loglines = []
302 with open(templogfile) as f:
303 loglines = f.readlines()
304 for i in range(len(loglines)):
305 loglines[i] = loglines[i].expandtabs()
306 casalog.post(''.join(loglines))
308 # prepare variables for parsing log lines to extract info table
309 hfound = False
310 hend = False
311 namel = []
312 wvrl = []
313 flagl = []
314 rmsl = []
315 discl = []
316 flagsd = {}
317 parsingok = True
319 for ll in loglines:
320 if hfound:
321 if "Expected performance" in ll:
322 hend = True
323 elif not hend:
324 vals = ll.split()
325 wvrv = False
326 flagv = False
327 rmsv = 0.
328 discv = 0.
329 if(len(vals)!=6):
330 casalog.post('Error parsing wvrgcal info table.line: '+ll,'WARN')
331 parsingok=False
332 else:
333 if vals[2]=='Yes':
334 wvrv=True
335 else:
336 wvrv=False
337 if vals[3]=='Yes':
338 flagv=True
339 else:
340 flagv=False
341 try:
342 rmsv = float(vals[4])
343 except:
344 casalog.post('Error parsing RMS value in info table line: '+ll,'WARN')
345 rmsv = -1.
346 parsingok=False
347 try:
348 discv = float(vals[5])
349 except:
350 casalog.post('Error parsing Disc. value in info table line: '+ll,'WARN')
351 discv = -1.
352 parsingok=False
354 namel.append(vals[1])
355 wvrl.append(wvrv)
356 flagl.append(flagv)
357 rmsl.append(rmsv)
358 discl.append(discv)
361 elif (rval==0) and (not hend) and ("Disc (um)" in ll):
362 hfound = True
363 elif 'WVR data points for antenna' in ll: # take note of antennas flagged because of too few good WVR data
364 token = ll.split('for antenna ')[1].split()
365 antennaID = int(token[0])
366 if 'All WVR' in ll:
367 flagsd[antennaID] = 0.
368 else:
369 unflagged = int(token[2])
370 total = float(token[5]) # ends in a period
371 if total>0:
372 flagsd[antennaID] = unflagged / total
373 else:
374 casalog.post('Error: zero datapoints reported for antenna id '+str(antennaID)+' in info table line: '+ll,'WARN')
375 parsingok=False
377 # end for ll
379 # create list of flagging fractions for each antenna
380 unflagfracl = list(np.ones(len(namel)))
381 for myid in flagsd.keys():
382 if myid >= len(unflagfracl):
383 casalog.post('Error: flagged antenna id '+str(myid)+' > max known antenna id '+str(len(unflagfracl)-1) ,'WARN')
384 parsingok=False
385 else:
386 unflagfracl[myid] = flagsd[myid]
389 os.system('rm -rf '+templogfile)
391 taskrval = { 'Name': namel,
392 'WVR': wvrl,
393 'Flag': flagl,
394 'Frac_unflagged': unflagfracl,
395 'RMS_um': rmsl,
396 'Disc_um': discl,
397 'rval': rval,
398 'success': False}
400 for k in range(len(namel)):
401 if(flagl[k] and rmsl[k]==0. and discl[k]==0.):
402 casalog.post('Solution for flagged antenna '+namel[k]
403 +' could not be interpolated due to insufficient number of near antennas. Was set to unity.',
404 'WARN')
406 if (rval==0) and parsingok:
407 taskrval['success'] = True
410 if(rval == 0):
411 if (smoothing>0):
412 mycb = calibrater()
413 mycb.open(filename=vis, compress=False, addcorr=False, addmodel=False)
414 mycb.smooth(tablein=caltable+'_unsmoothed', tableout=caltable,
415 smoothtype='mean', smoothtime=smoothing)
416 mycb.close()
417 return taskrval
418 else:
419 if(rval == 255):
420 casalog.post('wvr.gcal terminated with exit status '+str(rval),'SEVERE')
421 return taskrval
422 elif(rval == 134 or rval==1):
423 casalog.post('wvr.gcal terminated with exit status '+str(rval),'WARN')
424 casalog.post("No useful input data.",'SEVERE')
425 return taskrval
426 else:
427 casalog.post('wvrgcal terminated with exit status '+str(rval),'WARN')
428 return taskrval
430 except Exception as instance:
431 print('*** Error *** ', instance)
432 raise Exception from instance