`io.fits` does not properly handle VLA columns that have entries with different dimensions
- Dominant language
- Python
- Stars
- 5.3k
- Forks
- 2.2k
- Avg merge
- 1d 18h
- Merged PRs (30d)
- 82
Description
### Description
While the FITS Standard does not support variable-length array (VLA) columns that contain data of different shapes, `io.fits` does not raise warnings/errors when trying to do so. In some scenarios, this means `io.fits` can create FITS files that cannot be fully read back.
This is a follow-up to the great work done by @cmarmo on #13417, which was a big step forward in fully supporting VLAs in `io.fits`.
### Expected behavior
I described the expected behavior of VLAs in much greater detail in https://github.com/astropy/astropy/issues/12860#issuecomment-1074628472, but in short: 1) when trying to create/write a VLA column with data of mismatched shapes, a warning/exception should be thrown and 2) when reading a file that contains impossible dimension constraints, TDIM should be ignored and a warning should be thrown. See reproduction steps for more detail.
### How to Reproduce
```python
from astropy.io import fits
import numpy as np
cols = []
# Let's go with dim='(5, 4)' and always use the correct shape for the first entry in the table
# Keep in mind that TDIM is reversed from the NumPy shape
# Here the second entry has the same size but different shape
array = [np.arange(20).reshape((4, 5)), np.arange(20).reshape((10, 2))]
cols.append(fits.Column(name="col1", format="PI", array=array, dim='(5, 4)'))
# Here the second entry has more elements than what TDIM implies
array = [np.arange(20).reshape((4, 5)), np.arange(30).reshape((10, 3))]
cols.append(fits.Column(name="col2", format="PI", array=array, dim='(5, 4)'))
# Here the second entry has more elements and is three dimensional
array = [np.arange(20).reshape((4, 5)), np.arange(30).reshape((5, 3, 2))]
cols.append(fits.Column(name="col3", format="PI", array=array, dim='(5, 4)'))
# Here the second entry has less elements than what is necessary for the defined TDIM
array = [np.arange(20).reshape((4, 5)), np.arange(4).reshape((2, 2))]
cols.append(fits.Column(name="col4", format="PI", array=array, dim='(5, 4)'))
# Here the second entry is empty
array = [np.arange(20).reshape((4, 5)), np.array([])]
cols.append(fits.Column(name="col5", format="PI", array=array, dim='(5, 4)'))
fits.BinTableHDU.from_columns(name='tab', columns=cols).writeto('test.fits', overwrite=True)
tab = fits.getdata('test.fits')
```
Although we broke a lot of the rules set in the FITS standard, these lines of code actually run without raising any warning or exception. I'll break it down column by column:
___
**`"col1"`**
In my opinion, on writing `col1` Astropy should've warned that the dimensions of the arrays don't match. On reading it back we get:
```python
>>> tab['col1'][1].shape
(4, 5)
```
Although one might've expected `(10, 2)`, given that FITS cannot set different dimensions per row this result makes sense as we're using the dimension set on `TDIM`.
___
**`"col2"` and `"col3"`**
The standard says "If the number of elements in the array implied by the `TDIMn` is fewer than the allocated size of the array in the FITS file, then the unused trailing elements should be interpreted as containing undefined fill values." As I understand it, this essentially means that if the arrays have more elements than expected, then they should be ignored. So we would expect to still get data in the (4, 5) shape. Yet we get:
```python
>>> tab['col2'][1].shape
(6, 5)
>>> tab['col3'][1].shape
(6, 5)
```
For some reason `io.fits` is forcing the first dimension to match `TDIM` and adjusting whatever the other dimension is in order to fit the data. I don't think this should be the default behavior, at least not without any warning.
___
**`"col4"`**
The Standard says "The size [implied by `TDIM`] must be (...), in the case of columns that have a `’P’` or `’Q’` `TFORMn` data type, less than or equal to the array length specified in the variable-length array descriptor". This means it should never be possible to have a row in a VLA column that has less elements than what the `TDIM` shape requires. However, `io.fits` allows us to write this. Then if we try to read it back:
```python
>>> tab['col4'][1].shape
Traceback (most recent call last):
File "[...]\lib\site-packages\astropy\io\fits\fitsrec.py", line 513, in __getitem__
return self.field(key)
File "[...]\lib\site-packages\astropy\io\fits\fitsrec.py", line 723, in field
converted = self._convert_p(column, field, recformat)
File "[...]\lib\site-packages\astropy\io\fits\fitsrec.py", line 831, in _convert_p
dummy[idx] = dummy[idx].reshape((vla_first,) + vla_dim)
ValueError: cannot reshape array of size 4 into shape (0,5)
```
Very cryptic error message. Following, the "tolerant with input and strict with output" philosophy of Astropy, we should clearly not allow such a column to be written to a file as it would likely break most readers, However, when reading such a file we could simply throw a warning, ignore TDIM and assume the array is one-dimensional.
___
**`"col5"`**
The Standard says "In the special case where the variable-length array descriptor has a size of zero, then the `TDIMn` keyword is not applicable". This means that it's entirely legal to have an empty array in the column, and that array should not interact with `TDIM` in any way. We should just get back an empty array of shape `(0,)`, instead we get:
```python
>>>tab['col5'][1]
array([], shape=(0, 5), dtype=int16)
```
It is able to read back an empty array, which is great, as that is what is expected. However, the shape itself is wrong, it seems that once again `io.fits` is really forcing the first dimension to match `TDIM`.
### Versions
Windows-10-10.0.19044-SP0
Python 3.10.10 (tags/v3.10.10:aad5f6a, Feb 7 2023, 17:20:36) [MSC v.1929 64 bit (AMD64)]
astropy 5.3.1
Numpy 1.24.3
pyerfa 2.0.0.1
Scipy 1.10.0
Matplotlib 3.6.3
Contributor guide
Research direction
Start with astropy.io.fits, especially fits.Column, BinTableHDU.from_columns, and astropy/io/fits/fitsrec.py's _convert_p path; run the reproduction to observe each VLA case. Done means mismatched shapes produce a warning or exception on write, invalid TDIM constraints warn and fall back to one-dimensional data on read, and empty arrays retain shape (0,).
Written by the indexing model from the issue text.
Assessment
- Tech stack
- numpy, python
- Domain
- data
- Issue type
- Bug
- Difficulty
- 4/5
- Estimated time
- 3-5 days
- Activity status
- Stale
- Clarity
- Mostly clear
- Newbie friendliness
- 42/100