mhvk · Jan 21, 2021 · Jan 7, 2021 · Jan 7, 2021 · Jan 7, 2021 · Jan 7, 2021
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -138,6 +138,15 @@ astropy.io.ascii
 astropy.io.fits
 ^^^^^^^^^^^^^^^
 
+- For conversion between FITS tables and astropy ``Table``, the standard mask
+  values of ``NaN`` for float and null string for string are now properly
+  recognized, leading to a ``MaskedColumn`` with appropriately set mask
+  instead of a ``Column`` with those values exposed. Conversely, when writing
+  an astropy ``Table`` to a FITS tables, masked values are now consistently
+  converted to the standard FITS mask values of ``NaN`` for float and null
+  string for string (i.e., not just for tables with ``masked=True``, which no
+  longer is guaranteed to signal the presence of ``MaskedColumn``). [#11222]
+
 astropy.io.misc
 ^^^^^^^^^^^^^^^
 
@@ -165,6 +174,10 @@ astropy.table
   This makes accessing a column's underlying data array consistent with the
   ``.value`` attribute available for ``Time`` and ``Quantity`` objects. [#10962]
 
+- In reading from a FITS tables, the standard mask values of ``NaN`` for float
+  and null string for string are properly recognized, leading to a
+  ``MaskedColumn`` with appropriately set mask. [#11222]
+
 astropy.tests
 ^^^^^^^^^^^^^
 
@@ -258,6 +271,10 @@ astropy.table
   instance (which would lead to slow-down of, e.g., slicing afterwards).
   [#11077]
 
+- When writing to a FITS tables, masked values are again always converted to
+  the standard FITS mask values of ``NaN`` for float and null string
+  for string, not just for table with ``masked=True``. [#11222]
+
 astropy.tests
 ^^^^^^^^^^^^^
 

diff --git a/astropy/io/fits/connect.py b/astropy/io/fits/connect.py
@@ -6,6 +6,8 @@
 import warnings
 from copy import deepcopy
 
+import numpy as np
+
 from astropy.io import registry as io_registry
 from astropy import units as u
 from astropy.table import Table, serialize, meta, Column, MaskedColumn
@@ -227,25 +229,30 @@ def read_table_fits(input, hdu=None, astropy_native=False, memmap=False,
         finally:
             hdulist.close()
 
-    # Check if table is masked
-    masked = any(col.null is not None for col in table.columns)
-
-    # TODO: in future, it may make more sense to do this column-by-column,
-    # rather than via the structured array.
-
     # In the loop below we access the data using data[col.name] rather than
     # col.array to make sure that the data is scaled correctly if needed.
     data = table.data
 
     columns = []
     for col in data.columns:
-
-        # Set column data
-        if masked:
-            column = MaskedColumn(data=data[col.name], name=col.name, copy=False)
-            if col.null is not None:
-                column.set_fill_value(col.null)
-                column.mask[column.data == col.null] = True
+        # Check if column is masked. Here, we make a guess based on the
+        # presence of FITS mask values. For integer columns, this is simply
+        # the null header, for float and complex, the presence of NaN, and for
+        # string, empty strings.
+        masked = mask = False
+        if col.null is not None:
+            mask = col.array == col.null
+            # Return a MaskedColumn even if no elements are masked so
+            # we roundtrip better.
+            masked = True
+        elif issubclass(col.dtype.type, np.inexact):
+            mask = np.isnan(col.array)
+        elif issubclass(col.dtype.type, np.character):
+            mask = col.array == b''
+
+        if masked or np.any(mask):
+            column = MaskedColumn(data=data[col.name], name=col.name,
+                                  mask=mask, copy=False)
         else:
             column = Column(data=data[col.name], name=col.name, copy=False)
 
@@ -260,7 +267,7 @@ def read_table_fits(input, hdu=None, astropy_native=False, memmap=False,
         columns.append(column)
 
     # Create Table object
-    t = Table(columns, masked=masked, copy=False)
+    t = Table(columns, copy=False)
 
     # TODO: deal properly with unsigned integers
 

diff --git a/astropy/io/fits/convenience.py b/astropy/io/fits/convenience.py
@@ -473,33 +473,42 @@ def table_to_hdu(table, character_as_bytes=False):
             table, hdr = time_to_fits(table)
 
     # Create a new HDU object
-    if table.masked:
-        # float column's default mask value needs to be Nan
-        for column in table.columns.values():
-            fill_value = column.get_fill_value()
-            if column.dtype.kind == 'f' and np.allclose(fill_value, 1e20):
-                column.set_fill_value(np.nan)
+    tarray = table.as_array()
+    if isinstance(tarray, np.ma.MaskedArray):
+        # Fill masked values carefully:
+        # float column's default mask value needs to be Nan and
+        # string column's default mask should be an empty string.
+        # Note: getting the fill value for the structured array is
+        # more reliable than for individual columns for string entries.
+        # (no 'N/A' for a single-element string, where it should be 'N').
+        default_fill_value = np.ma.default_fill_value(tarray.dtype)
+        for colname, (coldtype, _) in tarray.dtype.fields.items():
+            if np.all(tarray.fill_value[colname] == default_fill_value[colname]):
+                if issubclass(coldtype.type, np.complexfloating):
+                    tarray.fill_value[colname] = complex(np.nan, np.nan)
+                elif issubclass(coldtype.type, np.inexact):
+                    tarray.fill_value[colname] = np.nan
+                elif issubclass(coldtype.type, np.character):
+                    tarray.fill_value[colname] = ''
 
         # TODO: it might be better to construct the FITS table directly from
         # the Table columns, rather than go via a structured array.
-        table_hdu = BinTableHDU.from_columns(np.array(table.filled()), header=hdr, character_as_bytes=True)
+        table_hdu = BinTableHDU.from_columns(tarray.filled(), header=hdr,
+                                             character_as_bytes=character_as_bytes)
         for col in table_hdu.columns:
             # Binary FITS tables support TNULL *only* for integer data columns
             # TODO: Determine a schema for handling non-integer masked columns
-            # in FITS (if at all possible)
+            # with non-default fill values in FITS (if at all possible).
             int_formats = ('B', 'I', 'J', 'K')
             if not (col.format in int_formats or
                     col.format.p_format in int_formats):
                 continue
 
-            # The astype is necessary because if the string column is less
-            # than one character, the fill value will be N/A by default which
-            # is too long, and so no values will get masked.
-            fill_value = table[col.name].get_fill_value()
-
-            col.null = fill_value.astype(table[col.name].dtype)
+            fill_value = tarray[col.name].fill_value
+            col.null = fill_value.astype(int)
     else:
-        table_hdu = BinTableHDU.from_columns(np.array(table.filled()), header=hdr, character_as_bytes=character_as_bytes)
+        table_hdu = BinTableHDU.from_columns(tarray, header=hdr,
+                                             character_as_bytes=character_as_bytes)
 
     # Set units and format display for output HDU
     for col in table_hdu.columns:

diff --git a/astropy/io/fits/tests/test_connect.py b/astropy/io/fits/tests/test_connect.py
@@ -1,12 +1,11 @@
-import os
 import gc
 import sys
 import pathlib
 import warnings
 
 import pytest
 import numpy as np
-from numpy.testing import assert_allclose
+from numpy.testing import assert_allclose, assert_array_equal
 
 from astropy.io.fits.column import (_parse_tdisp_format, _fortran_to_python_format,
                                     python_to_tdisp)
@@ -168,31 +167,35 @@ def test_masked(self, tmpdir):
         t1.mask['c'] = [0, 1, 1, 0]
         t1.write(filename, overwrite=True)
         t2 = Table.read(filename)
-        assert t2.masked
         assert equal_data(t1, t2)
         assert np.all(t1['a'].mask == t2['a'].mask)
-        # Disabled for now, as there is no obvious way to handle masking of
-        # non-integer columns in FITS
-        # TODO: Re-enable these tests if some workaround for this can be found
-        # assert np.all(t1['b'].mask == t2['b'].mask)
-        # assert np.all(t1['c'].mask == t2['c'].mask)
+        assert np.all(t1['b'].mask == t2['b'].mask)
+        assert np.all(t1['c'].mask == t2['c'].mask)
 
-    def test_masked_nan(self, tmpdir):
+    @pytest.mark.parametrize('masked', [True, False])
+    def test_masked_nan(self, masked, tmpdir):
         filename = str(tmpdir.join('test_masked_nan.fits'))
-        data = np.array(list(zip([5.2, 8.4, 3.9, 6.3],
-                                 [2.3, 4.5, 6.7, 8.9])),
-                        dtype=[('a', np.float64), ('b', np.float32)])
-        t1 = Table(data, masked=True)
-        t1.mask['a'] = [1, 0, 1, 0]
-        t1.mask['b'] = [1, 0, 0, 1]
+        a = np.ma.MaskedArray([5.25, 8.5, 3.75, 6.25], mask=[1, 0, 1, 0])
+        b = np.ma.MaskedArray([2.5, 4.5, 6.75, 8.875], mask=[1, 0, 0, 1], dtype='f4')
+        t1 = Table([a, b], names=['a', 'b'], masked=masked)
         t1.write(filename, overwrite=True)
         t2 = Table.read(filename)
-        np.testing.assert_array_almost_equal(t2['a'], [np.nan, 8.4, np.nan, 6.3])
-        np.testing.assert_array_almost_equal(t2['b'], [np.nan, 4.5, 6.7, np.nan])
-        # assert t2.masked
-        # t2.masked = false currently, as the only way to determine whether a table is masked
-        # while reading is to check whether col.null is present. For float columns, col.null
-        # is not initialized
+        assert_array_equal(t2['a'].data, [np.nan, 8.5, np.nan, 6.25])
+        assert_array_equal(t2['b'].data, [np.nan, 4.5, 6.75, np.nan])
+        assert np.all(t1['a'].mask == t2['a'].mask)
+        assert np.all(t1['b'].mask == t2['b'].mask)
+
+    def test_masked_serialize_data_mask(self, tmpdir):
+        filename = str(tmpdir.join('test_masked_nan.fits'))
+        a = np.ma.MaskedArray([5.25, 8.5, 3.75, 6.25], mask=[1, 0, 1, 0])
+        b = np.ma.MaskedArray([2.5, 4.5, 6.75, 8.875], mask=[1, 0, 0, 1])
+        t1 = Table([a, b], names=['a', 'b'])
+        t1.write(filename, overwrite=True, serialize_method='data_mask')
+        t2 = Table.read(filename)
+        assert_array_equal(t2['a'].data, [5.25, 8.5, 3.75, 6.25])
+        assert_array_equal(t2['b'].data, [2.5, 4.5, 6.75, 8.875])
+        assert np.all(t1['a'].mask == t2['a'].mask)
+        assert np.all(t1['b'].mask == t2['b'].mask)
 
     def test_read_from_fileobj(self, tmpdir):
         filename = str(tmpdir.join('test_read_from_fileobj.fits'))
@@ -493,9 +496,9 @@ def test_masking_regression_1795():
     """
     t = Table.read(get_pkg_data_filename('data/tb.fits'))
     assert np.all(t['c1'].mask == np.array([False, False]))
-    assert np.all(t['c2'].mask == np.array([False, False]))
-    assert np.all(t['c3'].mask == np.array([False, False]))
-    assert np.all(t['c4'].mask == np.array([False, False]))
+    assert not hasattr(t['c2'], 'mask')
+    assert not hasattr(t['c3'], 'mask')
+    assert not hasattr(t['c4'], 'mask')
     assert np.all(t['c1'].data == np.array([1, 2]))
     assert np.all(t['c2'].data == np.array([b'abc', b'xy ']))
     assert_allclose(t['c3'].data, np.array([3.70000007153, 6.6999997139]))

diff --git a/astropy/io/fits/tests/test_convenience.py b/astropy/io/fits/tests/test_convenience.py
@@ -6,6 +6,7 @@
 
 import pytest
 import numpy as np
+from numpy.testing import assert_array_equal
 
 from astropy.io import fits
 from astropy import units as u
@@ -61,6 +62,26 @@ def test_table_to_hdu(self):
         filename = self.temp('test_table_to_hdu.fits')
         hdu.writeto(filename, overwrite=True)
 
+    def test_masked_table_to_hdu(self):
+        i = np.ma.MaskedArray([1, 2, 3], mask=[True, False, False])
+        s = np.ma.MaskedArray(['a', 'b', 'c'], mask=[False, True, True])
+        c = np.ma.MaskedArray([2.3+1j, 4.5+0j, 6.7-1j], mask=[True, False, True])
+        f = np.ma.MaskedArray([2.3, 4.5, 6.7], mask=[True, False, True])
+        table = Table([i, s, c, f], names=['i', 's', 'c', 'f'])
+        # Check that FITS standard is used in replacing masked values.
+        hdu = fits.table_to_hdu(table)
+        assert isinstance(hdu, fits.BinTableHDU)
+        assert hdu.header['TNULL1'] == i.fill_value
+        assert_array_equal(hdu.data['i'], i.filled())
+        assert_array_equal(hdu.data['s'], s.filled(''))
+        assert_array_equal(hdu.data['c'], c.filled(np.nan))
+        assert_array_equal(hdu.data['c'].real, c.real.filled(np.nan))
+        assert_array_equal(hdu.data['c'].imag, c.imag.filled(np.nan))
+        assert_array_equal(hdu.data['c'], c.filled(complex(np.nan, np.nan)))
+        assert_array_equal(hdu.data['f'], f.filled(np.nan))
+        filename = self.temp('test_table_to_hdu.fits')
+        hdu.writeto(filename, overwrite=True)
+
     def test_table_non_stringifyable_unit_to_hdu(self):
         table = Table([[1, 2, 3], ['a', 'b', 'c'], [2.3, 4.5, 6.7]],
                       names=['a', 'b', 'c'], dtype=['i', 'U1', 'f'])

diff --git a/astropy/timeseries/io/kepler.py b/astropy/timeseries/io/kepler.py
@@ -4,7 +4,7 @@
 import numpy as np
 
 from astropy.io import registry, fits
-from astropy.table import Table
+from astropy.table import Table, MaskedColumn
 from astropy.time import Time, TimeDelta
 
 from astropy.timeseries.sampled import TimeSeries
@@ -62,10 +62,15 @@ def kepler_fits_reader(filename):
         tab.rename_column("T", "TIME")
 
     for colname in tab.colnames:
+        unit = tab[colname].unit
+        # Make masks nan for any column which will turn into a Quantity
+        # later. TODO: remove once we support Masked Quantities properly?
+        if unit and isinstance(tab[colname], MaskedColumn):
+            tab[colname] = tab[colname].filled(np.nan)
         # Fix units
-        if tab[colname].unit == 'e-/s':
+        if unit == 'e-/s':
             tab[colname].unit = 'electron/s'
-        if tab[colname].unit == 'pixels':
+        if unit == 'pixels':
             tab[colname].unit = 'pixel'
 
         # Rename columns to lowercase