Merge pull request #97 from lsst-camera-dh/LSSTTD-1549_CCD_signal_cor…

…relation_plots

Lssttd 1549 ccd signal correlation plots

data/BOT_jh_glob_patterns.ini

@@ -12,6 +12,7 @@ dark_current = dark_dark_*
 cte_high = sflat_flat_*H*
 cte_low = sflat_flat_*L*
 flat_pairs = flat_*_flat?_*
+raft_signal_correlations = sflat_flat_*H*
 ptc = flat_*_flat?_*
 overscan = flat_*_flat?_*
 brighter_fatter = flat_*_flat?_*

harnessed_jobs/flat_pairs_BOT/v0/producer_flat_pairs_BOT.py

@@ -3,11 +3,24 @@
 Producer script for BOT flat pairs (linearity and full-well) analysis.
 """
 import os
+from camera_components import camera_info
 from flat_pairs_jh_task import flat_pairs_jh_task
+from raft_jh_signal_correlations import raft_jh_signal_correlations
 from bot_eo_analyses import get_analysis_types, run_python_task_or_cl_script
 
 if 'linearity' in get_analysis_types():
     flat_pairs_script \
         = os.path.join(os.environ['EOANALYSISJOBSDIR'], 'harnessed_jobs',
                     'flat_pairs_BOT', 'v0', 'flat_pairs_jh_task.py')
     run_python_task_or_cl_script(flat_pairs_jh_task, flat_pairs_script)
+
+    raft_jh_signal_correlations_script \
+        =  os.path.join(os.environ['EOANALYSISJOBSDIR'], 'harnessed_jobs',
+                        'flat_pairs_BOT', 'v0',
+                        'raft_jh_signal_correlations.py')
+
+    # Run raft-level signal correlations just on science rafts for now.
+    installed_rafts = camera_info.installed_science_rafts
+    run_python_task_or_cl_script(raft_jh_signal_correlations,
+                                 raft_jh_signal_correlations_script,
+                                 device_names=installed_rafts)

harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env ipython
+"""
+Script to compute signal region correlations at the raft level.
+"""
+
+def raft_jh_signal_correlations(raft_name):
+    """JH version of raft-level signal-correlation analysis."""
+    import os
+    import logging
+    import numpy as np
+    import matplotlib.pyplot as plt
+    import siteUtils
+    from bot_eo_analyses import bias_filename, make_file_prefix, \
+        append_acq_run, glob_pattern
+    from signal_correlations import raft_level_signal_correlations
+
+    logger = logging.getLogger('raft_jh_noise_correlations')
+    logger.setLevel(logging.INFO)
+
+    # Use the high signal superflat files.
+    pattern = glob_pattern('raft_signal_correlations', f'{raft_name}_S??')
+    acq_jobname = siteUtils.getProcessName('BOT_acq')
+    sflat_files = siteUtils.dependency_glob(pattern, acq_jobname=acq_jobname)
+    folders = sorted(list(set([os.path.basename(os.path.dirname(_))
+                               for _ in sflat_files])))
+    logger.info(f'folders: {folders}')
+
+    flat1_files = dict()
+    flat2_files = dict()
+    for item in sflat_files:
+        folder = os.path.basename(os.path.dirname(item))
+        if folder not in folders[:2]:
+            continue
+        logger.info(f'item: {item}')
+        logger.info(f'folder: {folder}')
+        basename = os.path.basename(item)
+        logger.info(f'basename: {basename}')
+        slot = basename.split('_')[-1][:-len('.fits')]
+        if folder == folders[0]:
+            flat1_files[slot] = item
+        elif folder == folders[1]:
+            flat2_files[slot] = item
+    logger.info('flat pair files:')
+    for slot in flat1_files:
+        logger.info('  ' + flat1_files[slot])
+        logger.info('  ' + flat2_files[slot])
+
+    # Find the median bias files for the target raft.
+    run = siteUtils.getRunNumber()
+    bias_files = dict()
+    for slot in flat1_files.keys():
+        det_name = '_'.join((raft_name, slot))
+        bias_files[slot] = bias_filename(run, det_name)
+
+    file_prefix = make_file_prefix(run, raft_name)
+    title = append_acq_run("Imaging region correlations, "
+                           f"Run {run}, {raft_name}")
+    raft_level_signal_correlations(flat1_files, flat2_files, bias_files,
+                                   title=title)
+    plt.savefig('{}_imaging_region_correlations.png'.format(file_prefix))
+
+if __name__ == '__main__':
+    import sys
+    raft_name = sys.argv[1]
+    raft_jh_signal_correlations(raft_name)
+

python/bot_eo_validators.py

@@ -525,6 +525,14 @@ def validate_flat_pairs(results, det_names):
                                                    file_prefix,
                                                    metadata=metadata))
 
+    # Persist the raft-level imaging region correlation plots.
+    for raft in camera_info.get_installed_raft_names():
+        metadata = dict(TESTTYPE='FLAT', TEST_CATEGORY='EO', RAFT=raft, RUN=run)
+        file_prefix = make_file_prefix(run, raft)
+        filename = f'{file_prefix}_imaging_region_correlations.png'
+        results.extend(siteUtils.persist_png_files(filename, file_prefix,
+                                                   metadata=metadata))
+
     report_missing_data("validate_flat_pairs", missing_det_names)
 
     return results

python/signal_correlations.py

@@ -0,0 +1,107 @@
+import copy
+import itertools
+import numpy as np
+import matplotlib.pyplot as plt
+import astropy.visualization as viz
+from astropy.visualization.mpl_normalize import ImageNormalize
+import lsst.afw.math as afwMath
+import lsst.eotest.sensor as sensorTest
+from correlated_noise import set_ticks
+
+def diff_image_arrays(flat1_file, flat2_file, bias_frame, buffer=10):
+    """
+    Compute the difference images for each amp, using the Astier
+    weighting scheme to account for somewhat different exposure times,
+    and return a dict of the image arrays, keyed by amp.
+    """
+    image_arrays = dict()
+    ccd1 = sensorTest.MaskedCCD(flat1_file, bias_frame=bias_frame)
+    ccd2 = sensorTest.MaskedCCD(flat2_file, bias_frame=bias_frame)
+    imaging_bbox = ccd1.amp_geom.imaging
+    imaging_bbox.grow(-buffer)
+    for amp in ccd1:
+        image1 = ccd1.unbiased_and_trimmed_image(amp, imaging=imaging_bbox)
+        image2 = ccd2.unbiased_and_trimmed_image(amp, imaging=imaging_bbox)
+        mean1 = afwMath.makeStatistics(image1, afwMath.MEAN, ccd1.stat_ctrl)\
+                       .getValue()
+        mean2 = afwMath.makeStatistics(image2, afwMath.MEAN, ccd1.stat_ctrl)\
+                       .getValue()
+        fmean = (mean1 + mean2)/2.
+        image1 *= mean2/fmean
+        image2 *= mean1/fmean
+        image1 -= image2
+        image_arrays[amp] = copy.deepcopy(image1.getImage().getArray())
+    return image_arrays
+
+
+def raft_level_signal_correlations(flat1_files, flat2_files, bias_frames,
+                                   buffer=10, title='', vrange=None,
+                                   stretch=viz.LinearStretch, figsize=(8, 8)):
+
+    """
+    Compute the correlation coefficients between the imaging section
+    pixels for the difference images from a flat pair for the 144
+    amplifiers in raft.
+
+    Parameters
+    ----------
+    flat1_files: dict
+        Dictionary of flat1 image files, indexed by sensor slot id.
+        These should be from the same flat pair frame as the flat2_files.
+    flat2_files: dict
+        Dictionary of flat2 image files, indexed by sensor slot id.
+    bias_frames: dict
+        Dictionary of super bias frames, 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()
+    segments = []
+
+    ccd0 = sensorTest.MaskedCCD(list(flat1_files.values())[0])
+    bbox = ccd0.amp_geom.imaging
+    bbox.grow(-buffer)
+
+    for slot in slots:
+        if slot not in flat1_files:
+            for amp in ccd0:
+                segments.append(np.zeros((bbox.getHeight(), bbox.getWidth())))
+        else:
+            imarrs = diff_image_arrays(flat1_files[slot], flat2_files[slot],
+                                       bias_frame=bias_frames[slot],
+                                       buffer=buffer)
+            for amp in imarrs:
+                segments.append(imarrs[amp])
+    namps = len(segments)
+    data = np.array([np.corrcoef(segments[i[0]].ravel(),
+                                 segments[i[1]].ravel())[0, 1]
+                     for i in itertools.product(range(namps), range(namps))])
+    data = data.reshape((namps, namps))
+    fig = plt.figure(figsize=figsize)
+    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.ravel()))
+    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

python/stage_bot_data.py

@@ -34,6 +34,7 @@
                  'tearing_BOT': ('tearing',)}
 
 RAFT_DATA_KEYS = {'read_noise_BOT': ('raft_noise_correlations',),
+                  'flat_pairs_BOT': ('raft_signal_correlations',),
                   'tearing_BOT': ('divisadero_tearing',)}