fits.io fails to write numpy structured arrays with unicode columns

Using astropy/1.2.1, numpy/1.11.1, and python/3.5.2 or 2.7.11

import numpy as np
from astropy.io import fits

a = np.zeros(10, dtype=[('X', 'S5'), ('Y', 'U5')])
a['X'][0:5] = b'ELG'
a['X'][5:10] = b'LRG'
a['Y'][0:5] = 'blat'
a['Y'][5:10] = 'foo'
fits.writeto('test1.fits', a, clobber=True)
b, hdr = fits.getdata('test1.fits', header=True)

Under python 2.7.11, the fits.writeto line fails with:

AttributeError                            Traceback (most recent call last)
<ipython-input-8-1d967b090d27> in <module>()
----> 1 fits.writeto('test1.fits', a, clobber=True)

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/convenience.pyc in writeto(filename, data, header, output_verify, clobber, checksum)
    407         hdu = PrimaryHDU(data, header=header)
    408     hdu.writeto(filename, clobber=clobber, output_verify=output_verify,
--> 409                 checksum=checksum)
    410 
    411 

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/base.pyc in writeto(self, name, output_verify, clobber, checksum)
   1763         hdulist = HDUList([PrimaryHDU(), self])
   1764         hdulist.writeto(name, output_verify, clobber=clobber,
-> 1765                         checksum=checksum)
   1766 
   1767     def _verify(self, option='warn'):

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc in writeto(self, fileobj, output_verify, clobber, checksum)
    689             hdu._prewriteto(checksum=checksum)
    690             try:
--> 691                 hdu._writeto(hdulist._file)
    692             finally:
    693                 hdu._postwriteto()

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/base.pyc in _writeto(self, fileobj, inplace, copy)
    740         if not inplace or self._new:
    741             header_offset, _ = self._writeheader(fileobj)
--> 742             data_offset, data_size = self._writedata(fileobj)
    743 
    744             # Set the various data location attributes on newly-written HDUs

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/base.pyc in _writedata(self, fileobj)
    677         if self._data_loaded or self._data_needs_rescale:
    678             if self.data is not None:
--> 679                 size += self._writedata_internal(fileobj)
    680             # pad the FITS data block
    681             if size > 0:

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/table.pyc in _writedata_internal(self, fileobj)
    859                 # unicode columns directly to the file, so we have to switch
    860                 # to a slower row-by-row write
--> 861                 self._writedata_by_row(fileobj)
    862             else:
    863                 fileobj.writearray(data)

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/hdu/table.pyc in _writedata_by_row(self, fileobj)
    912                     item = np.char.encode(item, 'ascii')
    913 
--> 914                 fileobj.writearray(item)
    915 
    916     _tdump_file_format = textwrap.dedent("""

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/file.pyc in writearray(self, array)
    300 
    301         if hasattr(self._file, 'write'):
--> 302             _array_to_file(array, self._file)
    303 
    304     def flush(self):

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/io/fits/util.pyc in _array_to_file(arr, outfile)
    799     # https://github.com/numpy/numpy/issues/2256
    800 
--> 801     if (sys.platform == 'darwin' and arr.nbytes >= _OSX_WRITE_LIMIT + 1 and
    802             arr.nbytes % 4096 == 0):
    803         # chunksize is a count of elements in the array, not bytes

AttributeError: 'str' object has no attribute 'nbytes'

Under python 3.5.2, it silently writes a corrupted fits file, which can be seen by reading it back in:

In [2]: b, hdr = fits.getdata('test1.fits', header=True)
   ...: 
WARNING: File may have been truncated: actual file length (8605) is smaller than the expected size (8640) [astropy.io.fits.file]

In [3]: b
Out[3]: 
FITS_rec([('ELGbl', 'atELG'), ('blatE', 'LGbla'), ('tELGb', 'latEL'),
       ('Gblat', 'LRGfo'), ('oLRGf', 'ooLRG'), ('fooLR', 'GfooL'),
       ('RGfoo', ''), ('', ''), ('', ''), ('', '')], 
      dtype=(numpy.record, [('X', 'S5'), ('Y', 'S5')]))

In [4]: hdr
Out[4]: 
XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   10 / length of dimension 1                          
NAXIS2  =                   10 / length of dimension 2                          
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    2 / number of table fields                         
TTYPE1  = 'X       '                                                            
TFORM1  = '5A      '                                                            
TTYPE2  = 'Y       '                                                            
TFORM2  = '5A      '

Under both python 2.7.11 and 3.5.2, converting a to a Table first writes a valid fits file, but it has silently converted the unicode column into a byte column:

from astropy.table import Table
Table(a).write('test2.fits', overwrite=True)
b2, hdr2 = fits.getdata('test2.fits', header=True)

In [13]: b2
Out[13]: 
FITS_rec([('ELG', 'blat'), ('ELG', 'blat'), ('ELG', 'blat'), ('ELG', 'blat'),
       ('ELG', 'blat'), ('LRG', 'foo'), ('LRG', 'foo'), ('LRG', 'foo'),
       ('LRG', 'foo'), ('LRG', 'foo')], 
      dtype=(numpy.record, [('X', 'S5'), ('Y', 'S5')]))

In [14]: hdr2
Out[14]: 
XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   10 / length of dimension 1                          
NAXIS2  =                   10 / length of dimension 2                          
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    2 / number of table fields                         
TTYPE1  = 'X       '                                                            
TFORM1  = '5A      '                                                            
TTYPE2  = 'Y       '                                                            
TFORM2  = '5A      '

This appears to be a problem with the .writeto() or .getdata() methods. I ran the following with both Python 2.7.12/numpy 1.11.1/astropy 1.2.1 and Python 3.5.2/numpy 1.11.1/astropy 1.2.1:

import numpy as np
from astropy.io import fits
x = [b'ELG']*5 + [b'LRG']*5
y = [u'blat']*5 + [u'foo']*5  # I didn't use u'' in Python 3
x = np.array(x, dtype='S5')
y = np.array(y, dtype='U5')
col1 = fits.Column(name='X', format='5A', array=x)
col2 = fits.Column(name='Y', format='5A', array=y)
cols = fits.ColDefs([col1, col2])
tbhdu = fits.BinTableHDU.from_columns(cols)
tbhdu.writeto('test-py2.fits')  # Replace py2 with py3 for Python 3

Not only did I get two valid FITS files, the files were identical at the bit level. Reading the data back in also produces no errors.

@weaverba137

digging a bit further, here is another example that gives the same failures:

a = np.zeros(10, dtype=[('X', 'S5'), ('Y', 'U5')])
a['X'][0:5] = b'ELG'
a['X'][5:10] = b'LRG'
a['Y'][0:5] = 'blat'
a['Y'][5:10] = 'foo'
hdu1 = fits.PrimaryHDU()
hdu2 = fits.BinTableHDU(a)
hx = fits.HDUList([hdu1, hdu2])
hx.writeto('blat.fits')

Comparing the tbhdu created by hand from @weaverba137's example with the hdu2 created here:

In [22]: tbhdu.columns
Out[22]: 
ColDefs(
    name = 'X'; format = '5A'
    name = 'Y'; format = '5A'
)

In [23]: hdu2.columns
Out[23]: 
ColDefs(
    name = 'X'; format = '5A'
    name = 'Y'; format = '5A'
)

In [24]: tbhdu.data
Out[24]: 
FITS_rec([('ELG', 'blat'), ('ELG', 'blat'), ('ELG', 'blat'), ('ELG', 'blat'),
       ('ELG', 'blat'), ('LRG', 'foo'), ('LRG', 'foo'), ('LRG', 'foo'),
       ('LRG', 'foo'), ('LRG', 'foo')], 
      dtype=(numpy.record, [('X', 'S5'), ('Y', 'S5')]))

In [25]: hdu2.data
Out[25]: 
FITS_rec([('ELG', u'blat'), ('ELG', u'blat'), ('ELG', u'blat'),
       ('ELG', u'blat'), ('ELG', u'blat'), ('LRG', u'foo'),
       ('LRG', u'foo'), ('LRG', u'foo'), ('LRG', u'foo'), ('LRG', u'foo')], 
      dtype=(numpy.record, [('X', 'S5'), ('Y', '<U5')]))

i.e. they both got the FITS format right ( @weaverba137 set it explicitly; this example allows BinTableHDU to figure it out), but BinTableHDU isn't correctly downcasting the unicode column into an 'S5' array for output.

Someone can correct me if I am wrong... I think if you want FITS table here, you need to create the Column objects first, then construct a BinTableHDU using those columns, before writing it out. Not sure if there is a real bug here, so I am removing the "Bug" tag for now.

See http://docs.astropy.org/en/stable/io/fits/usage/table.html#creating-a-fits-table

It would be much easier to use the Table interface - I can include an example later.

So is astropy.io.fits.writeto() officially deprecated? If the only method for creating a fits table that is officially supported is the column-by-column method, then the methods that are "supposed" to be able to create a FITS table directly from a numpy array should be deprecated and removed.

@weaverba137

@weaverba137 , I never used that myself but my impression from the documentation of fits.writeto() convenience function is that you still have to create the FITS_rec yourself (it is not automagically creating a table for you). From my brief test, fits.writeto('somefile.fits', tbhdu.data) seems to work fine.

Ref: http://docs.astropy.org/en/latest/io/fits/api/files.html#astropy.io.fits.writeto

@astrofrog

@astrofrog See issue #5280 for some unit tests using the Table interface that fail under python 3.5 but pass under 2.7. Consider closing this ticket as a special case of the suite of failing tests in that ticket.

@pllim

@pllim the documentation for astropy.io.fits.writeto describes the data array as

    data : array, record array, or groups data object
        data to write to the new file

I would interpret that as meaning that a normal numpy array should work; not just an array as constructed via fits.BinTableHDU.from_columns(cols).

Note that the problem also occurs when using io.fits.BinTableHDU(data); i.e. the problem reported here isn't with the convenience function io.fits.writeto, but rather with the BinTableHDU(data) constructor (I think...)

My point is that you are not supposed to use BinTableHDU(data) but rather explicitly construct the columns, as documented. Perhaps this is merely a documentation issue rather than a bug.

@pllim

@pllim, if you are "not supposed to use it", why does it accept that input? Documentation or no, the code should reject invalid inputs!

Check this out:

Python 3.5:

>>> import numpy as np
>>> from astropy.io.fits.fitsrec import FITS_rec
>>> data = np.zeros((10,), dtype=[('X', 'S5'), ('Y', 'U5')])
>>> data['X'] = [b'ELG']*5 + [b'LRG']*5
>>> data['Y'] = ['foo']* 5 + ['blat']*5
>>> data
array([(b'ELG', 'foo'), (b'ELG', 'foo'), (b'ELG', 'foo'), (b'ELG', 'foo'),
       (b'ELG', 'foo'), (b'LRG', 'blat'), (b'LRG', 'blat'),
       (b'LRG', 'blat'), (b'LRG', 'blat'), (b'LRG', 'blat')], 
      dtype=[('X', 'S5'), ('Y', '<U5')])
>>> data2 = data.view(FITS_rec)
>>> data2
FITS_rec([('ELG', 'foo'), ('ELG', 'foo'), ('ELG', 'foo'), ('ELG', 'foo'),
       ('ELG', 'foo'), ('LRG', 'blat'), ('LRG', 'blat'), ('LRG', 'blat'),
       ('LRG', 'blat'), ('LRG', 'blat')], 
      dtype=[('X', 'S5'), ('Y', '<U5')])

What happened to the 'X' values? They were correctly displayed as bytes as a basic numpy array, but now are being displayed as str when viewed as FITS_rec.

Here's Python 2.7 for comparison:

>>> import numpy as np
>>> from astropy.io.fits.fitsrec import FITS_rec
>>> data = np.zeros((10,), dtype=[('X', 'S5'), ('Y', 'U5')])
>>> data['X'] = ['ELG']*5 + ['LRG']*5
>>> data['Y'] = [u'foo']* 5 + [u'blat']*5
>>> data
array([('ELG', u'foo'), ('ELG', u'foo'), ('ELG', u'foo'), ('ELG', u'foo'),
       ('ELG', u'foo'), ('LRG', u'blat'), ('LRG', u'blat'),
       ('LRG', u'blat'), ('LRG', u'blat'), ('LRG', u'blat')], 
      dtype=[('X', 'S5'), ('Y', '<U5')])
>>> data2 = data.view(FITS_rec)
>>> data2
FITS_rec([('ELG', u'foo'), ('ELG', u'foo'), ('ELG', u'foo'), ('ELG', u'foo'),
       ('ELG', u'foo'), ('LRG', u'blat'), ('LRG', u'blat'),
       ('LRG', u'blat'), ('LRG', u'blat'), ('LRG', u'blat')], 
      dtype=[('X', 'S5'), ('Y', '<U5')])

This looks correct.

BinTableHDU(data) works under 2.7, and it works under 3.5 for tables that don't have unicode strings, and the BinTableHDU.__init__ docstring implies that it should work with a data array input, and astropy.io.fits.convenience.py itself uses hdu = BinTableHDU(data) on vanilla numpy arrays (line 889). I.e. I think this is a real bug, not just a documentation issue.

@astrofrog

I am hesitant to invoke embray again today... But for now, just use from_columns as a workaround (or the example that @astrofrog promised).

@sbailey

As @sbailey suggested above, I am closing this one and we can continue the discussions on #5280. Thanks for bringing this up!

Navigation Menu

Repository navigation

fits.io fails to write numpy structured arrays with unicode columns #5268

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

fits.io fails to write numpy structured arrays with unicode columns #5268

Description

Activity

weaverba137 commented on Aug 29, 2016

sbailey commented on Aug 30, 2016

pllim commented on Aug 30, 2016

pllim commented on Aug 30, 2016

astrofrog commented on Aug 30, 2016

weaverba137 commented on Aug 30, 2016

pllim commented on Aug 30, 2016

sbailey commented on Aug 30, 2016

sbailey commented on Aug 30, 2016

pllim commented on Aug 30, 2016

weaverba137 commented on Aug 30, 2016

weaverba137 commented on Aug 30, 2016

sbailey commented on Aug 30, 2016

pllim commented on Aug 30, 2016

pllim commented on Aug 30, 2016

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions