<jupyter_start><jupyter_text>Test preproc file construction to avoid FITS errors  Replicate `desi_preproc` scriptExample call to `desi_preproc`:```desi_preproc -n 20210412 -e 84523 --cameras r0 \    --fibermap $DESI_ROOT/spectro/redux/f3/preproc/20210412/00084523/fibermap-00084523.fits \    --outfile preproc-r0-00084523.fits --model-variance```For details see [desihub/desidatamodel87](https://github.com/desihub/desidatamodel/issues/87).  Running `desi_preproc` in `f3`.* Log file: `/global/cfs/cdirs/desi/spectro/redux/f3/run/scripts/night/20210412/prestdstar-20210412-00084523-a0123456789-51794302.log`* Batch file: `/global/cfs/cdirs/desi/spectro/redux/f3/run/scripts/night/20210412/prestdstar-20210412-00084523-a0123456789.slurm`* srun command: `srun -N 3 -n 200 -c 4 --cpu-bind=cores desi_proc   --traceshift --nostdstarfit --nofluxcalib --cameras a0123456789 -n 20210412 -e 84523 --timingfile /global/cfs/cdirs/desi/spectro/redux/f3/run/scripts/night/20210412/prestdstar-20210412-00084523-a0123456789-timing-$SLURM_JOBID.json --mpi --nofluxcalib`  Set up options<jupyter_code>import os
import subprocess as sub
import time
import warnings
import fitsio
import numpy as np
from astropy.io import fits
from astropy.table import Table
from desiutil.depend import add_dependencies
from desiutil.log import get_logger
from desispec.calibfinder import parse_date_obs, CalibFinder
from desispec.io import findfile, read_raw, write_image, iotime, read_fibermap, read_xytraceset
from desispec.io.util import addkeys, fitsheader
from desispec.util import header2night
from desispec.preproc import preproc, parse_sec_keyword
import desispec.maskbits as maskbits

os.environ['DESI_LOGLEVEL'] = 'DEBUG'
night = 20210412
expid = 84523
camera = 'r0'
infile = findfile('raw', night, expid)
fibermap = os.path.join(os.environ['DESI_SPECTRO_REDUX'], os.environ['SPECPROD'], 'preproc', str(night), f'{expid:08d}', f'fibermap-{expid:08d}.fits')
outfile = os.path.join(os.environ['CSCRATCH'], f'preproc-{camera}-{expid:08d}.fits')
infile, fibermap, outfile<jupyter_output>DEBUG:meta.py:172:findfile: hpixdir = 'hpix'
DEBUG:meta.py:183:findfile: rawdata_dir = '/global/cfs/cdirs/desi/spectro/data'<jupyter_text>Read raw data`desi_preproc` executable corresponds to `desispec.scripts.preproc.main()`, which itself primarily coordinates multiprocessing based on the single-file function `desispec.scripts.preproc.preproc_file()`.  `preproc_file()` itself does very little except call `desispec.io.raw.read_raw()` and `desispec.io.image.write_image()`.<jupyter_code>#
# Interesting design choice: even though fill_header=None, it still fills from HDU 0.  It is triggered by the fact that the keyword exists, not by its value.  Very IDL.
#
img = read_raw(infile, camera, fibermapfile=fibermap, fill_header=None, model_variance=True)<jupyter_output><empty_output><jupyter_text>Write image<jupyter_code>write_image(outfile, img)<jupyter_output><empty_output><jupyter_text>Compare imageDifferences in `OSTEP[ABCD]` can be ignored.  This changed in PR [1543](https://github.com/desihub/desispec/pull/1543) relative to `f3`.<jupyter_code>command = ['fitsdiff', '-c', 'CHECKSUM,DATASUM', outfile, os.path.join(os.path.dirname(fibermap), os.path.basename(outfile))]
p = sub.Popen(command, stdout=sub.PIPE, stderr=sub.PIPE)
out, err = p.communicate()
if out:
    print(out.decode('ascii'))
if err:
    print(out.decode('ascii'))<jupyter_output>fitsdiff: 4.0.1.post1
 a: /global/cscratch1/sd/bweaver/preproc-r0-00084523.fits
 b: /global/cfs/cdirs/desi/spectro/redux/f3/preproc/20210412/00084523/preproc-r0-00084523.fits
 Keyword(s) whose comments are not to be compared:
  CHECKSUM DATASUM
 Maximum number of different data values to be reported: 10
 Relative tolerance: 0.0, Absolute tolerance: 0.0

Primary HDU:

   Headers contain differences:
     Headers have different number of cards:
      a: 548
      b: 545
     Extra keyword 'COMMENT' in a: "  FITS (Flexible Image Transport System) format is defined in 'Astronomy"
     Extra keyword 'EQUINOX' in a: 2000.0
     Extra keyword 'LONGSTRN' in a: 'OGIP 1.0'
     Extra keyword 'EPOCH'  in b: 2000.0
     Keyword CHECKSUM has different values:
        a> KFdZL9ZZKEdZK9ZZ
        b> 9lMgAkJZ9kJfAkJZ
     Keyword DEPVER17 has different values:
        a> 0.48.1
        b> 0.47.1.dev6182
     Keyword OSTEPA   has different values:
        a> 0.725162289410946
        b> 0.62965874851[...]<jupyter_text>Break down `read_raw`<jupyter_code>def test_read_raw(filename, camera, fibermapfile=None, fill_header=None, **kwargs):
    '''Returns preprocessed raw data from `camera` extension of `filename`.
    
    Parameters
    ----------
    filename : :class:`str`
        Input FITS filename with DESI raw data.
    camera : :class:`str`
        Camera name (B0, R1, ... Z9) or FITS extension name.
    fibermapfile : :class:`str`, optional
        Read fibermap from this file; if ``None`` create blank fibermap.
    fill_header : :class:`list`, optional
        A list of HDU names or numbers.  The header cards from these HDUs
        will be added to the header of the camera HDU read from `filename`.

    Returns
    -------
    :class:`desispec.image.Image`
        Image object with member variables pix, ivar, mask, readnoise.
    
    Raises
    ------
    IOError
        If `camera` is not a HDU in `filename`.
    KeyError
        If ``EXPTIME`` is not present in any header in `filename`.
    ValueError
        If ``NIGHT`` in the primary header does not match ``NIGHT`` in the
        camera header, or if `fill_header` is not a :class:`list`.
    
    Notes
    -----
    Other keyword arguments are passed to :func:`desispec.preproc.preproc`,
    *e.g.* bias, pixflat, mask.  See :func:`~desispec.preproc.preproc`
    documentation for details.
    '''

    log = get_logger()

    t0 = time.time()
    fx = fits.open(filename, memmap=False)
    if camera.upper() not in fx:
        raise IOError(f'Camera {camera} not in {filename}!')

    rawimage = fx[camera.upper()].data
    header = fx[camera.upper()].header
    hdu = 0
    #
    # primary_header will typically represent HDU 1 ('SPEC') since
    # HDU 0 is empty.
    #
    while True:
        primary_header = fx[hdu].header
        if "EXPTIME" in primary_header: break

        if len(fx) > hdu + 1:
            if hdu > 0:
                log.warning("Did not find header keyword EXPTIME in HDU %d, moving to the next.", hdu)
            hdu += 1
        else:
            msg = "Did not find header keyword EXPTIME in any HDU of %s!"
            log.critical(msg, filename)
            raise KeyError(msg % filename)

    #- Check if NIGHT keyword is present and valid; fix if needed
    #- e.g. 20210105 have headers with NIGHT='None' instead of YEARMMDD
    try:
        tmp = int(primary_header['NIGHT'])
    except (KeyError, ValueError, TypeError):
        primary_header['NIGHT'] = (header2night(primary_header), 'Observing night')

    try:
        tmp = int(header['NIGHT'])
    except (KeyError, ValueError, TypeError):
        try:
            header['NIGHT'] = (header2night(header), 'Observing night')
        except (KeyError, ValueError, TypeError):
            #- early teststand data only have NIGHT/timestamps in primary hdr
            header['NIGHT'] = (primary_header['NIGHT'], 'Observing night')

    #- early data (e.g. 20200219/51053) had a mix of int vs. str NIGHT
    primary_header['NIGHT'] = (int(primary_header['NIGHT']), 'Observing night')
    header['NIGHT'] = (int(header['NIGHT']), 'Observing night')

    if primary_header['NIGHT'] != header['NIGHT']:
        msg = 'Primary header NIGHT=%d != camera header NIGHT=%d!'
        log.critical(msg, primary_header['NIGHT'], header['NIGHT'])
        raise ValueError(msg % (primary_header['NIGHT'], header['NIGHT']))

    #- early data have >8 char FIBERASSIGN key; rename to match current data
    if 'FIBERASSIGN' in primary_header:
        log.warning('Renaming long header keyword FIBERASSIGN -> FIBASSGN in primary_header.')
        primary_header.rename_keyword('FIBERASSIGN', 'FIBASSGN')

    if 'FIBERASSIGN' in header:
        log.warning('Renaming long header keyword FIBERASSIGN -> FIBASSGN in header.')
        header.rename_keyword('FIBERASSIGN', 'FIBASSGN')
    #
    # A lot of this inheritance stuff is moot because real data files
    # have an empty HDU 0 with no interesting headers.
    #
    inherited = False
    if 'INHERIT' in header and header['INHERIT']:
        inherited = True
        log.info('Camera header %s will INHERIT from HDU 0.', camera) 
        header.extend(fx[0].header, strip=True, unique=True)

    if fill_header is None:
        if inherited:
            hdus = []
        else:
            hdus = [0,]
        if 'PLC' in fx:
            hdus.append('PLC')
    elif isinstance(fill_header, list):
        hdus = fill_header
    else:
        msg = 'Unknown type for fill_header!'
        log.critical(msg)
        raise ValueError(msg)

    log.info('Will add header keywords from HDUs %s.', str(hdus))
    
    for hdu in hdus:
        if hdu in fx or int(hdu) in fx:
            header.extend(fx[hdu].header, strip=True, unique=True)
        else:
            log.warning("HDU %s is not in FITS file.", str(hdu))

    fx.close()
    duration = time.time() - t0
    log.info(iotime.format('read', filename, duration))
    #
    # Other cleanup of headers
    #
    longstrn = fits.Card('LONGSTRN', 'OGIP 1.0', 'The OGIP Long String Convention may be used.')
    if 'MODULE' in primary_header:
        log.debug("Inserting LONGSTRN keyword before MODULE.")
        primary_header.insert('MODULE', longstrn)
    else:
        log.debug("Inserting LONGSTRN keyword before EXTNAME.")
        primary_header.insert('EXTNAME', longstrn)
    log.debug("Renaming EPOCH to EQUINOX in primary_header.")
    primary_header.rename_keyword('EPOCH', 'EQUINOX')
    log.debug("Renaming EPOCH to EQUINOX in header.")
    header.rename_keyword('EPOCH', 'EQUINOX')
    #
    # Run preproc()
    #
    img = preproc(rawimage, header, primary_header, **kwargs)
    #
    # Load fibermap data.
    #
    if fibermapfile is not None and os.path.exists(fibermapfile):
        fibermap = read_fibermap(fibermapfile)
    else:
        log.warning('creating blank fibermap')
        fibermap = empty_fibermap(5000)

    #- Add image header keywords inherited from raw data to fibermap too
    # but I think this was already done in assemble_fibermap
    addkeys(fibermap.meta, img.meta)

    #- Augment the image header with some tile info from fibermap if needed
    for key in ['TILEID', 'TILERA', 'TILEDEC']:
        if key in fibermap.meta:
            if key not in img.meta:
                log.info('Updating header from fibermap {}={}'.format(
                    key, fibermap.meta[key]))
                img.meta[key] = fibermap.meta[key]
            elif img.meta[key] != fibermap.meta[key]:
                #- complain loudly, but don't crash and don't override
                log.error('Inconsistent {}: raw header {} != fibermap header {}'.format(key, img.meta[key], fibermap.meta[key]))


    #- Trim to matching camera based upon PETAL_LOC, but that requires
    #- a mapping prior to 20191211

    #- HACK HACK HACK
    #- TODO: replace this with a mapping from calibfinder, as soon as
    #- that is implemented in calibfinder / desi_spectro_calib
    #- HACK HACK HACK

    #- From DESI-5286v5 page 3 where sp=sm-1 and
    #- "spectro logical number" = petal_loc
    spec_to_petal = {4:2, 2:9, 3:0, 5:3, 1:8, 0:4, 6:6, 7:7, 8:5, 9:1}
    assert set(spec_to_petal.keys()) == set(range(10))
    assert set(spec_to_petal.values()) == set(range(10))

    #- Mapping only for dates < 20191211
    if "NIGHT" in primary_header:
        dateobs = int(primary_header["NIGHT"])
    elif "DATE-OBS" in primary_header:
        dateobs=parse_date_obs(primary_header["DATE-OBS"])
    else:
        msg = "Need either NIGHT or DATE-OBS in primary header"
        log.error(msg)
        raise KeyError(msg)
    if dateobs < 20191211 :
        petal_loc = spec_to_petal[int(camera[1])]
        log.warning('Prior to 20191211, mapping camera {} to PETAL_LOC={}'.format(camera, petal_loc))
    else :
        petal_loc = int(camera[1])
        log.debug('Since 20191211, camera {} is PETAL_LOC={}'.format(camera, petal_loc))

    if 'PETAL_LOC' in fibermap.dtype.names : # not the case in early teststand data
        ii = (fibermap['PETAL_LOC'] == petal_loc)
        fibermap = fibermap[ii]

    cfinder = None

    camname = camera.upper()[0]
    if camname == 'B':
        badamp_bit = maskbits.fibermask.BADAMPB
    elif camname == 'R':
        badamp_bit = maskbits.fibermask.BADAMPR
    else:
        badamp_bit = maskbits.fibermask.BADAMPZ


    if 'FIBER' in fibermap.dtype.names : # not the case in early teststand data

        ## Mask fibers
        cfinder = CalibFinder([header,primary_header])
        fibers  = fibermap['FIBER'].data
        for k in ["BROKENFIBERS","BADCOLUMNFIBERS","LOWTRANSMISSIONFIBERS"] :
            log.debug("{}={}".format(k,cfinder.badfibers([k])))

        ## Mask bad fibers
        fibermap['FIBERSTATUS'][np.in1d(fibers,cfinder.badfibers(["BROKENFIBERS"]))] |= maskbits.fibermask.BROKENFIBER
        fibermap['FIBERSTATUS'][np.in1d(fibers,cfinder.badfibers(["BADCOLUMNFIBERS"]))] |= maskbits.fibermask.BADCOLUMN
        fibermap['FIBERSTATUS'][np.in1d(fibers,cfinder.badfibers(["LOWTRANSMISSIONFIBERS"]))] |= maskbits.fibermask.LOWTRANSMISSION
        # Also, for backward compatibility
        fibermap['FIBERSTATUS'][np.in1d(fibers%500,cfinder.badfibers(["BROKENFIBERS"])%500)] |= maskbits.fibermask.BROKENFIBER
        fibermap['FIBERSTATUS'][np.in1d(fibers%500,cfinder.badfibers(["BADCOLUMNFIBERS"])%500)] |= maskbits.fibermask.BADCOLUMN
        fibermap['FIBERSTATUS'][np.in1d(fibers%500,cfinder.badfibers(["LOWTRANSMISSIONFIBERS"])%500)] |= maskbits.fibermask.LOWTRANSMISSION

        # Mask Fibers that are set to be excluded due to CCD/amp/readout issues
        fibermap['FIBERSTATUS'][np.in1d(fibers,cfinder.badfibers(["BADAMPFIBERS"]))] |= badamp_bit
        fibermap['FIBERSTATUS'][np.in1d(fibers,cfinder.badfibers(["EXCLUDEFIBERS"]))] |= badamp_bit # for backward compatibiliyu
        fibermap['FIBERSTATUS'][np.in1d(fibers%500,cfinder.badfibers(["BADAMPFIBERS"])%500)] |= badamp_bit
        fibermap['FIBERSTATUS'][np.in1d(fibers%500,cfinder.badfibers(["EXCLUDEFIBERS"])%500)] |= badamp_bit # for backward compatibiliyu
        if cfinder.haskey("EXCLUDEFIBERS") :
            log.warning("please use BADAMPFIBERS instead of EXCLUDEFIBERS")

    if np.sum(img.mask & maskbits.ccdmask.BADREADNOISE > 0) >= img.mask.size//4 :
        log.info("Propagate ccdmask.BADREADNOISE to fibermap FIBERSTATUS")

        if cfinder is None :
            cfinder = CalibFinder([header,primary_header])

        psf_filename = cfinder.findfile("PSF")
        tset = read_xytraceset(psf_filename)
        mean_wave =(tset.wavemin+tset.wavemax)/2.
        xfiber  = tset.x_vs_wave(np.arange(tset.nspec),mean_wave)
        amp_ids = desispec.preproc.get_amp_ids(header)

        for amp in amp_ids :
            kk  = parse_sec_keyword(header['CCDSEC'+amp])
            ntot = img.mask[kk].size
            nbad = np.sum((img.mask[kk] & maskbits.ccdmask.BADREADNOISE) > 0)
            if nbad / ntot > 0.5 :
                # not just nbad>0 b/c/ there are always pixels with low QE
                # that have increased readnoise after pixel flatfield
                log.info("Setting BADREADNOISE bit for fibers of amp {}".format(amp))
                badfibers = (xfiber>=kk[1].start-3)&(xfiber<kk[1].stop+3)
                fibermap["FIBERSTATUS"][badfibers] |= ( maskbits.fibermask.BADREADNOISE | badamp_bit )

    img.fibermap = fibermap

    return img

img = test_read_raw(infile, camera, fibermapfile=fibermap, fill_header=None, model_variance=True)<jupyter_output>INFO:<ipython-input-16-e6c4137c090e>:122:test_read_raw: Will add header keywords from HDUs [0].
INFO:<ipython-input-16-e6c4137c090e>:132:test_read_raw: iotime 0.496 sec to read desi-00084523.fits.fz at 2022-01-05T15:41:09.174763
DEBUG:<ipython-input-16-e6c4137c090e>:138:test_read_raw: Inserting LONGSTRN keyword before MODULE.
DEBUG:<ipython-input-16-e6c4137c090e>:143:test_read_raw: Renaming EPOCH to EQUINOX in primary_header.
DEBUG:<ipython-input-16-e6c4137c090e>:145:test_read_raw: Renaming EPOCH to EQUINOX in header.
DEBUG:calibfinder.py:182:__init__: header['CAMERA']=r0
DEBUG:calibfinder.py:186:__init__: header['SPECID']=4
DEBUG:calibfinder.py:236:__init__: Use spectrograph hardware identifier SMY
DEBUG:calibfinder.py:246:__init__: reading calib data in /global/cfs/cdirs/desi/spectro/desi_spectro_calib/trunk/spec/sm4/sm4-r.yaml
DEBUG:calibfinder.py:258:__init__: Found 6 data for camera sm4-r in filename /global/cfs/cdirs/desi/spectro/desi_spectro_calib/trunk/spec/sm4/sm4-r.yaml
DEBUG[...]<jupyter_text>Break down `write_image`<jupyter_code>def test_write_image(outfile, image, meta=None):
    """Writes image object to outfile

    Args:
        outfile : output file string
        image : desispec.image.Image object
            (or any object with 2D array attributes image, ivar, mask)

    Optional:
        meta : dict-like object with metadata key/values (e.g. FITS header)
    """

    log = get_logger()
    if meta is not None:
        hdr = fitsheader(meta)
    else:
        hdr = fitsheader(image.meta)

    add_dependencies(hdr)

    #- Work around fitsio>1.0 writing blank keywords, e.g. on 20191212
    for key in hdr.keys():
        if type(hdr[key]) == fits.card.Undefined:
            log.warning('Setting blank keyword {} to None'.format(key))
            hdr[key] = None

    outdir = os.path.dirname(os.path.abspath(outfile))
    if not os.path.isdir(outdir):
        os.makedirs(outdir, exist_ok=True)

    hx = fits.HDUList()
    hdu = fits.ImageHDU(image.pix.astype(np.float32), name='IMAGE', header=hdr)
    if 'CAMERA' not in hdu.header:
        hdu.header.append( ('CAMERA', image.camera.lower(), 'Spectrograph Camera') )

    if 'RDNOISE' not in hdu.header and np.isscalar(image.readnoise):
        hdu.header.append( ('RDNOISE', image.readnoise, 'Read noise [RMS electrons/pixel]'))

    hx.append(hdu)
    hx.append(fits.ImageHDU(image.ivar.astype(np.float32), name='IVAR'))
    hx.append(fits.CompImageHDU(image.mask.astype(np.int16), name='MASK'))
    if not np.isscalar(image.readnoise):
        hx.append(fits.ImageHDU(image.readnoise.astype(np.float32), name='READNOISE'))

    if hasattr(image, 'fibermap'):
        if isinstance(image.fibermap, Table):
            with warnings.catch_warnings():
                #- nanomaggies aren't an official IAU unit but don't complain
                warnings.filterwarnings('ignore', ".*nanomaggies.*")
                fmhdu = fits.convenience.table_to_hdu(image.fibermap)
            fmhdu.name = 'FIBERMAP'
        else:
            fmhdu = fits.BinTableHDU(image.fibermap, name='FIBERMAP')

        hx.append(fmhdu)

    t0 = time.time()
    hx.writeto(outfile+'.tmp', overwrite=True, checksum=True)
    os.rename(outfile+'.tmp', outfile)
    duration = time.time() - t0
    log.info(iotime.format('write', outfile, duration))

    return outfile

test_write_image(outfile, img)<jupyter_output>INFO:<ipython-input-24-6e9d0c1c1e03>:61:test_write_image: iotime 63.597 sec to write preproc-r0-00084523.fits at 2022-01-05T15:45:46.862792<jupyter_text>Test fitsio on fibermap files<jupyter_code># Huh, fitsio changes the order of comment cards.
fmap, fmap_hdr = fitsio.read(fibermap, ext='FIBERMAP', header=True)
fmap_hdr
fits.Header.fromstring(str(fmap_hdr).lstrip(), sep='\n')
fmap_table = Table(fmap)
fmap_table
fmap_table.meta
fibermap
img = fits.PrimaryHDU(np.random.uniform(size=(10,10)))
img.header['COMMENT'] = 'This is the first comment.'
img.header['COMMENT'] = 'This is the second comment.'
img.header['HISTORY'] = 'This is the first history.'
img.header['HISTORY'] = 'This is the second history.'
img.header['EXTNAME'] = 'TEST'
img.header['PLANCK'] = (6.62607015e-34, "[J s] Planck's Constant")
hdulist = fits.HDUList([img])
hdulist.writeto(os.path.join(os.environ['CSCRATCH'], 'fitsio_header_test.fits'), overwrite=True)
img.header
t, hdr = fitsio.read(os.path.join(os.environ['CSCRATCH'], 'fitsio_header_test.fits'), ext='TEST', header=True)
t
hdr
repr(hdr)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12415548810339763, 'fraction_numerical': 0.034906491726231274, 'mean_word_length': 3.42238579779173, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 49, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '643823', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7633, 'n_tokens_neox': 7093, 'n_words': 1617, 'path': 'doc/nb/PreprocConstructionTest.ipynb', 'repo_name': 'echaussidon/desispec', 'size': '20426'}