import numpy as np import shutil import os from io import StringIO import math, re, sys, logging from functools import reduce import time as tm import datetime as dt import casatools msmd = casatools.msmetadata() tb = casatools.table() desc = {'ANTENNA_ID': {'comment': 'ID of antenna in this array', 'dataManagerGroup': 'StandardStMan', 'dataManagerType': 'StandardStMan', 'keywords': {}, 'maxlen': 0, 'option': 0, 'valueType': 'int'}, 'FEED_ID': {'comment': 'Feed id', 'dataManagerGroup': 'StandardStMan', 'dataManagerType': 'StandardStMan', 'keywords': {}, 'maxlen': 0, 'option': 0, 'valueType': 'int'}, 'INTERVAL': {'comment': 'Interval for which this set of parameters is accurate', 'dataManagerGroup': 'StandardStMan', 'dataManagerType': 'StandardStMan', 'keywords': {'QuantumUnits': np.array(['s'], dtype=' 9 or '.' in token: res = ('value', token) else: res = ('key', token.upper()) return res def s_newline(scanner, token): scanner.line_no += 1 return None def s_quote(scanner, token): return ('quote', token[1:-1]) def s_number(scanner, token): return ('number', float(token)) def s_angle(scanner, token): # also time l = token.split(":") ## It's not python to use reduce. ## But I neither remember nor care what is. val = reduce(lambda acc, x: 60*acc+math.copysign(acc, float(x)), l, 0.0) return ('number', val) def s_value(scanner, token): return ('value', token) def s_string(scanner, token): return ('value', token) def s_misc(scanner, token): logging.debug("Misc token %s at line %d" % (str(token), scanner.line_no)) return ('misc', token) def s_comment(scanner, token): # was only used for debugging. scanner.line_no += 1 return ('comment', token) scanner = re.Scanner([ ("\!.*\n", s_newline), ("[ \t]+", None), ("\n", s_newline), ("\r\n", s_newline), ("=", lambda s, t: ('equal',)), ("'[^'\n]*'", s_quote), ("\"[^'\n]*\"", s_quote), # Sigh. Double quotes used in freq.dat ("/", lambda s, t: ('end_chunk',)), (",", lambda s, t: ('comma',)), ("[+-]?[0-9]+:[0-9]+:[0-9]+(.[0-9]*)?", s_angle), ("[+-]?[0-9]+:[0-9]+(.[0-9]*)?", s_angle), ("[+-]?[0-9]*\.[0-9]+([Ee][+-][0-9]{1,3})?(?![A-Za-z_0-9()])", s_number), ("[+-]?[0-9]+\.?(?![A-Za-z_0-9()])", s_number), ## apparently parens and unquoted minus signs are allowed in keywords? ("[A-Za-z.0-9]([()A-Za-z_0-9._+-]+)?", s_keyword), (".*", s_misc) ]) def p_all(): global tok chunks = [] try: while True: if tok=='EOF':break chunks.append(p_chunk()) except StopIteration: pass return chunks def p_chunk(): global tok entries = [] while tok[0] != 'end_chunk': entries.append(p_item()) logging.debug("p_chunk %s", str(tok)) tok = next(tokIt) return entries def p_item(): global tok lhs = p_key() #if tok[0] == 'key': #print(tok) if tok==('equal',): logging.debug("p_item %s", str(tok)) tok = next(tokIt) rhs = p_rhs() elif tok[0] in ['value', 'quote', 'number']: rhs = p_rhs() else: rhs = True # for unitary expressions. return (lhs, rhs) def p_key(): global tok logging.debug("p_key: %s", str(tok)) if tok[0] == 'key': res = tok tok = next(tokIt) else: raise RuntimeError("Expected key token, got %s" % str(tok)) return res[1] def p_rhs(): global tok val = p_value() rhs = [val] while tok == ('comma',): logging.debug("p_rhs: %s", str(tok)) tok = next(tokIt) rhs.append(p_value()) # p_value advances tok beyond the value. if len(rhs)==1: rhs = rhs[0] return rhs def p_value(): global tok if tok[0] not in ['value', 'quote', 'number', 'key']: raise RuntimeError("Unexpected RHS token %s" % str(tok)) val = tok logging.debug("p_value: %s", str(val)) tok = next(tokIt) return val[1] scanner.line_no = 0 ################################################################################################# # data table Classes class TSysTable: def __init__(self): self.time = [] self.time_interval = [] self.source_id = [] self.antenna_no = [] self.array = [] self.freqid = [] self.tsys_1 = [] self.tsys_2 = [] self.tant = [] return class GainCurveTable: def __init__(self): self.antenna_no = [] self.array = [] self.freqid = [] self.nterm = [] self.y_typ = [] self.gain = [] self.sens_1 = [] self.sens_2 = [] return ################################################################## def read_keyfile(f): global tok, tokIt scanner.line_no = 0 try: res = scanner.scan(f.read()) #print(res[0]) if res[1]!='': raise RuntimeError("Unparsed text: %s." % (res[1][:20])) tokIt = iter(res[0]+['EOF']) try: tok = next(tokIt) res = p_all() except StopIteration: # empty file res = '' except RuntimeError as txt: #print("line %d: %s" % (scanner.line_no, txt), file=sys.stderr) raise RuntimeError return res def update_map(pols, spws, spwmap, index): idx = 0 if not isinstance(index, (list, tuple)): index = [index] pass for labels in index: for label in labels.split('|'): pol = label[0] rng = label[1:].split(':') if pol != 'X': if not pol in pols: pols.append(pol) pass if len(rng) == 1: rng.append(rng[0]) pass rng = [int(x) - 1 for x in rng] for spw in range(rng[0], rng[1] + 1): if not spw in spws: spws.append(spw) pass spwmap[(pol, spw)] = idx continue pass continue idx += 1 continue spws = sorted(spws) return def find_antenna(keys, ignore): for key in keys[1:]: if not type(key[1]) is bool: continue if key[0] in ignore: continue return key[0] return None def skip_values(infp): for line in infp: if line.startswith('!'): continue if line.strip().endswith('/'): break continue return def get_antenna_index(antenna_name, ant_names): return ant_names.index(antenna_name) def mjd_to_date(mjd): # Method used in FITSDateUtil mjd_epoch = dt.datetime(1858, 11 , 17, 0, 0, 0) mjd_int = int(mjd / 86400) delta = dt.timedelta(days=mjd_int) return mjd_epoch + delta def mjd_seconds(yy, mm, dd, d=0): if (mm < 3): yy -= 1 mm += 12 dd += d b = 0 if (yy>1582 or (yy==1582 and (mm>10 or (mm==10 and dd >= 15)))): b = np.floor(yy/100.) b = 2 - b + int(b/4) val = np.floor(365.25*yy) + np.floor(30.6001*(mm+1)) + dd - 679006.0 + b return val def get_timetuple(ts): # ts as string with these possible formats: # hh.hh # hh:mm.mm # hh:mm:ss.ss # NOTE: Regexs below will match any number of decimals on the last quantity (e.g. 19.8222222 and 19.8 both work) if re.match(r"[0-9]{2}\.[0-9]+", ts): # hh.hh tm_hour = int(ts.split('.')[0]) tm_min = math.modf(60*float(ts.split('.')[1])) tm_sec = int(60 * tm_min[0]) tm_min = int(tm_min[1]) elif re.match(r"[0-9]{2}:[0-9]{2}\.[0-9]+", ts): # hh:mm.mm tm_hour = int(ts.split(':')[0]) tm_min = math.modf(float(ts.split(':')[1])) tm_sec = int(60 * tm_min[0]) tm_min = int(tm_min[1]) elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}$", ts): # hh:mm:ss tm_hour = int(ts.split(':')[0]) tm_min = int(ts.split(':')[1]) tm_sec = int(ts.split(':')[2]) elif re.match(r"[0-9]{2}:[0-9]{2}:[0-9]{2}\.[0-9]+", ts): # hh:mm:ss.ss tm_hour = int(ts.split(':')[0]) tm_min = int(ts.split(':')[1]) tm_sec = float(ts.split(':')[2]) return tm_hour, tm_min, tm_sec def process_tsys(infp, keys, pols, data, n_band, spws, first_time, last_time, ant_names): # Get the antenna name antenna_name = find_antenna(keys[0], ['SRC/SYS']) # Skip if no antenna name found if not antenna_name: print('ANTENNA missing from TSYS group') skip_values(infp) return try: antenna = get_antenna_index(antenna_name, ant_names) except: print('Antenna %s not present in FITS-IDI file' % antenna_name) skip_values(infp) return keys = dict(keys[0]) spws = [] spwmap = {} update_map(pols, spws, spwmap, keys['INDEX']) if 'INDEX2' in keys: update_map(pols, spws, spwmap, keys['INDEX2']) pass if len(spws) != n_band: print('INDEX for antenna %s does not match FITS-IDI file' % antenna_name, file=sys.stderr) pass spws = range(n_band) timeoff = 0 if 'TIMEOFF' in keys: timeoff = float(keys['TIMEOFF']) for line in infp: if line.startswith('!'): continue fields = line.split() if len(fields) > 1: tm_yday = int(fields[0]) # Get timestamp data depending on data format obs_date = mjd_to_date(first_time) tm_year = obs_date.year idi_time = tm.strptime(f"{obs_date.year}-{obs_date.month}-{obs_date.day}","%Y-%m-%d") idi_ref = tm.mktime(idi_time) # Get the start time in sec for MJD using the method from the table filler mjd_sec = mjd_seconds(int(obs_date.year), int(obs_date.month), int(obs_date.day)) * 86400 tm_hour, tm_min, tm_sec = get_timetuple(fields[1]) days = (tm_hour*3600 + tm_min*60 + tm_sec) / 86400 t = "%dy%03dd%02dh%02dm%02ds" % \ (tm_year, tm_yday, tm_hour, tm_min, tm_sec) t = tm.mktime(tm.strptime(t, "%Yy%jd%Hh%Mm%Ss")) days = (t + timeoff - idi_ref) / 86400 curr_time = mjd_sec + (days*86400) source = 0 if (days) > ((last_time-mjd_sec)/86400) or (days) < ((first_time-mjd_sec)/86400): source = -1 #print(deg - mjd_sec) values = fields[2:] tsys = {'R': [], 'L': []} for spw in spws: for pol in ['R', 'L']: try: value = float(values[spwmap[(pol, spw)]]) if value == 999.9: value = float(-999.9) pass except: value = float(-999.9) pass tsys[pol].append(value) continue continue if source != -1: time_val = mjd_sec + (days * 86400) data.time.append(time_val) data.time_interval.append(0.0) data.antenna_no.append(antenna) data.tsys_1.append(tsys['R']) data.tsys_2.append(tsys['L']) pass pass if line.strip().endswith('/'): break continue return gain_keys = [ 'EQUAT', 'ALTAZ', 'ELEV', 'GCNRAO', 'TABLE', 'RCP', 'LCP' ] def process_gc(infp, keys, pols, data, ant_names): antenna_name = find_antenna(keys[0], gain_keys) if not antenna_name: print('Antenna missing from GAIN group') skip_values(infp) return try: antenna = get_antenna_index(antenna_name, ant_names) except: print('Antenna %s not present in FITS-IDI file' % antenna_name) skip_values(infp) return keys = dict(keys[0]) dpfu = {} try: dpfu['R'] = keys['DPFU'][0] dpfu['L'] = keys['DPFU'][1] pols.append('R') pols.append('L') except: dpfu['R'] = dpfu['L'] = keys['DPFU'] pols.append('X') pass try: value = keys['POLY'][0] except: keys['POLY'] = [keys['POLY']] pass y_typ = 0 if 'ELEV' in keys: y_typ = 1 elif 'EQUAT' in keys: y_typ = 1 elif 'ALTAZ' in keys: y_typ = 2 else: print('Unknown gain curve type for antenna %s' % antenna_name) return poly = keys['POLY'] data.antenna_no.append(antenna) data.array.append(1) data.freqid.append(1) data.y_typ.append(y_typ) data.nterm.append(len(poly)) data.gain.append(poly) data.sens_1.append(dpfu['R']) data.sens_2.append(dpfu['L']) return def antab_interp(vis, outvis, antab, ant_names, n_band, spws, gc_interval, gc_time, first_time, last_time): pols = [] data = TSysTable() data_gc = GainCurveTable() keys = StringIO() fp = open(antab, 'r') for line in fp: if line.startswith('!'): continue keys.write(line) if line.strip().endswith('/'): keys.seek(0) try: tsys = read_keyfile(keys) except RuntimeError: pass if tsys and tsys[0] and tsys[0][0][0] == 'TSYS': process_tsys(fp, tsys, pols, data, n_band, spws, first_time, last_time, ant_names) pass elif tsys and tsys[0] and tsys[0][0][0] == 'GAIN': process_gc(fp, tsys, pols, data_gc, ant_names) keys = StringIO() continue continue # Create the subtable for syscal if selected in params if append_tsys: create_syscal_subtable(data, outvis, spws) if append_gc: create_gaincurve_subtable(data_gc, outvis, spws, gc_interval, gc_time) def create_gaincurve_subtable(data, outvis, spws, gc_interval, gc_time): # Create the subtable for gain curve tb.create(f'{outvis}/GAIN_CURVE', desc_gc, dminfo=dminfo_gc) tb.putkeyword('GAIN_CURVE', f'Table: {outvis}/GAIN_CURVE') # assemble columns ANTENNA_ID = [] FEED_ID = [] SPECTRAL_WINDOW_ID = [] TIME = [] INTERVAL = [] TYPE = [] NUM_POLY = [] GAIN = [] SENSITIVITY = [[],[]] # Fill the arrays for i in range(len(data.antenna_no)): for j in spws: TIME.append(gc_time) ANTENNA_ID.append(data.antenna_no[i]) FEED_ID.append(-1) INTERVAL.append(gc_interval) TYPE.append(str(data.y_typ[i])) NUM_POLY.append(data.nterm[i]) # ASK GEORGE ABOUT THIS. SHAPE OF EACH ELEMENT VARIES # Might be an existing bug... GAIN.append(data.gain[i]) SPECTRAL_WINDOW_ID.append(int(j)) SENSITIVITY[0].append(data.sens_1[i]) SENSITIVITY[1].append(data.sens_2[i]) # Add the columns to the table tb.open(f'{outvis}/GAIN_CURVE', nomodify=False) tb.addrows(len(ANTENNA_ID)) tb.putcol('TIME', TIME) tb.putcol('ANTENNA_ID', ANTENNA_ID) tb.putcol('FEED_ID', FEED_ID) tb.putcol('INTERVAL', INTERVAL) tb.putcol('TYPE', TYPE) tb.putcol('NUM_POLY', NUM_POLY) #Fill gain row by row for i in range(len(GAIN)): # If the gain col shape doesn't match what the table filler does if len(np.shape(GAIN[i])): tb.putcell('GAIN', i, [GAIN[i],GAIN[i]]) else: tb.putcell('GAIN', i, GAIN[i]) tb.putcol('SPECTRAL_WINDOW_ID', SPECTRAL_WINDOW_ID) tb.putcol('SENSITIVITY', SENSITIVITY) #tb.flush() tb.close() def create_syscal_subtable(data, outvis, spws): # Create the subtable for syscal tb.create(f'{outvis}/SYSCAL', desc, dminfo=dminfo) tb.putkeyword('SYSCAL', f'Table: {outvis}/SYSCAL') # Assemble columns ANTENNA_ID = [] FEED_ID = [] INTERVAL = [] SPW_ID = [] TIME = [] TSYS = [[],[]] # Fill the arrays for i in range(len(data.time)): for j in spws: ANTENNA_ID.append(data.antenna_no[i]) FEED_ID.append(0) INTERVAL.append(data.time_interval[i]) SPW_ID.append(int(j)) TIME.append(data.time[i]) #TIME.append(ref_time + data.time[i]*86400) # Still need to figure out tsys... idx = int(i/len(spws)) TSYS[0].append(data.tsys_1[idx][j]) TSYS[1].append(data.tsys_2[idx][j]) # Add the columns to the table tb.open(f'{outvis}/SYSCAL', nomodify=False) tb.addrows(len(TIME)) tb.putcol('TIME', TIME) tb.putcol('ANTENNA_ID', ANTENNA_ID) tb.putcol('INTERVAL', INTERVAL) tb.putcol('SPECTRAL_WINDOW_ID', SPW_ID) tb.putcol('TSYS', TSYS) #tb.flush() tb.close()