Merge pull request #13 from lsst-camera-dh/LSSTTD-1271_poc_request

Add plot of overscan correlation coefficients for 144 amps in raft

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))

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()

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()