Update cube and cutout unittests

astrocut/tests/test_cube_cut.py

@@ -243,6 +243,42 @@ def test_header_keywords_quality(cube_file, ffi_type, tmp_path):
         assert ext1_header['BJDREFI'] == ext2_header['BJDREFI'] == primary_header['BJDREFI']
 
 
+@pytest.mark.parametrize("ffi_type", ["SPOC", "TICA"])
+def test_header_keywords_diffs(cube_file, ffi_type, tmp_path):
+    """Test known header keywords differences"""
+
+    tmpdir = str(tmp_path)
+
+    cutout_maker = CutoutFactory()
+    coord = "256.88 6.38"
+    cutout_size = [5, 3]
+    out_file = cutout_maker.cube_cut(
+        cube_file, coord, cutout_size, ffi_type, output_path=path.join(tmpdir, "out_dir"), verbose=False
+    )
+
+    with fits.open(out_file) as hdulist:
+
+        # Get units from BinTableHDU
+        cols = hdulist[1].columns.info('name, unit', output=False)
+        cols_dict = dict(zip(*cols.values()))
+
+        ncols = 12 if ffi_type == 'SPOC' else 11
+        assert len(cols_dict) == ncols
+
+        # Checking for known keyword differences
+        if ffi_type == 'SPOC':
+            assert hdulist[0].header['TIMEREF'] == 'SOLARSYSTEM', 'TIMEREF keyword does not match expected'
+            assert hdulist[0].header['TASSIGN'] == 'SPACECRAFT', 'TASSIGN keyword does not match expected'
+            assert cols_dict['FLUX'] == 'e-/s', f'Expected `FLUX` units of "e-/s", got units of "{cols_dict["FLUX"]}"'
+            assert hdulist[1].data.field('CADENCENO').all() == 0.0
+
+        if ffi_type == 'TICA':
+            assert hdulist[0].header['TIMEREF'] is None, 'TIMEREF keyword does not match expected'
+            assert hdulist[0].header['TASSIGN'] is None, 'TASSIGN keyword does not match expected'
+            assert cols_dict['FLUX'] == 'e-', f'Expected `FLUX` units of "e-", got units of "{cols_dict["FLUX"]}"'
+            assert hdulist[1].data.field('CADENCENO').all() != 0.0
+
+
 @pytest.mark.parametrize("ffi_type", ["SPOC", "TICA"])
 def test_get_cutout_limits(cube_file, ffi_type, tmp_path):
     """Test _get_cutout_limits"""
@@ -390,6 +426,97 @@ def test_target_pixel_file(cube_file, ffi_type, tmp_path):
     tpf.close()
 
 
+def test_tica_cutout_error(tmp_path):
+    """
+    Test cutouts created from existing TICA cubes (i.e., those with
+    error arrays)
+
+    Write test to check that cutouts can be made from *both* the TICA
+    cubes that still have the error array, as well as cubes that
+    *have* been updated and no longer have the error array. Needed to
+    verify that CutoutFactory will work on the cubes that have not 
+    been remade yet.
+    """
+
+    tmpdir = str(tmp_path)
+
+    img_sz = 10
+    num_im = 100
+
+    # Create a mock TICA cube that still has the error array
+    # Cube with TICA headers to copy to cube with error array
+    spoc_cube_error = CubeFactory() 
+    spoc_ffi_files = create_test_ffis(img_sz, num_im, dir_name=tmpdir)
+    cube_error_file = spoc_cube_error.make_cube(spoc_ffi_files, 
+                                                path.join(tmpdir, "out_dir", "test_error_cube.fits"), 
+                                                verbose=False)
+
+    # Cube with error array
+    tica_cube_no_error = TicaCubeFactory()
+    tica_ffi_files = create_test_ffis(img_sz, num_im, dir_name=tmpdir, product='TICA')
+    tica_cube_no_error_file = tica_cube_no_error.make_cube(tica_ffi_files, 
+                                                           path.join(tmpdir, "out_dir", "test_add_error_cube.fits"), 
+                                                           verbose=False)
+
+    # Take ImageHDU from the cube with the error array, write into cube with TICA headers
+    with fits.open(tica_cube_no_error_file, mode='update') as hdulist:
+
+        tica_hdr = hdulist[1].header
+        error_hdr = fits.getheader(cube_error_file, 1)
+
+        # Overwrite header keywords
+        for h in error_hdr.keys():
+            tica_hdr[h] = error_hdr[h]
+
+        # Overwrite header data
+        hdulist[1].data = fits.getdata(cube_error_file, 1)
+
+        # Close file, write all changes
+        hdulist.close()
+
+    # Update variable name to reflect addition of error array to TICA cube
+    tica_cube_error_file = tica_cube_no_error_file
+
+    # Verify dimensions of new TICA cube with added error array
+    assert np.shape(fits.getdata(tica_cube_error_file, 1)) == (img_sz, img_sz, num_im, 2)
+
+    ffi_type = 'TICA'
+    # Test cutouts from TICA cube with error array
+    test_cube_cutout(tica_cube_error_file, tica_ffi_files, ffi_type, tmp_path)
+    test_header_keywords_quality(tica_cube_error_file, ffi_type, tmp_path)
+    test_header_keywords_diffs(tica_cube_error_file, ffi_type, tmp_path)
+    test_target_pixel_file(tica_cube_error_file, ffi_type, tmp_path)
+
+
+@pytest.mark.parametrize("ffi_type", ["SPOC", "TICA"])
+def test_ffi_cube_header(cube_file, ffi_files, ffi_type):
+    """Test FFI headers versus cube headers"""
+
+    ffi_header_keys = list(fits.getheader(ffi_files[0], 0).keys())
+    cube_header_keys = list(fits.getheader(cube_file, 0).keys())
+
+    # Lists of expected keys for each product type
+    if ffi_type == 'SPOC':
+        # CHECKSUM in FFI header, not in cube header
+        # SECTOR, DATE, ORIGIN in cube header, not in FFI header
+        effi_keys = ['CHECKSUM']
+        ecube_keys = ['SECTOR', 'DATE', 'ORIGIN']
+
+    else:
+        # CHECKSUM, and image dimensions (NAXIS1, NAXIS2) in FFI header, not in cube header
+        # SECTOR, DATE, EXTNAME, ORIGIN in cube header, not in FFI header
+        effi_keys = ['CHECKSUM', 'NAXIS1', 'NAXIS2']
+        ecube_keys = ['SECTOR', 'DATE', 'EXTNAME', 'ORIGIN']
+
+    for k in effi_keys:
+        assert k in ffi_header_keys
+        assert k not in cube_header_keys
+
+    for k in ecube_keys:
+        assert k not in ffi_header_keys
+        assert k in cube_header_keys
+
+
 @pytest.mark.parametrize('ffi_type', ['SPOC', 'TICA'])
 def test_exceptions(cube_file, ffi_type):
     """
@@ -538,6 +665,51 @@ def test_s3_cube_cut(tmp_path: Path):
     hdulist.close()
 
 
+@pytest.mark.xfail
+def test_s3_tica_cube_cut(tmp_path: Path):
+    """Does using an S3-hosted TESS cube yield correct results?
+
+    This test implements a spot check which verifies whether a cutout
+    for Sector 27 obtained from an S3-hosted cube
+    file yields results that are identical to those returned
+    by the Tesscut service.
+
+    To speed up the test and avoid adding astroquery as a dependency,
+    the test uses hard-coded reference values which were obtained
+    as follows:
+
+    >>> from astroquery.mast import Tesscut  # doctest: +SKIP
+    >>> crd = SkyCoord(299.27269, -67.14491, unit="deg")  # doctest: +SKIP
+    >>> cut = Tesscut.get_cutouts(coordinates=crd, size=3, sector=27, product='TICA')  # doctest: +SKIP
+    >>> cut[0][1].data.shape  # doctest: +SKIP
+    (3360,)
+    >>> cut[0][1].data['TIME'][0]  # doctest: +SKIP
+    2036.281565948625
+    >>> cut[0][1].data['FLUX'][500][0, 0]  # doctest: +SKIP
+    59313.523
+    >>> np.sum(cut[0][1].data['FLUX_ERR'])  # doctest: +SKIP
+    0.0
+    >>> cut[0][0].header['CAMERA']  # doctest: +SKIP
+    2
+    >>> cut[0][0].header['CCD']  # doctest: +SKIP
+    1
+
+    TODO: Remove test `test_tica_cutout_error` when this test starts working
+    """
+    # Test case: RR Lyrae star in Sector 27 (Camera 2, CCD 1)
+    coord = SkyCoord(299.27269, -67.14491, unit="deg", frame="icrs")
+    cube_file = "s3://stpubdata/tess/public/mast/tica/tica-s0027-2-1-cube.fits"
+    cutout_file = CutoutFactory().cube_cut(cube_file, coord, 3, output_path=str(tmp_path))
+    hdulist = fits.open(cutout_file)
+    assert hdulist[1].data.shape == (3360,)
+    assert np.isclose(hdulist[1].data["TIME"][0], 2036.281565948625)
+    assert np.isclose(hdulist[1].data["FLUX"][500][0, 0], 59313.523)
+    assert np.sum(hdulist[1].data["FLUX_ERR"]) == 0.0
+    assert hdulist[0].header["CAMERA"] == 2
+    assert hdulist[0].header["CCD"] == 1
+    hdulist.close()
+
+
 def test_multithreading(tmp_path):
     tmpdir = str(tmp_path)