Merge pull request #2217 from desihub/psf-with-bad-amp

Support psf and fiberflat calibrations of CCDs with a bad amplifier

py/desispec/ccdcalib.py

@@ -415,7 +415,6 @@ def _find_zeros(night, cameras, nzeros=25, nskip=2):
     select_zeros=exptable['OBSTYPE']=='zero'
     bad = select_zeros & (exptable['LASTSTEP']!='all')
     badcam = select_zeros & (exptable['BADCAMWORD']!='')
-    badamp = select_zeros & (exptable['BADAMPS']!='')
     notallcams = select_zeros & (exptable['CAMWORD']!='a0123456789')
 
     expdicts = dict()
@@ -436,9 +435,9 @@ def _find_zeros(night, cameras, nzeros=25, nskip=2):
             log.info(f'Dropping {ndrop} bad {zerotype} ZEROs: {drop_expids}')
             expids = expids[~drop]
 
-        if np.any(badcam|badamp|notallcams):
+        if np.any(badcam|notallcams):
             #do the by spectrograph evaluation of bad spectra
-            drop = np.isin(expids, exptable['EXPID'][badcam|badamp|notallcams])
+            drop = np.isin(expids, exptable['EXPID'][badcam|notallcams])
             ndrop = np.sum(drop)
             drop_expids = expids[drop]
             expids = expids[~drop]
@@ -447,20 +446,19 @@ def _find_zeros(night, cameras, nzeros=25, nskip=2):
             if len(expids) >= nzeros:
                 #in this case we can just drop all partially bad exposures as we have enough that are good on all cams
                 log.info(f'Additionally dropped {ndrop} partially bad {zerotype} ZEROs for'
-                         + f' all cams because of BADCAM/BADAMP/CAMWORD: {drop_expids}')
+                         + f' all cams because of BADCAM/CAMWORD: {drop_expids}')
             else:
                 #in this case we want to recover as many as possible
                 log.info(f'additionally dropped {ndrop} bad {zerotype} ZEROs for some cams'
-                         + f'because of BADCAM/BADAMP/CAMWORD: {drop_expids}')
+                         + f'because of BADCAM/CAMWORD: {drop_expids}')
 
                 for expid in drop_expids:
                     select_exp=exptable['EXPID']==expid
-                    badampstring=exptable['BADAMPS'][select_exp][0]
                     goodcamword=difference_camwords(exptable['CAMWORD'][select_exp][0],
                                                     exptable['BADCAMWORD'][select_exp][0])
                     goodcamlist=decode_camword(goodcamword)
                     for camera in goodcamlist:
-                        if camera in cameras and camera not in badampstring:
+                        if camera in cameras:
                             expdict[camera].append(expid)
         else:
             expdict={f'{cam}': expids for cam in cameras}

py/desispec/fiberflat.py

@@ -12,7 +12,7 @@
 from desispec.linalg import cholesky_solve
 from desispec.linalg import cholesky_solve_and_invert
 from desispec.linalg import spline_fit
-from desispec.maskbits import specmask
+from desispec.maskbits import specmask, fibermask
 from desispec.maskedmedian import masked_median
 from desispec import util
 import scipy,scipy.sparse
@@ -323,7 +323,7 @@ def compute_fiberflat(frame, nsig_clipping=10., accuracy=5.e-4, minval=0.1, maxv
     # reset ivar
     ivar = frame.ivar.copy()
 
-    fiberflat_mask=12 # place holder for actual mask bit when defined
+    fiberflat_mask=specmask.BADFIBERFLAT
 
     nsig_for_mask=nsig_clipping # only mask out N sigma outliers
 
@@ -347,7 +347,7 @@ def compute_fiberflat(frame, nsig_clipping=10., accuracy=5.e-4, minval=0.1, maxv
                 smooth_fiberflat=spline_fit(wave,wave[w],fiberflat[fiber,w],smoothing_res,fiberflat_ivar[fiber,w])
             except (ValueError, TypeError) as e :
                 print("error in spline_fit")
-                mask[fiber] += fiberflat_mask
+                mask[fiber] |= fiberflat_mask
                 fiberflat_ivar[fiber] = 0.
                 break
 
@@ -360,7 +360,7 @@ def compute_fiberflat(frame, nsig_clipping=10., accuracy=5.e-4, minval=0.1, maxv
                     ii=np.argsort(chi2[bad])
                     bad=bad[ii[-nbadmax:]]
 
-                mask[fiber,bad] += fiberflat_mask
+                mask[fiber,bad] |= fiberflat_mask
                 fiberflat_ivar[fiber,bad] = 0.
                 nbad_tot += bad.size
             else :
@@ -445,6 +445,10 @@ def compute_fiberflat(frame, nsig_clipping=10., accuracy=5.e-4, minval=0.1, maxv
     log.info("add a systematic error of 0.0035 to fiberflat variance (calibrated on sims)")
     fiberflat_ivar = (fiberflat_ivar>0)/( 1./ (fiberflat_ivar+(fiberflat_ivar==0) ) + 0.0035**2)
 
+    # update mask based upon final ivar=0
+    bad = fiberflat_ivar == 0
+    mask[bad] |= fiberflat_mask
+
     fiberflat = FiberFlat(wave, fiberflat, fiberflat_ivar, mask, mean_spectrum,
                      chi2pdf=chi2pdf,header=frame.meta,fibermap=frame.fibermap)
 
@@ -655,11 +659,22 @@ def autocalib_fiberflat(fiberflats):
         ivar      = (var>0)/(var+(var==0))
         fiberflat /= ii.size
         meanspec /= ii.size
+
+        # set tiny ivars back to 0; leftover from not propagaging inf, 1/inf, nan...
+        ivar[ivar<1e-8] = 0.0
+        mask[ivar == 0.0] |= specmask.BADFIBERFLAT
+
+        # set a FIBERSTATUS mask bit for fibers with more than 25% of their pixels masked
+        fibermap = fflat0.fibermap.copy()
+        fraction_masked = np.sum(mask != 0, axis=1) / mask.shape[1]
+        bad = fraction_masked > 0.25
+        fibermap['FIBERSTATUS'][bad] |= fibermask.BADFLAT
+
         fflat = FiberFlat(fflat0.wave,fiberflat,ivar,mask,
                           meanspec=meanspec,
                           header=fflat0.header,
                           fibers=fflat0.fibers,
-                          fibermap=fflat0.fibermap,
+                          fibermap=fibermap,
                           spectrograph=fflat0.spectrograph)
         output_fiberflats[spec] = fflat
         mflat.append(np.median(fiberflat,axis=0))

py/desispec/io/fiberflat.py

@@ -21,6 +21,7 @@
 from .meta import findfile
 from .util import fitsheader, native_endian, makepath, checkgzip
 from .util import get_tempfilename
+from .fibermap import read_fibermap
 from . import iotime
 
 def write_fiberflat(outfile,fiberflat,header=None, fibermap=None):
@@ -112,7 +113,7 @@ def read_fiberflat(filename):
         meanspec  = native_endian(fx["MEANSPEC"].data.astype('f8'))
         wave      = native_endian(fx["WAVELENGTH"].data.astype('f8'))
         if 'FIBERMAP' in fx:
-            fibermap = fx['FIBERMAP'].data
+            fibermap = read_fibermap(fx)
         else:
             fibermap = None
 

py/desispec/io/util.py

@@ -6,6 +6,7 @@
 """
 import os
 import glob
+import re
 import time
 import datetime
 import subprocess
@@ -660,7 +661,7 @@ def camword_intersection(camwords, full_spectros_only=False):
     return camword
 
 
-def erow_to_goodcamword(erow, suppress_logging=False):
+def erow_to_goodcamword(erow, suppress_logging=False, exclude_badamps=False):
     """
     Takes a dictionary-like object with at minimum keys for
     OBSTYPE (type str), CAMWORD (type str), BADCAMWORD (type str),
@@ -670,7 +671,7 @@ def erow_to_goodcamword(erow, suppress_logging=False):
         erow (dict-like): the exposure table row with columns OBSTYPE, CAMWORD,
                           BADCAMWORD, and BADAMPS.
         suppress_logging (bool), whether to suppress logging messages. Default False.
-
+        exclude_badamps (bool), if True, discard cameras with a badamp set
 
     Returns:
         goodcamword (str): Camword for that observation given the obstype and
@@ -680,11 +681,12 @@ def erow_to_goodcamword(erow, suppress_logging=False):
                                   badcamword=erow['BADCAMWORD'],
                                   badamps=erow['BADAMPS'],
                                   obstype=erow['OBSTYPE'],
-                                  suppress_logging=suppress_logging)
+                                  suppress_logging=suppress_logging,
+                                  exclude_badamps=exclude_badamps)
 
 
 def columns_to_goodcamword(camword, badcamword, badamps=None, obstype=None,
-                           suppress_logging=False):
+                           suppress_logging=False, exclude_badamps=False):
     """
     Takes a camera information and obstype and returns the correct camword
     of good cameras for that type of observation.
@@ -695,13 +697,14 @@ def columns_to_goodcamword(camword, badcamword, badamps=None, obstype=None,
         badamps (str): comma seperated list of bad amplifiers [brz][0-9][A-D]
         obstype (str), obstype of exposure
         suppress_logging (bool), whether to suppress logging messages. Default False.
+        exclude_badamps (bool), if True, discard cameras with a badamp set
 
     Returns:
         goodcamword (str): Camword for that observation given the obstype and
                            input camera information.
     """
     log = get_logger()
-    if not suppress_logging:
+    if exclude_badamps and not suppress_logging:
         if badamps is None:
             log.info("No badamps given, proceeding without badamp removal")
         elif obstype is None:
@@ -711,14 +714,13 @@ def columns_to_goodcamword(camword, badcamword, badamps=None, obstype=None,
 
     goodcamword = difference_camwords(camword, badcamword,
                                       suppress_logging=suppress_logging)
-    if badcamword is not None and badamps != '':
-        if obstype is None or obstype != 'science':
-            badampcams = []
-            for (camera, petal, amplifier) in parse_badamps(badamps):
-                badampcams.append(f'{camera}{petal}')
-            badampcamword = create_camword(np.unique(badampcams))
-            goodcamword = difference_camwords(goodcamword, badampcamword,
-                                              suppress_logging=suppress_logging)
+    if exclude_badamps and badamps != '':
+        badampcams = []
+        for (camera, petal, amplifier) in parse_badamps(badamps):
+            badampcams.append(f'{camera}{petal}')
+        badampcamword = create_camword(np.unique(badampcams))
+        goodcamword = difference_camwords(goodcamword, badampcamword,
+                                          suppress_logging=suppress_logging)
 
     return goodcamword
 
@@ -853,6 +855,45 @@ def validate_badamps(badamps,joinsymb=','):
         log.info(f'Badamps given as: {badamps} verified to work with modifications to: {newbadamps}')
     return newbadamps
 
+def get_amps_for_camera(camera_amps, camera=None):
+    """
+    Return just the amplifier IDs for a given camera
+
+    Parameters
+    ----------
+    camera_amps : str
+        Either comma-separted camera+amps like "b7A,z3B,b2C"
+        or just the amplifier letters like "A" or "CD"
+    camera : str, optional
+        Camera like b7 or z3.  Must be specified if camera_amps
+        is full comma-separated camera+amps list.
+
+    Returns
+    -------
+    amps : str
+        String of amplifiers like 'A' or 'BC'.
+        Unique, sorted, uppercase string of amps like 'A' or 'BC'
+    """
+    log = get_logger()
+    #- if no number in camera_amps, interpret just as ABCD not r7A,b2C,z1D
+    if re.match(r'.*\d.*', camera_amps) is None:
+        result = ''.join(sorted(camera_amps.strip().replace(',','')))
+    elif camera is None:
+        msg = f'Must specify camera to filter {camera_amps=}'
+        log.error(msg)
+        raise ValueError(msg)
+    else:
+        camera = camera.lower()
+        amps = list()
+        for (cam, pet, amp) in parse_badamps(camera_amps):
+            thiscam = f'{cam}{pet}'.lower()
+            if camera == thiscam:
+                amps.append(amp)
+
+        result = ''.join(sorted(set(amps)))
+
+    return result.upper()
+
 def all_impacted_cameras(badcamword, badamps):
     """
     Checks badcamword string and badamps string (badamps joined by ',') and

py/desispec/maskbits.py

@@ -49,6 +49,7 @@
     - [PIXFLATLOW,      6, "pixflat < 0.1"]
     - [HIGHVAR,         7, "High variability in pixel value"]
     - [BADREADNOISE,    8, "Very high CCD amplifier read noise"]
+    - [BADAMP,          9, "Bad CCD amplifier"]
 
 #- Mask bits that apply to an entire fiber
 #- Note: some of these are informational and aren't necessarily bad
@@ -103,7 +104,7 @@
     specmask = BitMask('specmask', _bitdefs)
     ccdmask = BitMask('ccdmask', _bitdefs)
     fibermask = BitMask('fibermask', _bitdefs)
-    extractmaskval = ccdmask.mask("BAD|HOT|DEAD|SATURATED|COSMIC|PIXFLATZERO|PIXFLATLOW|HIGHVAR")
+    extractmaskval = ccdmask.mask("BAD|HOT|DEAD|SATURATED|COSMIC|PIXFLATZERO|PIXFLATLOW|HIGHVAR|BADAMP")
 
 except TypeError:
     #

py/desispec/preproc.py

@@ -704,7 +704,8 @@ def preproc(rawimage, header, primary_header, bias=True, dark=True, pixflat=True
             bias_img=None,model_variance=False,no_traceshift=False,bkgsub_science=False,
             keep_overscan_cols=False,no_overscan_per_row=False,no_ccd_region_mask=False,
             no_cte_corr=False,cte_params_filename=None,
-            fallback_on_dark_not_found=False):
+            fallback_on_dark_not_found=False,
+            badamps=''):
     '''
     preprocess image using metadata in header
 
@@ -756,6 +757,8 @@ def preproc(rawimage, header, primary_header, bias=True, dark=True, pixflat=True
 
     Optional disabling of CTE correction if no_cte_corr=True
 
+    Optional `badamps` masks those amplifiers
+
     Returns Image object with member variables:
         pix : 2D preprocessed image in units of electrons per pixel
         ivar : 2D inverse variance of image
@@ -818,6 +821,15 @@ def preproc(rawimage, header, primary_header, bias=True, dark=True, pixflat=True
     #- Check if this file uses amp names 1,2,3,4 (old) or A,B,C,D (new)
     amp_ids = get_amp_ids(header)
 
+    #- Check badamps input and add keyword
+    if badamps != '':
+        if not np.all(np.isin(list(badamps), amp_ids)):
+            msg = f'{badamps=} not in list of {amp_ids=}'
+            log.critical(msg)
+            raise ValueError(msg)
+
+        header['BADAMPS'] = badamps
+
     #- if CAMERA is missing, this will raise an exception in a few lines,
     #- but allows CAMERA logging in the meantime if it is present
     try:
@@ -1121,6 +1133,12 @@ def preproc(rawimage, header, primary_header, bias=True, dark=True, pixflat=True
             row_subtracted_overscan_col = raw_overscan_col - overscan_col[:,None]
             o,r = calc_overscan(row_subtracted_overscan_col)
             rdnoise  = np.repeat(r,nrows)
+
+        #- Mask bad amplifiers
+        if amp.upper() in badamps.upper():
+            log.warning(f'Masking {camera} amplifier {amp} as BADAMP')
+            mask[kk] |= ccdmask.BADAMP
+
         if bias is not False :
             # the master bias noise is already in the raw data
             # (because we already subtracted the bias)

py/desispec/scripts/preproc.py

@@ -12,7 +12,7 @@
 import multiprocessing as mp
 import numpy as np
 from desispec import io
-from desispec.io.util import get_tempfilename
+from desispec.io.util import get_tempfilename, get_amps_for_camera
 from desispec.parallel import default_nproc
 from desiutil.log import get_logger
 log = get_logger()
@@ -43,6 +43,10 @@ def parse(options=None):
                         help = 'output preprocessed image file')
     parser.add_argument('--cameras', type = str, default = None, required=False,
                         help = 'comma separated list of cameras')
+    parser.add_argument('--badamps', type=str, default='', required=False,
+                        help = 'comma separated list of bad CCD amps. ' +
+                               'If processing multiple cameras give full b8C,r7A.'
+                               'If processing a single camera can be just amp names (ABCD).')
     parser.add_argument('--bias', type = str, default = None, required=False,
                         help = 'bias image calibration file')
     parser.add_argument('--dark', type = str, default = None, required=False,
@@ -205,7 +209,8 @@ def main(args=None):
                 no_overscan_per_row=args.no_overscan_per_row,
                 no_cte_corr=args.no_cte_correction,
                 cte_params_filename=args.cte_params,
-                fallback_on_dark_not_found=args.fallback_on_dark_not_found
+                fallback_on_dark_not_found=args.fallback_on_dark_not_found,
+                badamps=get_amps_for_camera(args.badamps, camera)
         )
         opts_array.append(opts)
 

py/desispec/scripts/proc.py

@@ -430,6 +430,8 @@ def main(args=None, comm=None):
             if not args.obstype in ['ARC'] : # never model variance for arcs
                 if not args.no_model_pixel_variance and args.obstype != 'DARK' :
                     cmd += " --model-variance"
+            if args.badamps is not None:
+                cmd += f" --badamps {args.badamps}"
 
             inputs = [args.input]
 

py/desispec/scripts/proc_joint_fit.py

@@ -291,8 +291,8 @@ def main(args=None, comm=None):
             if len(inflats) < 3 * 4 * len(args.cameras):
                 log.critical("Fewer than 3 exposures with 4 lamps were available. Can't perform joint fit. Exiting...")
                 enough_inputs = False
-            elif len(args.cameras) < 6:
-                log.critical("Fewer than 6 cameras were available, so couldn't perform joint fit. Exiting ...")
+            elif len(args.cameras) < 4:
+                log.critical("Fewer than 4 cameras were available, so couldn't perform joint fit. Exiting ...")
                 enough_inputs = False
 
         if comm is not None: