Various work on distortion modelling

Change min_line_length to be a fraction of the slit length
Move GNIRS determineDistortion() wrapper to GNIRSSpect from GNIRSLongslit and handle XD
Put distortion model *after* rectification model and calculate it correctly
Allow peaks to be found in the non-linear regime

Author: chris-simpson

geminidr/core/parameters_spect.py

@@ -145,8 +145,8 @@ class determineDistortionConfig(config.Config):
     max_shift = config.RangeField("Maximum shift per pixel in line position",
                                   float, 0.05, min=0.001, max=0.1)
     max_missed = config.RangeField("Maximum number of steps to miss before a line is lost", int, 5, min=0)
-    min_line_length = config.RangeField("Exclude line traces shorter than this fraction of spatial dimension",
-                                        float, 0., min=0., max=1.)
+    min_line_length = config.RangeField("Exclude line traces shorter than this fraction of slit length",
+                                        float, 0.8, min=0., max=1.)
     debug_reject_bad = config.Field("Reject lines with suspiciously high SNR (e.g. bad columns)?", bool, True)
     debug = config.Field("Display line traces on image display?", bool, False)
 

geminidr/core/primitives_spect.py

@@ -1473,67 +1473,105 @@ def determineDistortion(self, adinputs=None, **params):
                     initial_peaks, _ = peak_finding.find_wavelet_peaks(
                         data, widths=widths, mask=mask & DQ.not_signal,
                         variance=variance, min_snr=min_snr, reject_bad=debug_reject_bad)
-                    log.stdinfo(f"Found {len(initial_peaks)} peaks")
                 # The coordinates are always returned as (x-coords, y-coords)
                 rwidth = 0.42466 * fwidth
 
-                # Straight slits, such as in longslit, can have all the lines
-                # traced simultaneously since they all have the same starting
-                # point. "Curved" slits need to be handled one-by-one. This is
-                # quite a bit slower, so this block of code does the line
-                # tracing based on the slit involved.
-                if constant_slit:
-                    traces = tracing.trace_lines(
-                        # Only need a single `start` value for all lines.
-                        ext, axis=1 - dispaxis,
-                        start=start, initial=initial_peaks,
-                        rwidth=rwidth, cwidth=max(int(fwidth), 5), step=step,
-                        nsum=nsum, max_missed=max_missed,
-                        max_shift=max_shift * ybin / xbin,
-                        viewer=self.viewer if debug else None,
-                        min_line_length=min_line_length)
-
-                else:
-                    traces = []
-                    for peak in initial_peaks:
-                        # Need to start midway along the slit, which varies
-                        # along the dispersion axis. `extract_info` here is the
-                        # polynomial describing that midway line.
-                        start = extract_info(peak)
-                        traces.extend(tracing.trace_lines(
+                if len(initial_peaks):
+                    # The slit length may be smaller than the width of the slice,
+                    # so we need to estimate the slit length and
+                    # "min_line_length" is the fraction of that, not the fraction
+                    # of the slice width.
+                    if ext.mask is not None:
+                        loc = int(np.median(initial_peaks))
+                        if dispaxis == 0:
+                            _slice = ext.mask[loc] & DQ.unilluminated
+                        else:
+                            _slice = ext.mask[:, loc] & DQ.unilluminated
+                        slit_length_frac = 1 - ((_slice.argmin() +
+                                                 _slice[::-1].argmin()) / _slice.size)
+                    else:
+                        try:
+                            slit_length_frac = ad.MDF['slitlength_pixels'] / ext.shape[1 - dispaxis]
+                        except (AttributeError, KeyError):
+                            slit_length_frac = 1
+
+                    # Straight slits, such as in longslit, can have all the lines
+                    # traced simultaneously since they all have the same starting
+                    # point. "Curved" slits need to be handled one-by-one. This is
+                    # quite a bit slower, so this block of code does the line
+                    # tracing based on the slit involved.
+                    if constant_slit:
+                        traces = tracing.trace_lines(
+                            # Only need a single `start` value for all lines.
                             ext, axis=1 - dispaxis,
-                            start=start, initial=[peak],
+                            start=start, initial=initial_peaks,
                             rwidth=rwidth, cwidth=max(int(fwidth), 5), step=step,
                             nsum=nsum, max_missed=max_missed,
                             max_shift=max_shift * ybin / xbin,
                             viewer=self.viewer if debug else None,
-                            min_line_length=0.1))
+                            min_line_length=min_line_length*slit_length_frac)
 
-                # List of traced peak positions
-                in_coords = np.array([coord for trace in traces for
-                                      coord in trace.input_coordinates()]).T
+                    else:
+                        traces = []
+                        for peak in initial_peaks:
+                            # Need to start midway along the slit, which varies
+                            # along the dispersion axis. `extract_info` here is the
+                            # polynomial describing that midway line.
+                            start = extract_info(peak)
+                            traces.extend(tracing.trace_lines(
+                                ext, axis=1 - dispaxis,
+                                start=start, initial=[peak],
+                                rwidth=rwidth, cwidth=max(int(fwidth), 5), step=step,
+                                nsum=nsum, max_missed=max_missed,
+                                max_shift=max_shift * ybin / xbin,
+                                viewer=self.viewer if debug else None,
+                                min_line_length=min_line_length*slit_length_frac))
+
+                    log.stdinfo(f"Traced {len(traces)} lines from "
+                                f"{len(initial_peaks)} peaks")
+
+                    # List of traced peak positions
+                    in_coords = np.array([coord for trace in traces for
+                                          coord in trace.input_coordinates()]).T
+
+                    # We can't do anything if we have no coordinates
+                    if in_coords.size == 0:
+                        log.warning("Failed to trace any lines for "
+                                    f"{ad.filename}:{ext.id}")
+                        continue
 
-                # We can't do anything if we have no coordinates
-                if in_coords.size == 0:
-                    log.warning("Failed to trace any lines for "
-                                f"{ad.filename}:{ext.id}")
+                else:
+                    log.warning("Failed to find any peaks in "
+                                "f{ad.filename}:{ext.id}")
                     continue
 
-                # If there's a "rectified" frame, we want to use the pixel
-                # coordinates in *that* frame as input so that the pixels
-                # -> rectified -> distortion_corrected transform works
-                # correctly.
-                try:
-                    t = ext.wcs.get_transform(ext.wcs.input_frame, 'rectified')
-                except CoordinateFrameError:
-                    pass
-                else:
-                    in_coords = np.array(t(*in_coords))
                 # List of "reference" positions (i.e., the coordinate
                 # perpendicular to the line remains constant at its initial value
                 ref_coords = np.array([coord for trace in traces for
                                        coord in trace.reference_coordinates()]).T
 
+                # If the frame has a rectification model, then we want to
+                # calculate the distortion transform *after* applying this
+                # model. This is important because, if we only have one line
+                # from which to determine the distortion, the model will only
+                # be a function of X and so we need a vertical spectrum.
+                try:
+                    rect_model = ext.wcs.get_transform(ext.wcs.input_frame,
+                                                       "rectified")
+                except CoordinateFrameError:
+                    has_rect_model = False
+                else:
+                    has_rect_model = True
+                    in_coords = rect_model(*in_coords)
+                    # For a vertically-dispersed spectrum, the rectification
+                    # model alters X as a function of Y. Therefore, because
+                    # incoords and ref_coords have different Y for the same X,
+                    # putting them *both* through the rectification model will
+                    # produce different X values, which is not what we want.
+                    # S0 replace the X values in ref_coords with the ones in
+                    # in_coords. Coords are *always* (x, y)
+                    ref_coords[dispaxis] = in_coords[dispaxis]
+
                 # The model is computed entirely in the pixel coordinate frame
                 # of the data, so it could be used as a gWCS object
                 m_init = models.Chebyshev2D(x_degree=orders[1 - dispaxis],
@@ -1573,8 +1611,25 @@ def determineDistortion(self, adinputs=None, **params):
                     ext.wcs = gWCS([(cf.Frame2D(name="pixels"), model),
                                     (cf.Frame2D(name="world"), None)])
                 else:
-                    ext.wcs.insert_frame(ext.wcs.input_frame, model,
-                                         cf.Frame2D(name="distortion_corrected"))
+                    try:
+                        frame_index = ext.wcs.available_frames.index("distortion_corrected")
+                    except ValueError:
+                        pass
+                    else:
+                        log.warning("Deleting existing distortion model in "
+                                    f"{ad.filename}:{ext.id}")
+                        ext.wcs = ext.wcs.__class__(
+                            ext.wcs.pipeline[:frame_index-1] +
+                            [(ext.wcs.pipeline[frame_index-1].frame,
+                              ext.wcs.pipeline[frame_index].transform)] +
+                            ext.wcs.pipeline[frame_index+1:]
+                        )
+                    if has_rect_model:
+                        ext.wcs.insert_frame("rectified", model,
+                                             cf.Frame2D(name="distortion_corrected"))
+                    else:
+                        ext.wcs.insert_frame(ext.wcs.input_frame, model,
+                                             cf.Frame2D(name="distortion_corrected"))
 
                 nsuccess += 1
 

geminidr/f2/tests/longslit/test_determine_distortion.py

@@ -37,7 +37,7 @@
     "nsum": 10,
     "spatial_order": 3,
     "spectral_order": 3,
-    "min_line_length": 0.3,
+    "min_line_length": 0.8,
     "debug_reject_bad": False
 }
 
@@ -186,7 +186,6 @@
 # Tests Definitions ------------------------------------------------------------
 @pytest.mark.f2ls
 @pytest.mark.preprocessed_data
-@pytest.mark.regression
 @pytest.mark.parametrize("ad,params", input_pars, indirect=['ad'])
 def test_regression_for_determine_distortion_using_wcs(
         ad, params, change_working_dir, ref_ad_factory):
@@ -214,18 +213,24 @@ def test_regression_for_determine_distortion_using_wcs(
         distortion_determined_ad = p.writeOutputs().pop()
 
     ref_ad = ref_ad_factory(distortion_determined_ad.filename)
-    model = distortion_determined_ad[0].wcs.get_transform(
-        "pixels", "distortion_corrected")[2]
-    ref_model = ref_ad[0].wcs.get_transform("pixels", "distortion_corrected")[2]
+
+    # Confirm that the distortion model is placed after the rectification model
+    assert (distortion_determined_ad[0].wcs.available_frames.index("distortion_corrected") >
+            distortion_determined_ad[0].wcs.available_frames.index("rectified"))
+    # assert (ref_ad[0].wcs.available_frames.index("distortion_corrected") >
+    #         ref_ad[0].wcs.available_frames.index("rectified"))
+
+    model = distortion_determined_ad[0].wcs.get_transform("pixels", "distortion_corrected")
+    ref_model = ref_ad[0].wcs.get_transform("pixels", "distortion_corrected")
 
     # Otherwise we're doing something wrong!
-    assert model.__class__.__name__ == ref_model.__class__.__name__ == "Chebyshev2D"
+    assert model[-1].__class__.__name__ == ref_model[-1].__class__.__name__ == "Chebyshev2D"
 
-    X, Y = np.mgrid[:ad[0].shape[0], :ad[0].shape[1]]
+    Y, X = np.mgrid[:ad[0].shape[0], :ad[0].shape[1]]
 
-    # Increasing atol to 0.07 due to S20180114S0104_flatCorrected.fits producing
-    # slightly different results on Jenkins vs. on MacOS. DB 20240820
-    np.testing.assert_allclose(model(X, Y), ref_model(X, Y), atol=0.07)
+    xx, yy = X[ad[0].mask == 0], Y[ad[0].mask == 0]
+    diffs = model(xx, yy)[1] - ref_model(xx, yy)[1]  # 1 is y-axis in astropy
+    np.testing.assert_allclose(diffs, 0, atol=1)
 
 
 @pytest.mark.f2ls
@@ -252,8 +257,9 @@ def test_fitcoord_table_and_gwcs_match(ad, params, change_working_dir):
         p.determineDistortion(**fixed_parameters_for_determine_distortion)
         distortion_determined_ad = p.writeOutputs().pop()
 
-    model = distortion_determined_ad[0].wcs.get_transform(
-        "pixels", "distortion_corrected")
+    model = distortion_determined_ad[0].wcs.pipeline[
+        distortion_determined_ad[0].wcs.available_frames.index(
+            "distortion_corrected") - 1].transform
 
     fitcoord = distortion_determined_ad[0].FITCOORD
     fitcoord_model = am.table_to_model(fitcoord[0])

geminidr/gnirs/parameters_gnirs_crossdispersed.py

@@ -4,6 +4,7 @@
 from geminidr.core import parameters_crossdispersed
 from geminidr.core import parameters_standardize
 from geminidr.core.parameters_standardize import addIllumMaskToDQConfig
+from . import parameters_gnirs_spect
 from gempy.library import config
 
 
@@ -17,6 +18,12 @@ def setDefaults(self):
         self.add_illum_mask = True
 
 
+class determineDistortionConfig(parameters_gnirs_spect.determineDistortionConfig):
+    def setDefaults(self):
+        self.spatial_order = 2
+        self.min_line_length = 0.5  # because some orders go off the edge
+
+
 class determineSlitEdgesConfig(parameters_spect.determineSlitEdgesConfig):
     # GNIRS XD has narrow slits with more curvature than the longslit flats
     # the default values were calibrated to, so adjust some values.

geminidr/gnirs/parameters_gnirs_longslit.py

@@ -11,16 +11,6 @@ def list_of_ints_check(value):
     [int(x) for x in str(value).split(',')]
     return True
 
-class determineDistortionConfig(parameters_spect.determineDistortionConfig):
-    spectral_order = config.RangeField("Fitting order in spectral direction", int, None, min=1, optional=True)
-    min_line_length = config.RangeField("Exclude line traces shorter than this fraction of spatial dimension",
-                                        float, None, min=0., max=1., optional=True)
-    max_missed = config.RangeField("Maximum number of steps to miss before a line is lost",
-                               int, None, min=0, optional=True)
-    def setDefaults(self):
-        self.min_snr = 10
-        self.debug_reject_bad = False
-
 class determineWavelengthSolutionConfig(parameters_spect.determineWavelengthSolutionConfig):
     order = config.RangeField("Order of fitting function", int, None, min=1,
                               optional=True)

geminidr/gnirs/parameters_gnirs_spect.py

@@ -2,7 +2,7 @@
 # in the primitives_gnirs_spect.py file, in alphabetical order.
 from astrodata import AstroData
 from gempy.library import config
-from geminidr.core import parameters_preprocess
+from geminidr.core import parameters_preprocess, parameters_spect
 from . import parameters_gnirs
 
 
@@ -16,6 +16,15 @@ def setDefaults(self):
         self.min_skies = 2
 
 
+class determineDistortionConfig(parameters_spect.determineDistortionConfig):
+    spectral_order = config.RangeField("Fitting order in spectral direction", int, None, min=1, optional=True)
+    max_missed = config.RangeField("Maximum number of steps to miss before a line is lost",
+                               int, None, min=0, optional=True)
+    def setDefaults(self):
+        self.min_snr = 10
+        self.debug_reject_bad = False
+
+
 class skyCorrectConfig(parameters_preprocess.skyCorrectConfig):
     def setDefaults(self):
         # self.scale_sky = False #MS: IF for whatever reason the exposure times are different between frames being subtracted, that case may require a special treatment