Source code for pipeline.hsd.tasks.skycal.display

import os
import matplotlib.pyplot as plt
import numpy
import traceback

import pipeline.infrastructure.logging as logging
import pipeline.infrastructure.renderer.logger as logger
from pipeline.h.tasks.common.displays import common as common
from pipeline.h.tasks.common.displays import bandpass as bandpass
from ..common import display as sd_display
from pipeline.infrastructure import casa_tasks
from pipeline.infrastructure import casa_tools
from pipeline.infrastructure.displays.plotstyle import casa5style_plot

LOG = logging.get_logger(__name__)


[docs]class SingleDishSkyCalDisplayBase(object):
[docs] def init_with_field(self, context, result, field): vis = self._vis ms = context.observing_run.get_ms(vis) num_fields = len(ms.fields) if field.isdigit() and int(field) < num_fields: self.field_id = int(field) self.field_name = ms.get_fields(self.field_id)[0].clean_name else: self.field_name = field fields = ms.get_fields(name=field) assert len(fields) == 1 self.field_id = fields[0].id LOG.debug('field: ID %s Name \'%s\''%(self.field_id, self.field_name)) old_prefix = self._figroot.replace('.png', '') self._figroot = self._figroot.replace('.png', '-%s.png' % (self.field_name)) new_prefix = self._figroot.replace('.png', '') self._update_figfile(old_prefix, new_prefix) if 'field' not in self._kwargs: self._kwargs['field'] = self.field_id
[docs] def add_field_identifier(self, plots): for plot in plots: if 'field' not in plot.parameters: plot.parameters['field'] = self.field_name
def _update_figfile(self): raise NotImplementedError()
[docs]class SingleDishSkyCalAmpVsFreqSummaryChart(common.PlotbandpassDetailBase, SingleDishSkyCalDisplayBase): def __init__(self, context, result, field): super(SingleDishSkyCalAmpVsFreqSummaryChart, self).__init__(context, result, 'freq', 'amp', showatm=True, overlay='antenna', solutionTimeThresholdSeconds=3600.) # self._figfile structure: {spw_id: {antenna_id: filename}} self.spw_ids = list(self._figfile.keys()) # take any value from the dictionary self._figfile = dict((spw_id, list(self._figfile[spw_id].values())[0]) for spw_id in self.spw_ids) self.init_with_field(context, result, field)
[docs] def plot(self): missing = [spw_id for spw_id in self.spw_ids if not os.path.exists(self._figfile[spw_id])] if missing: LOG.trace('Executing new plotbandpass job for missing figures') for spw_id in missing: # PIPE-110: show image sideband for DSB receivers. showimage = self._rxmap.get(spw_id, "") == "DSB" try: task = self.create_task(spw_id, '', showimage=showimage) task.execute(dry_run=False) except Exception as ex: LOG.error('Could not create plotbandpass summary plots') LOG.exception(ex) return None wrappers = [] for spw_id, figfile in self._figfile.items(): # PIPE-110: show image sideband for DSB receivers. showimage = self._rxmap.get(spw_id, "") == "DSB" if os.path.exists(figfile): task = self.create_task(spw_id, '', showimage=showimage) wrapper = logger.Plot(figfile, x_axis=self._xaxis, y_axis=self._yaxis, parameters={'vis': self._vis_basename, 'spw': spw_id, 'field': self.field_name}, command=str(task)) wrappers.append(wrapper) else: LOG.trace('No plotbandpass summary plot found for spw ' '%s' % spw_id) return wrappers
def _update_figfile(self, old_prefix, new_prefix): for spw_id, figfile in self._figfile.items(): self._figfile[spw_id] = figfile.replace(old_prefix, new_prefix) spw_indicator = 'spw{}'.format(spw_id) pieces = self._figfile[spw_id].split('.') try: spw_index = pieces.index(spw_indicator) except: spw_index = -3 # remove antenna name from the filename pieces.pop(spw_index - 1) self._figfile[spw_id] = '.'.join(pieces)
[docs]class SingleDishSkyCalAmpVsFreqDetailChart(bandpass.BandpassDetailChart, SingleDishSkyCalDisplayBase): def __init__(self, context, result, field): super(SingleDishSkyCalAmpVsFreqDetailChart, self).__init__( context, result, xaxis='freq', yaxis='amp', showatm=True, overlay='time') self.init_with_field(context, result, field)
[docs] def plot(self): wrappers = super(SingleDishSkyCalAmpVsFreqDetailChart, self).plot() self.add_field_identifier(wrappers) return wrappers
def _update_figfile(self, old_prefix, new_prefix): for spw_id in self._figfile: for antenna_id, figfile in self._figfile[spw_id].items(): new_figfile = figfile.replace(old_prefix, new_prefix) self._figfile[spw_id][antenna_id] = new_figfile
[docs]class SingleDishPlotmsLeaf(object): """ Class to execute plotms and return a plot wrapper. Task arguments for plotms is customized for single dish usecase. """ def __init__(self, context, result, calapp, xaxis, yaxis, spw='', ant='', coloraxis='', **kwargs): LOG.debug('__init__(caltable={caltable}, spw={spw}, ant={ant})'.format(caltable=calapp.gaintable, spw=spw, ant=ant)) self.xaxis = xaxis self.yaxis = yaxis self.field = calapp.gainfield self.caltable = calapp.gaintable self.vis = calapp.vis self.spw = str(spw) self.antenna = str(ant) self.coloraxis = coloraxis ms = context.observing_run.get_ms(self.vis) num_fields = len(ms.fields) if self.field.isdigit() and int(self.field) < num_fields: self.field_id = int(self.field) self.field_name = ms.get_fields(self.field_id)[0].clean_name else: self.field_name = self.field fields = ms.get_fields(name=self.field) assert len(fields) == 1 self.field_id = fields[0].id LOG.debug('field: ID %s Name \'%s\'' % (self.field_id, self.field_name)) self.antmap = dict((a.id, a.name) for a in ms.antennas) if len(self.antenna) == 0: self.antenna_selection = 'summary' else: self.antenna_selection = list(self.antmap.values())[int(self.antenna)] LOG.info('antenna: ID %s Name \'%s\'' % (self.antenna, self.antenna_selection)) # self.antenna_selection = '*&&&' self._figroot = os.path.join(context.report_dir, 'stage%s' % result.stage_number)
[docs] def plot(self): prefix = '{caltable}-{y}_vs_{x}-{field}-{ant}-spw{spw}'.format( caltable=os.path.basename(self.caltable), y=self.yaxis, x=self.xaxis, field=self.field_name, ant=self.antenna_selection, spw=self.spw) title = '{caltable} \nField "{field}" Antenna {ant} Spw {spw} \ncoloraxis={caxis}'.format( caltable=os.path.basename(self.caltable), field=self.field_name, ant=self.antenna_selection, spw=self.spw, caxis=self.coloraxis) figfile = os.path.join(self._figroot, '{prefix}.png'.format(prefix=prefix)) task = self._create_task(title, figfile) try: if os.path.exists(figfile): LOG.debug('Returning existing plot') else: task.execute() plot_objects = [self._get_plot_object(figfile, task)] except Exception as e: LOG.error(str(e)) LOG.debug(traceback.format_exc()) LOG.error('Failed to generate plot "{}"'.format(figfile)) plot_objects = [] return plot_objects
def _create_task(self, title, figfile): task_args = {'vis': self.caltable, 'xaxis': self.xaxis, 'yaxis': self.yaxis, 'plotfile': figfile, 'coloraxis': self.coloraxis, 'showgui': False, 'field': self.field, 'spw': self.spw, 'antenna': self.antenna, 'title': title, 'showlegend': True, 'averagedata': True, 'avgchannel': '1e8'} return casa_tasks.plotms(**task_args) def _get_plot_object(self, figfile, task): parameters = {'vis': os.path.basename(self.vis), 'ant': self.antenna_selection, 'spw': self.spw, 'field': self.field_name} return logger.Plot(figfile, x_axis='Time', y_axis='Amplitude', parameters=parameters, command=str(task))
[docs]class SingleDishPlotmsAntComposite(common.AntComposite): leaf_class = SingleDishPlotmsLeaf
[docs]class SingleDishPlotmsSpwComposite(common.SpwComposite): leaf_class = SingleDishPlotmsLeaf
[docs]class SingleDishPlotmsAntSpwComposite(common.AntSpwComposite): leaf_class = SingleDishPlotmsSpwComposite
[docs]class SingleDishSkyCalAmpVsTimeSummaryChart(SingleDishPlotmsSpwComposite): def __init__(self, context, result, calapp): super(SingleDishSkyCalAmpVsTimeSummaryChart, self).__init__(context, result, calapp, xaxis='time', yaxis='amp', coloraxis='ant1')
[docs]class SingleDishSkyCalAmpVsTimeDetailChart(SingleDishPlotmsAntSpwComposite): def __init__(self, context, result, calapp): super(SingleDishSkyCalAmpVsTimeDetailChart, self).__init__(context, result, calapp, xaxis='time', yaxis='amp', coloraxis='corr')
[docs]class SingleDishSkyCalIntervalVsTimeDisplay(common.PlotbandpassDetailBase, SingleDishSkyCalDisplayBase): """ Class to execute pyplot and return a plot (figure) of Interval vs. Time. If figtype='summary', the first spw is used, while all spw are used if figtype='detail'. """ def __init__(self, context, result, calapp, figtype=''): self.context = context self.result = result self.calapp = calapp self.figtype = str(figtype) LOG.info('figtype = {0}'.format(self.figtype)) @casa5style_plot def plot(self): context = self.context result = self.result calapp = self.calapp figtype = self.figtype vis = os.path.basename(result.inputs['vis']) ms = context.observing_run.get_ms(vis) antennas = ms.antennas fields = ms.fields science_spw = ms.get_spectral_windows(science_windows_only=True) if figtype == 'detail': spw_ids = [spw.id for spw in science_spw] elif figtype == 'summary': spw_ids = [science_spw[0].id] else: pass plot = None plots = [] antenna_ids = [ant.id for ant in antennas] field_strategy = ms.calibration_strategy['field_strategy'] for spw in spw_ids: spw_id = spw LOG.debug('spw_id={0}'.format(spw_id)) for antenna_id in antenna_ids: LOG.debug('antenna_id={0}'.format(antenna_id)) n = 0 for field_id_target, field_id_reference in field_strategy.items(): LOG.debug('field_id_target = {0}, field_id_reference = {1}'.format(field_id_target, field_id_reference)) field = fields[field_id_target] # make plots for the interval ratio (off-source/on-source) vs time; with casa_tools.TableReader(calapp.gaintable) as tb: t = tb.query('SPECTRAL_WINDOW_ID=={}&&ANTENNA1=={}&&FIELD_ID=={}'.format(spw_id, antenna_id, field_id_reference), sortlist='TIME', columns='TIME, SPECTRAL_WINDOW_ID, INTERVAL') mjd = t.getcol('TIME') ms_target = ms.get_scans(scan_intent='TARGET', field=field_id_target) ms_target0 = ms_target[0] interval_unit = ms_target0.mean_interval(spw_id=spw_id) interval_unit = interval_unit.total_seconds() interval = t.getcol('INTERVAL') / interval_unit t.close() if len(mjd) == 0: pass else: mjd_list = mjd.tolist() mjd_list = [v/86400.0 for v in mjd_list] mjd_list = sd_display.mjd_to_plotval(mjd_list) start_time = numpy.min([numpy.min(x) for x in mjd_list if len(mjd_list) > 0]) end_time = numpy.max([numpy.max(x) for x in mjd_list if len(mjd_list) > 0]) fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.xaxis.set_major_locator(sd_display.utc_locator(start_time=start_time, end_time=end_time)) ax.xaxis.set_major_formatter(sd_display.utc_formatter()) ax.tick_params( axis='both', labelsize=10 ) antenna_name = antennas[antenna_id].name field_name = field.clean_name plt.title('Interval vs. Time Plot\n{} Field:{} Antenna:{} Spw:{}'.format(vis, field_name, antenna_name, spw_id), fontsize=12) plt.ylabel('Interval of Off-Source / Interval of On-Source', fontsize=10) plt.xlabel("UTC", fontsize=10) ax.plot(mjd_list, interval, linestyle='None', marker=".", label="Interval of Off-Source\nUnit: {} seconds (Interval of On-Source)".format(interval_unit)) min_interval = numpy.min(interval) max_interval = numpy.max(interval) ax.set_ylim([min_interval-3.0, max_interval+3.0]) plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=10) if figtype == "summary": prefix = vis + "_" + '{}'.format(antenna_name)+ "_" + '{}'.format(field_name) + "_summary_hsd_skycal_offinterval" else: prefix = vis + "_" + '{}'.format(antenna_name)+ "_" + '{}'.format(field_name) + "_spw" + '{}'.format(spw_id) + "_hsd_skycal_offinterval" figroot = os.path.join(context.report_dir, 'stage%s' % result.stage_number) figpath = os.path.join(figroot, '{prefix}.png'.format(prefix=prefix)) LOG.info('Plot of Interval vs Time: figpath = {0}'.format(figpath)) plt.savefig(figpath) plt.close() if os.path.exists(figpath): parameters = { 'spw': spw_id, 'ant': antenna_name, 'vis': vis, 'type': 'Plot of Interval vs. Time', 'file': vis} plot = logger.Plot(figpath, x_axis='Time', y_axis='Off-Source Interval / On-Source Interval', field=field_name, parameters=parameters) plots.append(plot) n += 1 return plots
@casa5style_plot def plot_elevation_difference(context, result, eldiff, threshold=3.0): """ context result eldiff -- dictionary whose value is ElevationDifference named tuple instance that holds 'timeon': timestamp for ON-SOURCE pointings 'elon': ON-SOURCE elevation 'timecal': timestamp for OFF-SOURCE pointings 'elcal': OFF-SOURCE elevation 'time0': timestamp for preceding OFF-SOURCE pointings 'eldiff0': elevation difference between ON-SOURCE and preceding OFF-SOURCE 'time1': timestamp for subsequent OFF-SOURCE pointings 'eldiff1': elevation difference between ON-SOURCE and subsequent OFF-SOURCE eldiff is a nested dictionary whose first key is FIELD_ID for target, the second one is ANTENNA_ID, and the third one is SPW_ID. threhshold -- Elevation threshold for QA (default 3deg) """ calapp = result.final[0] vis = calapp.calto.vis ms = context.observing_run.get_ms(vis) figroot = os.path.join(context.report_dir, 'stage%s' % result.stage_number) figure0 = 'PERANTENNA_PLOT' figure1 = 'ALLANTENNA_PLOT' start_time = numpy.min([numpy.min(x.timeon) for z in eldiff.values() for y in z.values() for x in y.values() if len(x.timeon) > 0]) end_time = numpy.max([numpy.max(x.timeon) for z in eldiff.values() for y in z.values() for x in y.values() if len(x.timeon) > 0]) def init_figure(figure_id): plt.figure(figure_id) plt.clf() a0 = plt.axes([0.125, 0.51, 0.775, 0.39]) a0.xaxis.set_major_locator(sd_display.utc_locator(start_time=start_time, end_time=end_time)) a0.xaxis.set_major_formatter(plt.NullFormatter()) a0.tick_params( axis='both', labelsize=10 ) plt.ylabel('Elevation [deg]', fontsize=10) a1 = plt.axes([0.125, 0.11, 0.775, 0.39]) a1.xaxis.set_major_locator(sd_display.utc_locator(start_time=start_time, end_time=end_time)) a1.xaxis.set_major_formatter(sd_display.utc_formatter()) a1.tick_params( axis='both', labelsize=10 ) plt.ylabel('Elevation Difference [deg]', fontsize=10) plt.xlabel('UTC', fontsize=10) return a0, a1 def finalize_figure(figure_id, vis, field_name, antenna_name): figure = plt.figure(figure_id) axes = figure.axes assert len(axes) == 2 a0 = axes[0] a1 = axes[1] plt.gcf().sca(a1) ymin, ymax = plt.ylim() dy = ymax - ymin plt.ylim([0, max(ymax + 0.1 * dy, threshold + 0.1)]) plt.axhline(threshold, color='red') xmin, xmax = plt.xlim() dx = xmax - xmin x = xmax - 0.01 * dx y = threshold - 0.05 plt.text(x, y, 'QA threshold', ha='right', va='top', color='red', size='small') plt.gcf().sca(a0) labelon = False labeloff = False for l in a0.lines: if (labelon is False) and (l.get_color() == 'black'): l.set_label('ON') labelon = True if (labeloff is False) and (l.get_color() == 'blue'): l.set_label('OFF') labeloff = True if labelon and labeloff: break plt.legend(loc='best', numpoints=1, prop={'size': 'small'}) plt.title('Elevation Difference between ON and OFF\n{} Field {} Antenna {}'.format(vis, field_name, antenna_name), fontsize=12) def generate_plot(figure_id, vis, field_name, antenna_name): plt.figure(figure_id) vis_prefix = '.'.join(vis.split('.')[:-1]) figfile = 'elevation_difference_{}_{}_{}.png'.format(vis_prefix, field_name, antenna_name) figpath = os.path.join(figroot, figfile) #LOG.info('figpath={}'.format(figpath)) plt.savefig(figpath) plot = None if os.path.exists(figpath): parameters = {'intent': 'TARGET', 'spw': '', 'pol': '', 'ant': antenna_name, 'vis': vis, 'type': 'Elevation Difference vs. Time', 'file': vis} plot = logger.Plot(figpath, x_axis='Time', y_axis='Elevation Difference', field=field_name, parameters=parameters) return plot def close_figure( figure_id ): plt.close(figure_id) plots = [] for field_id, eldiff_field in eldiff.items(): # figure for summary plot a2, a3 = init_figure(figure1) field = ms.fields[field_id] field_name = field.clean_name plots_per_field = [] for antenna_id, eldant in eldiff_field.items(): # figure for per-antenna plots a0, a1 = init_figure(figure0) antenna_name = ms.antennas[antenna_id].name for spw_id, eld in eldant.items(): if len(eld.timeon) == 0 or len(eld.timecal) == 0: continue # Elevation vs. Time xon = sd_display.mjd_to_plotval(eld.timeon) xoff = sd_display.mjd_to_plotval(eld.timecal) for a in [a0, a2]: a.plot(xon, eld.elon, '.', color='black', mew=0) a.plot(xoff, eld.elcal, '.-', color='blue', mew=0) # Elevation Difference vs. Time time0 = eld.time0 eldiff0 = eld.eldiff0 time1 = eld.time1 eldiff1 = eld.eldiff1 io0 = numpy.where(numpy.abs(eldiff0) < threshold)[0] ix0 = numpy.where(numpy.abs(eldiff0) >= threshold)[0] io1 = numpy.where(numpy.abs(eldiff1) < threshold)[0] ix1 = numpy.where(numpy.abs(eldiff1) >= threshold)[0] index_list = [io0, io1, ix0, ix1] time_list = [time0, time1, time0, time1] eldiff_list = [eldiff0, eldiff1, eldiff0, eldiff1] style_list = ['g.', 'g.', 'rs', 'rs'] for idx, t, ed, ls in zip(index_list, time_list, eldiff_list, style_list): if len(idx) > 0: x = sd_display.mjd_to_plotval(t[idx]) y = numpy.abs(ed[idx]) for a in [a1, a3]: a.plot(x, y, ls, mew=0) # finalize figure for per-antenna plot finalize_figure(figure0, ms.basename, field_name, antenna_name) # generate plot object plot = generate_plot(figure0, ms.basename, field_name, antenna_name) if plot is not None: plots_per_field.append(plot) # finalize figure for summary plot finalize_figure(figure1, ms.basename, field_name, 'ALL') # generate plot object plot = generate_plot(figure1, ms.basename, field_name, '') if plot is not None: plots_per_field.append(plot) plots.extend(plots_per_field) close_figure( figure0 ) close_figure( figure1 ) return plots