Merge branch 'develop' into binospec_overscan

Author: tepickering

pypeit/core/coadd.py

@@ -2199,7 +2199,7 @@ def combspec(waves, fluxes, ivars, gpms, sn_smooth_npix,
 
     return wave_grid_mid, wave_stack, flux_stack, ivar_stack, gpm_stack
 
-def multi_combspec(waves, fluxes, ivars, masks, sn_smooth_npix=None,
+def multi_combspec(waves, fluxes, ivars, gpms, sn_smooth_npix=None,
                    wave_method='linear', dwave=None, dv=None, dloglam=None, spec_samp_fact=1.0, wave_grid_min=None,
                    wave_grid_max=None, ref_percentile=70.0, maxiter_scale=5,
                    sigrej_scale=3.0, scale_method='auto', hand_scale=None, sn_min_polyscale=2.0, sn_min_medscale=0.5,
@@ -2213,17 +2213,17 @@ def multi_combspec(waves, fluxes, ivars, masks, sn_smooth_npix=None,
     Parameters
     ----------
     waves : list
-        List of `numpy.ndarray`_ wavelength arrays with shape (nspec_i,) with
+        List of `numpy.ndarray`_ arrays with shape (nspec_i,) with
         the wavelength arrays of the spectra to be coadded.
     fluxes : list
-        List of `numpy.ndarray`_ wavelength arrays with shape (nspec_i,) with
+        List of `numpy.ndarray`_ arrays with shape (nspec_i,) with
         the flux arrays of the spectra to be coadded.
     ivars : list
-        List of `numpy.ndarray`_ wavelength arrays with shape (nspec_i,) with
+        List of `numpy.ndarray`_ arrays with shape (nspec_i,) with
         the ivar arrays of the spectra to be coadded.
-    masks : list
-        List of `numpy.ndarray`_ wavelength arrays with shape (nspec_i,) with
-        the mask arrays of the spectra to be coadded.
+    gpms : list
+        List of `numpy.ndarray`_ arrays with shape (nspec_i,) with
+        the GPM arrays of the spectra to be coadded (True=good).
     sn_smooth_npix : int, optional
         Number of pixels to median filter by when computing S/N used to decide
         how to scale and weight spectra. If set to None, the code will determine
@@ -2319,7 +2319,7 @@ def multi_combspec(waves, fluxes, ivars, masks, sn_smooth_npix=None,
         core.wavecal.wvutils.py for more.
     wave_stack : `numpy.ndarray`_, (ngrid,)
         Wavelength grid for stacked spectrum. As discussed above, this is the
-        weighted average of the wavelengths of each spectrum that contriuted to
+        weighted average of the wavelengths of each spectrum that contributed to
         a bin in the input wave_grid wavelength grid. It thus has ngrid
         elements, whereas wave_grid has ngrid+1 elements to specify the ngrid
         total number of bins. Note that wave_stack is NOT simply the wave_grid
@@ -2341,7 +2341,7 @@ def multi_combspec(waves, fluxes, ivars, masks, sn_smooth_npix=None,
         log.info('Using a sn_smooth_npix={:d} to decide how to scale and weight your spectra'.format(sn_smooth_npix))
 
     wave_grid_mid, wave_stack, flux_stack, ivar_stack, mask_stack = combspec(
-        waves, fluxes,ivars, masks, wave_method=wave_method, dwave=dwave, dv=dv, dloglam=dloglam,
+        waves, fluxes,ivars, gpms, wave_method=wave_method, dwave=dwave, dv=dv, dloglam=dloglam,
         spec_samp_fact=spec_samp_fact, wave_grid_min=wave_grid_min, wave_grid_max=wave_grid_max, ref_percentile=ref_percentile,
         maxiter_scale=maxiter_scale, sigrej_scale=sigrej_scale, scale_method=scale_method, hand_scale=hand_scale,
         sn_min_medscale=sn_min_medscale, sn_min_polyscale=sn_min_polyscale, sn_smooth_npix=sn_smooth_npix,

pypeit/core/wavecal/autoid.py

@@ -1100,6 +1100,8 @@ def full_template(spec, lamps, par, ok_mask, det, binspectral, nsnippet=2, slit_
         if par['reid_arxiv'] is None:
             raise PypeItError('WavelengthSolutionPar parameter `reid_arxiv` not '
                        'specified for "full_template" method.')
+        wvstr = "the entire wavelength range" if par['wvrng_arxiv'] is None else f"wavelength range {par['wvrng_arxiv']} Angstroms"
+        log.info(f"Loading template from {par['reid_arxiv']} for detector {det} and {wvstr}")
         temp_wv_og, temp_spec_og, temp_bin, order, lines_pix, lines_wav, lines_fit_ord = \
             waveio.load_template(par['reid_arxiv'], det, wvrng=par['wvrng_arxiv'])
     else:

pypeit/core/wavecal/templates.py

@@ -56,7 +56,7 @@
 
 
 def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
-                   normalize=False, subtract_conti=False, wvspec=None,
+                   scalespec=False, scalevals=None, subtract_conti=False, wvspec=None,
                    lowredux=False, ifiles=None, det_cut=None, chk=False,
                    miny=None, overwrite=True, ascii_tbl=False, in_vac=True,
                    shift_wave=False, binning=None, micron=False,
@@ -93,9 +93,16 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
             Show a plot or two
         miny (float):
             Impose a minimum value
-        normalize (bool, optional):
-            If provided multiple in_files, normalize each
-            snippet to have the same maximum amplitude.
+        scalespec (bool, optional):
+            If True, scale the spectrum by the provided scalevals. If
+            scalevals is None, the spectrum will be scaled to the
+            normalization value (default is 10000). Default is False, which will not scale.
+        scalevals (list, optional):
+            If scalespec is True, the values to scale each spectrum by.
+            Default is None, which will not scale. If provided, the
+            length of scalevals must match the number of in_files. If
+            scalespec is True and scalevals is None, the spectrum will
+            be scaled to the normalization value (default is 10000).
         subtract_conti (bool, optional):
             Subtract the continuum for the final archive
         ascii_tbl (bool, optional):
@@ -117,6 +124,12 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
         ifiles = np.arange(len(in_files))
     if binning is None:
         binning = [None]*len(ifiles)
+    if scalespec and scalevals is None:
+        scalevals = [None]*len(ifiles)
+    # Prepare an axis if we are debugging
+    if chk:
+        plt.clf()
+        ax = plt.gca()
     # Load xidl file
     # Grab it
     # Load and splice
@@ -128,6 +141,7 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
     # Loop on the files
     for kk, slit in enumerate(slits):
         # Load up
+        pypeitFit = None  # Default value, in case a PypeItFit is not available
         if wvspec is None:
             in_file = in_files[ifiles[kk]]
             if lowredux:
@@ -146,9 +160,24 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
         log.info("wvmin, wvmax of {}: {}, {}".format(in_file, wv_vac.min(), wv_vac.max()))
         # Cut
         if len(slits) > 1:
+            # Plot the inputs if we are debugging
+            if chk:
+                ax.plot(wv_vac, spec, drawstyle='steps-mid')
+            # Rebin spec if binning is different
+            if binning is not None and binning[kk] != binspec:
+                npix_orig = spec.size
+                x_orig = np.arange(npix_orig) / float(npix_orig - 1)
+                npix = int(npix_orig * binning[kk] / binspec)
+                x = np.arange(npix) / float(npix - 1)
+                # Interpolate the wavelengths and spectrum onto the new grid
+                spec = (interp1d(x_orig, spec, axis=0,
+                                 bounds_error=False, fill_value='extrapolate'))(x)
+                wv_vac = (interp1d(x_orig, wv_vac, axis=0,
+                                 bounds_error=False, fill_value='extrapolate'))(x)
+            # Default good pixels
             wvmin, wvmax = grab_wvlim(kk, wv_cuts, len(slits))
-            # Default
             gdi = (wv_vac > wvmin) & (wv_vac < wvmax)
+            npix = wv_vac.size
             if shift_wave:
                 if len(lvals) > 0:
                     # Find pixel closet to end
@@ -161,16 +190,10 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
                     # Delta pix -- approximate but should be pretty good
                     dpix = dwv_specs / dwv_snipp[ipix]
                     # Calculate new wavelengths
-                    npix = wv_vac.size
-                    # Rebin spec?
-                    if binning is not None and binning[kk] != binspec:
-                        npix_orig = spec.size
-                        x_orig = np.arange(npix_orig) / float(npix_orig - 1)
-                        x = np.arange(npix) / float(npix - 1)
-                        spec = (interp1d(x_orig, spec, axis=0,
-                                         bounds_error=False, fill_value='extrapolate'))(x)
-                    # Evaluate
-                    new_wave = pypeitFit.eval((-dpix + np.arange(npix)) / (npix - 1))
+                    if pypeitFit is None:
+                        new_wave = -dwv_specs + wv_vac
+                    else:
+                        new_wave = pypeitFit.eval((-dpix + np.arange(npix)) / (npix - 1))
                     # Range
                     iend = np.argmin(np.abs(new_wave - wvmax))
                     # Interpolate
@@ -189,15 +212,16 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
         for kk,spec in enumerate(yvals):
             _, _, _, _, spec_cont_sub = wvutils.arc_lines_from_spec(spec)
             yvals[kk] = spec_cont_sub
-    # Normalize?
-    if normalize:
-        norm_val = 10000.
-        # Max values
-        maxs = []
-        for kk,spec in enumerate(yvals):
-            mx = np.max(spec)
-            spec = spec * norm_val / mx
-            yvals[kk] = spec
+    # Scale the input spectra by the user-provided scalevals or by the maximum value of each snippet
+    if scalespec:
+        for kk, spec in enumerate(yvals):
+            if scalevals[kk] is not None:
+                scale_factor = scalevals[kk]
+            else:
+                norm_val = 10000.
+                # Use the normalization value if scalevals is not provided
+                scale_factor = norm_val / np.max(spec)
+            yvals[kk] = spec * scale_factor
     # Concatenate
     nwspec = np.concatenate(yvals)
     nwwv = np.concatenate(lvals)
@@ -206,9 +230,7 @@ def build_template(in_files, slits, wv_cuts, binspec, outroot, outdir=None,
         nwspec = np.maximum(nwspec, miny)
     # Check
     if chk:
-        plt.clf()
-        ax = plt.gca()
-        ax.plot(nwwv, nwspec)
+        ax.plot(nwwv, nwspec, 'k-', drawstyle='steps-mid')
         plt.show()
     # Generate the table
     wvutils.write_template(nwwv, nwspec, binspec, outdir, outroot, det_cut=det_cut, overwrite=overwrite)
@@ -603,7 +625,7 @@ def main(flg):
         lcut = [7840.]
         wfile = template_path / 'Keck_LRIS' / 'R600_7500' / 'MasterWaveCalib_I_1_01.json'
         build_template(wfile, slits, lcut, binspec, outroot, lowredux=False,
-                       chk=True, normalize=True, subtract_conti=True)
+                       chk=True, scalespec=True, subtract_conti=True)
 
     # ##################################
     # Magellan/MagE
@@ -877,6 +899,21 @@ def main(flg):
         outfile = dataPaths.reid_arxiv.path / iout
         tbl.write(outfile, overwrite=True)
         print("Wrote: {}".format(outfile))
+    # Keck KCWI
+    if flg & (2**35):
+        pass
+    # Keck KCRM
+    if flg & (2**36):
+        dirc = template_path / 'KCWI' / 'RM1'
+        infiles = [dirc / 'keck_kcrm_rm1_lcen6230.fits',dirc / 'keck_kcrm_rm1_lcen7010.fits']
+        outroot = 'keck_kcrm_RM1.fits'
+        slits = [0, 0]
+        binspec = 1 # Desired binning
+        binning = [1, 2] # Spectral binning of each infile
+        scalevals = [1.0, 3.4] # Scale the second one down by 2x to match the first
+        lcut = [6910.0]
+        build_template(infiles, slits, lcut, binspec, outroot, scalespec=True, scalevals=scalevals,
+                       binning=binning, reid_files=True, shift_wave=True, chk=True)
 
 
 # Command line execution
@@ -939,5 +976,11 @@ def main(flg):
     # P200 Triplespec
     #flg += 2**34
 
+    # Keck KCWI
+    #flg += 2**35
+
+    # Keck KCRM
+    flg += 2**36
+
     main(flg)
 

pypeit/data/arc_lines/reid_arxiv/README

@@ -9,6 +9,11 @@ DEIMOS
 
 Generated with the scripts in pypeit.core.wavecal.templates.py
 
+Keck/KCWI + Keck/KCRM
+---------------------
+
+Generated with the scripts in pypeit.core.wavecal.templates.py
+
 Kastb
 -----
 

pypeit/data/arc_lines/reid_arxiv/keck_kcrm_RM1.fits

@@ -1,8 +1,8 @@
 # This directory contains data files for the blackbody stars
-# presented in Suzuki & Fukugita, 2018, AJ, 156, 219:
+# presented in Suzuki & Fukugita (2018), AJ, 156, 219:
 # https://ui.adsabs.harvard.edu/abs/2018AJ....156..219S/abstract
-# with updates from Cooke, Suzuki, Prochaska 2025:
-# <insert abs link here>
+# with updates from Cooke, Suzuki, Prochaska (2026), MNRAS, 546, 2246:
+# https://ui.adsabs.harvard.edu/abs/2026MNRAS.546f2246C/abstract
 # The flux distribution over the range 3000 – 25000A can be
 # generated using the "a" and "T" parameters in Table below,
 # combined with Equation 1 from the above paper.

pypeit/data/standards/blackbody/blackbody_info.txt

@@ -667,7 +667,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False,
                 arcfitter.store_solution(final_fit[str(slit_idx)], "", self.binspectral,
                                          # JTS: Commented out as it seems the
                                          # store solution function does not
-                                         # take these arguments any more.
+                                         # take these arguments anymore.
                                          # specname=self.spectrograph.name,
                                          # gratname="UNKNOWN", dispangl="UNKNOWN"
                                          )
@@ -739,7 +739,7 @@ def build_wv_calib(self, arccen, method, skip_QA=False,
             raise PypeItError('Unrecognized wavelength calibration method: {:}'.format(method))
 
         # Build the DataContainer
-        if self.par['redo_slits'] is not None:
+        if self.par['redo_slits'] is not None and prev_wvcalib is not None:
             # If we are only redoing slits, we start from the
             #  previous wv_calib and update only the (good) redone slits
             self.wv_calib = prev_wvcalib