Coverage for /home/casatest/venv/lib/python3.12/site-packages/casatasks/private/task_appendantab.py: 9%
658 statements
« prev ^ index » next coverage.py v7.10.4, created at 2025-08-21 07:43 +0000
1import numpy as np
2import shutil
3import os
4from io import StringIO
5import math, re, sys, logging
6from functools import reduce
7import time as tm
8import datetime as dt
9import casatools
10msmd = casatools.msmetadata()
11tb = casatools.table()
13desc = {'ANTENNA_ID': {'comment': 'ID of antenna in this array',
14 'dataManagerGroup': 'StandardStMan',
15 'dataManagerType': 'StandardStMan',
16 'keywords': {},
17 'maxlen': 0,
18 'option': 0,
19 'valueType': 'int'},
20 'FEED_ID': {'comment': 'Feed id',
21 'dataManagerGroup': 'StandardStMan',
22 'dataManagerType': 'StandardStMan',
23 'keywords': {},
24 'maxlen': 0,
25 'option': 0,
26 'valueType': 'int'},
27 'INTERVAL': {'comment': 'Interval for which this set of parameters is accurate',
28 'dataManagerGroup': 'StandardStMan',
29 'dataManagerType': 'StandardStMan',
30 'keywords': {'QuantumUnits': np.array(['s'], dtype='<U1')},
31 'maxlen': 0,
32 'option': 0,
33 'valueType': 'double'},
34 'SPECTRAL_WINDOW_ID': {'comment': 'ID for this spectral window setup',
35 'dataManagerGroup': 'StandardStMan',
36 'dataManagerType': 'StandardStMan',
37 'keywords': {},
38 'maxlen': 0,
39 'option': 0,
40 'valueType': 'int'},
41 'TIME': {'comment': 'Midpoint of time for which this set of parameters is accurate',
42 'dataManagerGroup': 'StandardStMan',
43 'dataManagerType': 'StandardStMan',
44 'keywords': {'MEASINFO': {'Ref': 'UTC', 'type': 'epoch'},
45 'QuantumUnits': np.array(['s'], dtype='<U1')},
46 'maxlen': 0,
47 'option': 0,
48 'valueType': 'double'},
49 'TSYS': {'comment': 'System temp. for each of the two receptors',
50 'dataManagerGroup': 'StandardStMan',
51 'dataManagerType': 'StandardStMan',
52 'keywords': {'QuantumUnits': np.array(['K'], dtype='<U1')},
53 'maxlen': 0,
54 'ndim': -1,
55 'option': 0,
56 'valueType': 'float'},
57 '_define_hypercolumn_': {},
58 '_keywords_': {},
59 '_private_keywords_': {}}
60dminfo = {'*1': {'COLUMNS': np.array(['ANTENNA_ID', 'FEED_ID', 'INTERVAL', 'SPECTRAL_WINDOW_ID', 'TIME',
61 'TSYS'], dtype='<U18'),
62 'NAME': 'StandardStMan',
63 'SEQNR': 0,
64 'SPEC': {'BUCKETSIZE': 1152,
65 'IndexLength': 7486,
66 'MaxCacheSize': 2,
67 'PERSCACHESIZE': 2},
68 'TYPE': 'StandardStMan'}}
69desc_gc = {'ANTENNA_ID': {'comment': 'Antenna identifier',
70 'dataManagerGroup': 'StandardStMan',
71 'dataManagerType': 'StandardStMan',
72 'keywords': {},
73 'maxlen': 0,
74 'option': 0,
75 'valueType': 'int'},
76 'FEED_ID': {'comment': 'Feed identifier',
77 'dataManagerGroup': 'StandardStMan',
78 'dataManagerType': 'StandardStMan',
79 'keywords': {},
80 'maxlen': 0,
81 'option': 0,
82 'valueType': 'int'},
83 'GAIN': {'comment': 'Gain polynomial',
84 'dataManagerGroup': 'StandardStMan',
85 'dataManagerType': 'StandardStMan',
86 'keywords': {},
87 'maxlen': 0,
88 'ndim': -1,
89 'option': 0,
90 'valueType': 'float'},
91 'INTERVAL': {'comment': 'Interval for which this set of parameters is accurate',
92 'dataManagerGroup': 'StandardStMan',
93 'dataManagerType': 'StandardStMan',
94 'keywords': {'QuantumUnits': np.array(['s'], dtype='<U1')},
95 'maxlen': 0,
96 'option': 0,
97 'valueType': 'double'},
98 'NUM_POLY': {'comment': 'Number of terms in polynomial',
99 'dataManagerGroup': 'StandardStMan',
100 'dataManagerType': 'StandardStMan',
101 'keywords': {},
102 'maxlen': 0,
103 'option': 0,
104 'valueType': 'int'},
105 'SENSITIVITY': {'comment': 'Antenna sensitivity',
106 'dataManagerGroup': 'StandardStMan',
107 'dataManagerType': 'StandardStMan',
108 'keywords': {'QuantumUnits': np.array(['K/Jy'], dtype='<U4')},
109 'maxlen': 0,
110 'ndim': -1,
111 'option': 0,
112 'valueType': 'float'},
113 'SPECTRAL_WINDOW_ID': {'comment': 'Spectral window identifier',
114 'dataManagerGroup': 'StandardStMan',
115 'dataManagerType': 'StandardStMan',
116 'keywords': {},
117 'maxlen': 0,
118 'option': 0,
119 'valueType': 'int'},
120 'TIME': {'comment': 'Midpoint of time for which this set of parameters is accurate',
121 'dataManagerGroup': 'StandardStMan',
122 'dataManagerType': 'StandardStMan',
123 'keywords': {'MEASINFO': {'Ref': 'UTC', 'type': 'epoch'},
124 'QuantumUnits': np.array(['s'], dtype='<U1')},
125 'maxlen': 0,
126 'option': 0,
127 'valueType': 'double'},
128 'TYPE': {'comment': 'Gain curve type',
129 'dataManagerGroup': 'StandardStMan',
130 'dataManagerType': 'StandardStMan',
131 'keywords': {},
132 'maxlen': 0,
133 'option': 0,
134 'valueType': 'string'},
135 '_define_hypercolumn_': {},
136 '_keywords_': {},
137 '_private_keywords_': {}}
138dminfo_gc = {'*1': {'COLUMNS': np.array(['ANTENNA_ID', 'FEED_ID', 'GAIN', 'INTERVAL', 'NUM_POLY',
139 'SENSITIVITY', 'SPECTRAL_WINDOW_ID', 'TIME', 'TYPE'], dtype='<U18'),
140 'NAME': 'StandardStMan',
141 'SEQNR': 0,
142 'SPEC': {'BUCKETSIZE': 1920,
143 'IndexLength': 142,
144 'MaxCacheSize': 2,
145 'PERSCACHESIZE': 2},
146 'TYPE': 'StandardStMan'}}
148#################################################################################################
149# Main task
150def appendantab(vis=None, outvis=None, antab=None,
151 overwrite=False, append_tsys=True, append_gc=True):
153 # disallow writing to the input vis
154 if vis == outvis:
155 print("Please provide a path for outvis different to the input vis")
156 return
158 # Require an outvis to be specified
159 if not outvis:
160 print('Please provide an outvis name to write to')
161 return
163 # get info from vis
164 msmd.open(vis)
165 ant_names = msmd.antennastations()
166 n_band = len(msmd.bandwidths())
167 spws = msmd.spwfordatadesc()
168 msmd.close()
170 # get start time and end time from ms
171 tb.open(vis)
172 times = tb.getcol('TIME')
173 first_time = times[0]
174 last_time = times[-1]
175 tb.close()
177 # Check if the outvis already exists
178 if os.path.exists(outvis):
179 print("File ", outvis, " already exists, change outvis name")
180 return
181 # run the antab parsing and table filling
182 shutil.copytree(vis, outvis)
183 interpreter = AntabInterp(vis, outvis, antab, ant_names, n_band, spws,
184 first_time, last_time,
185 append_tsys, append_gc, overwrite)
186 #antab_interp(vis, outvis, antab, ant_names, n_band, spws,
187 # first_time, last_time,
188 # append_tsys, append_gc, overwrite)
189 interpreter.antab_interp()
193#################################################################################################
194# scanner regex from casavlbitools
195def s_err(scanner, error):
196 raise RuntimeError("line: %d" % scanner.line_no)
198def s_keyword(scanner, token):
199 if len(token) > 9 or '.' in token:
200 res = ('value', token)
201 else:
202 res = ('key', token.upper())
203 return res
205def s_newline(scanner, token):
206 scanner.line_no += 1
207 return None
209def s_quote(scanner, token):
210 return ('quote', token[1:-1])
212def s_number(scanner, token):
213 return ('number', float(token))
215def s_angle(scanner, token): # also time
216 l = token.split(":")
217 ## It's not python to use reduce.
218 ## But I neither remember nor care what is.
219 val = reduce(lambda acc, x: 60*acc+math.copysign(acc, float(x)), l, 0.0)
220 return ('number', val)
222def s_value(scanner, token):
223 return ('value', token)
225def s_string(scanner, token):
226 return ('value', token)
228def s_misc(scanner, token):
229 logging.debug("Misc token %s at line %d" % (str(token), scanner.line_no))
230 return ('misc', token)
232def s_comment(scanner, token): # was only used for debugging.
233 scanner.line_no += 1
234 return ('comment', token)
236scanner = re.Scanner([
237 ("\!.*\n", s_newline),
238 ("[ \t]+", None),
239 ("\n", s_newline),
240 ("\r\n", s_newline),
241 ("=", lambda s, t: ('equal',)),
242 ("'[^'\n]*'", s_quote),
243 ("\"[^'\n]*\"", s_quote), # Sigh. Double quotes used in freq.dat
244 ("/", lambda s, t: ('end_chunk',)),
245 (",", lambda s, t: ('comma',)),
246 ("[+-]?[0-9]+:[0-9]+:[0-9]+(.[0-9]*)?", s_angle),
247 ("[+-]?[0-9]+:[0-9]+(.[0-9]*)?", s_angle),
248 ("[+-]?[0-9]*\.[0-9]+([Ee][+-][0-9]{1,3})?(?![A-Za-z_0-9()])", s_number),
249 ("[+-]?[0-9]+\.?(?![A-Za-z_0-9()])", s_number),
250 ## apparently parens and unquoted minus signs are allowed in keywords?
251 ("[A-Za-z.0-9]([()A-Za-z_0-9._+-]+)?", s_keyword),
252 (".*", s_misc)
253])
255def p_all():
256 global tok
257 chunks = []
258 try:
259 while True:
260 if tok=='EOF':break
261 chunks.append(p_chunk())
262 except StopIteration:
263 pass
264 return chunks
266def p_chunk():
267 global tok
268 entries = []
269 while tok[0] != 'end_chunk':
270 entries.append(p_item())
271 logging.debug("p_chunk %s", str(tok))
272 tok = next(tokIt)
273 return entries
275def p_item():
276 global tok
277 lhs = p_key()
279 if tok==('equal',):
280 logging.debug("p_item %s", str(tok))
281 tok = next(tokIt)
282 rhs = p_rhs()
283 elif tok[0] in ['value', 'quote', 'number']:
284 rhs = p_rhs()
285 else:
286 rhs = True # for unitary expressions.
287 return (lhs, rhs)
289def p_key():
290 global tok
291 logging.debug("p_key: %s", str(tok))
292 if tok[0] == 'key':
293 res = tok
294 tok = next(tokIt)
295 else:
296 raise RuntimeError("Expected key token, got %s" % str(tok))
297 return res[1]
299def p_rhs():
300 global tok
301 val = p_value()
302 rhs = [val]
303 while tok == ('comma',):
304 logging.debug("p_rhs: %s", str(tok))
305 tok = next(tokIt)
306 rhs.append(p_value()) # p_value advances tok beyond the value.
307 if len(rhs)==1:
308 rhs = rhs[0]
309 return rhs
311def p_value():
312 global tok
313 if tok[0] not in ['value', 'quote', 'number', 'key']:
314 raise RuntimeError("Unexpected RHS token %s" % str(tok))
315 val = tok
316 logging.debug("p_value: %s", str(val))
317 tok = next(tokIt)
318 return val[1]
320scanner.line_no = 0
321#################################################################################################
322# data table Classes
324class TSysTable:
325 def __init__(self):
326 self.time = []
327 self.time_interval = []
328 self.source_id = []
329 self.antenna_no = []
330 self.array = []
331 self.freqid = []
332 self.tsys_1 = []
333 self.tsys_2 = []
334 self.tant = []
335 return
337class GainCurveTable:
338 def __init__(self):
339 self.antenna_no = []
340 self.array = []
341 self.freqid = []
342 self.spwid = []
343 self.nterm = []
344 self.y_typ = []
345 self.gain = []
346 self.sens_1 = []
347 self.sens_2 = []
348 return
350##################################################################
353def read_keyfile(f):
354 global tok, tokIt
355 scanner.line_no = 0
357 try:
358 res = scanner.scan(f.read())
359 if res[1]!='':
360 raise RuntimeError("Unparsed text: %s." % (res[1][:20]))
361 tokIt = iter(res[0]+['EOF'])
362 try:
363 tok = next(tokIt)
364 res = p_all()
365 except StopIteration: # empty file
366 res = ''
367 except RuntimeError as txt:
368 #print("line %d: %s" % (scanner.line_no, txt), file=sys.stderr)
369 raise RuntimeError
370 return res
372def update_map(pols, spws, spwmap, index):
373 idx = 0
374 if not isinstance(index, (list, tuple)):
375 index = [index]
376 pass
377 for labels in index:
378 for label in labels.split('|'):
379 pol = label[0]
380 rng = label[1:].split(':')
381 if pol != 'X':
382 if not pol in pols:
383 pols.append(pol)
384 pass
385 if len(rng) == 1:
386 rng.append(rng[0])
387 pass
388 rng = [int(x) - 1 for x in rng]
389 for spw in range(rng[0], rng[1] + 1):
390 if not spw in spws:
391 spws.append(spw)
392 pass
393 spwmap[(pol, spw)] = idx
394 continue
395 pass
396 continue
397 idx += 1
398 continue
399 spws = sorted(spws)
400 return
402def find_antenna(keys, ignore):
403 for key in keys[1:]:
404 if not type(key[1]) is bool:
405 continue
406 if key[0] in ignore:
407 continue
408 return key[0]
409 return None
411def skip_values(infp):
412 for line in infp:
413 if line.startswith('!'):
414 continue
415 if line.strip().endswith('/'):
416 break
417 continue
418 return
420def get_antenna_index(antenna_name, ant_names):
421 return ant_names.index(antenna_name)
423def mjd_to_date(mjd):
424 # Method used in FITSDateUtil
425 mjd_epoch = dt.datetime(1858, 11 , 17, 0, 0, 0)
427 mjd_int = int(mjd / 86400)
428 delta = dt.timedelta(days=mjd_int)
429 return mjd_epoch + delta
431def mjd_seconds(yy, mm, dd, d=0):
432 if (mm < 3):
433 yy -= 1
434 mm += 12
436 dd += d
437 b = 0
439 if (yy>1582 or (yy==1582 and (mm>10 or (mm==10 and dd >= 15)))):
440 b = np.floor(yy/100.)
441 b = 2 - b + int(b/4)
443 val = np.floor(365.25*yy) + np.floor(30.6001*(mm+1)) + dd - 679006.0 + b
444 return val
447def get_timetuple(ts):
448 # ts as string with these possible formats:
449 # hh.hh
450 # hh:mm.mm
451 # hh:mm:ss.ss
452 # NOTE: Regexs below will match any number of decimals on the last quantity (e.g. 19.8222222 and 19.8 both work)
453 if re.match(r"[0-9]{2}\.[0-9]+", ts):
454 # hh.hh
455 tm_hour = int(ts.split('.')[0])
456 tm_min = math.modf(60*float(ts.split('.')[1]))
457 tm_sec = int(60 * tm_min[0])
458 tm_min = int(tm_min[1])
459 elif re.match(r"[0-9]{2}:[0-9]{2}\.[0-9]+", ts):
460 # hh:mm.mm
461 tm_hour = int(ts.split(':')[0])
462 tm_min = math.modf(float(ts.split(':')[1]))
463 tm_sec = int(60 * tm_min[0])
464 tm_min = int(tm_min[1])
465 elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}$", ts):
466 # hh:mm:ss
467 tm_hour = int(ts.split(':')[0])
468 tm_min = int(ts.split(':')[1])
469 tm_sec = int(ts.split(':')[2])
470 elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}\.[0-9]+", ts):
471 # hh:mm:ss.ss
472 tm_hour = int(ts.split(':')[0])
473 tm_min = int(ts.split(':')[1])
474 tm_sec = float(ts.split(':')[2])
475 return tm_hour, tm_min, tm_sec
477############## interp class ###################
479class AntabInterp:
480 def __init__(self, vis, outvis, antab,
481 ant_names, n_band, spws,
482 first_time, last_time,
483 append_tsys, append_gc,
484 overwrite):
486 # Passed to class when constructing
487 self.vis = vis
488 self.outvis = outvis
489 self.antab = antab
490 self.ant_names = ant_names
491 self.n_band = n_band
492 self.first_time = first_time
493 self.last_time = last_time
494 self.spws = spws
495 self.append_tsys = append_tsys
496 self.append_gc = append_gc
497 self.overwrite = overwrite
499 # used during procesing
500 self.keys = None
501 self.data = None
502 self.data_gc = None
503 self.pols = None
504 self.spw_ranges = {}
505 self.gain_keys = [ 'EQUAT', 'ALTAZ', 'ELEV', 'GCNRAO', 'TABLE', 'RCP', 'LCP' ]
506 self.replacements = dict()
509 # Core antab interp function
510 def antab_interp(self):
511 self.pols = []
512 self.data = TSysTable()
513 self.data_gc = GainCurveTable()
514 keys = StringIO()
515 fp = open(self.antab, 'r')
517 for line in fp:
518 #print(line)
519 if line.startswith('!'):
520 continue
521 keys.write(line)
522 if line.strip().endswith('/'):
523 keys.seek(0)
524 try:
525 tsys = read_keyfile(keys)
526 #print('TSYS: ', tsys)
527 except RuntimeError:
528 print("\n", keys.getvalue(), file=sys.stderr)
529 raise RuntimeError('error parsing ANTAB file')
530 if tsys and tsys[0] and tsys[0][0][0] == 'TSYS':
531 self.process_tsys(fp, tsys)
532 pass
533 elif tsys and tsys[0] and tsys[0][0][0] == 'GAIN':
534 #print('ENTER PROCESS GC')
535 self.process_gc(fp, tsys)
536 keys = StringIO()
537 continue
538 continue
540 # Create the subtable for syscal if selected in params
541 if self.append_tsys:
542 self.create_syscal_subtable()
543 if self.append_gc:
544 self.create_gaincurve_subtable()
546 def get_spw_ranges(self):
547 # Get ranges for each spectral window
548 try:
549 tb.open(f'{self.vis}/SPECTRAL_WINDOW')
550 for spw in self.spws:
551 chan_freq = tb.getcol('CHAN_FREQ')[:,spw]
552 min_freq = min(chan_freq) / (1e6)
553 max_freq = max(chan_freq) / (1e6)
554 self.spw_ranges[spw] = (min_freq, max_freq)
555 tb.close()
556 except Exception as e:
557 print("ERROR FAILED TO GET SPW FREQ RANGES")
559 def check_which_spw(self, freq):
560 in_spw = []
561 for spw, ranges in self.spw_ranges.items():
562 if ranges[0] >= freq[0] and ranges[1] <= freq[1]:
563 in_spw.append(spw)
564 # If there is no spw that this belongs to default to what we were setting this value to before
565 if not in_spw:
566 print("WARNING: NO CORRESPONDING SPW FOUND FOR FREQ {}".format(freq))
567 return [-1]
568 else:
569 return in_spw
571 def process_tsys(self, infp, keys):
572 # Get the antenna name
573 antenna_name = find_antenna(keys[0], ['SRC/SYS'])
574 # Skip if no antenna name found
575 if not antenna_name:
576 print('ANTENNA missing from TSYS group')
577 skip_values(infp)
578 return
579 try:
580 antenna = get_antenna_index(antenna_name, self.ant_names)
581 except:
582 print('Antenna {0} not present in the Measurement Set. {0} values will be ignored.'.format(antenna_name))
583 skip_values(infp)
584 return
586 try:
587 self.get_spw_ranges()
588 spw_ids = self.check_which_spw(keys['FREQ'])
589 except Exception as e:
590 spw_ids = [-1]
591 print('ANTAB HAS NO FREQ INFORMATION DEFAULTING TO FREQ ID -1')
593 keys = dict(keys[0])
594 spws = []
595 spwmap = {}
597 if 'INDEX' in keys:
598 update_map(self.pols, spws, spwmap, keys['INDEX'])
599 if 'INDEX2' in keys:
600 update_map(self.pols, spws, spwmap, keys['INDEX2'])
601 pass
602 if len(spws) != self.n_band:
603 print('INDEX for antenna %s does not match FITS-IDI file'
604 % antenna_name, file=sys.stderr)
605 pass
607 spws = range(self.n_band)
608 timeoff = 0
609 if 'TIMEOFF' in keys:
610 timeoff = float(keys['TIMEOFF'])
612 for line in infp:
613 if line.startswith('!'):
614 continue
616 fields = line.split()
618 if len(fields) > 1:
619 tm_yday = int(fields[0])
620 # Get timestamp data depending on data format
621 obs_date = mjd_to_date(self.first_time)
622 tm_year = obs_date.year
623 idi_time = tm.strptime(f"{obs_date.year}-{obs_date.month}-{obs_date.day}","%Y-%m-%d")
624 idi_ref = tm.mktime(idi_time)
626 # Get the start time in sec for MJD using the method from the table filler
627 mjd_sec = mjd_seconds(int(obs_date.year), int(obs_date.month), int(obs_date.day)) * 86400
629 tm_hour, tm_min, tm_sec = get_timetuple(fields[1])
630 days = (tm_hour*3600 + tm_min*60 + tm_sec) / 86400
631 t = "%dy%03dd%02dh%02dm%02ds" % \
632 (tm_year, tm_yday, tm_hour, tm_min, tm_sec)
633 t = tm.mktime(tm.strptime(t, "%Yy%jd%Hh%Mm%Ss"))
634 days = (t + timeoff - idi_ref) / 86400
636 curr_time = mjd_sec + (days*86400)
638 source = 0
639 if (days) > ((self.last_time-mjd_sec)/86400) or (days) < ((self.first_time-mjd_sec)/86400):
640 source = -1
642 #print(deg - mjd_sec)
643 values = fields[2:]
644 tsys = {'R': [], 'L': []}
645 for spw in spws:
646 for pol in ['R', 'L']:
647 try:
648 value = float(values[spwmap[(pol, spw)]])
649 if value == 999.9:
650 value = float(-999.9)
651 pass
652 except:
653 value = float(-999.9)
654 pass
655 tsys[pol].append(value)
656 continue
657 continue
658 if source != -1:
659 time_val = mjd_sec + (days * 86400)
660 self.data.time.append(time_val)
661 self.data.time_interval.append(0.0)
662 self.data.antenna_no.append(antenna)
663 self.data.tsys_1.append(tsys['R'])
664 self.data.tsys_2.append(tsys['L'])
665 self.data.freqid.append(spw_ids)
667 pass
668 pass
670 if line.strip().endswith('/'):
671 break
672 continue
673 return
675 def process_gc(self, infp, keys):
676 antenna_name = find_antenna(keys[0], self.gain_keys)
677 if not antenna_name:
678 print('Antenna missing from GAIN group')
679 skip_values(infp)
680 return
681 try:
682 antenna = get_antenna_index(antenna_name, self.ant_names)
683 except:
684 print('Antenna {0} not present in the Measurement Set. {0} values will be ignored.'.format(antenna_name))
685 skip_values(infp)
686 return
687 keys = dict(keys[0])
689 try:
690 freq = keys['FREQ']
691 self.get_spw_ranges()
692 spw_ids = self.check_which_spw(keys['FREQ'])
693 except Exception as e:
694 spw_ids = [-1]
695 freq = [0, 0]
696 print('ANTAB HAS NO FREQ INFORMATION DEFAULTING TO SPW ID -1 and freq to [0, 0]')
698 dpfu = {}
699 try:
700 dpfu['R'] = keys['DPFU'][0]
701 dpfu['L'] = keys['DPFU'][1]
702 self.pols.append('R')
703 self.pols.append('L')
704 except:
705 dpfu['R'] = dpfu['L'] = keys['DPFU']
706 self.pols.append('X')
707 pass
708 try:
709 value = keys['POLY'][0]
710 except:
711 keys['POLY'] = [keys['POLY']]
712 pass
714 y_typ = 0
715 if 'ELEV' in keys:
716 y_typ = 1
717 elif 'EQUAT' in keys:
718 y_typ = 1
719 elif 'ALTAZ' in keys:
720 y_typ = 2
721 else:
722 print('Unknown gain curve type for antenna %s' % antenna_name)
723 return
725 poly = keys['POLY']
726 self.data_gc.antenna_no.append(antenna)
727 self.data_gc.array.append(1)
728 self.data_gc.freqid.append(freq)
729 self.data_gc.spwid.append(spw_ids)
730 self.data_gc.y_typ.append(y_typ)
731 self.data_gc.nterm.append(len(poly))
732 self.data_gc.gain.append(poly)
733 self.data_gc.sens_1.append(dpfu['R'])
734 self.data_gc.sens_2.append(dpfu['L'])
735 return
737 def create_existing_data_table(self):
738 existing_data = dict()
739 self.replacements = dict()
740 if os.path.exists(f'{self.outvis}/GAIN_CURVE'):
741 tb.open(f'{self.outvis}/GAIN_CURVE')
742 antennas = tb.getcol('ANTENNA_ID')
743 times = tb.getcol('TIME')
744 spws = tb.getcol('SPECTRAL_WINDOW_ID')
746 # add a set of times for each antenna
747 # add a list of rows and spws for each antenna
748 for i in range(len(antennas)):
749 if antennas[i] not in existing_data:
750 existing_data[antennas[i]] = dict()
751 existing_data[antennas[i]]['times'] = set()
752 existing_data[antennas[i]]['rows'] = []
753 existing_data[antennas[i]]['spws'] = []
755 existing_data[antennas[i]]['times'].add(times[i])
756 existing_data[antennas[i]]['rows'].append(i)
757 existing_data[antennas[i]]['spws'].append(spws[i])
758 else:
759 existing_data[antennas[i]]['times'].add(times[i])
760 existing_data[antennas[i]]['rows'].append(i)
761 existing_data[antennas[i]]['spws'].append(spws[i])
763 tb.close()
765 return existing_data
767 def create_gaincurve_subtable(self):
768 # If no subtable exists add one
769 if not os.path.exists(f'{self.outvis}/GAIN_CURVE'):
770 # Create the subtable for gain curve
771 tb.create(f'{self.outvis}/GAIN_CURVE', desc_gc, dminfo=dminfo_gc)
772 tb.putkeyword('GAIN_CURVE', f'Table: {self.outvis}/GAIN_CURVE')
774 # assemble columns
775 ANTENNA_ID = []
776 FEED_ID = []
777 SPECTRAL_WINDOW_ID = []
778 TIME = []
779 INTERVAL = []
780 TYPE = []
781 NUM_POLY = []
782 GAIN = []
783 SENSITIVITY = [[],[]]
785 # If a subtable exists then read contents into assembed data
786 # Create data structure for antenna an time stamps to ignore if not overwriting
787 existing_data = self.create_existing_data_table()
789 # try getting channel freq information
790 tb.open(f'{self.vis}/SPECTRAL_WINDOW')
791 chan_freqs = tb.getcol('CHAN_FREQ')
792 tb.close()
794 # Get interval information from the observation table instead of meta-data
795 tb.open(f'{self.outvis}/OBSERVATION')
796 time_range = tb.getcol("TIME_RANGE")
797 gc_interval = (time_range[1] - time_range[0])[0]
798 gc_time = np.mean(time_range)
799 tb.close()
801 # New Array filler
802 for row in range(len(self.data_gc.gain)):
803 # iterate over spw range for non-overlapping case
804 for j in self.spws:
805 # get antenna and freqency range from antab for this row
806 ant = self.data_gc.antenna_no[row]
807 freq_range = self.data_gc.freqid[row]
809 # If there is already data for that antenna and time and overwrite = False then skip
810 if ant in existing_data and gc_time in existing_data[ant]['times'] and not self.overwrite:
811 print("antenna: {} has data present for the time {}. The antab data for this value will be ignored".format(ant, gc_time))
812 continue
814 # If there is data present and overwrite is true then fill the replacement table
815 elif ant in existing_data and gc_time in existing_data[ant]['times'] and self.overwrite:
816 # Get all the relevant data to put in replacements
817 rows = existing_data[ant]['rows']
818 spws = existing_data[ant]['spws']
819 # For each row and get chan freq of corrisponding spw
820 for i in range(len(rows)):
821 # r is the row that the data exists in the original GAIN_CURVE table
822 r = rows[i]
823 spw = spws[i]
824 freq = chan_freqs[row, spw] / (1e6)
825 # get the sensitivity from the proper freqency range
826 if freq >= freq_range[0] and freq <= freq_range[1]:
827 sens = [self.data_gc.sens_1[row], self.data_gc.sens_2[row]]
828 else:
829 continue
831 # Make an entry for row r in the replacments table
832 # r is the row in the gaincurve table that will be replaced
833 # row is the row in the antab where the data is present
834 self.replacements[r] = {'ANTENNA_ID':ant, 'FEED_ID':-1, 'INTERVAL':gc_interval,
835 'TYPE':self.data_gc.y_typ[row], 'NUM_POLY':self.data_gc.nterm[row], 'GAIN':self.data_gc.gain[row],
836 'TIME':gc_time, 'SPECTRAL_WINDOW_ID':spw, 'SENSITIVITY':sens}
837 continue
839 # If the data doesn't overlap with existing data append as normal
840 else:
841 TIME.append(gc_time)
842 ANTENNA_ID.append(self.data_gc.antenna_no[row])
843 FEED_ID.append(-1)
844 INTERVAL.append(gc_interval)
845 TYPE.append(str(self.data_gc.y_typ[row]))
846 NUM_POLY.append(self.data_gc.nterm[row])
847 GAIN.append(self.data_gc.gain[row])
848 SPECTRAL_WINDOW_ID.append(int(j))
849 SENSITIVITY[0].append(self.data_gc.sens_1[row])
850 SENSITIVITY[1].append(self.data_gc.sens_2[row])
852 # If appending get existing cols and append values to it
853 tb.open(f'{self.outvis}/GAIN_CURVE', nomodify=False)
854 if not self.overwrite and tb.nrows() > 0:
855 #tb.addrows(len(TIME))
856 TIME = np.append(tb.getcol('TIME'), TIME)
857 ANTENNA_ID = np.append(tb.getcol('ANTENNA_ID'), ANTENNA_ID)
858 FEED_ID = np.append(tb.getcol('FEED_ID'), FEED_ID)
859 INTERVAL = np.append(tb.getcol('INTERVAL'), INTERVAL)
860 TYPE = np.append(tb.getcol('TYPE'), TYPE)
861 NUM_POLY = np.append(tb.getcol('NUM_POLY'), NUM_POLY)
863 # Can't call getcol on Gain since dimensions can vary
864 tmp_GAIN = []
865 # TODO ISSUE: Should only apply the gains for certain frequency ranges?
866 for row in range(tb.nrows()):
867 tmp_GAIN.append(tb.getcell('GAIN', row))
868 # re-assembled the old and new gain col
869 GAIN = tmp_GAIN + GAIN
870 SPECTRAL_WINDOW_ID = np.append(tb.getcol('SPECTRAL_WINDOW_ID'), SPECTRAL_WINDOW_ID)
871 tmp_0 = np.append(tb.getcol('SENSITIVITY')[0], SENSITIVITY[0])
872 tmp_1 = np.append(tb.getcol('SENSITIVITY')[1], SENSITIVITY[1])
873 SENSITIVITY = np.asarray([tmp_0, tmp_1])
875 # If overwriting then go over the existing rows to overwrite
876 if self.overwrite:
877 TIME = np.append(tb.getcol('TIME'), TIME)
878 ANTENNA_ID = np.append(tb.getcol('ANTENNA_ID'), ANTENNA_ID)
879 FEED_ID = np.append(tb.getcol('FEED_ID'), FEED_ID)
880 INTERVAL = np.append(tb.getcol('INTERVAL'), INTERVAL)
881 TYPE = np.append(tb.getcol('TYPE'), TYPE)
882 NUM_POLY = np.append(tb.getcol('NUM_POLY'), NUM_POLY)
883 # Can't call getcol on Gain since dimensions can vary
884 tmp_GAIN = []
885 # TODO ISSUE: Should only apply the gains for certain frequency ranges?
886 for row in range(tb.nrows()):
887 tmp_GAIN.append(tb.getcell('GAIN', row))
888 # re-assembled the old and new gain col
889 GAIN = tmp_GAIN + GAIN
890 SPECTRAL_WINDOW_ID = np.append(tb.getcol('SPECTRAL_WINDOW_ID'), SPECTRAL_WINDOW_ID)
891 tmp_0 = np.append(tb.getcol('SENSITIVITY')[0], SENSITIVITY[0])
892 tmp_1 = np.append(tb.getcol('SENSITIVITY')[1], SENSITIVITY[1])
893 SENSITIVITY = np.asarray([tmp_0, tmp_1])
895 for row, data in self.replacements.items():
896 TIME[row] = data['TIME']
897 ANTENNA_ID[row] = data['ANTENNA_ID']
898 FEED_ID[row] = data['FEED_ID']
899 INTERVAL[row] = data['INTERVAL']
900 TYPE[row] = data['TYPE']
901 NUM_POLY[row] = data['NUM_POLY']
902 GAIN[row] = data['GAIN']
903 SPECTRAL_WINDOW_ID[row] = data['SPECTRAL_WINDOW_ID']
904 SENSITIVITY[:, row] = [data['SENSITIVITY'][0], data['SENSITIVITY'][1]]
906 tb.close()
908 # Add the columns to the table
909 tb.open(f'{self.outvis}/GAIN_CURVE', nomodify=False)
910 to_add = len(TIME) - tb.nrows()
911 tb.addrows(to_add)
912 tb.putcol('TIME', TIME)
913 tb.putcol('ANTENNA_ID', ANTENNA_ID)
914 tb.putcol('FEED_ID', FEED_ID)
915 tb.putcol('SPECTRAL_WINDOW_ID', SPECTRAL_WINDOW_ID)
916 tb.putcol('INTERVAL', INTERVAL)
917 # Type 1 POWER(EL) Type 2 POWER(ZA)
918 tb.putcol('NUM_POLY', NUM_POLY)
920 #Fill gain row by row
921 for i in range(len(GAIN)):
922 # If the gain col shape doesn't match what the table filler does
923 if len(np.shape(GAIN[i])) < 2:
924 tb.putcell('GAIN', i, [GAIN[i],GAIN[i]])
925 else:
926 tb.putcell('GAIN', i, GAIN[i])
928 # Convert Type to proper string
929 if TYPE[i] == '1':
930 TYPE[i] = 'POWER(EL)'
931 elif TYPE[i] == '2':
932 TYPE[i] = 'POWER(ZA)'
934 tb.putcol('TYPE', TYPE)
935 tb.putcol('SENSITIVITY', SENSITIVITY)
936 #tb.flush()
937 tb.close()
940 def create_syscal_subtable(self):
941 # If no subtable exists
942 if not os.path.exists(f'{self.outvis}/SYSCAL'):
943 # Create the subtable for syscal
944 tb.create(f'{self.outvis}/SYSCAL', desc, dminfo=dminfo)
945 tb.putkeyword('SYSCAL', f'Table: {self.outvis}/SYSCAL')
947 # Assemble columns
948 ANTENNA_ID = []
949 FEED_ID = []
950 INTERVAL = []
951 SPW_ID = []
952 TIME = []
953 TSYS = [[],[]]
954 self.data.tsys_1 = np.asarray(self.data.tsys_1)
955 self.data.tsys_2 = np.asarray(self.data.tsys_2)
957 # If a subtable exists then read contents into assembed data
958 # Create data structure for antenna an time stamps to ignore if not overwriting
959 existing_data = self.create_existing_data_table()
961 # Fill the arrays
962 to_skip = []
963 self.replacements = dict()
964 for i in range(len(self.data.time)):
965 # if the antenna has been visited already with the same time stamp then ignore
966 if self.data.antenna_no[i] in existing_data and self.data.time[i] in existing_data[self.data.antenna_no[i]]['times'] and not self.overwrite:
967 # This data is already in the subtable so ignore
968 print("antenna: {} has data present for the time {}. The antab data for this value will be ignored".format(self.data.antenna_no[i], self.data.time[i]))
969 to_skip.append(True)
970 continue
971 else:
972 to_skip.append(False)
974 for j in self.spws:
975 # get antenna from antab for this row
976 ant = self.data.antenna_no[i]
978 # Check that this antenna is in this spectral window
979 if j not in self.data.freqid[i] and self.data.freqid[i] != [-1]:
980 continue
981 if ant in existing_data and self.data.time[i] in existing_data[ant]['times'] and self.overwrite:
982 # Get the row and information to overwrite the row from the ms
983 # Get all the relevant data to put in replacements
984 rows = existing_data[ant]['rows']
985 spws = existing_data[ant]['spws']
987 for r in range(len(rows)):
988 row = rows[r]
989 spw = spws[r]
991 self.replacements[row] = {'ANTENNA':self.data.antenna_no[i], 'FEED_ID':0, 'INTERVAL':self.data.time_interval[i],
992 'SPW_ID':spw, 'TIME':self.data.time[i], 'TSYS':[[],[]]}
993 print(row, spw, int(j))
994 #row = existing_data[ant]['row']
995 to_skip.append(True)
997 #self.replacements[row] = {'ANTENNA':self.data.antenna_no[i], 'FEED_ID':0, 'INTERVAL':self.data.time_interval[i],
998 # 'SPW_ID':int(j), 'TIME':self.data.time[i], 'TSYS':[[],[]]}
1000 continue
1002 ANTENNA_ID.append(self.data.antenna_no[i])
1003 FEED_ID.append(0)
1004 INTERVAL.append(self.data.time_interval[i])
1005 SPW_ID.append(int(j))
1006 TIME.append(self.data.time[i])
1008 for i in range(len(self.data.tsys_1)):
1009 if to_skip[i] and not self.overwrite:
1010 continue
1011 elif to_skip[i] and self.overwrite:
1012 self.replacements[i]['TSYS'][0].extend(self.data.tsys_1[i])
1013 self.replacements[i]['TSYS'][1].extend(self.data.tsys_2[i])
1014 continue
1015 TSYS[0].extend(self.data.tsys_1[i])
1016 TSYS[1].extend(self.data.tsys_2[i])
1018 # If appending get existing cols and append values to it also check if empty
1019 tb.open(f'{self.outvis}/SYSCAL', nomodify=False)
1020 if not self.overwrite and tb.nrows() > 0:
1021 TIME = np.append(tb.getcol('TIME'), TIME)
1022 ANTENNA_ID = np.append(tb.getcol('ANTENNA_ID'), ANTENNA_ID)
1023 INTERVAL = np.append(tb.getcol('INTERVAL'), INTERVAL)
1024 SPW_ID = np.append(tb.getcol('SPECTRAL_WINDOW_ID'), SPW_ID)
1025 tmp_0 = np.append(tb.getcol('TSYS')[0], TSYS[0])
1026 tmp_1 = np.append(tb.getcol('TSYS')[1], TSYS[1])
1027 TSYS = np.asarray([tmp_0, tmp_1])
1029 # If overwriting then go over the existing rows to overwrite
1030 if self.overwrite:
1031 for row, data in self.replacements.items():
1032 TIME[row] = data['TIME']
1033 ANTENNA_ID[row] = data['ANTENNA']
1034 INTERVAL[row] = data['INTERVAL']
1035 SPW_ID[row] = data['SPW_ID']
1036 #print(TSYS[0][row], TSYS[1][row])
1037 TSYS[0][row] = data['TSYS'][0]
1038 TSYS[1][row] = data['TSYS'][1]
1039 tb.close()
1041 # Add the columns to the table
1042 tb.open(f'{self.outvis}/SYSCAL', nomodify=False)
1043 to_add = len(TIME) - tb.nrows()
1044 tb.addrows(to_add)
1045 tb.putcol('TIME', TIME)
1046 tb.putcol('ANTENNA_ID', ANTENNA_ID)
1047 tb.putcol('INTERVAL', INTERVAL)
1048 tb.putcol('SPECTRAL_WINDOW_ID', SPW_ID)
1049 tb.putcol('TSYS', TSYS)
1050 #tb.flush()
1051 tb.close()