import argparse
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import math
import mwalib
matplotlib.use("GTK3Agg")
dpi = 100
MODE_RANGE = "RANGE"
MODE_BASELINE = "BASELINE"
class ViewFITSArgs:
def __init__(self, passed_args):
self.filename = passed_args["filename"]
self.metafits_filename = passed_args["metafits"]
self.time_step1 = passed_args["timestep1"]
self.time_step2 = passed_args["timestep2"]
self.tile1 = passed_args["ant1"]
self.tile2 = passed_args["ant2"]
self.autos_only = passed_args["autosonly"]
self.channel1 = passed_args["channel1"]
self.channel2 = passed_args["channel2"]
self.ppd_plot = passed_args["ppdplot"]
self.ppd_plot2 = passed_args["ppdplot2"]
self.grid_plot = passed_args["gridplot"]
self.grid_plot2 = passed_args["gridplot2"]
self.grid_pol = passed_args["gridpol"]
self.phase_plot_all = passed_args["phaseplot_all"]
self.phase_plot_one = passed_args["phaseplot_one"]
self.mode = passed_args["mode"]
self.dumpraw = passed_args["dump_raw"]
self.dumpplot = passed_args["dump_plot"]
self.any_plotting = (
self.ppd_plot or self.ppd_plot2 or self.grid_plot or self.phase_plot_all or self.phase_plot_one
)
if not self.any_plotting and self.dumpplot:
print("You must be plotting to dump out the plot data to file!")
exit(-1)
self.values = 2
self.validate_params()
def validate_params(self): print(f"Opening with mwalib using metafits file {self.metafits_filename} and data file {self.filename}...")
self.context = mwalib.CorrelatorContext(
self.metafits_filename,
[
self.filename,
],
)
for p in self.context.provided_coarse_chan_indices:
print(self.context.coarse_chans[p])
self.correlator_version: mwalib.MWAVersion = self.context.mwa_version
self.pols = self.context.metafits_context.num_visibility_pols
self.fits_tiles = len(self.context.metafits_context.antennas)
self.chan_x_pols_x_vals = self.context.metafits_context.num_corr_fine_chans_per_coarse * self.pols * self.values
self.fits_channels = self.context.metafits_context.num_corr_fine_chans_per_coarse
self.fits_has_weights = True
if self.mode != MODE_RANGE and self.mode != MODE_BASELINE:
print(f"Error, mode should be {MODE_RANGE} or {MODE_BASELINE}!")
exit(-1)
if self.tile1 == -1:
self.tile1 = 0
if self.tile2 == -1:
self.tile2 = self.fits_tiles - 1
if self.tile1 > self.fits_tiles - 1:
print("Error ant1 is more than the last tile index!")
exit(-1)
if self.tile2 > self.fits_tiles - 1:
print("Error ant2 is more than the last tile index!")
exit(-1)
if self.tile1 > self.tile2:
print("Error ant1 is more than the ant2!")
exit(-1)
self.fits_time_steps = self.context.num_timesteps
if self.time_step1 == -1:
self.time_step1 = 1
if self.time_step2 == -1:
self.time_step2 = self.fits_time_steps
if self.time_step1 > self.fits_time_steps:
print("Error t1 is more than the max time step!")
exit(-1)
if self.time_step2 > self.fits_time_steps:
print("Error t2 is more than the max time step!")
exit(-1)
if self.time_step1 > self.time_step2:
print("Error t1 is more than the t2!")
exit(-1)
if self.channel1 == -1:
self.channel1 = 0
if self.channel2 == -1:
self.channel2 = self.fits_channels - 1
if self.channel1 > self.fits_channels - 1:
print("Error c1 is more than the number of the last fine channel!")
exit(-1)
if self.channel2 > self.fits_channels - 1:
print("Error c2 is more than the number of the last fine channel!")
exit(-1)
if self.channel1 > self.channel2:
print("Error c1 is more than the c2!")
exit(-1)
self.tile_count = self.tile2 - self.tile1 + 1
if self.mode == MODE_BASELINE:
self.baseline_count = 1
else:
self.baseline_count = int((self.tile_count * (self.tile_count + 1)) / 2)
self.fits_baseline_count = int((self.fits_tiles * (self.fits_tiles + 1)) / 2)
self.channel_count = self.channel2 - self.channel1 + 1
self.time_step_count = self.time_step2 - self.time_step1 + 1
self.hdu_time1 = self.time_step1 - 1
self.unix_time1 = self.context.timesteps[self.hdu_time1].unix_time_ms / 1000.0
self.hdu_time2 = self.time_step2 - 1
self.unix_time2 = self.context.timesteps[self.hdu_time2].unix_time_ms / 1000.0
if self.grid_plot2:
self.param_string = (
f"[{self.grid_pol}] {self.filename} "
f"t={self.time_step1}-{self.time_step2} "
f"t_hdu={self.hdu_time1}-{self.hdu_time2} \n"
f"t_unix={self.unix_time1}-{self.unix_time2} "
f"tile={self.tile1}-{self.tile2} "
f"ch={self.channel1}-{self.channel2} "
f"autosonly?={self.autos_only} "
f"{self.tile_count}t/{self.fits_tiles}t"
)
elif self.phase_plot_one:
self.param_string = (
f"[{self.grid_pol}] {self.filename} "
f"t={self.time_step1}-{self.time_step2} "
f"t_hdu={self.hdu_time1}-{self.hdu_time2} \n"
f"t_unix={self.unix_time1}-{self.unix_time2} "
f"baseline={self.tile1}-{self.tile2} "
f"chs={self.channel1}-{self.channel2} "
)
else:
self.param_string = (
f"{self.filename} t={self.time_step1}-{self.time_step2} "
f"t_hdu={self.hdu_time1}-{self.hdu_time2} \n"
f"t_unix={self.unix_time1}-{self.unix_time2} "
f"tile={self.tile1}-{self.tile2} "
f"ch={self.channel1}-{self.channel2} "
f"autosonly?={self.autos_only} "
f"{self.tile_count}t/{self.fits_tiles}t"
)
print(self.param_string)
def meets_criteria(i, j, a1, a2, mode):
if mode == MODE_RANGE:
return i >= a1 and j <= a2
elif mode == MODE_BASELINE:
return i == a1 and j == a2
else:
return None
def peek_fits(program_args: ViewFITSArgs): print("Initialising data structures...")
plot_ppd_data_x = None
plot_ppd_data_y = None
plot_ppd2_data_x = None
plot_ppd2_data_y = None
plot_grid_data = None
plot_grid2_data = None
plot_phase_data_x = None
plot_phase_data_y = None
if program_args.ppd_plot:
plot_ppd_data_x = np.empty(shape=(program_args.channel_count, program_args.time_step_count))
plot_ppd_data_x.fill(0)
plot_ppd_data_y = np.empty(shape=(program_args.channel_count, program_args.time_step_count))
plot_ppd_data_y.fill(0)
if program_args.ppd_plot2:
plot_ppd2_data_x = np.empty(
shape=(
program_args.time_step_count,
program_args.baseline_count,
program_args.channel_count,
)
)
plot_ppd2_data_x.fill(0)
plot_ppd2_data_y = np.empty(
shape=(
program_args.time_step_count,
program_args.baseline_count,
program_args.channel_count,
)
)
plot_ppd2_data_y.fill(0)
if program_args.grid_plot:
plot_grid_data = np.empty(
shape=(
program_args.time_step_count,
program_args.tile_count,
program_args.tile_count,
)
)
plot_grid_data.fill(0)
if program_args.grid_plot2:
plot_grid2_data = np.empty(
shape=(
program_args.time_step_count,
program_args.tile_count,
program_args.tile_count,
)
)
plot_grid2_data.fill(0)
if program_args.phase_plot_all or program_args.phase_plot_one:
plot_phase_data_x = np.empty(
shape=(
program_args.time_step_count,
program_args.baseline_count,
program_args.channel_count,
)
)
plot_phase_data_x.fill(0)
plot_phase_data_y = np.empty(
shape=(
program_args.time_step_count,
program_args.baseline_count,
program_args.channel_count,
)
)
plot_phase_data_y.fill(0)
raw_dump_file = None
plot_dump_file = None
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = f"{program_args.context.metafits_context.obs_id}_mwax.csv"
else:
filename = f"{program_args.context.metafits_context.obs_id}_mwa.csv"
if program_args.dumpplot:
plot_dump_filename = f"plot_dump_{filename}"
print(f"Openning {plot_dump_filename} for writing")
plot_dump_file = open(plot_dump_filename, "w")
if program_args.ppd_plot:
plot_dump_file.write("time_index, fine_chan, x, y\n")
elif program_args.ppd_plot2:
plot_dump_file.write("plot_number, time_index, fine_chan, x, y\n")
elif program_args.grid_plot:
plot_dump_file.write("unix_time, tile1, tile2, log10_scaled_power\n")
elif program_args.grid_plot2:
plot_dump_file.write("unix_time, tile1, tile2, log10_scaled_power\n")
elif program_args.phase_plot_one:
plot_dump_file.write("time_index, baseline, x, y\n")
elif program_args.phase_plot_all:
plot_dump_file.write("time_index, baseline, x, y\n")
if program_args.dumpraw:
raw_dump_filename = f"raw_dump_{filename}"
print(f"Openning {raw_dump_filename} for writing")
raw_dump_file = open(raw_dump_filename, "w")
raw_dump_file.write(
"unix_time, baseline, ant1, ant2, fine_ch, "
"xx_r, xx_i, xy_r, xy_i, yx_r, yx_i, yy_r, yy_i, "
"xx_pow, xy_pow, yx_pow, yy_pow, x_phase_deg, y_phase_deg\n"
)
for timestep_index, timestep in enumerate(program_args.context.timesteps):
time_index = timestep_index - program_args.hdu_time1
if timestep_index < program_args.hdu_time1 or timestep_index > program_args.hdu_time2:
continue
else:
print(f"Processing timestep: {timestep_index} (time index: {time_index})...")
data = program_args.context.read_by_baseline(
timestep_index,
program_args.context.provided_coarse_chan_indices[0],
)
data = data.reshape(
program_args.context.metafits_context.num_baselines,
program_args.chan_x_pols_x_vals,
)
baseline = 0
selected_baseline = 0
if time_index == 0:
print(f"QUAKTIME:{program_args.context.metafits_context.quack_time_duration_ms / 1000.0} s\n")
print("\nUnflagged tiles:")
print("================")
for i in range(0, len(program_args.context.metafits_context.antennas)):
if (
program_args.context.metafits_context.antennas[i].rfinput_x.flagged is False
and program_args.context.metafits_context.antennas[i].rfinput_y.flagged is False
):
print(
f"Index {i}, TileID: {program_args.context.metafits_context.antennas[i].tile_id} "
f"{program_args.context.metafits_context.antennas[i].tile_name}"
f" (rec:{program_args.context.metafits_context.antennas[i].rfinput_x.rec_number},"
f"slot:{program_args.context.metafits_context.antennas[i].rfinput_x.rec_slot_number})"
)
print("\nFlagged tiles:")
print("================")
for i in range(0, len(program_args.context.metafits_context.antennas)):
if (
program_args.context.metafits_context.antennas[i].rfinput_x.flagged
or program_args.context.metafits_context.antennas[i].rfinput_y.flagged
):
print(
f"Index {i}, TileID: {program_args.context.metafits_context.antennas[i].tile_id} "
f"{program_args.context.metafits_context.antennas[i].tile_name}"
f" (rec:{program_args.context.metafits_context.antennas[i].rfinput_x.rec_number},"
f"slot:{program_args.context.metafits_context.antennas[i].rfinput_x.rec_slot_number})"
)
for i in range(0, program_args.tile2 + 1):
for j in range(i, program_args.fits_tiles):
if ((i == j and program_args.autos_only is True) or (program_args.autos_only is False)) and (
meets_criteria(
i,
j,
program_args.tile1,
program_args.tile2,
program_args.mode,
)
):
for chan in range(program_args.channel1, program_args.channel2 + 1):
index = chan * (program_args.pols * program_args.values)
if (
program_args.context.metafits_context.antennas[i].rfinput_x.flagged
or program_args.context.metafits_context.antennas[i].rfinput_y.flagged
or program_args.context.metafits_context.antennas[j].rfinput_x.flagged
or program_args.context.metafits_context.antennas[j].rfinput_y.flagged
) and 1 == 0:
xx_r = 0
xx_i = 0
xy_r = 0
xy_i = 0
yx_r = 0
yx_i = 0
yy_r = 0
yy_i = 0
power_xx = 0
power_xy = 0
power_yx = 0
power_yy = 0
phase_x = 0
phase_y = 0
else:
xx_r = data[baseline][index]
xx_i = data[baseline][index + 1]
xy_r = data[baseline][index + 2]
xy_i = data[baseline][index + 3]
yx_r = data[baseline][index + 4]
yx_i = data[baseline][index + 5]
yy_r = data[baseline][index + 6]
yy_i = data[baseline][index + 7]
power_xx = xx_i * xx_i + xx_r * xx_r
power_xy = xy_i * xy_i + xy_r * xy_r
power_yx = yx_i * yx_i + yx_r * yx_r
power_yy = yy_i * yy_i + yy_r * yy_r
phase_x = math.degrees(math.atan2(xx_i, xx_r))
phase_y = math.degrees(math.atan2(yy_i, yy_r))
if program_args.ppd_plot:
plot_ppd_data_x[chan][time_index] = plot_ppd_data_x[chan][time_index] + (
power_xx / program_args.baseline_count
)
plot_ppd_data_y[chan][time_index] = plot_ppd_data_y[chan][time_index] + (
power_yy / program_args.baseline_count
)
elif program_args.ppd_plot2:
plot_ppd2_data_x[time_index][selected_baseline][chan] = power_xx
plot_ppd2_data_y[time_index][selected_baseline][chan] = power_yy
elif program_args.grid_plot:
plot_grid_data[time_index][j][i] = plot_grid_data[time_index][j][i] + power_xx + power_yy
elif program_args.grid_plot2:
if program_args.grid_pol == "XX":
plot_grid2_data[time_index][j][i] += power_xx
elif program_args.grid_pol == "XY":
plot_grid2_data[time_index][j][i] += power_xy
elif program_args.grid_pol == "YX":
plot_grid2_data[time_index][j][i] += power_yx
elif program_args.grid_pol == "YY":
plot_grid2_data[time_index][j][i] += power_yy
else:
print(
"Grid Plot requires you to specify "
"-gp (--gridpol) and "
"takes XX,XY,YX,YY as parameters"
)
exit(-1)
elif program_args.phase_plot_all or program_args.phase_plot_one:
plot_phase_data_x[time_index][selected_baseline][chan] = phase_x
plot_phase_data_y[time_index][selected_baseline][chan] = phase_y
if program_args.dumpraw:
raw_dump_file.write(
f"{timestep.unix_time_ms / 1000.0},"
f"{baseline},"
f"{i},"
f"{j},"
f"{chan},"
f"{xx_r},"
f"{xx_i},"
f"{xy_r},"
f"{xy_i},"
f"{yx_r},"
f"{yx_i},"
f"{yy_r},"
f"{yy_i},"
f"{power_xx},"
f"{power_xy},"
f"{power_yx},"
f"{power_yy},"
f"{phase_x},"
f"{phase_y}\n"
)
selected_baseline = selected_baseline + 1
baseline = baseline + 1
print("Processing of data done!")
if program_args.dumpraw:
if raw_dump_file:
raw_dump_file.close()
if program_args.ppd_plot:
convert_to_db = False
do_ppd_plot(
program_args.param_string,
program_args,
plot_ppd_data_x,
plot_ppd_data_y,
convert_to_db,
)
if program_args.ppd_plot2:
convert_to_db = False
do_ppd_plot2(
program_args.param_string,
program_args,
plot_ppd2_data_x,
plot_ppd2_data_y,
convert_to_db,
)
if program_args.grid_plot:
do_grid_plot(program_args.param_string, program_args, plot_grid_data)
if program_args.grid_plot2:
do_grid_plot2(program_args.param_string, program_args, plot_grid2_data)
if program_args.phase_plot_all:
do_phase_plot(
program_args.param_string,
program_args,
plot_phase_data_x,
plot_phase_data_y,
)
if program_args.phase_plot_one:
do_phase_plot1(
program_args.param_string,
program_args,
plot_phase_data_x,
plot_phase_data_y,
)
if program_args.dumpplot:
if plot_dump_file:
plot_dump_file.close()
print("view_fits Done!\n")
def do_ppd_plot(
title,
program_args: ViewFITSArgs,
plot_ppd_data_x,
plot_ppd_data_y,
convert_to_db,
):
print("Preparing ppd plot...")
min_db = 0
if convert_to_db:
for t in range(0, program_args.time_step_count):
for c in range(0, program_args.channel_count):
plot_ppd_data_x[c][t] = math.log10(plot_ppd_data_x[c][t] + 1) * 10
plot_ppd_data_y[c][t] = math.log10(plot_ppd_data_y[c][t] + 1) * 10
min_db = 0
fig, ax = plt.subplots(
figsize=(20, 10),
nrows=1,
ncols=1,
squeeze=False,
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
for t in range(0, program_args.time_step_count):
print(f"Adding data points for time ({t})...")
for c in range(0, program_args.channel_count):
plot_ppd_data_x[c][t] = plot_ppd_data_x[c][t] - min_db
plot_ppd_data_y[c][t] = plot_ppd_data_y[c][t] - min_db
plot = ax[0][0]
plot.plot(plot_ppd_data_x, "o", color="blue")
plot.plot(plot_ppd_data_y, "o", color="green")
if convert_to_db:
plot.set_ylabel("dB", size=6)
else:
plot.set_ylabel("Raw value", size=6)
plot.set_xlabel("fine channel", size=6)
plot.set_title(f"t={program_args.unix_time1}", size=6)
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = "ppd_plot_mwax.png"
else:
filename = "ppd_plot_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
def do_ppd_plot2(
title,
program_args: ViewFITSArgs,
plot_ppd_data_x,
plot_ppd_data_y,
convert_to_db,
): print("Preparing ppd plot2...")
plots = program_args.time_step_count
plot_rows = math.floor(math.sqrt(plots))
plot_cols = math.ceil(plots / plot_rows)
plot_row = 0
plot_col = 0
min_db = 0
if convert_to_db:
for t in range(0, program_args.time_step_count):
for c in range(0, program_args.channel_count):
for b in range(0, program_args.baseline_count):
plot_ppd_data_x[t][b][c] = math.log10(plot_ppd_data_x[t][b][c] + 1) * 10
plot_ppd_data_y[t][b][c] = math.log10(plot_ppd_data_y[t][b][c] + 1) * 10
min_db = min(
np.array(plot_ppd_data_x).min(initial=0),
np.array(plot_ppd_data_y).min(initial=0),
)
fig, ax = plt.subplots(
figsize=(20, 10),
nrows=plot_rows,
ncols=plot_cols,
squeeze=False,
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
for t in range(0, program_args.time_step_count):
if min_db != 0:
for c in range(0, program_args.channel_count):
for b in range(0, program_args.baseline_count):
plot_ppd_data_x[t][b][c] = plot_ppd_data_x[t][b][c] - min_db
plot_ppd_data_y[t][b][c] = plot_ppd_data_y[t][b][c] - min_db
plot = ax[plot_row][plot_col]
for b in range(0, program_args.baseline_count):
plot.plot(plot_ppd_data_x[t][b], "o", markersize=1, color="blue")
plot.plot(plot_ppd_data_y[t][b], "o", markersize=1, color="green")
if convert_to_db:
plot.set_ylabel("dB", size=6)
else:
plot.set_ylabel("Raw value", size=6)
plot.set_xlabel("fine channel", size=6)
plot.set_title(f"t={t + program_args.time_step1}", size=6)
if plot_col < plot_cols - 1:
plot_col = plot_col + 1
else:
plot_row = plot_row + 1
plot_col = 0
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = "ppd_plot2_mwax.png"
else:
filename = "ppd_plot2_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
def do_grid_plot(title, program_args: ViewFITSArgs, plot_grid_data): print("Preparing grid plot...")
plots = program_args.time_step_count
plot_rows = math.floor(math.sqrt(plots))
plot_cols = math.ceil(plots / plot_rows)
plot_row = 0
plot_col = 0
if 1 == 1:
for time_index in range(0, program_args.time_step_count):
for t1 in range(0, program_args.tile_count):
for t2 in range(t1, program_args.tile_count):
plot_grid_data[time_index][t2][t1] = math.log10(plot_grid_data[time_index][t2][t1] + 1) * 10
fig, ax = plt.subplots(
figsize=(30, 30),
nrows=plot_rows,
ncols=plot_cols,
squeeze=False,
sharex="all",
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
n_step = math.ceil(plots / 1.25)
if n_step % 2 != 0:
n_step = n_step + 1
if n_step % 4 != 0:
n_step = n_step - 2
if n_step <= 1:
n_step = 1
else:
n_step = n_step * 2
if n_step > 16:
n_step = 16
for time_index in range(0, program_args.time_step_count):
for t1 in range(0, program_args.tile_count):
for t2 in range(t1, program_args.tile_count):
print(
time_index,
program_args.tile1 + t1,
program_args.tile1 + t2,
plot_grid_data[time_index][t2][t1],
)
print(f"Adding data points for plot({time_index})...")
plot = ax[plot_row][plot_col]
plot.imshow(plot_grid_data[time_index], cmap="inferno", interpolation="None")
plot.set_title(f"t={time_index + program_args.time_step1}", size=6)
plot.set_xticks(np.arange(program_args.tile_count, step=n_step))
plot.set_yticks(np.arange(program_args.tile_count, step=n_step))
plot.set_xticklabels(np.arange(program_args.tile1, program_args.tile2 + 1, step=n_step))
plot.set_yticklabels(np.arange(program_args.tile1, program_args.tile2 + 1, step=n_step))
if plot_col == 0:
plot.set_ylabel("ant2", size=6)
if plot_row == plot_rows - 1:
plot.set_xlabel("ant1", size=6)
plt.setp(
plot.get_xticklabels(),
rotation=90,
ha="right",
va="center",
rotation_mode="anchor",
)
if plot_col < plot_cols - 1:
plot_col = plot_col + 1
else:
plot_row = plot_row + 1
plot_col = 0
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = "grid_plot_mwax.png"
else:
filename = "grid_plot_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
def do_grid_plot2(title, program_args: ViewFITSArgs, plot_grid_data): print("Preparing grid plot2...")
plots = program_args.time_step_count
plot_rows = math.floor(math.sqrt(plots))
plot_cols = math.ceil(plots / plot_rows)
plot_row = 0
plot_col = 0
scaling_value: float = np.max(plot_grid_data)
for time_index in range(0, program_args.time_step_count):
for t1 in range(0, program_args.tile_count):
for t2 in range(t1, program_args.tile_count):
value: float = plot_grid_data[time_index][t2][t1] / scaling_value
try:
if abs(value) < 0.0000001:
plot_grid_data[time_index][t2][t1] = 0
else:
plot_grid_data[time_index][t2][t1] = math.log10(value * 10000000)
except Exception as e:
print(f"Exception trying math.log10({value * 10000000} / {scaling_value})")
raise e
print(
f"scaling value was: {scaling_value}; now scaled and "
f"log10'd-> min = {np.min(plot_grid_data)}; "
f"max = {np.max(plot_grid_data)}"
)
fig, ax = plt.subplots(
figsize=(30, 30),
nrows=plot_rows,
ncols=plot_cols,
squeeze=False,
sharex="all",
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
n_step = 1
for time_index in range(0, program_args.time_step_count):
print(f"Adding data points for plot({time_index})...")
plot = ax[plot_row][plot_col]
plot.imshow(plot_grid_data[time_index], cmap="inferno", interpolation="None")
plot.set_title(f"t={time_index + program_args.time_step1}", size=6)
plot.set_xticks(np.arange(program_args.tile_count, step=n_step))
plot.set_yticks(np.arange(program_args.tile_count, step=n_step))
plot.set_xticklabels(np.arange(program_args.tile1, program_args.tile2 + 1, step=n_step))
plot.set_yticklabels(np.arange(program_args.tile1, program_args.tile2 + 1, step=n_step))
if plot_col == 0:
plot.set_ylabel("ant2", size=6)
if plot_row == plot_rows - 1:
plot.set_xlabel("ant1", size=6)
plt.setp(
plot.get_xticklabels(),
rotation=90,
ha="right",
va="center",
rotation_mode="anchor",
)
if plot_col < plot_cols - 1:
plot_col = plot_col + 1
else:
plot_row = plot_row + 1
plot_col = 0
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = f"grid_plot2_{program_args.grid_pol}_mwax.png"
else:
filename = f"grid_plot2_{program_args.grid_pol}_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
def do_phase_plot(title, program_args: ViewFITSArgs, plot_phase_data_x, plot_phase_data_y): print("Preparing phase plot...")
plots = program_args.baseline_count
print(
f"Timesteps: {program_args.time_step_count}, "
f"baselines: {program_args.baseline_count}, "
f"tiles: {program_args.tile_count}, "
f"channels: {program_args.channel_count}, "
f"plots: {plots}"
)
plot_cols = program_args.tile_count
plot_rows = program_args.tile_count
plot_row = 0
plot_col = 0
baseline = 0
fig, ax = plt.subplots(
figsize=(11.6, 8.3),
nrows=plot_rows,
ncols=plot_cols,
squeeze=False,
sharex="all",
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
for i in range(0, program_args.tile_count):
for j in range(i, program_args.tile_count):
if program_args.phase_plot_one:
if not (i == 0 and j == (program_args.tile_count - 1)):
print(f"{i} vs {j} skip")
baseline = baseline + 1
continue
print(f"Adding data points for plot({i},{j})...")
channel_list = range(program_args.channel1, program_args.channel2 + 1)
plot = ax[plot_row][plot_col]
for t in range(0, program_args.time_step_count):
plot.plot(
channel_list,
plot_phase_data_x[t][baseline],
"o",
markersize=1,
color="blue",
)
plot.plot(
channel_list,
plot_phase_data_y[t][baseline],
"o",
markersize=1,
color="green",
)
if plot_col == 0:
plot.set_ylabel("phase (deg)", size=6)
if plot_row == plot_rows - 1:
plot.set_xlabel("fine channel", size=6)
plot.set_ylim([-180, 180])
tile1_name = program_args.context.metafits_context.antennas[i + program_args.tile1].tile_name
tile2_name = program_args.context.metafits_context.antennas[j + program_args.tile1].tile_name
plot.set_title(f"{tile1_name} v {tile2_name}", size=6, pad=2)
if plot_col < plot_cols - 1:
plot_col = plot_col + 1
else:
plot_row = plot_row + 1
plot_col = plot_row
baseline = baseline + 1
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = "phase_plot_mwax.png"
else:
filename = "phase_plot_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
def do_phase_plot1(title, program_args: ViewFITSArgs, plot_phase_data_x, plot_phase_data_y): print("Preparing phase plot1...")
plots = program_args.channel_count
sqrt_of_plots: int = math.floor(math.sqrt(plots))
plot_cols = sqrt_of_plots
plot_rows = program_args.channel_count // sqrt_of_plots
print(plot_rows, plot_cols)
print(
f"Timesteps: {program_args.time_step_count}, "
f"baselines: {program_args.baseline_count}, "
f"tiles: {program_args.tile_count}, "
f"channels: {program_args.channel_count}, "
f"plots: {plots}"
)
plot_row = 0
plot_col = 0
baseline = 0
fine_chan = 0
fig, ax = plt.subplots(
figsize=(11.6, 8.3),
nrows=plot_rows,
ncols=plot_cols,
squeeze=False,
sharex="all",
sharey="all",
dpi=dpi,
)
fig.suptitle(title)
print("Adding data points for plot...")
timestep_list = range(program_args.time_step1, program_args.time_step2 + 1)
for plot_row in range(0, plot_rows):
for plot_col in range(0, plot_cols):
plot_data_x = plot_phase_data_x[0 : program_args.time_step_count, baseline, fine_chan]
plot_data_y = plot_phase_data_y[0 : program_args.time_step_count, baseline, fine_chan]
plot = ax[plot_row][plot_col]
plot.plot(
timestep_list,
plot_data_x,
"o",
markersize=1,
color="blue",
)
plot.plot(
timestep_list,
plot_data_y,
"o",
markersize=1,
color="green",
)
plot.set_ylabel("phase (deg)", size=6)
plot.set_xlabel("timestep", size=6)
plot.set_ylim([-180, 180])
plot.set_title(
f"Fine ch: {program_args.context.metafits_context.metafits_fine_chan_freqs_hz[fine_chan] / (1000 * 1000):.3f} MHz",
size=6,
pad=2,
)
fine_chan += 1
print("Saving figure...")
if program_args.correlator_version == mwalib.MWAVersion.CorrMWAXv2:
filename = "phase_plot1_mwax.png"
else:
filename = "phase_plot1_mwa.png"
plt.savefig(filename, bbox_inches="tight", dpi=dpi)
print(f"saved {filename}")
plt.show()
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("filename", help="fits filename")
parser.add_argument("-m", "--metafits", required=True, help="Path to the metafits file.")
parser.add_argument(
"-t1",
"--timestep1",
required=False,
help="timestep start (1 based index)",
default=1,
type=int,
)
parser.add_argument(
"-t2",
"--timestep2",
required=False,
help="timestep end (defaults to last index)",
default=-1,
type=int,
)
parser.add_argument(
"-a1",
"--ant1",
required=False,
help="antenna (start)",
default=-1,
type=int,
)
parser.add_argument(
"-a2",
"--ant2",
required=False,
help="antenna (end)",
default=-1,
type=int,
)
parser.add_argument(
"-c1",
"--channel1",
required=False,
help="fine channel number (start)",
default=-1,
type=int,
)
parser.add_argument(
"-c2",
"--channel2",
required=False,
help="fine channel number (end)",
default=-1,
type=int,
)
parser.add_argument(
"-a",
"--autosonly",
required=False,
help="Only output the auto correlations",
action="store_true",
)
parser.add_argument(
"-p",
"--ppdplot",
required=False,
help="Create a ppd plot",
action="store_true",
)
parser.add_argument(
"-p2",
"--ppdplot2",
required=False,
help="Create a ppd plot that does not sum across all baselines. ie it plots all baselines",
action="store_true",
)
parser.add_argument(
"-g",
"--gridplot",
required=False,
help="Create a grid / baseline plot",
action="store_true",
)
parser.add_argument(
"-g2",
"--gridplot2",
required=False,
help="Create a grid / baseline plot but show a single pol (XX,XY,YX,YY) for each tile. Use gridpol to specify",
action="store_true",
)
parser.add_argument(
"-gp",
"--gridpol",
required=False,
help="If gridplot2 used, use this to specify the pol. Default is 'XX'",
default="XX",
)
parser.add_argument(
"-ph",
"--phaseplot_all",
required=False,
help="Will do a phase plot for all baselines for given antennas and timesteps",
action="store_true",
)
parser.add_argument(
"-ph1",
"--phaseplot_one",
required=False,
help="Will do a phase plot for given baseline and timesteps",
action="store_true",
)
parser.add_argument(
"-o",
"--mode",
required=True,
help="How to interpret a1 and a2: RANGE or BASELINE",
)
parser.add_argument(
"-dr",
"--dump-raw",
required=False,
help="Dump the raw data",
action="store_true",
)
parser.add_argument(
"-dp",
"--dump-plot",
required=False,
help="Dump the plot data",
action="store_true",
)
args = vars(parser.parse_args())
parsed_args = ViewFITSArgs(args)
peek_fits(parsed_args)