Merge branch 'enh/air_vacuum'

astrodata/wcs.py

@@ -14,6 +14,13 @@
 
 AffineMatrices = namedtuple("AffineMatrices", "matrix offset")
 
+FrameMapping = namedtuple("FrameMapping", "cls description")
+
+# Type of CoordinateFrame to construct for a FITS keyword and
+# readable name so user knows what's going on
+frame_mapping = {'WAVE': FrameMapping(cf.SpectralFrame, "Wavelength in vacuo"),
+                 'AWAV': FrameMapping(cf.SpectralFrame, "Wavelength in air")}
+
 #-----------------------------------------------------------------------------
 # FITS-WCS -> gWCS
 #-----------------------------------------------------------------------------
@@ -24,8 +31,6 @@ def fitswcs_to_gwcs(hdr):
     Create and return a gWCS object from a FITS header. If it can't
     construct one, it should quietly return None.
     """
-    # Type of CoordinateFrame to construct for a FITS keyword
-    frame_mapping = {'WAVE': cf.SpectralFrame}
     # coordinate names for CelestialFrame
     coordinate_outputs = {'alpha_C', 'delta_C'}
 
@@ -51,14 +56,19 @@ def fitswcs_to_gwcs(hdr):
         except TypeError:
             unit = None
         try:
-            frame = frame_mapping[output[:4].upper()](axes_order=(i,), unit=unit,
-                                          axes_names=(output,), name=output)
+            frame_type = output[:4].upper()
+            frame_info = frame_mapping[frame_type]
         except KeyError:
             if output in coordinate_outputs:
                 continue
             frame = cf.CoordinateFrame(naxes=1, axes_type=("SPATIAL",),
                                        axes_order=(i,), unit=unit,
                                        axes_names=(output,), name=output)
+        else:
+            frame = frame_info.cls(axes_order=(i,), unit=unit,
+                                   axes_names=(frame_type,),
+                                   name=frame_info.description)
+
         out_frames.append(frame)
 
     if coordinate_outputs.issubset(outputs):
@@ -175,7 +185,7 @@ def gwcs_to_fits(ndd, hdr=None):
             continue
         if axis_type == "SPECTRAL":
             wcs_dict[f'CRVAL{i}'] = hdr.get('CENTWAVE', wcs_center[i-1] if nworld_axes > 1 else wcs_center)
-            wcs_dict[f'CTYPE{i}'] = 'WAVE'
+            wcs_dict[f'CTYPE{i}'] = wcs.output_frame.axes_names[i-1]  # AWAV/WAVE
         else:  # Just something
             wcs_dict[f'CRVAL{i}'] = 0
 

geminidr/core/parameters_spect.py

@@ -92,6 +92,7 @@ class determineWavelengthSolutionConfig(config.Config):
                                            min=300., max=25000., optional=True)
     dispersion = config.Field("Estimated dispersion (nm/pixel)", float, None, optional=True)
     linelist = config.Field("Filename of arc line list", str, None, optional=True)
+    in_vacuo = config.Field("Use vacuum wavelength scale (rather than air)?", bool, False)
     alternative_centers = config.Field("Try alternative wavelength centers?", bool, False)
     debug = config.Field("Make diagnostic plots?", bool, False)
 

geminidr/gmos/lookups/CuAr_GMOS.dat

@@ -1,5 +1,6 @@
 # CuAr Line List: 3094A-10470A
-# units Angstroms
+# units Angstrom
+# wavelengths in AIR
 #3093.4019
 #3139.0176
 #3161.3726

geminidr/gmos/lookups/CuAr_GMOS_mixord.dat

@@ -1,5 +1,6 @@
 # CuAr Line List: 3094A-10470A
-# units Angstroms
+# units Angstrom
+# wavelengths in AIR
 #3093.4019
 #3139.0176
 #3161.3726

geminidr/gmos/primitives_gmos_spect.py

@@ -449,11 +449,11 @@ def findAcquisitionSlits(self, adinputs=None, **params):
             ad.update_filename(suffix=params["suffix"], strip=True)
         return adinputs
 
-    def _get_arc_linelist(self, ext, w1=None, w2=None, dw=None):
+    def _get_arc_linelist(self, ext, w1=None, w2=None, dw=None, in_vacuo=False):
         use_second_order = w2 > 1000 and abs(dw) < 0.2
         use_second_order = False
-        lookup_dir = os.path.dirname(import_module('.__init__', self.inst_lookups).__file__)
+        lookup_dir = os.path.dirname(import_module('.__init__',
+                                                   self.inst_lookups).__file__)
         filename = os.path.join(lookup_dir,
                                 'CuAr_GMOS{}.dat'.format('_mixord' if use_second_order else ''))
-        wavelengths = np.loadtxt(filename, usecols=[0])
-        return wavelengths, None
+        return self._read_and_convert_linelist(filename, w2=w2, in_vacuo=in_vacuo)

geminidr/gmos/tests/spect/test_determine_wavelength_solution.py

@@ -33,10 +33,13 @@
 import os
 import tarfile
 import logging
+from copy import deepcopy
 
 import numpy as np
 import pytest
 from matplotlib import pyplot as plt
+from astropy import units as u
+from specutils.utils.wcs_utils import air_to_vac
 
 import astrodata
 import geminidr
@@ -211,6 +214,30 @@ def test_regression_determine_wavelength_solution(
         do_plots(wcalibrated_ad)
 
 
+# We only need to test this with one input
+@pytest.mark.gmosls
+@pytest.mark.preprocessed_data
+@pytest.mark.parametrize("ad, fwidth, order, min_snr", input_pars[:1],
+                         indirect=True)
+def test_consistent_air_and_vacuum_solutions(ad, fwidth, order, min_snr):
+    p = GMOSLongslit([])
+    p.viewer = geminidr.dormantViewer(p, None)
+
+    ad_air = p.determineWavelengthSolution(
+        [deepcopy(ad)], order=order, min_snr=min_snr, fwidth=fwidth,
+        in_vacuo=False, **determine_wavelength_solution_parameters).pop()
+    ad_vac = p.determineWavelengthSolution(
+        [ad], order=order, min_snr=min_snr, fwidth=fwidth,
+        in_vacuo=True, **determine_wavelength_solution_parameters).pop()
+    wave_air = am.get_named_submodel(ad_air[0].wcs.forward_transform, "WAVE")
+    wave_vac = am.get_named_submodel(ad_vac[0].wcs.forward_transform, "WAVE")
+    x = np.arange(ad_air[0].shape[1])
+    wair = wave_air(x)
+    wvac = air_to_vac(wair * u.nm).to(u.nm).value
+    dw = wvac - wave_vac(x)
+    assert abs(dw).max() < 0.001
+
+
 # Local Fixtures and Helper Functions ------------------------------------------
 @pytest.fixture(scope='function')
 def ad(path_to_inputs, request):

gempy/library/transform.py

@@ -1439,8 +1439,8 @@ def add_longslit_wcs(ad):
                   for kw in ('reference_frame', 'unit', 'axes_names',
                              'axis_physical_types')}
         sky_frame = cf.CelestialFrame(axes_order=(1,2), name='sky', **kwargs)
-        spectral_frame = cf.SpectralFrame(name='wavelength', unit=u.nm,
-                                          axes_names='WAVE')
+        spectral_frame = cf.SpectralFrame(name='Wavelength in air', unit=u.nm,
+                                          axes_names='AWAV')
         output_frame = cf.CompositeFrame([spectral_frame, sky_frame], name='world')
 
         transform = ext.wcs.forward_transform