From 7d628d8930d61fdae8b4fd0417b90e1258469cc5 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 16 Dec 2020 11:20:48 -0800 Subject: [PATCH 1/4] code to compute raft-level correlations of imaging regions of amps --- python/signal_correlations.py | 105 ++++++++++++++++++++++++++++++++++ 1 file changed, 105 insertions(+) create mode 100644 python/signal_correlations.py diff --git a/python/signal_correlations.py b/python/signal_correlations.py new file mode 100644 index 0000000..14199c1 --- /dev/null +++ b/python/signal_correlations.py @@ -0,0 +1,105 @@ +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.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 + image_arrays[amp] = copy.deepcopy((image1 - image2).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 From 6e5342ecf44d1447df774ea57f7ab510728c4e20 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 16 Dec 2020 21:03:31 -0800 Subject: [PATCH 2/4] import fixes --- python/signal_correlations.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/python/signal_correlations.py b/python/signal_correlations.py index 14199c1..ff995d1 100644 --- a/python/signal_correlations.py +++ b/python/signal_correlations.py @@ -4,8 +4,9 @@ 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 +from correlated_noise import set_ticks def diff_image_arrays(flat1_file, flat2_file, bias_frame, buffer=10): """ @@ -28,8 +29,8 @@ def diff_image_arrays(flat1_file, flat2_file, bias_frame, buffer=10): fmean = (mean1 + mean2)/2. image1 *= mean2/fmean image2 *= mean1/fmean - image_arrays[amp] = copy.deepcopy((image1 - image2).getImage() - .getArray()) + image1 -= image2 + image_arrays[amp] = copy.deepcopy(image1.getImage().getArray()) return image_arrays @@ -38,8 +39,9 @@ def raft_level_signal_correlations(flat1_files, flat2_files, bias_frames, 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. + Compute the correlation coefficients between the imaging section + pixels for the difference images from a flat pair for the 144 + amplifiers in raft. Parameters ---------- From 5d310a26307727ef77a153d4c348b386ecda1ce2 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Thu, 17 Dec 2020 12:57:15 -0800 Subject: [PATCH 3/4] add code to flat_pairs_BOT harnessed job to produce raft-level imaging region correlation plots --- .../v0/producer_flat_pairs_BOT.py | 13 ++++ .../v0/raft_jh_signal_correlations.py | 76 +++++++++++++++++++ python/bot_eo_validators.py | 8 ++ 3 files changed, 97 insertions(+) create mode 100644 harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py diff --git a/harnessed_jobs/flat_pairs_BOT/v0/producer_flat_pairs_BOT.py b/harnessed_jobs/flat_pairs_BOT/v0/producer_flat_pairs_BOT.py index b9207e0..1cc6ac6 100755 --- a/harnessed_jobs/flat_pairs_BOT/v0/producer_flat_pairs_BOT.py +++ b/harnessed_jobs/flat_pairs_BOT/v0/producer_flat_pairs_BOT.py @@ -3,7 +3,9 @@ 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(): @@ -11,3 +13,14 @@ = 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) diff --git a/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py b/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py new file mode 100644 index 0000000..487bff9 --- /dev/null +++ b/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py @@ -0,0 +1,76 @@ +#!/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) + + # Find all of the flat files, and select the pair with signal + # levels closest to 50k e-/pixel. + pattern = glob_pattern('flat_pairs', f'{raft_name}_*') + acq_jobname = siteUtils.getProcessName('BOT_acq') + flat_files = siteUtils.dependency_glob(pattern, acq_jobname=acq_jobname) + + target_signal = 5e4 + selected = None + min_diff = None + for item in flat_files: + frame = os.path.basename(os.path.dirname(item)) + signal = frame.split('_')[3] + diff = np.abs(target_signal - float(signal)) + if selected is None or diff < min_diff : + selected = signal + min_diff = diff + flat1_files = dict() + flat2_files = dict() + for item in flat_files: + folder = os.path.basename(os.path.dirname(item)) + basename = os.path.basename(item) + if '_' + signal + '_' in folder: + if raft_name in basename: + slot = basename.split('_')[-1][:-len('.fits')] + if 'flat0' in folder: + flat1_files[slot] = item + elif 'flat1' in folder: + flat2_files[slot] = item + else: + raise RuntimeError('unmatched flat pair file') + logger.info('flat pair files:') + for slot in flat1_files: + logger.info(' ' + flat1_files[slot]) + logger.info(' ' + flat2_files[slot]) + logger.info('') + + # Find the median bias files for the target raft. + run = siteUtils.getRunNumber() + bias_files = dict() + logger.info('median bias files:') + for slot in flat1_files.keys(): + det_name = '_'.join((raft_name, slot)) + bias_files[slot] = bias_filename(run, det_name) + logger.info(' ' + bias_files[slot]) + + 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) + diff --git a/python/bot_eo_validators.py b/python/bot_eo_validators.py index d1f208b..9eb72c4 100644 --- a/python/bot_eo_validators.py +++ b/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 From 5838e57a9d971656712fa5b46f746114b5f6de80 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Thu, 17 Dec 2020 17:57:21 -0800 Subject: [PATCH 4/4] use high signal superflat files instead of flat pairs for imaging region correlation calculations --- data/BOT_jh_glob_patterns.ini | 1 + .../v0/raft_jh_signal_correlations.py | 48 ++++++++----------- python/stage_bot_data.py | 1 + 3 files changed, 21 insertions(+), 29 deletions(-) diff --git a/data/BOT_jh_glob_patterns.ini b/data/BOT_jh_glob_patterns.ini index 355a9f0..802ad3c 100644 --- a/data/BOT_jh_glob_patterns.ini +++ b/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?_* diff --git a/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py b/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py index 487bff9..09e6950 100644 --- a/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py +++ b/harnessed_jobs/flat_pairs_BOT/v0/raft_jh_signal_correlations.py @@ -17,50 +17,40 @@ def raft_jh_signal_correlations(raft_name): logger = logging.getLogger('raft_jh_noise_correlations') logger.setLevel(logging.INFO) - # Find all of the flat files, and select the pair with signal - # levels closest to 50k e-/pixel. - pattern = glob_pattern('flat_pairs', f'{raft_name}_*') + # Use the high signal superflat files. + pattern = glob_pattern('raft_signal_correlations', f'{raft_name}_S??') acq_jobname = siteUtils.getProcessName('BOT_acq') - flat_files = siteUtils.dependency_glob(pattern, acq_jobname=acq_jobname) + 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}') - target_signal = 5e4 - selected = None - min_diff = None - for item in flat_files: - frame = os.path.basename(os.path.dirname(item)) - signal = frame.split('_')[3] - diff = np.abs(target_signal - float(signal)) - if selected is None or diff < min_diff : - selected = signal - min_diff = diff - flat1_files = dict() - flat2_files = dict() - for item in flat_files: + 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) - if '_' + signal + '_' in folder: - if raft_name in basename: - slot = basename.split('_')[-1][:-len('.fits')] - if 'flat0' in folder: - flat1_files[slot] = item - elif 'flat1' in folder: - flat2_files[slot] = item - else: - raise RuntimeError('unmatched flat pair file') + 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]) - logger.info('') # Find the median bias files for the target raft. run = siteUtils.getRunNumber() bias_files = dict() - logger.info('median bias files:') for slot in flat1_files.keys(): det_name = '_'.join((raft_name, slot)) bias_files[slot] = bias_filename(run, det_name) - logger.info(' ' + bias_files[slot]) file_prefix = make_file_prefix(run, raft_name) title = append_acq_run("Imaging region correlations, " diff --git a/python/stage_bot_data.py b/python/stage_bot_data.py index 02d3a17..e8967d0 100644 --- a/python/stage_bot_data.py +++ b/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',)}