Merge branch 'release/4.x' into enh/gnirs_nonlinearity

Author: chris-simpson

geminidr/core/__init__.py

@@ -10,6 +10,7 @@
 from .primitives_standardize import Standardize
 from .primitives_telluric import Telluric
 from .primitives_visualize import Visualize
+from .primitives_stats import Stats
 
 # These are currently inherited by Image
 #from geminidr.core.primitives_register import Register

geminidr/core/parameters_spect.py

@@ -563,7 +563,7 @@ class resampleToCommonFrameConfig(config.Config):
     output_wave_scale = config.ChoiceField("Output wavelength scale", str,
                                            allowed={"reference": "Reference input",
                                                     "linear": "Linear",
-                                                    #"loglinear": "Log-linear",
+                                                    "loglinear": "Log-linear",
                                                     },
                                            default="linear", optional=False)
     dq_threshold = config.RangeField("Fraction from DQ-flagged pixel to count as 'bad'",

geminidr/core/parameters_stats.py

@@ -0,0 +1,6 @@
+from gempy.library import config
+
+
+class statsConfig(config.Config):
+    prefix = config.Field('4-character header prefix', str, None, optional=True)
+

geminidr/core/primitives_spect.py

@@ -435,7 +435,7 @@ def slit_pa(ext):
                         else:
                             log.warning(f"{ad.filename}:{ext.id} Cross-correlation failed")
 
-                        if debug_plots:
+                        if debug_plots:  # pragma: no cover
                             fig, ax = plt.subplots()
                             print(f"Comparing {ad.filename} to reference {refad.filename} "
                                   f"(expected peak at {expected_peak - corr.size // 2})")
@@ -1545,7 +1545,7 @@ def determineDistortion(self, adinputs=None, **params):
                 # TODO: Some logging about quality of fit
                 # print(np.min(diff), np.max(diff), np.std(diff))
 
-                if debug:
+                if debug:  # pragma: no cover
                     self.viewer.color = "red"
                     spatial_coords = np.linspace(ref_coords[dispaxis].min(), ref_coords[dispaxis].max(),
                                                 ext.shape[1 - dispaxis] // (step * 10))
@@ -1924,7 +1924,7 @@ def find_slits(mdf, ystep=50):
                     log.warning(f"Did not find expected number of {slit_name}s "
                                 f"(found {nfound}, expected {len(slit_lengths)}).")
 
-                if debug_plots:
+                if debug_plots:  # pragma: no cover
                     for pos in edges_1:
                         if pos:
                             plt.axvline(pos, color='Blue', alpha=0.5,
@@ -3151,7 +3151,7 @@ def flagCosmicRays(self, adinputs=None, **params):
                 del crmask, skyfit_input
                 gc.collect()
 
-            if debug:
+            if debug:  # pragma: no cover
                 fig, axes = plt.subplots(5, 3, sharex=True, sharey=True,
                                          tight_layout=True)
                 for i, ext in enumerate(ad_tiled):
@@ -3387,10 +3387,18 @@ def fluxCalibrate(self, adinputs=None, **params):
                 # Reconstruct the spline and evaluate it at every wavelength
                 sens_factor = sensfunc(waves.to(std_wave_unit).value) * std_flux_unit
                 try:  # conversion from magnitude/logarithmic units
-                    sens_factor = sens_factor.physical
+                    # See comment below
+                    with warnings.catch_warnings(category=RuntimeWarning,
+                                                 action="ignore"):
+                        sens_factor = sens_factor.physical
                 except AttributeError:
                     pass
 
+                # This avoids extrapolative blow-ups when flux-calibrating XD
+                # data that covers the full wavelength range but the standard
+                # only covers a small part.
+                sens_factor[(ext.mask & DQ.no_data) > 0] = 0
+
                 # Apply airmass correction. If none is needed/possible, we
                 # don't need to try to do this
                 if delta_airmass:
@@ -3942,6 +3950,8 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                                        r + ref_pixels_dict[j][0]).T.astype(int)
                                        for ad, r in zip(adinputs,
                                                         ref_pixels_dict[j])]
+        else:
+            dispaxis = 0
 
         # Gather information from all the spectra (Chebyshev1D model,
         # w1, w2, dw, npix), and compute the final bounds (w1out, w2out)
@@ -3963,9 +3973,13 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                 adinfo.append(model_info)
                 w1_arr[i, iext] = model_info['w1']
                 w2_arr[i, iext] = model_info['w2']
-
             info.append(adinfo)
 
+        # Are we combining multiple spectra with different wavelength settings
+        # into a single spectrum? This is important for later.
+        combining_multiple_wavelengths = (single_spectral and
+                                          len(set(w1_arr.ravel())) > 1)
+
         # Compute the output wavelength range for each extension. We can
         # calculate the overall output range if we're combining to a single
         # extension later
@@ -3994,15 +4008,16 @@ def resampleToCommonFrame(self, adinputs=None, **params):
             # parameters as the 4th is then calculable. First, we copy the
             # start and end wavelengths if those aren't specified. If neither
             # dw nor npix are specified, the behaviour depends on whether we
-            # are resampling to a single wavelength scale: if so, then we want
-            # to preserve the dispersion to avoid undersampling but, if not,
-            # then we want to preserve the number of pixels per extension.
+            # are resampling multiple spectra to a single wavelength scale:
+            # if so, then we want to preserve the dispersion to avoid
+            # undersampling but, if not,  then we want to preserve the number
+            # of pixels per extension.
             while nparams < 3:
                 if w1 is None:
                     w1 = wave_min
                 elif w2 is None:
                     w2 = wave_max
-                elif single_spectral and dw is None:
+                elif combining_multiple_wavelengths and dw is None:
                     w1 = np.full_like(w1, np.nanmin(w1))
                     w2 = np.full_like(w2, np.nanmax(w2))
                     if output_spectral == "linear":
@@ -4011,7 +4026,7 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                     else:
                         # dw has been calculated assuming the spectrum is
                         # linear, so we repeat that assumption
-                        dw = np.array([extinfo['dw'] / extinfo['w2'] - 1
+                        dw = np.array([extinfo['dw'] / extinfo['w2']
                                        for adinfo in info for extinfo in adinfo])
                     dw = np.full_like(w1, dw.min())
                 elif npix is None:
@@ -4025,12 +4040,18 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                     npix = np.ceil((w2 - w1) / dw).astype(int) + 1
                     w2 = w1 + (npix - 1) * dw
                 else:  # loglinear
-                    npix = np.ceil(np.log(w2 / w1) / np.log(1 + dw) - 1)
+                    npix = np.ceil(np.log(w2 / w1) / np.log(1 + dw) - 1).astype(int) + 1
                     w2 = w1 * (1 + dw) ** (npix - 1)
             elif w1 is None:
-                w1 = w2 - (npix - 1) * dw
+                if output_spectral == "linear":
+                    w1 = w2 - (npix - 1) * dw
+                else:  # loglinear
+                    w1 = w2 / (1 + dw) ** (npix - 1)
             elif w2 is None:
-                w2 = w1 + (npix - 1) * dw
+                if output_spectral == "linear":
+                    w2 = w1 + (npix - 1) * dw
+                else:  # loglinear
+                    w2 = w1 * (1 + dw) ** (npix - 1)
             elif output_spectral == "linear":  # dw is None
                 dw = (w2 - w1) / (npix - 1)
             else:  # dw is None and we're loglinearizing
@@ -4084,7 +4105,6 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                 actual_limits = new_wave_model([0, this_npix - 1])
                 w1[iext] = actual_limits.min()
                 w2[iext] = actual_limits.max()
-                yy = new_wave_model([this_npix-3,this_npix-2,this_npix-1])
 
             # Calculation for all extensions
             dw = (w2 - w1) / (npix - 1)
@@ -4164,6 +4184,8 @@ def resampleToCommonFrame(self, adinputs=None, **params):
                 if i == 0 and not new_wave_scale:
                     log.fullinfo(f"{ad.filename}: No interpolation")
                 msg = "Resampling"
+                if this_conserve:
+                    msg += " (with flux conservation)"
                 if new_wave_scale:
                     msg += f" and {output_spectral}izing"
                 dwstr = (f"{dw[iext]:.6f}" if output_spectral == "loglinear"
@@ -5220,7 +5242,8 @@ def _get_spectrophotometry(self, filename, in_vacuo=False):
             if isinstance(unit, u.UnrecognizedUnit):
                 # Try chopping off the trailing 's'
                 try:
-                    unit = u.Unit(re.sub(r's$', '', col.unit.name.lower()))
+                    unit = u.Unit(re.sub(r's$', '',
+                                         col.unit.name.lower()))
                 except:
                     unit = None
             if unit is None:
@@ -5232,7 +5255,8 @@ def _get_spectrophotometry(self, filename, in_vacuo=False):
                 else:
                     if orig_colname == 'FNU':
                         unit = u.Unit("erg cm-2 s-1 Hz-1")
-                    elif orig_colname in ('FLAM', 'FLUX') or np.median(col.data) < 1:
+                    elif (orig_colname in ('FLAM', 'FLUX') or
+                          np.median(col.data) < 1):
                         unit = u.Unit("erg cm-2 s-1 AA-1")
                     else:
                         unit = u.mag
@@ -5241,15 +5265,21 @@ def _get_spectrophotometry(self, filename, in_vacuo=False):
             # We've created a column called "MAGNITUDE" but it might be a flux
             if col.name == 'MAGNITUDE':
                 try:
-                    unit.to(u.W / u.m ** 3, equivalencies=u.spectral_density(1. * u.m))
+                    unit.to(u.W / u.m ** 3,
+                            equivalencies=u.spectral_density(1. * u.m))
                 except:
                     pass
                 else:
                     col.name = 'FLUX'
 
-        wavecol = spec_table["WAVELENGTH"].quantity
         if in_vacuo is None:
-            in_vacuo = min(wavecol) < 300 * u.nm
+            in_vacuo = min(spec_table["WAVELENGTH"].quantity) < 300 * u.nm
+
+        # The default (and best) specutils vacuum/air conversion has a
+        # singularity in the FUV, so we cut the wavelength scale.
+        # See https://github.com/astropy/specutils/issues/1162
+        spec_table = spec_table[spec_table["WAVELENGTH"] >= 300 * u.nm]
+        wavecol = spec_table["WAVELENGTH"].quantity
 
         if in_vacuo:
             spec_table["WAVELENGTH_VACUUM"] = spec_table["WAVELENGTH"]

geminidr/core/primitives_stats.py

@@ -0,0 +1,63 @@
+from copy import deepcopy
+
+import numpy as np
+
+from geminidr import PrimitivesBASE
+from gempy.gemini import gemini_tools as gt
+from recipe_system.utils.decorators import parameter_override, capture_provenance
+
+from . import parameters_stats
+
+@parameter_override
+class Stats(PrimitivesBASE):
+
+    def _initialize(self, adinputs, **kwargs):
+        super()._initialize(adinputs, **kwargs)
+        self._param_update(parameters_stats)
+
+    def stats(self, adinputs=None, **params):
+        """
+        Adds headers to the AD object giving some statistics of the unmasked
+        pixel values
+        By default, adds:
+        PIXMEAN - the arithmetic mean of the pixel values
+        PIXSTDV - the standard deviation of the pixel values
+        PIXMED - the median of the pixel values.
+
+        Parameters
+        ----------
+
+        adinputs: list of :class:`~astrodata.AstroData`
+
+        prefix: Prefix for header keywords. Maximum of 4 characters, defaults
+                to PIX.
+        """
+        log = self.log
+        log.debug(gt.log_message("primitive", self.myself(), "starting"))
+
+        pre = params['prefix']
+
+        if pre is None:
+            pre = 'PIX'
+        else:
+            pre = pre[:4]
+
+        for ad in adinputs:
+            for ext in ad:
+                try:
+                    ext.hdr[pre+'MEAN'] = (np.mean(ext.data[ext.mask==0]),
+                                          "Mean of unmasked pixel values")
+                except ValueError:
+                    # Things like NaN can't be written to FITS headers
+                    pass
+                try:
+                    ext.hdr[pre+'STDV'] = (np.std(ext.data[ext.mask==0]),
+                                           "Standard Deviation of pixel values")
+                except ValueError:
+                    pass
+                try:
+                    ext.hdr[pre+'MED'] = (np.median(ext.data[ext.mask == 0]),
+                                         "Median of unmasked of pixel values")
+                except ValueError:
+                    pass
+        return adinputs

geminidr/core/tests/test_spect.py

@@ -38,12 +38,14 @@
 from astropy.modeling import models
 from scipy import optimize
 
+from gwcs import coordinate_frames as cf
+from gwcs.wcs import WCS as gWCS
+
 from specutils.utils.wcs_utils import air_to_vac
 
 import astrodata, gemini_instruments
-from astrodata.testing import ad_compare
+from astrodata.testing import ad_compare, assert_most_close
 from gempy.library import astromodels as am
-from gempy.library.wavecal import LineList
 from gempy.library.config.config import FieldValidationError
 from geminidr.core import primitives_spect
 from geminidr.f2.primitives_f2_longslit import F2Longslit
@@ -624,6 +626,40 @@ def test_slit_rectification(filename, instrument, change_working_dir,
         np.testing.assert_allclose(ad_out[0].SLITEDGE[coeff], 0, atol=0.25)
 
 
+@pytest.mark.parametrize("gnirs1d", [np.ones((1000,), dtype=np.float32),
+                                     np.arange(1000, dtype=np.float32)],
+                         indirect=True)
+@pytest.mark.parametrize("wavescale", ["linear", "loglinear"])
+def test_resample1d_conserve(gnirs1d, wavescale):
+    """
+    Simple test to resample a synthetic 1D spectrum with a linear wavelength
+    solution to linear and loglinear and check that the flux is conserved.
+    """
+    gnirs1d[0].hdr['BUNIT'] = 'electron'
+    p = GNIRSLongslit([gnirs1d])
+    sum_before = gnirs1d[0].data.sum()
+    adout = p.resampleToCommonFrame(output_wave_scale=wavescale).pop()
+    sum_after = adout[0].data.sum()
+    # Tolerance allows for edge effects
+    assert sum_after == pytest.approx(sum_before, rel=0.002)
+
+
+@pytest.mark.parametrize("gnirs1d", [np.arange(1000, 2000, dtype=np.float32)],
+                         indirect=True)
+@pytest.mark.parametrize("wavescale", ["linear", "loglinear"])
+def test_resample1d_interpolate(gnirs1d, wavescale):
+    """
+    Simple test to resample a synthetic 1D spectrum with a linear wavelength
+    solution to linear and loglinear and check that the data are interpolated.
+    The input spectrum has signal=wavelength so it's easy to check the result.
+    """
+    gnirs1d[0].hdr['BUNIT'] = 'W / (m2 AA)'  # will interpolate
+    p = GNIRSLongslit([gnirs1d])
+    adout = p.resampleToCommonFrame(output_wave_scale=wavescale).pop()
+    waves = adout[0].wcs(np.arange(adout[0].data.size))
+    assert_most_close(waves, adout[0].data, max_miss=2, rtol=1e-5)
+
+
 def test_trace_apertures():
     # Input parameters ----------------
     width = 400
@@ -788,6 +824,22 @@ def test_transfer_distortion_model(change_working_dir, path_to_inputs, path_to_r
 
 # --- Fixtures and helper functions -------------------------------------------
 
+@pytest.fixture
+def gnirs1d(request, astrofaker):
+    ad = astrofaker.create('GNIRS', mode="SPECT")
+    ad.init_default_extensions()
+    data = request.param
+    mask = np.zeros_like(data, dtype=np.uint16)
+    variance = np.ones_like(data, dtype=np.float32)
+    ad[0].reset(data=data, mask=mask, variance=variance)
+    wave_model = models.Shift(1000) | models.Scale(1)
+    input_frame = astrodata.wcs.pixel_frame(1)
+    output_frame = cf.SpectralFrame(axes_order=(0,), unit=u.nm,
+                                    axes_names=("WAVE",))
+    ad[0].wcs = gWCS([(input_frame, wave_model),
+                      (output_frame, None)])
+    return ad
+
 
 def create_zero_filled_fake_astrodata(height, width):
     """

geminidr/gemini/primitives_gemini.py

@@ -17,7 +17,7 @@
 from gempy.library import astromodels as am, astrotools as at, transform
 
 from geminidr.core import Bookkeeping, CalibDB, Preprocess
-from geminidr.core import Visualize, Standardize, Stack
+from geminidr.core import Visualize, Standardize, Stack, Stats
 
 from .primitives_qa import QA
 from . import parameters_gemini
@@ -29,7 +29,7 @@
 @parameter_override
 @capture_provenance
 class Gemini(Standardize, Bookkeeping, Preprocess, Visualize, Stack, QA,
-             CalibDB):
+             CalibDB, Stats):
     """
     This is the class containing the generic Gemini primitives.