fix WCS issues and write combineOrders() test for GHOST

astrodata/wcs.py

@@ -203,17 +203,17 @@ def gwcs_to_fits(ndd, hdr=None):
         # Remove projection parts so we can calculate the CD matrix
         if projcode:
             nat2cel.name = 'nat2cel'
-            transform_inverse = transform.inverse.copy()
+            transform_inverse = transform.inverse
             for m in transform_inverse:
                 if isinstance(m, models.RotateCelestial2Native):
                     m.name = 'cel2nat'
                 elif isinstance(m, models.Sky2PixProjection):
                     m.name = 'sky2pix'
-            #transform_inverse = transform_inverse.replace_submodel('cel2nat', models.Identity(2))
-            #transform_inverse = transform_inverse.replace_submodel('sky2pix', models.Identity(2))
+            transform_inverse = transform_inverse.replace_submodel('sky2pix', models.Identity(2))
+            transform_inverse = transform_inverse.replace_submodel('cel2nat', models.Identity(2))
             transform = transform.replace_submodel('pix2sky', models.Identity(2))
             transform = transform.replace_submodel('nat2cel', models.Identity(2))
-            #transform.inverse = transform_inverse
+            transform.inverse = transform_inverse
 
     # Replace a log-linear axis with a linear axis representing the log
     # and a Tabular axis with Identity to ensure the affinity check is passed
@@ -292,15 +292,11 @@ def gwcs_to_fits(ndd, hdr=None):
                      if f'CTYPE{i}' not in wcs_dict})
 
     crval = [wcs_dict[f'CRVAL{i+1}'] for i, _ in enumerate(world_axes)]
-    #try:
-    #    crval[lon_axis] = 0
-    #    crval[lat_axis] = 0
-    #except NameError:
-    #    pass
-
-    # This (commented) line fails for un-invertable Tabular2D
-    crpix = np.array(wcs.backward_transform(*crval)) + 1
-    #crpix = np.array(transform.inverse(*crval)) + 1
+    try:
+        crval[lon_axis] = 0
+        crval[lat_axis] = 0
+    except NameError:
+        pass
 
     # Find any world axes that we previous logarithmed and fix the CDij
     # matrix -- we follow FITS-III (Greisen et al. 2006; A&A 446, 747)
@@ -315,6 +311,10 @@ def gwcs_to_fits(ndd, hdr=None):
                 wcs_dict[f'CTYPE{world_axis}'] = wcs_dict[f'CTYPE{world_axis}'][:4] + "-LOG"
             crval[world_axis-1] = np.log(crval[world_axis-1])
 
+    # This (commented) line fails for un-invertable Tabular2D
+    #crpix = np.array(wcs.backward_transform(*crval)) + 1
+    crpix = np.array(transform.inverse(*crval)) + 1
+
     # Cope with a situation where the sky projection center is not in the slit
     # We may be able to fix this in future, but FITS doesn't handle it well.
     if len(ndd.shape) > 1:

geminidr/ghost/tests/spect/test_combine_orders.py

@@ -0,0 +1,66 @@
+import pytest
+import os
+
+import numpy as np
+
+import astrodata, gemini_instruments
+
+from geminidr.ghost.primitives_ghost_spect import GHOSTSpect, make_wavelength_table
+
+
+FILENAMES = ["S20230130S0104_2x8_blue001_wavelengthSolutionAttached.fits",
+             "S20230130S0104_2x8_red001_wavelengthSolutionAttached.fits"]
+
+
+@pytest.mark.ghostspect
+@pytest.mark.parametrize("filename", FILENAMES)
+def test_combine_orders_single_file(change_working_dir, path_to_inputs, filename):
+    """Check that combineOrders() works on a single file"""
+    ad = astrodata.open(os.path.join(path_to_inputs, filename))
+    orig_wavl = make_wavelength_table(ad[0])
+    p = GHOSTSpect([ad])
+    ad_out = p.combineOrders().pop()
+    pixels = np.arange(ad_out[0].data.size)
+    x = ad_out[0].wcs(pixels[1:]) / ad_out[0].wcs(pixels[:-1])
+    assert x.std() < 1e-12  # all wavelength ratios are the same
+    assert ad_out[0].wcs(0) == pytest.approx(orig_wavl.min())
+
+    # Check that we can write and read back the _ordersCombined file
+    with change_working_dir():
+        ad_out.write("test.fits", overwrite=True)
+        ad2 = astrodata.open("test.fits")
+        np.testing.assert_allclose(ad_out[0].wcs(pixels), ad2[0].wcs(pixels))
+
+
+@pytest.mark.ghostspect
+def test_combine_orders_red_and_blue_no_stacking(path_to_inputs):
+    """Check that combineOrders() works on two files without stacking"""
+    adinputs = [astrodata.open(os.path.join(path_to_inputs, filename))
+                for filename in FILENAMES]
+    p = GHOSTSpect(adinputs)
+    adoutputs = p.combineOrders(stacking_mode="none")
+    assert len(adoutputs) == len(adinputs)  # no stacking
+
+    # Check that the wavelength solutions haven't merged in some way
+    for adin, adout in zip(adinputs, adoutputs):
+        orig_wavl = make_wavelength_table(adin[0])
+        assert adout[0].wcs(0) == pytest.approx(orig_wavl.min())
+
+
+@pytest.mark.ghostspect
+@pytest.mark.parametrize("stacking_mode", ("scaled", "unscaled"))
+def test_combine_orders_red_and_blue_stacking(path_to_inputs, stacking_mode):
+    """Check that combineOrders() works on two files without stacking"""
+    adinputs = [astrodata.open(os.path.join(path_to_inputs, filename))
+                for filename in FILENAMES]
+    orig_wavl = np.array([make_wavelength_table(ad[0]) for ad in adinputs])
+    p = GHOSTSpect(adinputs)
+    adoutputs = p.combineOrders(stacking_mode=stacking_mode)
+    assert len(adoutputs) == 1  # stacking
+    ad_out = adoutputs[0]
+
+    # Check that this combined arc makes sense
+    pixels = np.arange(ad_out[0].data.size)
+    x = ad_out[0].wcs(pixels[1:]) / ad_out[0].wcs(pixels[:-1])
+    assert x.std() < 1e-12  # all wavelength ratios are the same
+    assert ad_out[0].wcs(0) == pytest.approx(orig_wavl.min())