From bf630d88fb960de3b7fd2f67158482d7b8f6054d Mon Sep 17 00:00:00 2001 From: chris-simpson Date: Thu, 1 Feb 2024 15:43:05 -1000 Subject: [PATCH] fix WCS issues and write combineOrders() test for GHOST --- astrodata/wcs.py | 26 ++++---- .../ghost/tests/spect/test_combine_orders.py | 66 +++++++++++++++++++ 2 files changed, 79 insertions(+), 13 deletions(-) create mode 100644 geminidr/ghost/tests/spect/test_combine_orders.py diff --git a/astrodata/wcs.py b/astrodata/wcs.py index a72f35c10..a94e01d00 100644 --- a/astrodata/wcs.py +++ b/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: diff --git a/geminidr/ghost/tests/spect/test_combine_orders.py b/geminidr/ghost/tests/spect/test_combine_orders.py new file mode 100644 index 000000000..488519e8b --- /dev/null +++ b/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())