Merge pull request #1830 from pypeit/uves_popler

UVES Popler support for coadding data

deprecated/sensfunc.py

@@ -0,0 +1,56 @@
+def compute_blaze(self, wave, trace_spec, trace_spat, flatfile, box_radius=10.0,
+                  min_blaze_value=1e-3, debug=False):
+    """
+    Compute the blaze function from a flat field image.
+
+    Args:
+        wave (`numpy.ndarray`_):
+            Wavelength array. Shape = (nspec, norddet)
+        trace_spec (`numpy.ndarray`_):
+            Spectral pixels for the trace of the spectrum. Shape = (nspec, norddet)
+        trace_spat (`numpy.ndarray`_):
+            Spatial pixels for the trace of the spectrum. Shape = (nspec, norddet)
+        flatfile (:obj:`str`):
+            Filename for the flat field calibration image
+        box_radius (:obj:`float`, optional):
+            Radius of the boxcar extraction region used to extract the blaze function in pixels
+        min_blaze_value (:obj:`float`, optional):
+            Minimum value of the blaze function. Values below this are clipped and set to this value. Default=1e-3
+        debug (:obj:`bool`, optional):
+            Show plots useful for debugging. Default=False
+
+    Returns:
+        `numpy.ndarray`_: The log10 blaze function. Shape = (nspec, norddet)
+        if norddet > 1, else shape = (nspec,)
+    """
+    flatImages = flatfield.FlatImages.from_file(flatfile, chk_version=self.chk_version)
+
+    pixelflat_raw = flatImages.pixelflat_raw
+    pixelflat_norm = flatImages.pixelflat_norm
+    pixelflat_proc, flat_bpm = flat.flatfield(pixelflat_raw, pixelflat_norm)
+
+    flux_box = moment1d(pixelflat_proc * np.logical_not(flat_bpm), trace_spat, 2 * box_radius, row=trace_spec)[0]
+
+    pixtot = moment1d(pixelflat_proc * 0 + 1.0, trace_spat, 2 * box_radius, row=trace_spec)[0]
+    pixmsk = moment1d(flat_bpm, trace_spat, 2 * box_radius, row=trace_spec)[0]
+
+    mask_box = (pixmsk != pixtot) & np.isfinite(wave) & (wave > 0.0)
+
+    # TODO This is ugly and redundant with spec_atleast_2d, but the order of operations compels me to do it this way
+    blaze_function = (np.clip(flux_box * mask_box, 1e-3, 1e9)).reshape(-1, 1) \
+        if flux_box.ndim == 1 else flux_box * mask_box
+    wave_debug = wave.reshape(-1, 1) if wave.ndim == 1 else wave
+    log10_blaze_function = np.zeros_like(blaze_function)
+    norddet = log10_blaze_function.shape[1]
+    for iorddet in range(norddet):
+        blaze_function_smooth = utils.fast_running_median(blaze_function[:, iorddet], 5)
+        blaze_function_norm = blaze_function_smooth / blaze_function_smooth.max()
+        log10_blaze_function[:, iorddet] = np.log10(np.clip(blaze_function_norm, min_blaze_value, None))
+        if debug:
+            plt.plot(wave_debug[:, iorddet], log10_blaze_function[:, iorddet])
+    if debug:
+        plt.show()
+
+    # TODO It would probably better to just return an array of shape (nspec, norddet) even if norddet = 1, i.e.
+    # to get rid of this .squeeze()
+    return log10_blaze_function.squeeze()

doc/coadd1d.rst

@@ -303,3 +303,40 @@ and launches a GUI from the `linetools`_ package. e.g.:
 
     lt_xspec J1217p3905_coadd.fits
 
+UVES_popler coaddition
+======================
+
+If you prefer to use a GUI for the coaddition (to manually remove
+bad pixels, ghosts, cosmic rays etc.), then you can use the UVES_popler tool.
+This tool is developed by Michael Murphy and is available at
+`this link <https://github.com/MTMurphy77/UVES_popler/>`__.
+Here is an example of a coadded spectrum using UVES_popler:
+
+.. image:: ../figures/uves_popler.png
+   :scale: 60%
+
+UVES_popler was originally written to coadd ESO/UVES echelle spectra
+that were reduced by the ESO pipeline, and it has been recently modified
+to support the reduction of PypeIt longslit and echelle data.
+For details on how to use the tool, please refer to the
+`UVES_popler documentation <https://astronomy.swin.edu.au/~mmurphy/UVES_popler/>`__.
+To get you started with reading in PypeIt :doc:`out_spec1D` files,
+you need to generate a text file that lists the absolute paths to the
+:doc:`out_spec1D` files. Here is an example of how to generate this file:
+
+.. code-block:: console
+
+    ls -1 /path/to/your/pypeit_output/Science/spec1d/*.fits > /path/to/your/pypeit_output/pypeit_spec1d_files.txt
+
+Then you can use this file as input to UVES_popler, by using the following command:
+
+.. code-block:: console
+
+    cd /path/to/your/pypeit_output/
+    UVES_popler -disp 50 -filetype 11 pypeit_spec1d_files.txt
+
+This will launch the GUI, where you can interactively coadd your spectra. The
+``-disp 50`` option is used to set the pixel sampling of the spectra to 50 km/s,
+and the ``-filetype 11`` option is used to specify that the input files are PypeIt
+:doc:`out_spec1D` files. For more information on the options available, you can
+specify the ``-h`` option.

doc/releases/1.16.1dev.rst

@@ -92,6 +92,8 @@ Script Changes
 Datamodel Changes
 -----------------
 
+- Adjusted spec1d datamodel to enable use with UVES_popler GUI tool
+
 Under-the-hood Improvements
 ---------------------------
 

pypeit/coadd1d.py

@@ -235,7 +235,7 @@ def load(self):
                 msgs.error("Error in spec1d file for exposure {:d}: "
                            "More than one object was identified with the OBJID={:s} in file={:s}".format(
                     iexp, self.objids[iexp], self.spec1dfiles[iexp]))
-            wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, _, _, _, header = \
+            wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, blaze_iexp, _, header = \
                 sobjs[indx].unpack_object(ret_flam=self.par['flux_value'], extract_type=self.par['ex_value'])
             waves.append(wave_iexp)
             fluxes.append(flux_iexp)
@@ -486,7 +486,7 @@ def load_ech_arrays(self, spec1dfiles, objids, sensfuncfiles):
             indx = sobjs.name_indices(objids[iexp])
             if not np.any(indx):
                 msgs.error("No matching objects for {:s}.  Odds are you input the wrong OBJID".format(objids[iexp]))
-            wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, _, _, meta_spec, header = \
+            wave_iexp, flux_iexp, ivar_iexp, gpm_iexp, blaze_iexp, meta_spec, header = \
                     sobjs[indx].unpack_object(ret_flam=self.par['flux_value'], extract_type=self.par['ex_value'])
             # This np.atleast2d hack deals with the situation where we are wave_iexp is actually Multislit data, i.e. we are treating
             # it like an echelle spectrograph with a single order. This usage case arises when we want to use the

pypeit/core/coadd.py

@@ -688,12 +688,12 @@ def interp_spec(wave_new, waves, fluxes, ivars, gpms, log10_blaze_function=None,
             for ii in range(nexp):
                 fluxes_inter[:,ii], ivars_inter[:,ii], gpms_inter[:,ii], log10_blazes_inter[:,ii] \
                     = interp_oned(wave_new, waves[:,ii], fluxes[:,ii], ivars[:,ii], gpms[:,ii],
-                                  log10_blaze_function = log10_blaze_function[:, ii], sensfunc=sensfunc, kind=kind)
+                                  log10_blaze_function=log10_blaze_function[:, ii], sensfunc=sensfunc, kind=kind)
         else:
             for ii in range(nexp):
                 fluxes_inter[:,ii], ivars_inter[:,ii], gpms_inter[:,ii], _ \
                     = interp_oned(wave_new, waves[:,ii], fluxes[:,ii], ivars[:,ii], gpms[:,ii], sensfunc=sensfunc, kind=kind)
-            log10_blazes_inter=None
+            log10_blazes_inter = None
 
         return fluxes_inter, ivars_inter, gpms_inter, log10_blazes_inter
 

pypeit/core/extract.py

@@ -24,7 +24,8 @@
 
 
 def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
-                    spec, min_frac_use=0.05, fwhmimg=None, base_var=None, count_scale=None, noise_floor=None):
+                    spec, min_frac_use=0.05, fwhmimg=None, flatimg=None,
+                    base_var=None, count_scale=None, noise_floor=None):
 
     r"""
     Perform optimal extraction `(Horne 1986)
@@ -42,6 +43,7 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
       - spec.OPT_COUNTS_NIVAR  -->  Optimally extracted noise variance (sky + read noise) only
       - spec.OPT_MASK  -->   Mask for optimally extracted flux
       - spec.OPT_FWHM  -->   Spectral FWHM (in A) for optimally extracted flux
+      - spec.OPT_FLAT  -->   Flat field spectrum, normalised at the peak value, for the optimally extracted flux
       - spec.OPT_COUNTS_SKY  -->  Optimally extracted sky
       - spec.OPT_COUNTS_SIG_DET  -->  Square root of optimally extracted read noise squared
       - spec.OPT_FRAC_USE  -->  Fraction of pixels in the object profile subimage used for this extraction
@@ -88,6 +90,10 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
     fwhmimg : `numpy.ndarray`_, None, optional:
         Floating-point image containing the modeled spectral FWHM (in pixels) at every pixel location.
         Must have the same shape as ``sciimg``, :math:`(N_{\rm spec}, N_{\rm spat})`.
+    flatimg : `numpy.ndarray`_, None, optional:
+        Floating-point image containing the unnormalized flat-field image. This image
+        is used to extract the blaze function. Must have the same shape as ``sciimg``,
+        :math:`(N_{\rm spec}, N_{\rm spat})`.
     base_var : `numpy.ndarray`_, optional
         Floating-point "base-level" variance image set by the detector properties and
         the image processing steps. See :func:`~pypeit.core.procimg.base_variance`.
@@ -145,6 +151,8 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
     oprof_sub = oprof[:,mincol:maxcol]
     if fwhmimg is not None:
         fwhmimg_sub = fwhmimg[:,mincol:maxcol]
+    if flatimg is not None:
+        flatimg_sub = flatimg[:,mincol:maxcol]
     # enforce normalization and positivity of object profiles
     norm = np.nansum(oprof_sub,axis = 1)
     norm_oprof = np.outer(norm, np.ones(nsub))
@@ -190,6 +198,9 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
     fwhm_opt = None
     if fwhmimg is not None:
         fwhm_opt = np.nansum(mask_sub*ivar_sub*fwhmimg_sub*oprof_sub, axis=1) * utils.inverse(tot_weight)
+    blaze_opt = None
+    if flatimg is not None:
+        blaze_opt = np.nansum(mask_sub*ivar_sub*flatimg_sub*oprof_sub, axis=1) * utils.inverse(tot_weight)
     # Interpolate wavelengths over masked pixels
     badwvs = (mivar_num <= 0) | np.invert(np.isfinite(wave_opt)) | (wave_opt <= 0.0)
     if badwvs.any():
@@ -221,6 +232,9 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
     # Calculate the Angstroms/pixel and Spectral FWHM
     if fwhm_opt is not None:
         fwhm_opt *= np.gradient(wave_opt)  # Convert pixel FWHM to Angstroms
+    # Normalize the blaze function
+    if blaze_opt is not None:
+        blaze_opt /= np.nanmax(blaze_opt)
     # Fill in the optimally extraction tags
     spec.OPT_WAVE = wave_opt    # Optimally extracted wavelengths
     spec.OPT_COUNTS = flux_opt    # Optimally extracted flux
@@ -229,6 +243,8 @@ def extract_optimal(imgminsky, ivar, mask, waveimg, skyimg, thismask, oprof,
     spec.OPT_COUNTS_NIVAR = None if nivar_opt is None else nivar_opt*np.logical_not(badwvs)  # Optimally extracted noise variance (sky + read noise) only
     spec.OPT_MASK = mask_opt*np.logical_not(badwvs)     # Mask for optimally extracted flux
     spec.OPT_FWHM = fwhm_opt  # Spectral FWHM (in Angstroms) for the optimally extracted spectrum
+    if blaze_opt is not None:
+        spec.OPT_FLAT = blaze_opt   # Flat field spectrum, normalised to the peak value
     spec.OPT_COUNTS_SKY = sky_opt      # Optimally extracted sky
     spec.OPT_COUNTS_SIG_DET = base_opt      # Square root of optimally extracted read noise squared
     spec.OPT_FRAC_USE = frac_use    # Fraction of pixels in the object profile subimage used for this extraction
@@ -307,7 +323,7 @@ def extract_asym_boxcar(sciimg, left_trace, righ_trace, gpm=None, ivar=None):
         return flux_out, gpm_box, box_npix, ivar_out
 
 
-def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, base_var=None,
+def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, flatimg=None, base_var=None,
                    count_scale=None, noise_floor=None):
     r"""
     Perform boxcar extraction for a single :class:`~pypeit.specobj.SpecObj`.
@@ -325,6 +341,7 @@ def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, b
       - spec.BOX_COUNTS_NIVAR -->  Box car extracted noise variance
       - spec.BOX_MASK -->  Box car extracted mask
       - spec.BOX_FWHM -->  Box car extracted spectral FWHM
+      - spec.BOX_FLAT -->  Box car extracted flatfield spectrum function (normalized to peak value)
       - spec.BOX_COUNTS_SKY -->  Box car extracted sky
       - spec.BOX_COUNTS_SIG_DET -->  Box car extracted read noise
       - spec.BOX_NPIX  -->  Number of pixels used in boxcar sum
@@ -354,9 +371,12 @@ def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, b
         Container that holds object, trace, and extraction
         information for the object in question. **This object is altered in place!**
         Note that this routine operates one object at a time.
-    fwhmimg : `numpy.ndarray`_, None, optional:
+    fwhmimg : `numpy.ndarray`_, None, optional
         Floating-point image containing the modeled spectral FWHM (in pixels) at every pixel location.
         Must have the same shape as ``sciimg``, :math:`(N_{\rm spec}, N_{\rm spat})`.
+    flatimg : `numpy.ndarray`_, None, optional
+        Floating-point image containing the normalized flat-field. Must have the same shape as
+        ``sciimg``, :math:`(N_{\rm spec}, N_{\rm spat})`.
     base_var : `numpy.ndarray`_, optional
         Floating-point "base-level" variance image set by the detector properties and
         the image processing steps. See :func:`~pypeit.core.procimg.base_variance`.
@@ -407,6 +427,9 @@ def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, b
     fwhm_box = None
     if fwhmimg is not None:
         fwhm_box = moment1d(fwhmimg*mask, spec.TRACE_SPAT, 2*box_radius, row=spec.trace_spec)[0]
+    blaze_box = None
+    if flatimg is not None:
+        blaze_box = moment1d(flatimg*mask, spec.TRACE_SPAT, 2*box_radius, row=spec.trace_spec)[0]
     varimg = 1.0/(ivar + (ivar == 0.0))
     var_box = moment1d(varimg*mask, spec.TRACE_SPAT, 2*box_radius, row=spec.trace_spec)[0]
     nvar_box = None if var_no is None \
@@ -438,7 +461,11 @@ def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, b
     # Calculate the Angstroms/pixel and the final spectral FWHM value
     if fwhm_box is not None:
         ang_per_pix = np.gradient(wave_box)
-        fwhm_box *= ang_per_pix / (pixtot - pixmsk)  # Need to divide by total number of unmasked pixels
+        fwhm_box *= ang_per_pix * utils.inverse(pixtot - pixmsk)  # Need to divide by total number of unmasked pixels
+    # Normalize the blaze function
+    if blaze_box is not None:
+        blaze_box *= utils.inverse(pixtot - pixmsk)  # Need to divide by total number of unmasked pixels
+        blaze_box *= utils.inverse(np.nanmax(blaze_box[mask_box]))  # Now normalize to the peak value
     # Fill em up!
     spec.BOX_WAVE = wave_box
     spec.BOX_COUNTS = flux_box*mask_box
@@ -447,6 +474,8 @@ def extract_boxcar(imgminsky, ivar, mask, waveimg, skyimg, spec, fwhmimg=None, b
     spec.BOX_COUNTS_NIVAR = None if nivar_box is None else nivar_box*mask_box*np.logical_not(bad_box)
     spec.BOX_MASK = mask_box*np.logical_not(bad_box)
     spec.BOX_FWHM = fwhm_box  # Spectral FWHM (in Angstroms) for the boxcar extracted spectrum
+    if blaze_box is not None:
+        spec.BOX_FLAT = blaze_box  # Flat field spectrum, normalised to the peak value
     spec.BOX_COUNTS_SKY = sky_box
     spec.BOX_COUNTS_SIG_DET = base_box
     # TODO - Confirm this should be float, not int