Merge/release4.x

doc/DRAGONS/changes.rst

@@ -8,6 +8,49 @@
 Change Logs
 ***********
 
+4.2.2
+=====
+
+Bug fixes
+---------
+** geminidr.gmos **
+
+* Prevent ``addIllumMaskToDQ`` from crashing when the GMOS longslit
+  bridges cannot be located in the spatial profile, and instead output a
+  warning.
+
+Compatibility
+-------------
+A change to the behavior of Tabs in bokeh v3.9.0 caused the residual/ratio
+plot in the interactive tools to be corrupted. This release ensures
+compatibility with v3.9.0 as well as earlier versions.
+
+
+4.2.1-goawavecal
+================
+
+Improvements
+------------
+
+This internal (used in GOA) release includes the following improvements to
+the determineWavelengthSolution primitive:
+
+* Check the final fit a-posteriori to verify that it is within a factor of two
+  of the initial fit boundary constraints and rejects it if not
+* In Science-Quality mode, if no arc lines are found for all orders being
+  processed, an exception is raised and processing stops. In Quick-Look mode,
+  the previous behavior of adding a linear wavelength solution based on the
+  headers is unchanged
+* If only one arc line is found, a linear fit based on the headers is shifted
+  to match the one arc line, a warning is generated, and that fit is returned as
+  the wavelength solution
+
+This release also includes functionality in gempy to embed auxillary files
+(intended for eg PDF plot figures) in binary tables in astrodata instances, and
+hence FITS files. This is used for example to capture the PDFs of the arc line
+identifications into the reduced arc FITS files, and a corresponding script
+in gempy to extract them.
+
 4.2.1
 =====
 

geminidr/core/parameters_spect.py

@@ -579,6 +579,8 @@ class maskBeyondSlitConfig(config.Config):
         "Minimum fraction of pixel that must be illuminated to not be masked",
         float, 0.9, min=0., max=1., inclusiveMax=True)
 
+class monitorWavelengthSolutionConfig(config.Config):
+    suffix = config.Field("Filename suffix", str, "_wavelengthSolutionMonitoring", optional=True)
 
 class normalizeFlatConfig(config.core_1Dfitting_config):
     suffix = config.Field("Filename suffix", str, "_normalized", optional=True)

geminidr/core/primitives_ccd.py

@@ -225,7 +225,7 @@ def subtractOverscan(self, adinputs=None, **params):
                         # need to match asec location and size
                         pixels = np.arange(asec.y1, asec.y2)
                         data = np.mean(ext.data[asec.y1:asec.y2, x1:x2], axis=axis)
-                        stds = np.std(ext.data[asec.y1:asec.y2, x1:x2], axis=axis)
+                        stds = np.std(ext.data[asec.y1:asec.y2, x1:x2], axis=axis, ddof=1)
                     else:
                         if y1 > asec.y1:  # Bias on top
                             y1 += nbiascontam
@@ -240,17 +240,29 @@ def subtractOverscan(self, adinputs=None, **params):
                         data = np.mean(ext.data[y1:y2, asec.x1:asec.x2], axis=axis)
                         stds = np.std(ext.data[y1:y2, asec.x1:asec.x2], axis=axis)
 
-                    # We have 1 sample (of x2-x1 values) from each of N
-                    # populations each of which has a different population mean,
-                    # (ie the bias leve for that row) but we assume they all
-                    # have the same standard deviation (ie the read noise) and
-                    # we want to estimate that population standard deviation.
-                    # We've calculated the mean (data) and standard deviation
-                    # (stds) of each of the samples, and want to estimate the
-                    # population standard deviation. Because the means are
-                    # different, we can't just calculate the std of all the
-                    # samples.
-                    overstd = np.sqrt(np.mean(stds * stds))
+                    # We have (x2-x1) samples from each of N populations each
+                    # of which has a different population mean, (ie the bias
+                    # level for that row) but we assume they all have the same
+                    # standard deviation (ie the read noise) and we want to
+                    # estimate that population standard deviation.
+                    # We've calculated the mean (data) and estimates of the
+                    # population standard deviation[*] (stds) of each of the set
+                    # of samples, and want to estimate the population standard
+                    # deviation. Because the means are different, we can't just
+                    # calculate the std of all the samples.
+                    # [*] [by which I mean to say we set ddof=1 in the call to
+                    # np.stds to apply the (N/N-1) factor to estimate the
+                    # population standard deviation from the calculaiton of the
+                    # sample standard deviation.
+                    # We use the mean of these per-row population variance
+                    # estimates as our estimate of the global population
+                    # variance. This seems to work well, but I'm not sure if it
+                    # is truly an optimal estimator.
+                    # It turns out that the first and last row of the readout
+                    # often have pixels with massively outlying values, so we
+                    # disregard those samples
+                    trimmedstds = stds[1:-1]
+                    overstd = np.sqrt(np.mean(trimmedstds * trimmedstds))
 
                     # Readnoise descriptor always returns value in electrons
                     if ext.is_in_adu():
@@ -309,7 +321,10 @@ def subtractOverscan(self, adinputs=None, **params):
                     ext.hdr.set('OVERSCAN', previous_overscan + bias_level,
                                 self.keyword_comments['OVERSCAN'])
                     ext.hdr.set('OVERRMS', sigma, self.keyword_comments['OVERRMS'])
-                    ext.hdr.set('OVERRDNS', overstd, self.keyword_comments['OVERRDNS'])
+                    # By convention, "readnoise" is always in electrons
+                    overrdns = overstd * ext.gain() if ext.is_in_adu() else overstd
+                    ext.hdr.set('OVERRDNS', overrdns, self.keyword_comments['OVERRDNS'])
+
                     for desc in ('saturation_level', 'non_linear_level'):
                         with suppress(AttributeError, KeyError):
                             ext.hdr[ad._keyword_for(desc)] -= bias_level

geminidr/core/primitives_spect.py

@@ -16,6 +16,7 @@
 from importlib import import_module
 
 import matplotlib
+import numpy
 import numpy as np
 from astropy import units as u
 from astropy.io.ascii.core import InconsistentTableError
@@ -1218,6 +1219,7 @@ def _get_fit1d_input_data(ext, exptime, spec_table):
                                                    primitive_name="calculateSensitivity",
                                                    filename_info=filename_info,
                                                    help_text=CALCULATE_SENSITIVITY_HELP_TEXT,
+                                                   plot_ratios=False,
                                                    ui_params=uiparams)
                 geminidr.interactive.server.interactive_fitter(visualizer)
 
@@ -2909,22 +2911,44 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                     fit1d.image = image
                     wavecal.update_wcs_with_solution(ext, fit1d, other, config)
             else:
+                bad_solutions = []
+                single_line_solutions = []
                 for ext in ad:
-                    if len(ad) > 1:
-                        log.stdinfo(f"Determining solution for extension {ext.id}")
-
-                    input_data, fit1d, acceptable_fit = wavecal.get_automated_fit(
+                    input_data, fit1d = wavecal.get_automated_fit(
                         ext, uiparams, p=self, linelist=linelist, bad_bits=DQ.not_signal,
                         absorption=absorption)
                     wavecal.update_wcs_with_solution(ext, fit1d, input_data, config)
-                    if not acceptable_fit:
-                        log.warning("No acceptable wavelength solution found")
-                    else:
+                    if fit1d.image.size:
                         figures.append(wavecal.create_pdf_plot(
                             input_data, fit1d.points[~fit1d.mask],
                             fit1d.image[~fit1d.mask],
                             title=f"{ad.filename}:{ext.id}",
                             absorption=absorption))
+                        if fit1d.image.size == 1:
+                            single_line_solutions.append(ext.id)
+                    else:  # no line matches so not an acceptable solution
+                        bad_solutions.append(ext.id)
+
+                if bad_solutions:
+                    msg = f"{ad.filename}: failed to find an acceptable wavelength solution"
+                    if len(ad) > 1:
+                        if len(ad) > len(bad_solutions):
+                            msg += "in extensions " + ", ".join(str(i) for i in bad_solutions)
+                        else:
+                            msg += " in any extensions"
+                    if self.mode == 'sq' and len(ad) == len(bad_solutions):
+                        raise RuntimeError(msg)
+                    log.warning(msg)
+
+                if single_line_solutions:
+                    msg = f"{ad.filename}: a single-line solution has been applied"
+                    if len(ad) > 1:
+                        if len(ad) > len(single_line_solutions):
+                            msg += "in extensions " + ", ".join(str(i) for i in single_line_solutions)
+                        else:
+                            msg += " in all extensions"
+                    log.warning(msg)
+
 
             ad.update_filename(suffix=sfx, strip=True)
             if figures:
@@ -3995,6 +4019,63 @@ def maskBeyondSlit(self, adinputs=None, **params):
 
         return adinputs
 
+    def monitorWavelengthSolution(self, adinputs=None, **params):
+        """
+        This captures some info from the wavelength solution for the GOA
+        instrument monitoring system into various MWS_ headers for the instrument
+        monitoring system to pick up.
+
+        Values Captured include the polynomial coefficients and RMS from the
+        .WAVECAL table, the central wavelength from both the WCS and central_wavlength()
+        descriptor and the difference between them.
+        :return:
+        """
+        log = self.log
+        log.debug(gt.log_message("primitive", self.myself(), "starting"))
+
+        tocapture = ('c0', 'c1', 'c2', 'c3', 'rms', 'inv_rms', 'fwidth')
+
+        for ad in adinputs:
+            for ext in ad:
+                wavecal = getattr(ext, 'WAVECAL', None)
+                if wavecal is None:
+                    continue
+
+                for row in wavecal:
+                    # Capture values from WAVEVAL table into monitoring headers
+                    if row['name'] in tocapture:
+                        keyword = 'MWS_' + row['name'].upper()[:4]
+                        if row['name'] == 'inv_rms':
+                            keyword = 'MWS_IRMS'
+                        ext.hdr[keyword] = (row['coefficients'],
+                                            f'WAVECAL {row['name']} coefficient')
+
+                # Get the model and evaluate it at the mean of its domain
+                model = am.get_named_submodel(ext.wcs.forward_transform, "WAVE")
+                if model:
+                    # Central Wavelength check
+                    center = np.mean(model.domain)
+                    model_central_wlen = model(center)
+                    # We want actual_central_wavelength rather than central_wavelength
+                    header_central_wlen = ext.actual_central_wavelength(asNanometers=True)
+
+                    ext.hdr['MWS_DCWL'] = (model_central_wlen - header_central_wlen,
+                                           'Delta between model and header central_wavelength')
+
+                    # Dispersion check
+                    model_dispersion = model(center+0.5) - model(center-0.5)
+                    header_dispersion = ext.dispersion(asNanometers=True)
+                    ext.hdr['MWS_DDIS'] = (model_dispersion - header_dispersion,
+                                           'Delta between model and header dispersion')
+
+
+                # Count the number of arc lines in the wavecal table
+                ext.hdr['MWS_NUML'] = (numpy.count_nonzero(wavecal['peaks']),
+                                       'Number of non-zero peaks in the WAVECAL table')
+
+            ad.update_filename(suffix=params["suffix"], strip=True)
+
+        return adinputs
 
     def normalizeFlat(self, adinputs=None, **params):
         """

geminidr/core/tests/test_spect.py

@@ -44,11 +44,12 @@
 from specutils.utils.wcs_utils import air_to_vac
 
 import astrodata, gemini_instruments
-from astrodata.testing import ad_compare, assert_most_close
+from astrodata.testing import ad_compare, assert_most_close, download_from_archive
 from gempy.library import astromodels as am
 from gempy.library.config.config import FieldValidationError
 from geminidr.core import primitives_spect
 from geminidr.f2.primitives_f2_longslit import F2Longslit
+from geminidr.gmos.primitives_gmos_longslit import GMOSLongslit
 from geminidr.niri.primitives_niri_image import NIRIImage
 from geminidr.niri.primitives_niri_longslit import NIRILongslit
 from geminidr.gnirs import primitives_gnirs_longslit
@@ -574,6 +575,26 @@ def test_adjust_wavelength_zero_point_controlled(filename, center, shift,
     np.testing.assert_allclose(waves, new_waves, atol=0.1*dw)
 
 
+@pytest.mark.preprocessed_data
+def test_determine_wavelength_solution_exits_with_no_solution_in_sq(path_to_inputs):
+    """Primitive should crash in SQ mode without a solution"""
+    ad = astrodata.open(os.path.join(path_to_inputs, "S20260331S0149_mosaic.fits"))
+    p = GMOSLongslit([ad])
+    p.mode = "sq"  # it's the default, but just to be sure
+    with pytest.raises(RuntimeError):
+        p.determineWavelengthSolution()
+
+
+@pytest.mark.preprocessed_data
+def test_determine_wavelength_solution_continues_with_no_solution_in_qa(path_to_inputs, caplog):
+    """Primitive should crash in SQ mode without a solution"""
+    ad = astrodata.open(os.path.join(path_to_inputs, "S20260331S0149_mosaic.fits"))
+    p = GMOSLongslit([ad])
+    p.mode = "qa"
+    p.determineWavelengthSolution()
+    assert any(record.levelname == 'WARNING' for record in caplog.records)
+
+
 @pytest.mark.preprocessed_data
 @pytest.mark.parametrize('in_file,instrument',
                          [# GNIRS 111/mm LongBlue, off right edge of detector

geminidr/gemini/lookups/keyword_comments.py

@@ -48,7 +48,7 @@
     "OBSTYPE": "Observation type",
     "ORIGNAME": "Original filename prior to processing",
     "OVERRMS": "RMS of fit to overscan mean values",
-    "OVERRDNS": "Readnoise estimate derived from overscan",
+    "OVERRDNS": "Readnoise estimate from overscan [electron]",
     "OVERSCAN": "Median of Overscan mean values",
     "OVERSEC": "Section used for overscan calculation",
     "PIXSCALE": "Pixel scale [arcsec/pixel]",

geminidr/gmos/primitives_gmos_longslit.py

@@ -242,13 +242,20 @@ def addIllumMaskToDQ(self, adinputs=None, suffix=None, illum_mask=None,
                         plt.show()
                         plt.ion()
 
-                    yshift = mshift - maxima[0]
-                    if len(maxima) > 1 or abs(yshift) > mshift:
-                        log.warning(f"{ad.filename}: cross-correlation peak is"
-                                    " untrustworthy so not adding illumination "
-                                    "mask. Please re-run with a specified shift.")
+                    if not maxima:
+                        log.warning(f"{ad.filename}: no peaks found in "
+                                    "cross-correlation, so not adding "
+                                    "illumination mask. Please re-run with a specified shift.")
                         yshift = None
                         log.stdinfo(slit_location_msg)
+                    else:
+                        yshift = mshift - maxima[0]
+                        if len(maxima) > 1 or abs(yshift) > mshift:
+                            log.warning(f"{ad.filename}: cross-correlation peak is"
+                                        " untrustworthy so not adding illumination "
+                                        "mask. Please re-run with a specified shift.")
+                            yshift = None
+                            log.stdinfo(slit_location_msg)
                 else:
                     yshift = shift
 

geminidr/gmos/recipes/qa/recipes_ARC_LS_SPECT.py

@@ -19,5 +19,42 @@ def makeProcessedArc(p):
     p.storeProcessedArc()
     p.writeOutputs()
 
+def checkProcessedArc(p):
+    """
+    Extracts some values from the .WAVECAL extension for instrument monitoring
+    :param p:
+    :return:
+    """
+
+    p.monitorWavelengthSolution()
+    p.writeOutputs()
+
+def checkArc(p):
+    """
+    Reduce an arc, extract some values from the .WAVECAL extension for
+    instrument monitoring
+    :param p:
+    :return:
+    """
+    # In QA mode, if determineWavelengthSolution cannot find a good arc fit it
+    # uses a linear model based on the headers, whereas in SQ mode it raises an
+    # exception. For the purposes of checking arcs, we want everything to
+    # behave as if it is in SQ mode, so we just set that here
+    p.mode = 'sq'
+
+    p.prepare(require_wcs=False)
+    p.addDQ()
+    p.addVAR(read_noise=True)
+    p.overscanCorrect()
+    p.ADUToElectrons()
+    p.addVAR(poisson_noise=True)
+    p.mosaicDetectors()
+    # We force center=None to ensure it always gets the central row.
+    # This is the default anyway, but this prevents it being overridden.
+    p.determineWavelengthSolution(center=None)
+    # We don't need determineDistortion to check Arcs, and it's slow.
+    # p.determineDistortion()
+    p.monitorWavelengthSolution()
+    p.writeOutputs()
 
 _default = makeProcessedArc

geminidr/gnirs/tests/crossdispersed/test_determine_wavelength_solution.py

@@ -328,9 +328,9 @@ def test_regression_determine_wavelength_solution(
             np.testing.assert_allclose(wavelength[indices], ref_wavelength[indices],
                                    atol=tolerance)
             print(f"Test passed for {ad.filename}, extension {wcalibrated_ext.id}")
-        except AssertionError:
+        except AssertionError as e:
             failed = True
-            raise
+            raise AssertionError(f"Test failed for {ad.filename}, extension {wcalibrated_ext.id}\n{repr(model)}\n{repr(ref_model)}")
         finally:
             if write_report:
                 do_report(wcalibrated_ext, ref_ext, failed=failed)

geminidr/gnirs/tests/longslit/test_determine_wavelength_solution.py

@@ -246,6 +246,7 @@ def test_regression_determine_wavelength_solution(
         pixel_scale = ad[0].pixel_scale()  # arcsec / px
         p.viewer = geminidr.dormantViewer(p, None)
 
+        p.mode = 'qa'
         p.determineWavelengthSolution(**{**determine_wavelength_solution_parameters,
                                          **params})