Repository: astropy/astropy

Issue:
Slicing FITS_rec objects produces unexpected dtypes (VLA columns)
Original conversation from https://github.com/astropy/astropy/issues/5490#issuecomment-2336789274 onwards.

> @peppedilillo , can you please give us a minimally reproducible code snippet that _you_ encountered? What numpy version? Python version, OS, distro, etc?

```python
from astropy.io import fits


if __name__ == "__main__":
    hdul = fits.open("mre.fits")
    fields = ["CHANNEL", "PHA", "CID", "SDDID"]
    for field in fields:
        print(f"slice, then access {field} field: ", hdul[2].data[:1][field])
    for field in fields:
        print(f"access field, then access {field} slice: ", hdul[2].data[field][:1])
    hdul.close()
```

This outputs:

```python
slice, then access CHANNEL field:  [array([  0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   2,   2,
          0,   0,   3, 228,   0,   0,   0,   0,   0,  76, 106, 162,  65,
        140,  52,  98,  12,   1, 165,  55,   0,  64, 215,  77,  89, 166,
        110,   0,   0,   0,   2,   0,   0,   0,   0,   0,   0,   0,   2,
          0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
          0,   0,   2,   0,   0,   0,   0,   0,  12,  61, 246,   0,   0,
          0,   0,   0,   0,   0,   2,   0,   0,   0,   0,   0,  36, 185,
        226,   0,   0,   0,   0,   0,   0,   0,   2,   0,   0,   0,   0,
          0,  48, 247, 216,   1,   1,  68,  37, 228, 119,   0,   0,   0,
        203], dtype=uint8)                                              ]
slice, then access PHA field:  [array([], dtype=int16)]
slice, then access CID field:  [array([], dtype=uint8)]
slice, then access SDDID field:  [array([], dtype=uint8)]

access CHANNEL field, then slice:  [array([13, 19], dtype=uint8)]
access PHA field, then slice:  [array([21786, 22649], dtype=int16)]
access CID field, then slice:  [array([7, 7], dtype=uint8)]
access SDDID field, then slice:  [array([1, 0], dtype=uint8)]
```


The first four prints are garbage. The last four are correct, and take a long time in each statement. You can download the input file and the script [here](https://drive.google.com/file/d/1x6vTx8cfwkklaiPAerTv2-2MYZS01T5W/view?usp=sharing). Note that the problem occurs also if `fields` is set to `[PHA]`, which is signed. Some of these columns are variable length.

This occurred to me because I am working on a FITS viewer application for the terminal. The application show FITS tables in short "pages". Each page consists of a slice of `FITS_rec`. I could not get the right output for certain input files because of this problem. I'm trying to work on slices because, if some scaling or other slow astropy operation are carried out over the whole table, that would result in annoying delays on the readout.

Just tested on Debian 12 bookworm, but observed the same behaviour also on Windows 11 and MacOS Ventura.
Yes i'm on numpy 2 and Python 3.12.5. The whole environment is the following, I'm on anaconda presently:

```
# packages in environment at /home/deppep/miniforge3/envs/misfits:
#
# Name                    Version                   Build  Channel
_libgcc_mutex             0.1                 conda_forge    conda-forge
_openmp_mutex             4.5                       2_gnu    conda-forge
astropy                   6.1.3                    pypi_0    pypi
astropy-iers-data         0.2024.9.2.0.33.23          pypi_0    pypi
bzip2                     1.0.8                h4bc722e_7    conda-forge
ca-certificates           2024.8.30            hbcca054_0    conda-forge
click                     8.1.7                    pypi_0    pypi
ld_impl_linux-64          2.40                 hf3520f5_7    conda-forge
libexpat                  2.6.2                h59595ed_0    conda-forge
libffi                    3.4.2                h7f98852_5    conda-forge
libgcc                    14.1.0               h77fa898_1    conda-forge
libgcc-ng                 14.1.0               h69a702a_1    conda-forge
libgomp                   14.1.0               h77fa898_1    conda-forge
libnsl                    2.0.1                hd590300_0    conda-forge
libsqlite                 3.46.0               hde9e2c9_0    conda-forge
libuuid                   2.38.1               h0b41bf4_0    conda-forge
libxcrypt                 4.4.36               hd590300_1    conda-forge
libzlib                   1.3.1                h4ab18f5_1    conda-forge
linkify-it-py             2.0.3                    pypi_0    pypi
markdown-it-py            3.0.0                    pypi_0    pypi
mdit-py-plugins           0.4.1                    pypi_0    pypi
mdurl                     0.1.2                    pypi_0    pypi
misfits                   0.0.1                    pypi_0    pypi
ncurses                   6.5                  he02047a_1    conda-forge
numpy                     2.1.0                    pypi_0    pypi
openssl                   3.3.1                hb9d3cd8_3    conda-forge
packaging                 24.1                     pypi_0    pypi
pandas                    2.2.2                    pypi_0    pypi
pip                       24.2               pyh8b19718_1    conda-forge
platformdirs              4.2.2                    pypi_0    pypi
pyerfa                    2.0.1.4                  pypi_0    pypi
pygments                  2.18.0                   pypi_0    pypi
python                    3.12.5          h2ad013b_0_cpython    conda-forge
python-dateutil           2.9.0.post0              pypi_0    pypi
pytz                      2024.1                   pypi_0    pypi
pyyaml                    6.0.2                    pypi_0    pypi
readline                  8.2                  h8228510_1    conda-forge
rich                      13.8.0                   pypi_0    pypi
setuptools                72.2.0             pyhd8ed1ab_0    conda-forge
six                       1.16.0                   pypi_0    pypi
textual                   0.79.1                   pypi_0    pypi
tk                        8.6.13          noxft_h4845f30_101    conda-forge
typing-extensions         4.12.2                   pypi_0    pypi

```

<hr/>

@saimn : @peppedilillo - your issue is a different one, it has to do with VLA columns. Interestingly if you first access the full columns and then slice it, it works.

Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/development_details.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/development_details.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/development_details.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
saimn: Seems to create problems with other tests, was a bit optimistic when running only VLA tests :grimacing: .

Patch:
diff --git a/astropy/io/fits/fitsrec.py b/astropy/io/fits/fitsrec.py
index 506579ecb7..01d4161433 100644
--- a/astropy/io/fits/fitsrec.py
+++ b/astropy/io/fits/fitsrec.py
@@ -208,6 +208,7 @@ class FITS_rec(np.recarray):
             "_converted",
             "_heapoffset",
             "_heapsize",
+            "_tbsize",
             "_nfields",
             "_gap",
             "_uint",
@@ -238,6 +239,7 @@ class FITS_rec(np.recarray):
             self._converted = obj._converted
             self._heapoffset = obj._heapoffset
             self._heapsize = obj._heapsize
+            self._tbsize = obj._tbsize
             self._col_weakrefs = obj._col_weakrefs
             self._coldefs = obj._coldefs
             self._nfields = obj._nfields
@@ -251,6 +253,7 @@ class FITS_rec(np.recarray):
 
             self._heapoffset = getattr(obj, "_heapoffset", 0)
             self._heapsize = getattr(obj, "_heapsize", 0)
+            self._tbsize = getattr(obj, "_tbsize", 0)
 
             self._gap = getattr(obj, "_gap", 0)
             self._uint = getattr(obj, "_uint", False)
@@ -273,6 +276,7 @@ class FITS_rec(np.recarray):
         self._converted = {}
         self._heapoffset = 0
         self._heapsize = 0
+        self._tbsize = 0
         self._col_weakrefs = weakref.WeakSet()
         self._coldefs = None
         self._gap = 0
@@ -517,6 +521,10 @@ class FITS_rec(np.recarray):
 
         # We got a view; change it back to our class, and add stuff
         out = out.view(type(self))
+        out._heapoffset = self._heapoffset
+        out._heapsize = self._heapsize
+        out._tbsize = self._tbsize
+        out._gap = self._gap
         out._uint = self._uint
         out._coldefs = ColDefs(self._coldefs)
         arrays = []
@@ -1032,7 +1040,7 @@ class FITS_rec(np.recarray):
         May return ``None`` if no array resembling the "raw data" according to
         the stated criteria can be found.
         """
-        raw_data_bytes = self.nbytes + self._heapsize
+        raw_data_bytes = self._tbsize + self._heapsize
         base = self
         while hasattr(base, "base") and base.base is not None:
             base = base.base
diff --git a/astropy/io/fits/hdu/table.py b/astropy/io/fits/hdu/table.py
index ece25e720a..f5365adb76 100644
--- a/astropy/io/fits/hdu/table.py
+++ b/astropy/io/fits/hdu/table.py
@@ -216,8 +216,8 @@ class _TableLikeHDU(_ValidHDU):
         # pass datLoc, for P format
         data._heapoffset = self._theap
         data._heapsize = self._header["PCOUNT"]
-        tbsize = self._header["NAXIS1"] * self._header["NAXIS2"]
-        data._gap = self._theap - tbsize
+        data._tbsize = self._header["NAXIS1"] * self._header["NAXIS2"]
+        data._gap = self._theap - data._tbsize
 
         # pass the attributes
         for idx, col in enumerate(columns):
@@ -515,6 +515,7 @@ class _TableBaseHDU(ExtensionHDU, _TableLikeHDU):
             # calculate PCOUNT, for variable length tables
             tbsize = self._header["NAXIS1"] * self._header["NAXIS2"]
             heapstart = self._header.get("THEAP", tbsize)
+            self.data._tbsize = tbsize
             self.data._gap = heapstart - tbsize
             pcount = self.data._heapsize + self.data._gap
             if pcount > 0:
diff --git a/docs/changes/io.fits/16996.bugfix.rst b/docs/changes/io.fits/16996.bugfix.rst
new file mode 100644
index 0000000000..f1e27100fb
--- /dev/null
+++ b/docs/changes/io.fits/16996.bugfix.rst
@@ -0,0 +1,1 @@
+Fix access to VLA columns after slicing ``.data``.


Test patch:
diff --git a/astropy/io/fits/hdu/compressed/tests/test_compressed.py b/astropy/io/fits/hdu/compressed/tests/test_compressed.py
index f4c8e8ef7e..3e7d1ec60a 100644
--- a/astropy/io/fits/hdu/compressed/tests/test_compressed.py
+++ b/astropy/io/fits/hdu/compressed/tests/test_compressed.py
@@ -1349,84 +1349,3 @@ def test_hdu_lazy_loading(tmp_path):
                 assert_equal(hdulist[idx + 1].data, idx)
                 assert hdulist[idx + 1].data.shape == (2**idx, 2**idx)
                 fileobj.seek(0)
-
-
-@pytest.mark.parametrize(
-    "init_kwargs, expected_zflags",
-    [
-        pytest.param(
-            dict(quantize_level=-32, data=np.array([0.0])),
-            {
-                "ZNAME1": "BLOCKSIZE",
-                "ZVAL1": 32,
-                "ZNAME2": "BYTEPIX",
-                "ZVAL2": 4,
-                "ZNAME3": "NOISEBIT",
-                "ZVAL3": -32,
-                "ZCMPTYPE": "RICE_1",
-                "ZQUANTIZ": "NO_DITHER",
-            },
-            id="quantize_level_w_data",
-        ),
-        pytest.param(
-            dict(quantize_level=-32),
-            {
-                "ZNAME1": "BLOCKSIZE",
-                "ZVAL1": 32,
-                "ZNAME2": "BYTEPIX",
-                "ZVAL2": 4,
-                "ZNAME3": "NOISEBIT",
-                "ZVAL3": -32,
-                "ZCMPTYPE": "RICE_1",
-                "ZQUANTIZ": "NO_DITHER",
-            },
-            id="quantize_level_wo_data",
-        ),
-        pytest.param(
-            dict(quantize_method=2, dither_seed=1, data=np.array([0.0])),
-            {
-                "ZNAME1": "BLOCKSIZE",
-                "ZVAL1": 32,
-                "ZNAME2": "BYTEPIX",
-                "ZVAL2": 4,
-                "ZNAME3": "NOISEBIT",
-                "ZVAL3": 16,
-                "ZCMPTYPE": "RICE_1",
-                "ZQUANTIZ": "SUBTRACTIVE_DITHER_2",
-            },
-            id="quantize_level_dither_seed_w_data",
-        ),
-        pytest.param(
-            dict(quantize_method=2, dither_seed=1),
-            {
-                "ZNAME1": "BLOCKSIZE",
-                "ZVAL1": 32,
-                "ZNAME2": "BYTEPIX",
-                "ZVAL2": 4,
-                "ZNAME3": "NOISEBIT",
-                "ZVAL3": 16,
-                "ZCMPTYPE": "RICE_1",
-                "ZQUANTIZ": "SUBTRACTIVE_DITHER_2",
-            },
-            id="quantize_level_dither_seed_wo_data",
-        ),
-    ],
-)
-def test_compression_options_with_mutated_data(
-    init_kwargs,
-    expected_zflags,
-    tmp_path,
-):
-    # see https://github.com/astropy/astropy/issues/12216
-    prng = np.random.default_rng(seed=0)
-    data = prng.random((50, 30)).astype("float32")
-
-    hdu = fits.CompImageHDU(**init_kwargs)
-    hdu.data = data
-
-    hdu.writeto(tmp_path / "t.fits")
-
-    hdr = fits.getheader(tmp_path / "t.fits", ext=1, disable_image_compression=True)
-
-    zflags = {k: hdr[k] for k in expected_zflags}
-    assert zflags == expected_zflags
diff --git a/astropy/io/fits/tests/test_table.py b/astropy/io/fits/tests/test_table.py
index a72eb19db4..a377123843 100644
--- a/astropy/io/fits/tests/test_table.py
+++ b/astropy/io/fits/tests/test_table.py
@@ -3349,6 +3349,40 @@ class TestVLATables(FitsTestCase):
         ):
             fits.BinTableHDU.from_columns([c1])
 
+    def test_vla_slice(self):
+        filename = self.data("variable_length_table.fits")
+
+        # slicing the VLA array
+        with fits.open(filename) as hdul:
+            hdu = hdul[1]
+            arr = hdu.data["var"][:1]
+            assert arr.tolist() == [[45, 56]]
+            assert arr[0].dtype == np.int16
+
+        # slicing .data
+        with fits.open(filename) as hdul:
+            hdu = hdul[1]
+            arr = hdu.data[:1]["var"]
+            assert arr.tolist() == [[45, 56]]
+            assert arr[0].dtype == np.int16
+
+    def test_vla_slice2(self):
+        filename = self.data("theap-gap.fits")
+
+        # slicing the VLA array
+        with fits.open(filename) as hdul:
+            hdu = hdul[1]
+            arr = hdu.data["arr"][:3]
+            assert arr.tolist() == [[], [0], [0, 1]]
+            assert arr[1].dtype == np.int32
+
+        # slicing .data
+        with fits.open(filename) as hdul:
+            hdu = hdul[1]
+            arr = hdu.data[:3]["arr"]
+            assert arr.tolist() == [[], [0], [0, 1]]
+            assert arr[1].dtype == np.int32
+
 
 # These are tests that solely test the Column and ColDefs interfaces and
 # related functionality without directly involving full tables; currently there