diff --git a/README.md b/README.md index 3b0ddfc..b8e9c51 100644 --- a/README.md +++ b/README.md @@ -281,7 +281,7 @@ Detailed instructions here: + Bonus RFI Workshop content --> + +#### `1060550888` channel 143 - Tile104 Narrow Swoosh + +```bash +python demo/04_ssins.py demo/data/1060550888/raw/1060550888{.metafits,_20130814212851_gpubox12_01.fits} --no-diff --sigchain --suffix '.ch143' --sel-ants Tile104 --sel-pols yy +``` + +![Tile104 Narrow Swoosh](demo/data/1060550888/raw/1060550888.auto.ch143.Tile104.yy.sigchain.png) + +#### `1087596040` channel 134 - Snake + +```bash +python demo/04_ssins.py demo/data/1087596040/raw/1087596040{.metafits,_20140623220027_gpubox21_00.fits} --no-diff --crosses --suffix '.ch134' --sel-pols yy +``` + +![ch134 Snake](demo/data/1087596040/raw/1087596040.cross.ch134.yy.spectrum.png) + +#### `1087596040` channel 143 - Tile108 Swoosh + +```bash +python demo/04_ssins.py demo/data/1087596040/raw/1087596040{.metafits,_20140623220027_gpubox12_00.fits} --suffix '.ch143' --sel-ants Tile108 --sel-pols xx +``` + +![ch143 Tile108 Swoosh](demo/data/1087596040/raw/1087596040.diff.auto.ch143.Tile108.xx.spectrum.png) + +#### `1088806248` channel 133 - Snake + +```bash +python demo/04_ssins.py demo/data/1088806248/raw/1088806248{.metafits,_20140707221035_gpubox22_00.fits} --crosses --suffix '.ch133' --sel-pols xx +``` + +![ch133 Snake](demo/data/1088806248/raw/1088806248.diff.cross.ch133.xx.spectrum.png) + +#### `1088806248` channel 145 - Tile055 Discontinuity + +```bash +python demo/04_ssins.py demo/data/1088806248/raw/1088806248{.metafits,_20140707221035_gpubox10_00.fits} --no-diff --sigchain --suffix '.ch145' --sel-ants Tile055 --sel-pols yy +``` + +![ch145 narrowband discontinuity on Tile55 yy](demo/data/1088806248/raw/1088806248.auto.ch145.Tile055.yy.sigchain.png) + +#### `1089238040` channel 133 - Snake + +```bash +python demo/04_ssins.py demo/data/1089238040/raw/1089238040{.metafits,_20140712220707_gpubox22_00.fits} --no-diff --suffix '.ch133' --sel-pols xx +``` + +![ch133 Snake](demo/data/1089238040/raw/1089238040.auto.ch133.xx.spectrum.png) + +#### `1090871744` channel 137 - Slow-moving TV + +```bash +python demo/04_ssins.py demo/data/1090871744/raw/1090871744{.metafits,_20140731195531_gpubox18_00.fits} --no-diff --crosses --suffix '.ch137' --sel-pols xx +``` + +![ch137 Slow TV](demo/data/1090871744/raw/1090871744.cross.ch137.xx.spectrum.png) + +#### `1094319712` channel 142 - Tile108 Swoosh + +```bash +python demo/04_ssins.py demo/data/1094319712/raw/1094319712{.metafits,_20140909174139_gpubox13_00.fits} --no-diff --sigchain --suffix '.ch142' --sel-pols xx +``` + +![ch142 Tile108 Swoosh](demo/data/1094319712/raw/1094319712.auto.ch142.xx.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1094319712/raw/1094319712{.metafits,_20140909174139_gpubox13_00.fits} --no-diff --crosses --suffix '.ch142' --skip-ants Tile108 --sel-pols xx +``` + +![ch142 without Tile108](demo/data/1094319712/raw/1094319712.cross.ch142.noTile108.xx.spectrum.png) + +#### `1252516448` channel 142 - Slow-moving TV + +```bash +python demo/04_ssins.py demo/data/1252516448/raw/1252516448{.metafits,_20190914171353_gpubox13_00.fits} --no-diff --sigchain --suffix '.ch142' --sel-pols xx +``` + +![ch142 slow TV](demo/data/1252516448/raw/1252516448.auto.ch142.xx.sigchain.png) + +#### `1255099440` channel 135 - Satellite Reflection + +```bash +python demo/04_ssins.py demo/data/1255099440/raw/1255099440{.metafits,_20191014144345_gpubox20_00.fits} --no-diff --crosses --suffix '.ch135' --sel-pols xx +``` + +![ch135 satellite reflection](demo/data/1255099440/raw/1255099440.cross.ch135.xx.spectrum.png) + +#### `1261482120` channel 177 - Tile027 Sparkle + +```bash +python demo/04_ssins.py demo/data/1261482120/raw/1261482120{.metafits,_20191227114144_gpubox16_00.fits} --sigchain --no-diff --suffix '.ch177' --sel-pol xx +``` + +![ch177 Tile027 Sparkle](demo/data/1261482120/raw/1261482120.auto.ch177.xx.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1261482120/raw/1261482120{.metafits,_20191227114144_gpubox16_00.fits} --no-diff --crosses --suffix '.ch177' --skip-ants Tile027 --sel-pol yy +``` + +![ch177 Blip](demo/data/1261482120/raw/1261482120.cross.ch177.noTile027.yy.spectrum.png) + +#### `1324134264` channel 142 - slow TV 181-188MHz + +```bash +python demo/04_ssins.py demo/data/1324134264/raw/1324134264{.metafits,_20211221150406_ch142_000.fits} --suffix '.ch142.xxyy' --sigchain --no-diff --sel-pols xx yy +``` + +![ch142 Tile065 + HexE24 issue](demo/data/1324134264/raw/1324134264.auto.ch142.xxyy.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1324134264/raw/1324134264{.metafits,_20211221150406_ch142_000.fits} --suffix '.ch142.noT65E23' --crosses --no-diff --skip-ants Tile065 HexE23 --sel-pol xx +``` + +![ch142 Slow TV](demo/data/1324134264/raw/1324134264.cross.ch142.noT65E23.xx.spectrum.png) + +#### `1344506888` channel 137 - Rx07 issue + +```bash +python demo/04_ssins.py demo/data/1344506888/raw/1344506888{.metafits,_20220814100750_ch137_000.fits} --suffix '.ch137' --sigchain --no-diff ``` +![ch137 Rx07 + LBA5 + LBB5 issue](demo/data/1344506888/raw/1344506888.auto.ch137.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1344506888/raw/1344506888{.metafits,_20220814100750_ch137_000.fits} --suffix '.ch137.noRx7LBAB5' --no-dif --crosses --skip-ants Tile07{1..8} LBA5 LBB5 --sel-pol xx +``` + +![ch137 streaks](demo/data/1344506888/raw/1344506888.cross.ch137.noRx7LBAB5.xx.spectrum.png) + +#### `1360791928` - LB E,F,G issue + +```bash +python demo/04_ssins.py demo/data/1360791928/raw/1360791928{.metafits,_20230218214510_ch137_000.fits} --suffix '.ch137' --sigchain --no-diff --sel-pol xx +``` + +![ch137 LBE issue](demo/data/1360791928/raw/1360791928.auto.ch137.xx.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1360791928/raw/1360791928{.metafits,_20230218214510_ch137_000.fits} --suffix '.ch137' --sigchain --no-diff --skip-ants LBE{1..8} --sel-pol yy +``` + +![ch137 LB F,G issue](demo/data/1360791928/raw/1360791928.auto.ch137.yy.sigchain.png) + +#### `1361310560` + +```bash +python demo/04_ssins.py demo/data/1361310560/raw/1361310560{.metafits,_20230224214902_ch137_000.fits} --suffix '.ch137' --sigchain --no-diff --sel-pol xx +``` + +![ch137 LBE issue](demo/data/1361310560/raw/1361310560.auto.ch137.xx.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1361310560/raw/1361310560{.metafits,_20230224214902_ch137_000.fits} --suffix '.ch137.noLBE' --sigchain --no-diff --skip-ants LBE{1..8} --sel-pol xx +``` + +![ch137 Rx 07, 10, 12 issue](demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.xx.sigchain.png) + +```bash +python demo/04_ssins.py demo/data/1361310560/raw/1361310560{.metafits,_20230224214902_ch137_000.fits} --suffix '.ch137.noLBE.noRx7_10_12' --sigchain --no-diff --skip-ants LBE{1..8} Tile07{1..8} Tile10{1..8} Tile12{1..8} --sel-pol yy +``` + +![ch137 LB F,G issue](demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.noRx7_10_12.yy.sigchain.png) + +#### `1361397000` + +- Tile151 issue +- Rx10 & Rx12 issue + +```bash +python demo/04_ssins.py demo/data/1361397000/raw/1361397000{.metafits,_20230225214942_ch137_000.fits} --suffix '.ch137' --sigchain --no-diff +``` + +![ch137 Tile151, Rx 10, 12 issue](demo/data/1361397000/raw/1361397000.auto.ch137.sigchain.png) + +#### `1362519024` + +```bash +python demo/04_ssins.py demo/data/1362519024/raw/1362519024{.metafits,_20230310213006_ch065_000.fits} --suffix '.ch065' --sigchain --no-diff --sel-pol yy +``` + +![ch065 LBA1 issue](demo/data/1362519024/raw/1362519024.auto.ch065.yy.sigchain.png) + +### Configurations + now let's look at the rest of the obsids ```bash @@ -370,6 +598,8 @@ done ![images of each main MWA configuration](imgs/config_images.png) +### Joint deconvolution + combine them all into a single image ```bash @@ -377,6 +607,8 @@ rm -rf ${outdir}/combined/img/ obsid="combined" cal_ms=$(ls -1d ${outdir}/13*/cal/hyp_cal_*.ms ) demo/07_img.sh ``` +### Cleanup + clean up outdir to start fresh ```bash diff --git a/demo/04_ssins.py b/demo/04_ssins.py index aeae7b2..9e136f3 100755 --- a/demo/04_ssins.py +++ b/demo/04_ssins.py @@ -1,4 +1,5 @@ #!/usr/bin/env python +from pyuvdata import UVData from SSINS import SS, INS, MF import os import matplotlib as mpl @@ -6,227 +7,371 @@ import numpy as np from astropy.time import Time import argparse -import pandas as pd - -parser = argparse.ArgumentParser() -parser.add_argument("--no-diff", default=False, action="store_true") -parser.add_argument("--no-cache", default=False, action="store_true") -group = parser.add_mutually_exclusive_group() -group.add_argument("--crosses", default=False, action="store_true", help="plot crosses instead of default: autos") -group.add_argument("--sigchain", default=False, action="store_true", help="plot sigchain instead of default: spectrum", ) - -parser.add_argument("--threshold", default=5, help="match filter significance threshold") -parser.add_argument("--fontsize", default=8, help="plot font size") -parser.add_argument("files", nargs="+") -args = parser.parse_args() - -if len(args.files) > 1 and args.files[-1].endswith(".fits"): - metafits = args.files[0] - # output name is basename of metafits, or uvfits if provided - base, _ = os.path.splitext(metafits) -elif len(args.files) == 1: - vis = args.files[-1] - base, _ = os.path.splitext(vis) -else: - parser.print_usage() - exit(1) - -# sky-subtract https://ssins.readthedocs.io/en/latest/sky_subtract.html -ss = SS() - -def get_unflagged_ants(ss, args): - # TODO: args to flag additional ants - return np.unique(ss.ant_1_array) - -suffix = "" if args.no_diff else ".diff" -spectrum_type = "cross" if args.crosses else "auto" -suffix += f".{spectrum_type}" -cache = f"{base}{suffix}.h5" -if not args.no_cache and os.path.exists(cache): - print(f"reading from {cache=}") - ss.read_uvh5(cache, read_data=True, use_future_array_shapes=True) - unflagged_ants = get_unflagged_ants(ss, args) -else: - print(f"reading from {args.files=}") - ss.read( - args.files, - read_data=True, - diff=(not args.no_diff), # difference timesteps - remove_coarse_band=False, # does not work with low freq res - correct_van_vleck=False, # slow - remove_flagged_ants=True, # remove flagged antennas + +from matplotlib.axis import Axis + + +def get_parser(): + parser = argparse.ArgumentParser() + # arguments for UVData.read() + group_uvd_read = parser.add_argument_group("UVData.read") + group_uvd_read.add_argument("files", nargs="+") + group_uvd_read.add_argument( + "--no-diff", + default=False, + action="store_true", + help="don't difference visibilities in time (sky-subtract)", + ) + group_uvd_read.add_argument( + "--no-flag-init", + default=False, + action="store_true", + help="skip flagging of edge channels, quack time", ) - unflagged_ants = get_unflagged_ants(ss, args) - def cross_filter(pair): - a, b = pair - return a != b and a in unflagged_ants and b in unflagged_ants - def auto_filter(pair): - a, b = pair - return a == b and a in unflagged_ants - bls = [*filter( {"cross": cross_filter, "auto": auto_filter}[spectrum_type], ss.get_antpairs() )] - ss = ss.select(bls=bls, inplace=False) - ss.apply_flags(flag_choice="original") - if not args.no_cache: - ss.write_uvh5(cache) - print(f"wrote ss to {cache=}") + # arguments for UVData.select() + group_uvd_sel = parser.add_argument_group("UVData.select") + group_mutex = group_uvd_sel.add_mutually_exclusive_group() + group_mutex.add_argument( + "--sel-ants", + default=[], + nargs="*", + type=str, + help="antenna names to select, default: all unflagged", + ) + group_mutex.add_argument( + "--skip-ants", + default=[], + nargs="*", + type=str, + help="antenna names to skip, default: none", + ) + + group_uvd_sel.add_argument( + "--sel-pols", + default=[], + nargs="*", + type=str, + help="polarizations to select, default: all", + ) + + # arguments for SSINS.INS + parser.add_argument( + "--spectrum-type", + default="auto", + choices=["auto", "cross"], + help="analyse auto-correlations or cross-correlations. default: auto", + ) + + parser.add_argument( + "--crosses", + action="store_const", + const="cross", + dest="spectrum_type", + help="shorthand for --spectrum-type=cross", + ) + + parser.add_argument( + "--sigchain", + default=False, + action="store_true", + help="analyse z-scores for each tile and sum", + ) + + # arguments for SSINS.MF + parser.add_argument( + "--threshold", + default=5, + type=float, + help="match filter significance threshold. 0 disables match filter", + ) + parser.add_argument( + "--no-narrow", + default=False, + help="Don't look for narroband RFI", + ) + + # other + + parser.add_argument( + "--suffix", + default="", + type=str, + help="additional text to add to filename", + ) + + parser.add_argument("--fontsize", default=8, help="plot font size") + return parser + + +def get_unflagged_ants(ss: UVData, args): + all_ant_nums = np.array(ss.antenna_numbers) + all_ant_names = np.array(ss.antenna_names) + present_ant_nums = np.unique(ss.ant_1_array) + present_ant_mask = np.where(np.isin(all_ant_nums, present_ant_nums))[0] + present_ant_names = np.array([*map(str.upper, all_ant_names[present_ant_mask])]) + assert len(present_ant_nums) == len(present_ant_names) + + if args.sel_ants: + sel_ants = np.array([*map(str.upper, args.sel_ants)]) + return present_ant_nums[np.where(np.isin(present_ant_names, sel_ants))[0]] + elif args.skip_ants: + skip_ants = np.array([*map(str.upper, args.skip_ants)]) + return present_ant_nums[np.where(~np.isin(present_ant_names, skip_ants))[0]] + + return present_ant_nums + + +def get_gps_times(uvd: UVData): + return [*Time(np.unique(uvd.time_array), format="jd").gps] -# incoherent noise spectrum https://ssins.readthedocs.io/en/latest/incoherent_noise_spectrum.html -ins = INS(ss, spectrum_type=spectrum_type) -# match filter https://ssins.readthedocs.io/en/latest/match_filter.html -mf = MF(ss.freq_array, args.threshold, streak=True, narrow=False) -mf.apply_match_test(ins) -ins.sig_array[~np.isfinite(ins.sig_array)] = 0 +def get_suffix(args): + suffix = args.suffix + suffix = f".{args.spectrum_type}{suffix}" + if not args.no_diff: + suffix = f".diff{suffix}" + if len(args.sel_ants) == 1: + suffix = f"{suffix}.{args.sel_ants[0]}" + elif len(args.skip_ants) == 1: + suffix = f"{suffix}.no{args.skip_ants[0]}" + if len(args.sel_pols) == 1: + suffix = f"{suffix}.{args.sel_pols[0]}" + return suffix + + +def get_match_filter(ss, args): + """ + https://ssins.readthedocs.io/en/latest/match_filter.html + """ + return MF(ss.freq_array, args.threshold, streak=True, narrow=(not args.no_narrow)) + # #### # # PLOT # # #### # -plt.style.use("dark_background") -cmap = mpl.colormaps.get_cmap("viridis") -cmap.set_bad(color="#00000000") - -pols = ss.get_pols() -gps_times = np.unique(Time(ss.time_array, format="jd").gps) -int_time = np.unique(ss.integration_time)[0] -obsid = round(gps_times[0] - int_time / 2) -title = f"{obsid}" -plt.title(title) - -def plot_sigchain(ss, args): + + +def plot_sigchain(ss, args, obsname, suffix, cmap): + mf = get_match_filter(ss, args) + + unflagged_ants = get_unflagged_ants(ss, args) ant_mask = np.where(np.isin(ss.antenna_numbers, unflagged_ants))[0] - antenna_numbers = np.array(ss.antenna_numbers)[ant_mask] - antenna_names = np.array(ss.antenna_names)[ant_mask] + ant_numbers = np.array(ss.antenna_numbers)[ant_mask] + ant_names = np.array(ss.antenna_names)[ant_mask] + pols = ss.get_pols() + gps_times = get_gps_times(ss) + freqs_mhz = (ss.freq_array) / 1e6 + + # pad names + name_len = max(len(name) for name in ant_names) + ant_labels = [f"{name: <{name_len}}" for name in ant_names] # build a scores array for each signal chain scores = np.zeros((len(unflagged_ants), ss.Ntimes, ss.Nfreqs, ss.Npols)) - for ant_idx, (ant_num, ant_name) in enumerate(zip(antenna_numbers, antenna_names)): - print(ant_idx, ant_num, ant_name) + for ant_idx, (ant_num, ant_name) in enumerate(zip(ant_numbers, ant_names)): if ant_num not in unflagged_ants: continue # select only the auto-correlation for this antenna ssa = ss.select(antenna_nums=[(ant_num, ant_num)], inplace=False) - ins = INS(ssa, spectrum_type="auto") + ins = INS(ssa, spectrum_type=args.spectrum_type) mf.apply_match_test(ins) ins.sig_array[~np.isfinite(ins.sig_array)] = 0 scores[ant_idx] = ins.sig_array - plt.style.use("dark_background") subplots = plt.subplots( 2, len(pols), height_ratios=[len(unflagged_ants), ss.Nfreqs], - )[1] + )[1].reshape((2, len(pols))) def slice(scores, axis): - pscore = np.sqrt(np.sum(scores**2, axis=-1)) - # make per-antenna anomalies stand out more - # pscore = pscore / np.nanstd(pscore, axis=1)[:, np.newaxis] - # subtract minimum value - # pscore -= np.nanmin(pscore) - return pscore - - # pad names - namelen = max(len(name) for name in antenna_names) - antnames = [f"{name: <{namelen}}" for name in antenna_names] - channames = [f"{ch: 8.4f}" for ch in ss.freq_array / 1e6] - - for i, pol in enumerate(ss.get_pols()): - ax_signal, ax_spectrum = subplots[:, i] - title = f"{obsid} Autos z-score magnitude {pol}" + return np.sqrt(np.sum(scores**2, axis=axis)) + for i, pol in enumerate(pols): # by signal chain: [ant, time] + ax_signal: Axis = subplots[0, i] + ax_signal.set_title(f"{obsname} zscore{suffix} {pol if len(pols) > 1 else ''}") + if i == 0: + ax_signal.set_ylabel("Antenna") + signal_pscore = slice(scores[..., i], axis=-1) - print(signal_pscore.shape) - signal_df = pd.DataFrame( - signal_pscore.transpose(), index=gps_times, columns=antnames - ) - signal_df.to_csv( - f"{base}.auto_sig.{pol}.tsv", - index_label="gps_time", - float_format="%.3f", - sep="\t", - ) - ax_signal.imshow(signal_pscore, aspect="auto", interpolation="none", cmap=cmap) - ax_signal.set_ylabel("Antenna") - ax_spectrum.set_xlabel("Timestep") + ax_signal.imshow( + signal_pscore, + aspect="auto", + interpolation="none", + cmap=cmap, + extent=[ + np.min(gps_times), + np.max(gps_times), + len(ant_labels) - 0.5, + -0.5, + ], + ) ax_signal.set_yticks(np.arange(len(unflagged_ants))) ax_signal.tick_params(pad=True) - ax_signal.set_yticklabels(antnames, fontsize=args.fontsize) - ax_signal.set_title(f"{obsid} Autos z-score row-normalized {pol}") + ax_signal.set_yticklabels( + ant_labels, fontsize=args.fontsize, fontfamily="monospace" + ) # by spectrum: [freq, time] + ax_spectrum: Axis = subplots[1, i] + ax_spectrum.set_xlabel("GPS Time [s]") + if i == 0: + ax_spectrum.set_ylabel("Frequency channel [MHz]") + spectrum_pscore = slice(scores[..., i].transpose(2, 1, 0), axis=-1) - spectrum_df = pd.DataFrame( - spectrum_pscore.transpose(), index=gps_times, columns=channames - ) - spectrum_df.to_csv( - f"{base}.auto_spx.{pol}.tsv", - index_label="gps_time", - float_format="%.3f", - sep="\t", - ) - ax_spectrum.set_xlabel("Timestep") - ax_spectrum.set_ylabel("Frequency channel") - ax_spectrum.set_yticks(np.arange(ss.Nfreqs)) - ax_spectrum.set_yticklabels(channames, fontsize=args.fontsize) + ax_spectrum.tick_params(pad=True) - ax_spectrum.imshow(spectrum_pscore, aspect="auto", interpolation="none", cmap=cmap) + ax_spectrum.imshow( + spectrum_pscore, + aspect="auto", + interpolation="none", + cmap=cmap, + extent=[ + np.min(gps_times), + np.max(gps_times), + np.max(freqs_mhz), + np.min(freqs_mhz), + ], + ) + + +def plot_spectrum(ss, args, obsname, suffix, cmap): + # incoherent noise spectrum https://ssins.readthedocs.io/en/latest/incoherent_noise_spectrum.html + ins = INS(ss, spectrum_type=args.spectrum_type) + + mf = get_match_filter(ss, args) + mf.apply_match_test(ins) + ins.sig_array[~np.isfinite(ins.sig_array)] = 0 -def plot_spectrum(ss, args): pols = ss.get_pols() - subplot_rows = 2 + gps_times = get_gps_times(ss) + freqs_mhz = (ss.freq_array) / 1e6 subplots = plt.subplots( - subplot_rows, + 2, len(pols), sharex=True, sharey=True, - )[1] - - time_labels = [*Time(np.unique(ss.time_array), format="jd").iso] - ntimes_visible = 16 * 72 / args.fontsize / 2.4 / subplot_rows - time_stride = int(max(1, len(time_labels) // ntimes_visible)) - chan_labels = [f"{ch: 8.3f}" for ch in ss.freq_array / 1e6] - nchans_visible = 16 * 72 / args.fontsize / 2.0 / len(pols) - chan_stride = int(max(1, len(chan_labels) // nchans_visible)) + )[ + 1 + ].reshape((2, len(pols))) for i, pol in enumerate(pols): - ax_met, ax_sig = subplots[:, i] + # axis for metric being plotted + ax_met: Axis = subplots[0, i] - ax_met.set_title( - f'{base.split("/")[-1]} ss vis amps{suffix} {pol}' - ) + ax_met.set_title(f"{obsname} vis amps{suffix} {pol if len(pols) > 1 else ''}") ax_met.imshow( - ins.metric_array[..., i], aspect="auto", interpolation="none", cmap=cmap + ins.metric_array[..., i], + aspect="auto", + interpolation="none", + cmap=cmap, + extent=[ + np.min(freqs_mhz), + np.max(freqs_mhz), + np.max(gps_times), + np.min(gps_times), + ], + ) + + # axis for significance + ax_sig: Axis = subplots[1, i] + ax_sig.set_title(f"{obsname} z-score{suffix} {pol if len(pols) > 1 else ''}") + ax_sig.imshow( + ins.sig_array[..., i], + aspect="auto", + interpolation="none", + cmap=cmap, + extent=[ + np.min(freqs_mhz), + np.max(freqs_mhz), + np.max(gps_times), + np.min(gps_times), + ], ) - ax_sig.set_title(f'{base.split("/")[-1]} z-score {pol}') - ax_sig.imshow(ins.sig_array[..., i], aspect="auto", interpolation="none", cmap=cmap) if i == 0: - ax_met.set_ylabel("Timestep [UTC]") - ax_met.set_yticks(np.arange(ss.Ntimes)[::time_stride]) - ax_met.set_yticklabels( - time_labels[::time_stride], fontsize=args.fontsize, fontfamily="monospace" - ) - ax_sig.set_ylabel("Timestep [UTC]") - ax_sig.set_yticks(np.arange(ss.Ntimes)[::time_stride]) - ax_sig.set_yticklabels( - time_labels[::time_stride], fontsize=args.fontsize, fontfamily="monospace" - ) + ax_met.set_ylabel("GPS Time [s]") + ax_sig.set_ylabel("GPS Time [s]") ax_sig.set_xlabel("Frequency channel [MHz]") - ax_sig.set_xticks(np.arange(ss.Nfreqs)[::chan_stride]) - ax_sig.set_xticklabels(chan_labels[::chan_stride], fontsize=args.fontsize, rotation=90) - -if args.sigchain: - plot_sigchain(ss, args) - plt.gcf().set_size_inches(8 * len(pols), (len(unflagged_ants) + ss.Ntimes) * args.fontsize / 72) - figname = f"{base}{suffix}.sigchain.png" -else: - plot_spectrum(ss, args) - plt.gcf().set_size_inches(8 * len(pols), 16) - figname = f"{base}{suffix}.spectrum.png" - -plt.savefig(figname, bbox_inches="tight") -print(f"wrote {figname}") \ No newline at end of file + + +def main(): + parser = get_parser() + args = parser.parse_args() + print(f"{args=}") + + if len(args.files) > 1 and args.files[-1].endswith(".fits"): + metafits = args.files[0] + # output name is basename of metafits, or uvfits if provided + base, _ = os.path.splitext(metafits) + elif len(args.files) == 1: + vis = args.files[-1] + base, _ = os.path.splitext(vis) + else: + parser.print_usage() + exit(1) + + print(f"reading from {args.files=}") + # sky-subtract https://ssins.readthedocs.io/en/latest/sky_subtract.html + ss = SS() + ss.read( + args.files, + read_data=True, + diff=(not args.no_diff), # difference timesteps + remove_coarse_band=False, # does not work with low freq res + correct_van_vleck=False, # slow + remove_flagged_ants=True, # remove flagged antennas + flag_init=(not args.no_flag_init), + ) + + unflagged_ants = get_unflagged_ants(ss, args) + + select_kwargs = { + "inplace": False, + } + if args.spectrum_type == "cross": + select_kwargs["bls"] = [ + (a, b) + for (a, b) in ss.get_antpairs() + if a != b and (b in unflagged_ants or a in unflagged_ants) + ] + else: + select_kwargs["bls"] = [ + (a, b) for (a, b) in ss.get_antpairs() if a == b and a in unflagged_ants + ] + if args.sel_pols: + select_kwargs["polarizations"] = args.sel_pols + + ss = ss.select(**select_kwargs) + ss.apply_flags(flag_choice="original") + + plt.style.use("dark_background") + cmap = mpl.colormaps.get_cmap("viridis") + cmap.set_bad(color="#00000000") + + suffix = get_suffix(args) + + pols = ss.get_pols() + obsname = base.split("/")[-1] + plt.suptitle(f"{obsname}{suffix}") + if args.sigchain: + plot_sigchain(ss, args, obsname, suffix, cmap) + plt.gcf().set_size_inches( + 8 * len(pols), (len(unflagged_ants) + ss.Nfreqs) * args.fontsize / 72 + ) + figname = f"{base}{suffix}.sigchain.png" + else: + plot_spectrum(ss, args, obsname, suffix, cmap) + plt.gcf().set_size_inches(8 * len(pols), 16) + figname = f"{base}{suffix}.spectrum.png" + + plt.savefig(figname, bbox_inches="tight") + print(f"wrote {figname}") + + +if __name__ == "__main__": + main() diff --git a/demo/data/1060550888/raw/1060550888.auto.ch143.Tile104.yy.sigchain.png b/demo/data/1060550888/raw/1060550888.auto.ch143.Tile104.yy.sigchain.png new file mode 100644 index 0000000..71b5770 Binary files /dev/null and b/demo/data/1060550888/raw/1060550888.auto.ch143.Tile104.yy.sigchain.png differ diff --git a/demo/data/1087596040/raw/1087596040.cross.ch134.yy.spectrum.png b/demo/data/1087596040/raw/1087596040.cross.ch134.yy.spectrum.png new file mode 100644 index 0000000..ec804fd Binary files /dev/null and b/demo/data/1087596040/raw/1087596040.cross.ch134.yy.spectrum.png differ diff --git a/demo/data/1087596040/raw/1087596040.diff.auto.ch143.Tile108.xx.spectrum.png b/demo/data/1087596040/raw/1087596040.diff.auto.ch143.Tile108.xx.spectrum.png new file mode 100644 index 0000000..b9cae75 Binary files /dev/null and b/demo/data/1087596040/raw/1087596040.diff.auto.ch143.Tile108.xx.spectrum.png differ diff --git a/demo/data/1088806248/raw/1088806248.auto.ch145.Tile055.yy.sigchain.png b/demo/data/1088806248/raw/1088806248.auto.ch145.Tile055.yy.sigchain.png new file mode 100644 index 0000000..3faba23 Binary files /dev/null and b/demo/data/1088806248/raw/1088806248.auto.ch145.Tile055.yy.sigchain.png differ diff --git a/demo/data/1088806248/raw/1088806248.diff.cross.ch133.xx.spectrum.png b/demo/data/1088806248/raw/1088806248.diff.cross.ch133.xx.spectrum.png new file mode 100644 index 0000000..0a3a88a Binary files /dev/null and b/demo/data/1088806248/raw/1088806248.diff.cross.ch133.xx.spectrum.png differ diff --git a/demo/data/1089238040/raw/1089238040.auto.ch133.xx.spectrum.png b/demo/data/1089238040/raw/1089238040.auto.ch133.xx.spectrum.png new file mode 100644 index 0000000..88c95bc Binary files /dev/null and b/demo/data/1089238040/raw/1089238040.auto.ch133.xx.spectrum.png differ diff --git a/demo/data/1090871744/raw/1090871744.cross.ch137.xx.spectrum.png b/demo/data/1090871744/raw/1090871744.cross.ch137.xx.spectrum.png new file mode 100644 index 0000000..f5fa0e7 Binary files /dev/null and b/demo/data/1090871744/raw/1090871744.cross.ch137.xx.spectrum.png differ diff --git a/demo/data/1094319712/raw/1094319712.auto.ch142.xx.sigchain.png b/demo/data/1094319712/raw/1094319712.auto.ch142.xx.sigchain.png new file mode 100644 index 0000000..2d67a12 Binary files /dev/null and b/demo/data/1094319712/raw/1094319712.auto.ch142.xx.sigchain.png differ diff --git a/demo/data/1094319712/raw/1094319712.cross.ch142.noTile108.xx.spectrum.png b/demo/data/1094319712/raw/1094319712.cross.ch142.noTile108.xx.spectrum.png new file mode 100644 index 0000000..bd99ab4 Binary files /dev/null and b/demo/data/1094319712/raw/1094319712.cross.ch142.noTile108.xx.spectrum.png differ diff --git a/demo/data/1252516448/raw/1252516448.auto.ch142.xx.sigchain.png b/demo/data/1252516448/raw/1252516448.auto.ch142.xx.sigchain.png new file mode 100644 index 0000000..b92baef Binary files /dev/null and b/demo/data/1252516448/raw/1252516448.auto.ch142.xx.sigchain.png differ diff --git a/demo/data/1255099440/raw/1255099440.cross.ch135.xx.spectrum.png b/demo/data/1255099440/raw/1255099440.cross.ch135.xx.spectrum.png new file mode 100644 index 0000000..74eb6b4 Binary files /dev/null and b/demo/data/1255099440/raw/1255099440.cross.ch135.xx.spectrum.png differ diff --git a/demo/data/1261482120/raw/1261482120.auto.ch177.xx.sigchain.png b/demo/data/1261482120/raw/1261482120.auto.ch177.xx.sigchain.png new file mode 100644 index 0000000..7cb00bd Binary files /dev/null and b/demo/data/1261482120/raw/1261482120.auto.ch177.xx.sigchain.png differ diff --git a/demo/data/1261482120/raw/1261482120.cross.ch177.noTile027.yy.spectrum.png b/demo/data/1261482120/raw/1261482120.cross.ch177.noTile027.yy.spectrum.png new file mode 100644 index 0000000..a32f4b2 Binary files /dev/null and b/demo/data/1261482120/raw/1261482120.cross.ch177.noTile027.yy.spectrum.png differ diff --git a/demo/data/1324134264/raw/1324134264.auto.ch142.xxyy.sigchain.png b/demo/data/1324134264/raw/1324134264.auto.ch142.xxyy.sigchain.png new file mode 100644 index 0000000..b3d68b6 Binary files /dev/null and b/demo/data/1324134264/raw/1324134264.auto.ch142.xxyy.sigchain.png differ diff --git a/demo/data/1324134264/raw/1324134264.cross.ch142.noT65E23.xx.spectrum.png b/demo/data/1324134264/raw/1324134264.cross.ch142.noT65E23.xx.spectrum.png new file mode 100644 index 0000000..4c03c1b Binary files /dev/null and b/demo/data/1324134264/raw/1324134264.cross.ch142.noT65E23.xx.spectrum.png differ diff --git a/demo/data/1341914000/raw/1341914000.auto.ch137.LBA.xx.sigchain.png b/demo/data/1341914000/raw/1341914000.auto.ch137.LBA.xx.sigchain.png new file mode 100644 index 0000000..d5535e4 Binary files /dev/null and b/demo/data/1341914000/raw/1341914000.auto.ch137.LBA.xx.sigchain.png differ diff --git a/demo/data/1341914000/raw/1341914000.auto.ch137.noLBA5.xx.sigchain.png b/demo/data/1341914000/raw/1341914000.auto.ch137.noLBA5.xx.sigchain.png new file mode 100644 index 0000000..8e234c6 Binary files /dev/null and b/demo/data/1341914000/raw/1341914000.auto.ch137.noLBA5.xx.sigchain.png differ diff --git a/demo/data/1341914000/raw/1341914000.diff.auto.ch137.noLBA5.yy.spectrum.png b/demo/data/1341914000/raw/1341914000.diff.auto.ch137.noLBA5.yy.spectrum.png new file mode 100644 index 0000000..70328e1 Binary files /dev/null and b/demo/data/1341914000/raw/1341914000.diff.auto.ch137.noLBA5.yy.spectrum.png differ diff --git a/demo/data/1341914000/raw/1341914000.diff.cross.ch137.noLBA5.yy.spectrum.png b/demo/data/1341914000/raw/1341914000.diff.cross.ch137.noLBA5.yy.spectrum.png new file mode 100644 index 0000000..3b8df82 Binary files /dev/null and b/demo/data/1341914000/raw/1341914000.diff.cross.ch137.noLBA5.yy.spectrum.png differ diff --git a/demo/data/1344506888/raw/1344506888.auto.ch137.sigchain.png b/demo/data/1344506888/raw/1344506888.auto.ch137.sigchain.png new file mode 100644 index 0000000..2ba0f0b Binary files /dev/null and b/demo/data/1344506888/raw/1344506888.auto.ch137.sigchain.png differ diff --git a/demo/data/1344506888/raw/1344506888.cross.ch137.noRx7LBAB5.xx.spectrum.png b/demo/data/1344506888/raw/1344506888.cross.ch137.noRx7LBAB5.xx.spectrum.png new file mode 100644 index 0000000..3344096 Binary files /dev/null and b/demo/data/1344506888/raw/1344506888.cross.ch137.noRx7LBAB5.xx.spectrum.png differ diff --git a/demo/data/1360791928/raw/1360791928.auto.ch137.xx.sigchain.png b/demo/data/1360791928/raw/1360791928.auto.ch137.xx.sigchain.png new file mode 100644 index 0000000..70808f0 Binary files /dev/null and b/demo/data/1360791928/raw/1360791928.auto.ch137.xx.sigchain.png differ diff --git a/demo/data/1360791928/raw/1360791928.auto.ch137.yy.sigchain.png b/demo/data/1360791928/raw/1360791928.auto.ch137.yy.sigchain.png new file mode 100644 index 0000000..732bf39 Binary files /dev/null and b/demo/data/1360791928/raw/1360791928.auto.ch137.yy.sigchain.png differ diff --git a/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.noRx7_10_12.yy.sigchain.png b/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.noRx7_10_12.yy.sigchain.png new file mode 100644 index 0000000..9996fd9 Binary files /dev/null and b/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.noRx7_10_12.yy.sigchain.png differ diff --git a/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.xx.sigchain.png b/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.xx.sigchain.png new file mode 100644 index 0000000..a6f49f3 Binary files /dev/null and b/demo/data/1361310560/raw/1361310560.auto.ch137.noLBE.xx.sigchain.png differ diff --git a/demo/data/1361310560/raw/1361310560.auto.ch137.xx.sigchain.png b/demo/data/1361310560/raw/1361310560.auto.ch137.xx.sigchain.png new file mode 100644 index 0000000..e7e9178 Binary files /dev/null and b/demo/data/1361310560/raw/1361310560.auto.ch137.xx.sigchain.png differ diff --git a/demo/data/1361397000/raw/1361397000.auto.ch137.sigchain.png b/demo/data/1361397000/raw/1361397000.auto.ch137.sigchain.png new file mode 100644 index 0000000..6a147bb Binary files /dev/null and b/demo/data/1361397000/raw/1361397000.auto.ch137.sigchain.png differ diff --git a/demo/data/1362519024/raw/1362519024.auto.ch065.yy.sigchain.png b/demo/data/1362519024/raw/1362519024.auto.ch065.yy.sigchain.png new file mode 100644 index 0000000..d4a6091 Binary files /dev/null and b/demo/data/1362519024/raw/1362519024.auto.ch065.yy.sigchain.png differ