diff --git a/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py b/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py index c091491..89ccd95 100755 --- a/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py +++ b/harnessed_jobs/read_noise_raft/v0/producer_read_noise_raft.py @@ -3,12 +3,14 @@ Producer script for raft-level read noise analysis. """ from __future__ import print_function +import os import matplotlib.pyplot as plt import lsst.eotest.sensor as sensorTest import siteUtils import eotestUtils -from correlated_noise import correlated_noise +from correlated_noise import correlated_noise, raft_level_oscan_correlations from multiprocessor_execution import sensor_analyses +import camera_components def run_read_noise_task(sensor_id): @@ -30,10 +32,32 @@ def run_read_noise_task(sensor_id): mask_files=mask_files, use_overscan=True) # Compute amp-amp correlated noise. - bias_stats, corr_fig, _ = correlated_noise(bias_files, target=0, - make_plots=True, title=sensor_id) + _, corr_fig, _ = correlated_noise(bias_files, target=0, + make_plots=True, title=sensor_id) plt.figure(corr_fig.number) plt.savefig('%s_correlated_noise.png' % file_prefix) + +def get_bias_files(raft_id=None): + """Get the bias files from the Fe55 acquisition.""" + if raft_id is None: + raft_id = os.environ['LCATR_UNIT_ID'] + raft = camera_components.Raft.create_from_etrav(raft_id) + bias_files = dict() + for slot, sensor_id in raft.items(): + bias_files[slot] \ + = siteUtils.dependency_glob('S*/%s_fe55_bias_*.fits' % sensor_id, + jobname=siteUtils.getProcessName('fe55_raft_acq'), + description='Bias files for noise correlations:')[0] + return bias_files + if __name__ == '__main__': sensor_analyses(run_read_noise_task) + + raft_id = os.environ['LCATR_UNIT_ID'] + run = siteUtils.getRunNumber() + bias_files = get_bias_files(raft_id) + title = 'Overscan correlations, {}, Run {}'.format(raft_id, run) + plt.rcParams['figure.figsize'] = (8, 8) + raft_level_oscan_correlations(bias_files, title=title) + plt.savefig('{}_{}_overscan_correlations.png'.format(raft_id, run)) diff --git a/harnessed_jobs/read_noise_raft/v0/validator_read_noise_raft.py b/harnessed_jobs/read_noise_raft/v0/validator_read_noise_raft.py index 7ff4e4f..18b0f22 100755 --- a/harnessed_jobs/read_noise_raft/v0/validator_read_noise_raft.py +++ b/harnessed_jobs/read_noise_raft/v0/validator_read_noise_raft.py @@ -52,6 +52,11 @@ sensor_id, folder=slot, metadata=metadata)) +# Persist the raft-level overscan correlation plot. +metadata = dict(LSST_NUM=raft_id, TESTTYPE='FE55', TEST_CATEGORY='EO') +results.extend(siteUtils.persist_png_files('%s*.png' % raft_id, raft_id, + metadata=metadata)) + results.extend(siteUtils.jobInfo()) lcatr.schema.write_file(results) lcatr.schema.validate_file() diff --git a/python/correlated_noise.py b/python/correlated_noise.py index 1f81a1c..d0c4fae 100644 --- a/python/correlated_noise.py +++ b/python/correlated_noise.py @@ -1,15 +1,20 @@ import glob from collections import namedtuple, defaultdict +import itertools import numpy as np import scipy import astropy.io.fits as fits +import astropy.visualization as viz +from astropy.visualization.mpl_normalize import ImageNormalize import matplotlib.pyplot as plt +import matplotlib.ticker as ticker import lsst.eotest.image_utils as imutils import lsst.eotest.sensor as sensorTest +import camera_components plt.rcParams['xtick.labelsize'] = 'x-small' plt.rcParams['ytick.labelsize'] = 'x-small' - +plt.rcParams['image.cmap'] = 'jet' def get_oscan_indices(target_file): "Return the pixel indices of the overscan region." @@ -118,7 +123,7 @@ def correlated_noise(bias_files, target=0, make_plots=False, plot_corr=True, reduced_mean_oscan = np.zeros(oscan_shape) num_oscan = 0 for oamp, oscan in bias_oscans.items(): - if oamp == amp or not (4. < np.std(oscan) < 25): + if oamp == amp: continue reduced_mean_oscan += (oscan - mean_oscans[oamp]) num_oscan += 1 @@ -155,6 +160,74 @@ def correlated_noise(bias_files, target=0, make_plots=False, plot_corr=True, return bias_stats, f1, f2 +def raft_level_oscan_correlations(bias_files, buffer=10, title='', + vrange=None, stretch=viz.LinearStretch): + """ + Compute the correlation coefficients between the overscan pixels + of the 144 amplifiers in raft. + + Parameters + ---------- + bias_files: dict + Dictionary of bias image files, indexed by sensor slot id. + buffer: int [10] + Buffer region around perimeter of serial overscan region to + avoid when computing the correlation coefficients. + title: str [''] + Plot title. + vrange: (float, float) [None] + Minimum and maximum values for color scale range. If None, then + the range of the central 98th percentile of the absolute value + of the data is used. + stretch: astropy.visualization.BaseStretch [LinearStretch] + Stretch to use for the color scale. + + Returns + ------- + (matplotlib.figure.Figure, np.array): The figure containing the plot and + the numpy array containing the correlation coefficients. + """ + slots = 'S00 S01 S02 S10 S11 S12 S20 S21 S22'.split() + bbox = None + overscans = [] + for slot in slots: + ccd = sensorTest.MaskedCCD(bias_files[slot]) + if bbox is None: + bbox = ccd.amp_geom.serial_overscan + bbox.grow(-buffer) + for amp in ccd: + image = ccd[amp].getImage() + overscans.append(image.Factory(image, bbox).getArray()) + namps = len(overscans) + data = np.array([np.corrcoef(overscans[i[0]].ravel(), + overscans[i[1]].ravel())[0, 1] + for i in itertools.product(range(namps), range(namps))]) + data = data.reshape((namps, namps)) + fig = plt.figure() + ax = fig.add_subplot(111) + ax.set_title(title, fontsize='medium') + + interval = viz.PercentileInterval(98.) + if vrange is None: + vrange = interval.get_limits(np.abs(data.flatten())) + norm = ImageNormalize(vmin=vrange[0], vmax=vrange[1], stretch=stretch()) + image = ax.imshow(data, interpolation='none', norm=norm) + plt.colorbar(image) + + set_ticks(ax, slots, amps=16) + + return fig, data + +def set_ticks(ax, slots, amps=16): + """Set the tick labels, centering the slot names between amps 1 and 16.""" + major_locs = [i*amps - 0.5 for i in range(len(slots) + 1)] + minor_locs = [amps//2 + i*amps for i in range(len(slots))] + for axis in (ax.xaxis, ax.yaxis): + axis.set_tick_params(which='minor', length=0) + axis.set_major_locator(ticker.FixedLocator(major_locs)) + axis.set_major_formatter(ticker.FixedFormatter(['']*len(major_locs))) + axis.set_minor_locator(ticker.FixedLocator(minor_locs)) + axis.set_minor_formatter(ticker.FixedFormatter(slots)) if __name__ == '__main__': plt.ion()