Merge branch 'enh/zero_point_adjustment'

geminidr/core/parameters_spect.py

@@ -28,6 +28,17 @@ def validate_regions_int(value, multiple=False):
     ranges = at.parse_user_regions(value, dtype=int, allow_step=True)
     return multiple or len(ranges) == 1
 
+class adjustWavelengthZeroPointConfig(config.Config):
+    suffix = config.Field("Filename suffix", str, "_wavelengthZeroPointAdjusted",
+                          optional=True)
+    center = config.RangeField("Central row/column to extract", int, None,
+                                min=1, optional=True)
+    shift = config.RangeField("Shift to apply in pixels (None: determine automatically)",
+                              float, 0, min=-2048, max=2048, optional=True)
+    verbose = config.Field("Print extra information", bool, False,
+                           optional=True)
+    debug_max_shift = config.RangeField("Maximum shift to allow (in pixels)",
+                                        float, 5, min=0)
 
 class adjustWCSToReferenceConfig(config.Config):
     suffix = config.Field("Filename suffix",

geminidr/core/primitives_preprocess.py

@@ -836,8 +836,8 @@ def flatCorrect(self, adinputs=None, suffix=None, flat=None, do_cal=True):
                 rectified = False
 
             if rectified:
-                log.stdinfo(f"{ad.filename}: applying slit rectification model "
-                            f"from the flat {flat.filename}")
+                log.stdinfo(f"{ad.filename}: adding slit rectification model "
+                            f"derived from {flat.filename} to WCS")
                 for ext in ad:
                     ext.wcs.insert_frame(ext.wcs.input_frame, rect_model,
                                          cf.Frame2D(name='rectified'))

geminidr/core/primitives_spect.py

@@ -82,6 +82,197 @@ def _initialize(self, adinputs, **kwargs):
         super()._initialize(adinputs, **kwargs)
         self._param_update(parameters_spect)
 
+    def adjustWavelengthZeroPoint(self, adinputs=None, **params):
+        """
+        Find sky lines and match them to a linelist in order to shift the
+        wavelength scale zero point slightly to account for flexure in the
+        telescope.
+
+        Parameters
+        ----------
+        adinputs : list of :class:`~astrodata.AstroData`
+            Wavelength calibrated 2D spectra.
+        suffix : str
+            Suffix to be added to output files
+        center : None or int
+            Central row/column for 1D extraction (None => use middle).
+        shift : float, Default : None
+            An optional shift to apply directly to the wavelength scale, in
+            pixels. If not given, the shift will be calculated from the sky
+            lines present in the image.
+        verbose : bool, Default : False
+            Print additional information on the fitting process.
+        """
+
+        def _add_shift_model_to_wcs(shift, dispaxis, ext):
+            """ Create a model to shift the wavelength scale by and add to WCS
+
+            This function creates a compound model of two Shift models in
+            parallel, with the one in the `dispaxis` direction getting the value
+            of `shift`. This allows the resulting compound model to applied to
+            a WCS without any special handling for different spectral axis
+            orientations.
+
+            Parameters
+            ----------
+            shift : float
+                The shift to apply, in pixels (will be converted to wavelength)
+            dispaxis : int, 0 or 1
+                The dispersion axis, in the Python sense
+            ext : Astrodata extension
+                The extension to apply the shift to.
+            """
+            dx, dy = 0, 0
+            if dispaxis == 0:
+                dy = shift
+            elif dispaxis == 1:
+                dx = shift
+            else:
+                raise ValueError("'dispaxis' must be 0 (vertical) or "
+                                 "1 (horizontal)")
+
+            # This should work for both orientations without having to code
+            # them separately.
+            model = models.Shift(dx) & models.Shift(dy)
+            model.name = 'FLEXCORR'
+            ext.wcs.insert_frame(ext.wcs.input_frame, model,
+                                 cf.Frame2D(name="wavelength_scale_adjusted"))
+
+        # Set up log
+        log = self.log
+        log.debug(gt.log_message("primitive", self.myself(), "starting"))
+        timestamp_key = self.timestamp_keys[self.myself()]
+        sfx = params["suffix"]
+        center = params["center"]
+        shift = params["shift"]
+        max_shift = params["debug_max_shift"]
+        verbose = params["verbose"]
+
+        # Check given shift, if there is one.
+        if shift and abs(shift) > max_shift:
+            raise ValueError("Provided shift is larger than parameter "
+                             f"'debug_max_shift': |{shift:.3g}| > "
+                             f"{max_shift:.3g}")
+
+        for ad in adinputs:
+            log.stdinfo(f"{ad.filename}:")
+
+            for ext in ad:
+                dispaxis = 2 - ext.dispersion_axis()  # Python sense
+                row_or_col = 'row' if dispaxis == 1 else 'column'
+
+                # If the user specifies a shift value, apply it and continue
+                # (case of shift == 0 caught and handled above)
+                if shift is not None:
+                    if shift == 0:
+                        msg = "    No wavelength shift from sky lines will be"\
+                              " performed since shift=0"
+                    else:
+                        msg = "    Shifted wavelength scale for extension "\
+                              f"{ext.id} by {shift:0.4g} pixels"
+                    _add_shift_model_to_wcs(shift, dispaxis, ext)
+                    log.stdinfo(msg)
+                    continue
+
+                # Otherwise, we'll need to automatically find the shift.
+                # Values (generally) taken from the defaults for
+                # determineWavelengthSolution, with some changes.
+                config_dict = {
+                        "order": 3,
+                        "niter": 3,
+                        "lsigma": 3.0,
+                        "hsigma": 3.0,
+                        "central_wavelength": None,
+                        "interactive": False,
+                        "center": center,
+                        "nsum": 10,
+                        "fwidth": None,
+                        "min_snr": 10,
+                        "min_sep": 2,
+                        "weighting": "local",
+                        "nbright": 0,
+                        "dispersion": None,
+                        "debug_min_lines": 15,
+                        "debug_alternative_centers": False,
+                    }
+
+                wave_scale = ext.wcs.output_frame.axes_names[0]
+                if wave_scale == 'WAVE':
+                    config_dict["in_vacuo"] = True
+                    log.fullinfo("Calibrated in vacuum.")
+                elif wave_scale == 'AWAV':
+                    log.fullinfo("Calibrated in air.")
+                    config_dict["in_vacuo"] = False
+                else:
+                    raise ValueError("Cannot interpret wavelength scale "
+                                     f"for {ext.filename}:{ext.id} "
+                                     f"(found '{wave_scale}'")
+
+                # get_all_input_data() outputs several lines of information
+                # which can be useful but confusing if many files are processed,
+                # so use the verbose parameter to allow users to control it.
+                # loglevel is the level at which the output should be logged,
+                # so higher levels (e.g. stdinfo) print more to the console.
+                loglevel = "stdinfo" if verbose else "fullinfo"
+                try:
+                    input_data = wavecal.get_all_input_data(
+                        ext, self, config_dict, linelist=None,
+                        bad_bits=DQ.not_signal, skylines=True,
+                        loglevel=loglevel)
+                except ValueError:
+                    raise ValueError("Something went wrong in finding sky "
+                                     "lines - check that the spectrum is being "
+                                     f"taken in a {row_or_col} free of the "
+                                     "object aperture, and change it with the "
+                                     "`center` parameter if necessary")
+
+                spectrum = input_data["spectrum"]
+                init_models = input_data["init_models"]
+                domain = init_models[0].meta["domain"]
+                peaks, weights = input_data["peaks"], input_data["weights"]
+                sky_lines = input_data["linelist"].wavelengths(
+                    in_vacuo=config_dict["in_vacuo"], units='nm')
+                sky_weights = input_data["linelist"].weights
+                if sky_weights is None:
+                    sky_weights = np.ones_like(sky_lines)
+                    log.debug("No weights were found for the reference linelist")
+
+                m_init = init_models[0]
+                # Fix all parameters in the model so that they don't change
+                # (only the Shift which will be added next).
+                for p in m_init.param_names:
+                    getattr(m_init, p).fixed = True
+
+                # Add a bounded Shift model in front of the wavelength solution
+                m_init = models.Shift(0, bounds={'offset': (-max_shift,
+                                                            max_shift)}) | m_init
+                m_init.meta["domain"] = domain
+
+                fwidth = input_data["fwidth"]
+                dw = np.diff(m_init(np.arange(spectrum.size))).mean()
+                kdsigma = fwidth * abs(dw)
+                k = 1 if kdsigma < 3 else 2
+
+                fit_it = KDTreeFitter(sigma=2 * abs(dw), maxsig=5,
+                                      k=k, method='Nelder-Mead')
+                m_final = fit_it(m_init, peaks, sky_lines,
+                                 in_weights=weights[config_dict["weighting"]],
+                                 ref_weights=sky_weights, matches=None,
+                                 options={'disp': True if verbose else False})
+
+                # Apply the shift to the wavelength scale
+                shift_final = m_final.offset_0.value
+                log.stdinfo(f"    Shifted wavelength scale for "
+                            f"extension {ext.id} by {shift_final:0.4g} "
+                            f"pixels ({shift_final * dw:0.4g} nm)")
+                _add_shift_model_to_wcs(shift_final, dispaxis, ext)
+
+            # Timestamp and update the filename
+            gt.mark_history(ad, primname=self.myself(), keyword=timestamp_key)
+            ad.update_filename(suffix=sfx, strip=True)
+
+        return adinputs
+
     def adjustWCSToReference(self, adinputs=None, **params):
         """
         Compute offsets along the slit by cross-correlation, or use offset
@@ -1386,7 +1577,7 @@ def determineSlitEdges(self, adinputs=None, **params):
 
                 if not edge_pairs:
                     log.warning("No edges could be determined for "
-                                f"{ad.filename} - no SLITEDGE table wil be "
+                                f"{ad.filename} - no SLITEDGE table will be "
                                 "attached. Consider setting positions of "
                                 "edges manually using the `edges1` and "
                                 "`edges2` parameters.")
@@ -2019,7 +2210,8 @@ def determineWavelengthSolution(self, adinputs=None, **params):
                         log.stdinfo(f"Determining solution for extension {ext.id}")
 
                     input_data, fit1d, acceptable_fit = wavecal.get_automated_fit(
-                        ext, uiparams, p=self, linelist=linelist, bad_bits=DQ.not_signal)
+                        ext, uiparams, p=self, linelist=linelist,
+                        bad_bits=DQ.not_signal)
                     if not acceptable_fit:
                         log.warning("No acceptable wavelength solution found "
                                     f"for {ext.id}")
@@ -2545,7 +2737,7 @@ def flagCosmicRays(self, adinputs=None, **params):
         log.fullinfo(f"  spatial_order: {y_order_in}")
         log.fullinfo(f"  bkgmodel: {bkgmodel}")
         log.fullinfo(f"  sigclip: {params['sigclip']}")
-        log.fullinfo(f"  sigfrac: {params['sigfrac']}\m")
+        log.fullinfo(f"  sigfrac: {params['sigfrac']}\n")
 
         for ad in adinputs:
             is_in_adu = ad[0].is_in_adu()

geminidr/core/tests/test_spect.py

@@ -338,6 +338,97 @@ def test_sky_correct_from_slit_with_multiple_sources():
 
     np.testing.assert_allclose(ad_out[0].data, source, atol=1e-3)
 
+@pytest.mark.preprocessed_data
+@pytest.mark.parametrize('in_shift', [0, -1.2, 2.75])
+def test_adjust_wavelength_zero_point_shift(in_shift, change_working_dir,
+                                            path_to_inputs):
+    with change_working_dir(path_to_inputs):
+        ad = astrodata.open('N20220706S0337_wavelengthSolutionAttached.fits')
+
+    p = GNIRSLongslit([ad])
+    ad_out = p.adjustWavelengthZeroPoint(shift=in_shift).pop()
+    transform = ad_out[0].wcs.get_transform('pixels',
+                                            'wavelength_scale_adjusted')
+    shift = getattr(transform, 'offset_1')
+    assert shift == pytest.approx(in_shift)
+
+@pytest.mark.preprocessed_data
+@pytest.mark.parametrize('in_shift', [-16, 7.7])
+def test_adjust_wavelength_zero_point_overlarge_shift(in_shift,
+                                                      change_working_dir,
+                                                      path_to_inputs):
+    with change_working_dir(path_to_inputs):
+        ad = astrodata.open('N20220706S0337_wavelengthSolutionAttached.fits')
+
+    p = GNIRSLongslit([ad])
+    with pytest.raises(ValueError):
+        p.adjustWavelengthZeroPoint(shift=in_shift).pop()
+
+@pytest.mark.preprocessed_data
+@pytest.mark.regression
+@pytest.mark.parametrize('filename,instrument',
+                         [('N20220706S0337', 'GNIRS'),
+                          ('N20110331S0400', 'GNIRS'),
+                          ('N20150511S0123', 'GNIRS'),
+                          ('N20220718S0140', 'GNIRS'),
+                          ('N20130827S0128', 'GNIRS'),
+                          ('S20140610S0077', 'F2'),
+                          ('S20210430S0138', 'F2'),
+                          ('S20150629S0230', 'F2'),
+                          ('S20210709S0035', 'F2'),
+                          ('S20170215S0111', 'F2'),
+                          ('S20180125S0028', 'F2'),
+                          ('N20050918S0135', 'NIRI'),
+                          ('N20050627S0040', 'NIRI'),
+                          ('N20070615S0118', 'NIRI'),
+                          ('N20061114S0193', 'NIRI'),
+                          ])
+def test_adjust_wavelength_zero_point_auto_shift(filename, instrument,
+                                                 change_working_dir,
+                                                 path_to_inputs):
+    # Dictionary of shift values (in pixels) for each file.
+    results = {'N20220706S0337': -0.0119375, # GNIRS 111/mm 0.10" LongBlue
+               'N20110331S0400': 0.1669375,  # GNIRS 111/mm 0.30" ShortBlue
+               'N20150511S0123': -0.0273750, # GNIRS 32/mm  0.45" ShortBlue
+               'N20220718S0140': 2.083875,   # GNIRS 32/mm  0.10" LongBlue
+               'N20130827S0128': -3.2903125, # GNIRS 10/mm  0.10" LongBlue
+               'S20140610S0077': -1.755625,  # F2    R3K 1pix-slit f/16
+               'S20210430S0138': 0.2556250,  # F2    R3K 2pix-slit f/16
+               'S20150629S0230': 0.3927500,  # F2    JH  3pix-slit f/16
+               'S20210709S0035': 0.3030625,  # F2    JH  4pix-slit f/16
+               'S20170215S0111': 0.0551250,  # F2    HK  6pix-slit f/16
+               'S20180125S0028': -0.046375,  # F2    JH  8pix-slit f/16
+               'N20050918S0135': 0.6130625,  # NIRI  Hgrism f6-6pix
+               'N20050627S0040': -0.059625,  # NIRI  Hgrism f6-6pix
+               'N20070615S0118': -0.029875,  # NIRI  Jgrism f6-6pix
+               'N20061114S0193': 0.1915000,  # NIRI  Kgrism f6-2pix
+               }
+
+    classes_dict = {'GNIRS': GNIRSLongslit,
+                    'F2': F2Longslit,
+                    'NIRI': NIRILongslit}
+
+    # In some files the aperture lies (at least partly) in the cental column/row
+    centers = {'N20130827S0128': 800,
+               'N20220718S0140': 300,
+               'S20140610S0077': 600}
+    try:
+        center = centers[filename]
+    except KeyError:
+        center = None
+
+    with change_working_dir(path_to_inputs):
+        ad = astrodata.open(filename + '_wavelengthSolutionAttached.fits')
+
+    p = classes_dict[instrument]([ad])
+    ad_out = p.adjustWavelengthZeroPoint(center=center, shift=None).pop()
+    transform = ad_out[0].wcs.get_transform('pixels',
+                                             'wavelength_scale_adjusted')
+    param = 'offset_0' if instrument == 'NIRI' else 'offset_1'
+    shift = getattr(transform, param)
+
+    assert shift == pytest.approx(results[filename])
+
 
 @pytest.mark.preprocessed_data
 @pytest.mark.parametrize('filename,instrument',

geminidr/f2/primitives_f2_spect.py

@@ -141,7 +141,7 @@ def _get_arc_linelist(self, waves=None, ad=None):
         else:
             linelist = 'nearIRsky.dat'
 
-        self.log.stdinfo(f"Using linelist {linelist}")
+        self.log.debug(f"Using linelist '{linelist}'")
         filename = os.path.join(lookup_dir, linelist)
         return wavecal.LineList(filename)
 

geminidr/f2/recipes/sq/recipes_LS_SPECT.py

@@ -6,7 +6,7 @@
 
 def reduceScience(p):
     """
-    To be updated as development continues: This recipe processes F2 longslit 
+    To be updated as development continues: This recipe processes F2 longslit
     spectroscopic data, currently up to basic extraction (no telluric correction).
 
     Parameters
@@ -22,6 +22,7 @@ def reduceScience(p):
     p.darkCorrect()
     p.flatCorrect()
     p.attachWavelengthSolution()
+    p.adjustWavelengthZeroPoint()
     p.separateSky()
     p.associateSky()
     p.skyCorrect()
@@ -59,5 +60,5 @@ def  makeWavelengthSolution(p):
     p.determineDistortion(debug=True)
     p.storeProcessedArc(force=True)
     p.writeOutputs()
-    
+
 _default = reduceScience

geminidr/gemini/lookups/timestamp_keywords.py

@@ -11,6 +11,7 @@
     "addObjectMaskToDQ": "ADOBJMSK",
     "addReferenceCatalog": "ADDRECAT",
     "addVAR": "ADDVAR",
+    "adjustWavelengthZeroPoint": "CORRWVZP",
     "adjustWCSToReference": "CORRWCS",
     "ADUToElectrons": "ADUTOELE",
     "applyDQPlane": "APLDQPLN",