# What version of the Stack are we using?
! echo $IMAGE_DESCRIPTION
! eups list -s | grep lsst_distrib

Release r29.2.0 (RSP Build 2244)
lsst_distrib          gc675d380bf+f75de59d28 	current o_latest v29_2_0 v29_2_0_rc1 setup

import numpy as np
import matplotlib.pyplot as plt
import lsst.afw.display as afwDisplay

from lsst.daf.butler import Butler

from lsst.ip.isr import isrFunctions
import lsst.afw.math as afwMath
from lsst.meas.algorithms import backgroundFlatContext

afwDisplay.setDefaultBackend("matplotlib")

# Always use this repo and collection in DP1
repo = 'dp1'
collections = 'LSSTComCam/DP1'

butler = Butler(repo, collections=collections)

dataId = {'instrument': 'LSSTComCam', 'detector': 8, 'visit': 2024113000204, 'band': 'u', 'day_obs': 20241130, 'physical_filter': 'u_02'}

Grab the configs used by ReprocessVisitImageTask and IsrTaskLSST

reprocessVisitImage_config = butler.get('reprocessVisitImage_config', dataId=dataId)
isr_config = butler.get('isr_config', dataId=dataId)

If this is True, an illuminationCorrection was applied to the visit_image.

reprocessVisitImage_config.do_apply_flat_background_ratio

True

illuminationCorrection = butler.get('illuminationCorrection', dataId=dataId)

If these two things are False, we don't need to worry about some other potential complications.

reprocessVisitImage_config.remove_initial_photo_calib, reprocessVisitImage_config.do_use_sky_corr

(False, False)

Grab the visit_image, and its background.

visit_image = butler.get('visit_image', dataId=dataId)

visit_image_background = butler.get('visit_image_background', dataId=dataId)

If this is "electron", then the relevant steps of ISR assume a gain of 1.0.

isElectrons = (visit_image.metadata["LSST ISR UNITS"] == "electron")
print('LSST ISR UNITS:', visit_image.metadata['LSST ISR UNITS'])

LSST ISR UNITS: electron

Items needed to invert the flat correction

maskedImage = visit_image.maskedImage

flat = butler.get('flat', dataId=dataId)
flatMaskedImage = flat.maskedImage

scalingType = isr_config.flatScalingType
userScale = isr_config.flatUserScale

Pathological pixels have known sources of disagreement between their variance and flux. We keep track of these using the mask plane.

maskPlaneDict = maskedImage.getMask().getMaskPlaneDict()

bad_mask_names = ['BAD', 'CR', 'CROSSTALK', 'INTRP', 'NO_DATA', 'SAT', 'SUSPECT', 'UNMASKEDNAN']

Illustrations of some of the things we're working with.

plt.imshow(np.arcsinh(maskedImage.image.array), origin='lower')
plt.show()

plt.imshow(np.arcsinh(maskedImage.variance.array), origin='lower')
plt.show()

plt.imshow(flat.image.array, origin='lower')
plt.colorbar()
plt.show()

plt.imshow(illuminationCorrection.image.array, origin='lower')
plt.colorbar()
plt.show()

plt.imshow(visit_image_background.getImage().array, origin='lower')
plt.colorbar()
plt.show()

Do these steps in the reverse order in which they're applied in ReprocessVisitImageTask.run(), so that we effectively work backward from a visit_image to a post_isr_image.

background_to_photometric_ratio = illuminationCorrection.image.clone()

# Restore the visit image background
with backgroundFlatContext(
    maskedImage,
    reprocessVisitImage_config.do_apply_flat_background_ratio,
    backgroundToPhotometricRatio=background_to_photometric_ratio,
):
    maskedImage += visit_image_background.getImage()

# Invert the illumination correction
if reprocessVisitImage_config.do_apply_flat_background_ratio:
    maskedImage *= background_to_photometric_ratio

Now look at IsrTaskLSST.run() and work backward to just after the point where the variance was computed.

# Invert the flat correction
isrFunctions.flatCorrection(
    maskedImage=maskedImage,
    flatMaskedImage=flatMaskedImage,
    scalingType=scalingType,
    userScale=userScale,
    invert=True
)

At this stage, the variance in each amplifier region should equal the image flux in each pixel (if everything is in electron units), plus a constant readnoise value. We can check the extent to which that's true by seeing if the variance vs. flux in each pixel obeys a straight-line relationship.

amps = visit_image.getDetector().getAmplifiers()
amp_bboxes = [amp.getBBox() for amp in amps]

gn = maskedImage.image.clone()
rn = maskedImage.image.clone()
gains, readnoises = [], []
fit_slopes, fit_slope_errors = [], []
fit_intercepts, fit_intercept_errors = [], []

for amp_bbox in amp_bboxes:
    amp_im_arr = maskedImage.image[amp_bbox].array
    amp_var_arr = maskedImage.variance[amp_bbox].array
    amp_mask_arr = maskedImage.mask[amp_bbox].array
    
    good = np.isfinite(amp_im_arr) & np.isfinite(amp_var_arr)
    for mask_name in bad_mask_names:
        if mask_name in maskPlaneDict:
            good &= ((amp_mask_arr & (1 << maskPlaneDict[mask_name])) == 0)
    
    fit, fit_cov = np.polyfit(amp_im_arr[good], amp_var_arr[good], deg=1, cov=True)
    gain = 1.0 / fit[0]
    readnoise = np.sqrt(gain * fit[1])
    gn[amp_bbox].array[:, :] = gain
    rn[amp_bbox].array[:, :] = readnoise

    plt.figure(figsize=(4,4))
    plt.scatter(amp_im_arr[good].flatten(), amp_var_arr[good].flatten())
    flux_range = np.array([amp_im_arr[good].min(), amp_im_arr[good].max()])
    var_range = fit[0] * flux_range + fit[1]
    plt.plot(flux_range, var_range, color='red', label='Least squares fit')
    plt.xlabel('Instrumental flux', size=16)
    plt.ylabel('Variance', size=16)
    plt.xticks(size=14)
    plt.yticks(size=14)
    plt.legend(loc='upper left')
    plt.show()

    gains.append(float(gain))
    readnoises.append(float(readnoise))
    fit_slopes.append(float(fit[0]))
    fit_slope_errors.append(float(np.sqrt(fit_cov[0][0])))
    fit_intercepts.append(float(fit[1]))
    fit_intercept_errors.append(float(np.sqrt(fit_cov[1][1])))

print('fit gains:', gains)
print('target gains:', [amp.getGain() if not isElectrons else 1.0 for amp in amps])
print('stored gains:', [amp.getGain() for amp in amps])
print('gains from metadata:', [visit_image.metadata[f"LSST ISR GAIN {amp.getName()}"] for amp in amps])
print('fit read noises (using fit gain):', readnoises)
print('stored readnoises:', [amp.getReadNoise() for amp in amps])
print('readnoises from metadata:', [visit_image.metadata[f"LSST ISR READNOISE {amp.getName()}"] for amp in amps])
print('fit_slopes:', fit_slopes)
print('fit_slope_errors:', fit_slope_errors)
print('fit_intercepts:', fit_intercepts)
print('fit_intercept_errors:', fit_intercept_errors)

plt.figure(figsize=(6,6))
plt.imshow(gn.array, origin='lower', cmap='gray')
plt.title('Effective gains [electrons/adu]')
plt.colorbar()
plt.show()

plt.figure(figsize=(6,6))
plt.imshow(rn.array, origin='lower', cmap='gray')
plt.title('Effective readnoises [electrons/pixel]')
plt.colorbar()
plt.show()

fit gains: [0.4012335058643878, 0.39134485805814817, 0.3953049250369523, 0.4041127281093877, 0.3787299135378056, 0.33498936892425557, 0.3979216444780112, 0.39818580711378043, 0.2802852982871276, 0.4011837604644625, 0.40288499547131484, 0.38712360275036845, 0.39791402278459176, 0.38306058351035105, 0.4029936515054312, 0.40182419795152696]
target gains: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
stored gains: [1.4044, 1.3889, 1.3973, 1.4016, 1.3937, 1.3369, 1.379, 1.3528, 1.3745, 1.3918, 1.3844, 1.3792, 1.4057, 1.4183, 1.4144, 1.4273]
gains from metadata: [1.70924427200861, 1.6998016925389, 1.70168603319942, 1.70044034400547, 1.69489085722979, 1.68396249983445, 1.67826381282722, 1.64244505984655, 1.64017217346943, 1.67084618273791, 1.67510440984028, 1.67580603497722, 1.67822074394158, 1.69733321199622, 1.69503179889932, 1.70818255134929]
fit read noises (using fit gain): [21.375420247192977, 22.440284452712852, 24.003877788480377, 24.941844524491454, 23.12314781698935, 20.819873338125404, 20.696712025705587, 20.538172160326734, 16.826695947268366, 22.467914467034767, 24.29162519247011, 24.881804599304306, 25.61992887517685, 25.08800715166701, 24.042497004167245, 22.496293678827325]
stored readnoises: [8.3, 9.5, 10.0, 10.4, 10.2, 9.3, 8.0, 7.7, 8.1, 8.8, 10.0, 10.8, 11.0, 10.9, 10.1, 9.0]
readnoises from metadata: [11.3885494940233, 12.3220056878356, 13.2911127234589, 13.6725500034007, 12.9861768926762, 12.444939028196, 10.8320087246049, 10.6212321297079, 11.0022978528847, 11.9689924588829, 13.1451619997689, 13.9867322288309, 14.2260228070891, 14.22095622502, 13.0727549547783, 12.0056058375806]
fit_slopes: [2.492314289270718, 2.5552910161181024, 2.529692742650555, 2.4745570491640487, 2.640404056438964, 2.985169359885298, 2.51305756768217, 2.511390366342849, 3.5677932667577448, 2.492623327629887, 2.4820978970193477, 2.5831543023865606, 2.513105703091403, 2.610553116261773, 2.481428668328595, 2.4886505220390744]
fit_slope_errors: [0.0003252776624748071, 0.0001611959110032728, 0.00020225970872419334, 0.0006960689207502442, 0.00013312933531125753, 0.00034763155867758726, 0.00026479413956254967, 0.0005364954981752878, 0.0006053309036035208, 0.000398546395154989, 0.0006686237249611295, 0.00013730004879079113, 0.000313003519083561, 0.00016141373268683905, 0.0006666129156430662, 0.0005007322530851704]
fit_intercepts: [1138.7598096020856, 1286.7586118733252, 1457.57391924839, 1539.4110727329323, 1411.7711483937653, 1293.9727825022303, 1076.4779815807601, 1059.3459338611917, 1010.1767671444829, 1258.2941540643612, 1464.6439086199773, 1599.2416781603013, 1649.554220220201, 1643.1032842748418, 1434.3691520599566, 1259.4642927529624]
fit_intercept_errors: [0.0330403862101013, 0.029226589815756272, 0.03195464397986767, 0.059604765318242645, 0.039169319531257935, 0.3055744007986836, 0.02808445235382666, 0.04683075961816629, 1.1583113503004363, 0.037422774664628865, 0.057043326563551405, 0.03568548415483605, 0.03893618308528594, 0.03689998425808659, 0.058917122615201385, 0.04736760709804708]