Ringing near patch boundaries

[dmargala]

The 2D extraction is sensitive to nsubbundles and nwavestep at the ~0.01 dflux/sigma level near the patch boundaries. For example, comparing nsubbundles=6 (default in desi_extract_spectra) with nsubbundles=5:

I checked that running with nsubbundles=6 reproduces the cascades result exactly. See below for details.

Code to reproduce using latest desi environment on cori

import os

import numpy as np
import matplotlib.pyplot as plt

import specter
import desispec.io
import desispec.scripts.extract

def extract_and_compare(ax=None, nsubbundles=6, nwavestep=50):
    specprod_dir = "/global/cfs/cdirs/desi/spectro/redux/cascades"

    night = 20201221
    expid = 69276
    camera = "z0"
    # first 300 wavebins of zcam wavelength range
    wavestr = "7520,7759.2,0.8"
    nspec = 25

    frame   = desispec.io.findfile('frame'  , night=night, expid=expid, camera=camera, specprod_dir=specprod_dir)
    preproc = desispec.io.findfile('preproc', night=night, expid=expid, camera=camera, specprod_dir=specprod_dir)
    psf     = desispec.io.findfile('psf'    , night=night, expid=expid, camera=camera, specprod_dir=specprod_dir)

    output = os.path.join(os.environ.get('SCRATCH'), 'tmp-frame.fits')

    cmdargs = (
        f" --input {preproc}"
        f" --psf {psf}"
        f" --output {output}"
        f" --wavelength {wavestr}"
        f" --nspec {nspec}"
        f" --nsubbundles {nsubbundles}"
        f" --nwavestep {nwavestep}"
        f" --barycentric-correction"
    )

    args = desispec.scripts.extract.parse(cmdargs.split())
    desispec.scripts.extract.main(args)

    cascades = desispec.io.read_frame(frame, nspec=nspec)
    compare = desispec.io.read_frame(output, nspec=nspec)

    waveslice = slice(0, len(compare.wave))

    combined_var = 1.0/cascades.ivar[:, waveslice] + 1.0/compare.ivar
    pull_flux = (cascades.flux[:, waveslice] - compare.flux) / np.sqrt(combined_var)

    if ax is None:
        ax = plt.gca()
    ax.plot(compare.wave, np.max(np.abs(pull_flux), axis=0), label='max')
    ax.plot(compare.wave, np.median(np.abs(pull_flux), axis=0), label='median')
    ax.legend()
    ax.set_xlabel("Wavelength [A]")
    ax.set_ylabel(r"$|\Delta f|/ \sigma$")

    title = f"{frame}"
    title += "\nvs\n"
    title += f"specter (nspec={nspec}, nsubbundles={nsubbundles}, nwavestep={nwavestep})"
    ax.set_title(title)

fig, ax = plt.subplots(figsize=(8, 4))
extract_and_compare(ax=ax, nsubbundles=6)

fig, ax = plt.subplots(figsize=(8, 4))
extract_and_compare(ax=ax, nsubbundles=5)

fig, ax = plt.subplots(figsize=(8, 4))
extract_and_compare(ax=ax, nsubbundles=6, nwavestep=100)

[dmargala]

Here is a look at the full wavelength range and all 500 spec for the same frame. I also added comparisons with blanc and gpu_specter (using preproc and psf inputs from cascades).