Broken FITS table round-trip with any nan in a float column.

Description

This is a follow-up to the long-standing issue #4708 (which is still open). I'm filing this as a "new" issue because the previous report does not cover all the cases, but I understand if this is somewhat of a duplicate.

The basic problem is that any float column that is serialized to fits that has a nan value gets transformed on read into a MaskedColumn that behaves in strange and inconsistent ways.

from astropy.table import Table
import numpy as np
import tempfile

# Round-trip a table correctly (no nans)
tab1 = Table(data=np.zeros(10, dtype=[("a", "f8")]))
f, name = tempfile.mkstemp(suffix=".fits")
tab1.write(name, overwrite=True)

tab2 = Table.read(name)
print(tab2["a"].__class__)
# <class 'astropy.table.column.Column'>

tab1_nan = tab1.copy()
tab1_nan["a"][0] = np.nan
f_nan, name_nan = tempfile.mkstemp(suffix=".fits")
tab1_nan.write(name_nan, overwrite=True)

tab2_nan_broken = Table.read(name_nan)
print(tab1_nan["a"].__class__)
# <class 'astropy.table.column.Column'>

print(tab2_nan_broken["a"].__class__)
# <class 'astropy.table.column.MaskedColumn'>

print(np.array(tab2_nan_broken["a"]))
# [nan  0.  0.  0.  0.  0.  0.  0.  0.  0.]

print(tab2_nan_broken["a"].tolist())
# [None, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

print(np.array(tab2_nan_broken["a"]).tolist())
# [nan, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

print(tab2_nan_broken["a"].filled())
#   a
# -----
# 1e+20
#   0.0
#   0.0
#   0.0
#   0.0
#   0.0
#   0.0
#   0.0
#   0.0
#   0.0

There are a number of surprising (and I believe incorrect) behaviors here.

First of all, the type of column on read is different than that used on write. Second of all, the __array__ access of the masked column seems to return the correct filled array, but the .filled() option doesn't (the masked array has the wrong fill value, which makes me wonder how __array__() even knows what the fill value should be). Third, converting to a list with tolist() gives inconsistent and different answers depending on the order of operations. Finally (not shown), using isnan() gives different results depending on the table is used before or after serialization because of the column type change. This led to a very subtle but critical index error in our codebase.

Fortunately, reading the table with mask_invalid=False does fix the issue (or, interestingly, with memmap=True, which oddly gives a different table format than without). However, it is not apparent from the documentation that this should be the solution to the problem. First, it is buried in the unified interface section https://docs.astropy.org/en/stable/io/unified.html#masked-columns and not about fits reading. Second, this section is talking about masked columns, but at no point did I intentionally create a masked column.

I would suggest that mask_invalid=True should be the default. I understand this is an API change, but the current behavior is fundamentally broken. If one wants to write masked columns, there is the option with mixins to store the mask information explicitly. At the very least, the tolist() and filled() methods of masked columns should be fixed.

Expected behavior

No response

How to Reproduce

Get package from '...'
Then run '...'
An error occurs.

# Put your Python code snippet here.

Versions

macOS-12.6.3-arm64-arm-64bit
Python 3.10.8 | packaged by conda-forge | (main, Nov 22 2022, 08:25:29) [Clang 14.0.6 ]
astropy 5.2
Numpy 1.21.6
pyerfa 2.0.0.1
Scipy 1.9.3
Matplotlib 3.6.2

First of all, the type of column on read is different than that used on write.

Yes, by default NaNs and null strings are masked, since #11222 (v4.3). The mask_invalid option was added later (#12544), so the masking by default was kept for compatibility. IMO masking should not be done by default, but changing that would be a breaking change.

Second of all, the array access of the masked column seems to return the correct filled array, but the .filled() option doesn't (the masked array has the wrong fill value, which makes me wonder how array() even knows what the fill value should be).

__array__ doesn't know the fill value, it just ignores the mask. It's the behavior from Numpy:

In [9]: np.ma.masked_invalid([np.nan, 0, 0])
Out[9]: 
masked_array(data=[--, 0.0, 0.0],
             mask=[ True, False, False],
       fill_value=1e+20)

In [10]: np.ma.masked_invalid([np.nan, 0, 0]).filled()
Out[10]: array([1.e+20, 0.e+00, 0.e+00])

In [11]: np.array(np.ma.masked_invalid([np.nan, 0, 0]))
Out[11]: array([nan,  0.,  0.])

Third, converting to a list with tolist() gives inconsistent and different answers depending on the order of operations.

Related to 2:

In [12]: np.ma.masked_invalid([np.nan, 0, 0]).tolist()
Out[12]: [None, 0.0, 0.0]

Fortunately, reading the table with mask_invalid=False does fix the issue (or, interestingly, with memmap=True, which oddly gives a different table format than without).

See discussion in #12544 for memmap, which is only useful when masking is not done so mask_invalid is set to False when memmap=True. This is explained in the docstring.

However, it is not apparent from the documentation that this should be the solution to the problem. First, it is buried in the unified interface section https://docs.astropy.org/en/stable/io/unified.html#masked-columns and not about fits reading.

Documentation can probably be improved but can you explain where you would expect to find this. Currently
Maybe https://docs.astropy.org/en/stable/io/unified.html#fits is linked in https://docs.astropy.org/en/stable/io/fits/index.html#working-with-table-data.

Thanks for your response.

That's very interesting that the swap between nan and 1e20 and None is a numpy mis-feature; very frustrating. Perhaps astropy can (partly) paper over this by setting the fill value to nan?

But my basic complaint is that I think that the default is incorrect (as you have suggested). I understand that it would be a breaking change, but right now reading a table changes data types from what was written, which is (imo) incorrect behavior.

I'm not sure where in the docs it should go, but what really threw me for a loop is that the discussion of this feature is in the section about masked columns -- but at no point did I ever create a masked column (and in fact they were not stored as masked columns on disk). Related, I also think it's a problem that you get a different data type if you read a table with or without memmap=True!

Hi,
After getting bitten by this issue several times,
I'll post here an example of code that shows how dangerous the current default behaviour is.
Simply saving a column with nan and then doing a histogram of that column after saving will return different result with current astropy/numpy.

import astropy.table as atpy
import numpy as np
import matplotlib.pyplot as plt

# create the array
x = np.array([1, np.nan])
hh, loc = np.histogram(np.log10(x), range=[-1, 5], bins=2)
print('Before', hh)

# writing
T = atpy.Table({'x': x})
T.write('/tmp/xx.fits', overwrite=True)

# reading
T = atpy.Table().read('/tmp/xx.fits')
hh, loc = np.histogram(np.log10(T['x']), range=[-1, 5], bins=2)
print('After', hh)

Output:

Before [1 0]
After [2 0]

Here we see interplay of

astropy converting things to masked array
then numpy ignoring/replacing fill_value when doing log10,
then np.histogram ignoring mask altogether.
The danger of this is this all happens silently and just leads to incorrect results.

In my opinion the right short term fix here should be make mask_invalid=False the default, as suggested above, or put a big warning.

Previously in a lot of my codes I had to do .filled(np.nan), otherwise I'd get weird results (because of this issue & I didn't know the mask_invalid option).

Thanks

PS everything here was tested with python 3.10.12, numpy 2.1.3, astropy 6.1.7

@saimn

Can we solve the issue partially by, when including the state of the column in the fits header when writing the file (as part of the astropy metadata), so that at least we roundtrip correctly? As @saimn noted, I don't think we can easily change the default for the reader, because that will affect any file that has nan in them (not sure whether the FITS standard has anything to say about how that should be interpreted; for integer columns, one has the ability to give the NULL value). But we can write our intent to the FITS file.

I think from user's perspective that certainly would be better.
From what I can see simple tables like I had in my example don't produce any astropy serialization metadata at all in the header. I don't remember when it is produced, but I would rather have that extra stuff in the header than having the current behaviour.
TBH I am not sure I fully follow the reasoning why it's impossible to change current behaviour wrt reading.
IMO I understand why columns with TNULL* would produce a masked array (i.e. integers, strings). The columns without TNULL should not (and apparently floating point TNULLs aren't allowed, which I interpret that the standard does not allow any kind of a mask for those at all).

@segasai

@segasai - I think it is more that for integers and strings, there is no natural way to indicate a masked value, and hence TNULL was introduced. For float, np.nan is a natural value for a mask, which in astropy is represented as such. As @saimn said, this is not something we can change now - this has been around for too long and likely lots of code relies on it. Furthermore, there is an option to get what you want -- even if it could use being more obvious...

That said, we can and should fix some things. The most obvious are to document this in a place where it would be easy to find (suggestions/PRs welcome!) and to set the fill value to np.nan. Next up would be to see whether we can serialize the column type. For this, let me also ping @taldcroft as it is about the behaviour of Table too. It might be possible to have this part available without having to pass on an option on the writing side (i.e., write out the column as serialized if it contains any nan).

Definitely documentation would be a great start. Adding some red warning boxes in a the key entry points for reading FITS tables would help.

About adding serialization information to FITS files written by astropy seems a little tricky. The default for astropy_native when reading is False, so then we would be adding serialization information that always gets read, regardless of astropy_native? In general this idea may be adding complexity under the hood (an astropy-specific FITS schema) that will make some users happy but end up surprising others.

Another consideration is that for users that write tables with astropy and then read them back, this would be a breaking change. Where previously they expected NaN to be converted to a masked element, now it comes back as a Column with NaNs. Again we may be boxed in by needing to be extremely conservative about any breaking changes in astropy FITS handling.

It's interesting that the FITS standard does not allow a NULL value for float columns. That pretty much forces astropy to guess about the intent of NaN in the column. Either way, someone will be surprised. Having mask_invalid=False be the default would mean a MaskedColumn would not roundtrip by default.

@taldcroft

FWIW, I've tried to repeat my example code from above for different formats and see which ones works

Success:
parquet
hdf5
ascii.ecsv
ascii

Fail:
pandas.csv
csv
votable
fits

To my surprise votable failed here as well. I think for internal use I'll probably switch to parquet as an alternative.

As a reply to @taldcroft . I understand the desire to avoid breaking changes. But to be honest the example above, is a case of clearly incorrect results ( partially caused by numpy's issues TBH). IMO that justifies breaking (maybe over the period of few releases with a big warning of the impending change)

@taldcroft

@taldcroft - I think I had misremembered what astropy_native actually does and hence made a misleading example. Now that I looked again, what it does is just converting FITS time-standard columns to Time directly. Which is not really relevant here, sorry!

What I tried to suggest is that we could use serialization to get things right. Currently, if I write a table with a SkyCoord to fits, it reads it back in correctly. So, it would be relatively trivial to set up writing with an option to serialize in the same way all columns, so that all column classes are explicitly stored in the file and that on reading even a Column with nan in it gets converted back to a Column. The question then would be on how to implement that, with the simplest being to always do that if the column (or quantity) contains NaN, but possibly requiring someone passing in an option to request that check. I think it would be nice to have at least the option for any astropy table to round-trip perfectly through FITS, which I think should be possible.

@segasai - there is really no way we can change reading of standard FITS files - a NaN in a float column will have to continue to become a MaskedColumn - the best we can do on our own is to ensure that if we write a regular Column which happens to have NaN values to FITS, then that is preserved (at least for astropy; for any other software, I fear we'd need a change in the FITS standard).

Anyway, as I wrote above, there is the other intermediate step which will at least help, which is to set fill value to nan. I find it rather unlikely that people are counting and/or happy with 1e20... I think a PR for that (or for documentation) would be welcome...

Not saying this is pretty, but as a proof of concept, with the following diff:

diff --git a/astropy/table/column.py b/astropy/table/column.py
index 1cf4a9b1a4..4910d99fed 100644
--- a/astropy/table/column.py
+++ b/astropy/table/column.py
@@ -373,12 +373,16 @@ class ColumnInfo(BaseColumnInfo):
         names = self._parent.dtype.names
         # For a regular column, we are done, but for a structured
         # column, we use a SerializedColumns to store the pieces.
-        if names is None:
+        if names is None and (self._parent.meta.get("__done__", False) or not any(np.isnan(self._parent))):
             return result
 
         from .serialize import SerializedColumn
 
         data = SerializedColumn()
+        if names is None:  # nan elements
+            result["data"] = self._parent.__class__(self._parent)  # copy meta
+            result["data"].meta["__done__"] = True
+            return result
         # If this column has a StructuredUnit, we split it and store
         # it on the corresponding part. Otherwise, we just store it
         # as an attribute below.  All other attributes we remove from
@@ -430,6 +434,10 @@ class ColumnInfo(BaseColumnInfo):
 
     def _construct_from_dict(self, map):
         if not isinstance(map.get("data"), dict):
+            if map.get("meta", {}).pop("__done__", False):
+                assert isinstance(map["data"], MaskedColumn)
+                map["data"] = map["data"].view(Column)
+
             return super()._construct_from_dict(map)

I get

import numpy as np
from astropy.table import Table
t = Table({"c": [1.0, np.nan]})
t.write("x.fits", overwrite=True)
Table.read("x.fits")
<Table length=2>
   c   
float64
-------
    1.0
    nan

Separately, if one uses Quantity and QTable this problem does not exist -- so if you have columns that have meaningful units (can be dimensionless), there is a work-around:

import numpy as np
import astropy.units as u
from astropy.table import Table, QTable
t = QTable({"c": [1.0, np.nan]*u.m})
t.write("x.fits", overwrite=True)
QTable.read("x.fits")
<QTable length=2>
   c   
   m   
float64
-------
    1.0
    nan

@mvhk

@mvhk - (I also misremembered astropy_native, sigh).

On the bigger point, your change will still potentially break workflows with a FITS file written by astropy (e.g. some pipeline product) and then later read using astropy. Users may have code that looks for a MaskedColumn and acts accordingly (instead of looking for NaNs).

Instead, if a broad consensus of the astropy community thinks that this issue is a big problem, we should change the default for mask_invalid. To my mind that could rise to the level of an APE. For instance in astropy 8.0, we could start emitting a warning if NaN is found when reading a FITS file where the code does not specify mask_invalid as False, True, or None, where the new default would be .... The warning would say to provide some value for mask_invalid and that the default will change in 9.0. Then in astropy 9.0 make the change.

That is a fair bit of churn, so the APE would need to weigh whether that is worth it relative to a documentation fix and accepting some risk of mistakes if people don't read the docs.

But to be honest the example above, is a case of clearly incorrect results ( partially caused by numpy's issues TBH).

For the record, I would consider those incorrect results to be entirely upstream issues of numpy treating masked arrays either incorrectly or inconsistently. Can those problems be fixed in numpy?

I would not recommend changing astropy FITS handling solely on the basis of those numpy issues. But I do have the feeling that the original choice for mask_invalid=True might merit reconsideration.

Just for reference, here are the relevant numpy issues, I believe
numpy/numpy#3762
numpy/numpy#10019
I got the impression that numpy masked array issues are non-trivial to fix either.

From my perspective (back to OP), I do think it is a mistake (surprising behavior) to automagically convert a floating point numpy array column to a masked column on read. I assume (with no statistics to back this up!) that the use of nans in floating point columns is much more common than the intentional use of masked arrays (partly because numpy masked arrays have surprising and slightly broken behaviors in some cases, so I try to avoid them unless necessary). If any additional metadata is to be stored in a fits file then I think it would be astropy saying "this was serialized as a masked array" and not the inverse. That's why I would prefer mask_invalid=False to be the default (understanding this is a breaking change). At the very least the documentation on the tables needs to be made more clear. As I posted originally the documentation is in the section on masked arrays and I had not thought of looking there because at no point did I explicitly create a masked array.

In the meantime, however, I think it would be helpful if the current behavior was improved in a couple of ways (that have already been suggested above):

Set the fill_value for the masked array to nan for floating-point columns. The numpy default of 1e20 is just crazy and not appropriate when we are using nan and masked interchangeably.
Make sure that the underlying masked values in the masked array are set to nan. Sometimes numpy grabs these instead of the fill value (blech) and these should reflect the actual value that we have serialized.

Navigation Menu

Repository navigation

Broken FITS table round-trip with any nan in a float column. #14558

Description

Expected behavior

How to Reproduce

Versions

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Broken FITS table round-trip with any nan in a float column. #14558

Description

Description

Expected behavior

How to Reproduce

Versions

Activity

saimn commented on Mar 23, 2023

erykoff commented on Mar 23, 2023

segasai commented on Nov 28, 2024

mhvk commented on Nov 29, 2024

segasai commented on Nov 29, 2024

mhvk commented on Nov 29, 2024

taldcroft commented on Dec 1, 2024

segasai commented on Dec 1, 2024

mhvk commented on Dec 1, 2024

mhvk commented on Dec 1, 2024

mhvk commented on Dec 1, 2024

taldcroft commented on Dec 1, 2024

taldcroft commented on Dec 1, 2024

segasai commented on Dec 1, 2024

erykoff commented on Dec 1, 2024

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions