Data mismatch when converting FITS random group data to binary table by embray · Pull Request #3561 · astropy/astropy

embray · 2015-03-28T00:03:27Z

Hi -- was trying to convert a FITS random group into a binary table, and I went a bit crazy because the data in wasn't the same as the data out.

Here's a script that will reproduce what I was seeing (it will wget a small file from the fits registry):

import os
from astropy.io import fits as pf
import numpy as np

filename_out = "test_bin.fits"

# Load random group data
# here's one from the fits registry that's pretty small:
# Shape of DATA is (7956, 1, 1, 1, 4, 3)
os.system("wget http://fits.gsfc.nasa.gov/samples/DDTSUVDATA.fits")
a = pf.open('DDTSUVDATA.fits')
d_in = a[0].data["DATA"]

# Convert to column
c_a =  pf.Column(name='TEST', array=d_in, dim="(3,4,1,1,1)", format="12E")
hdu = pf.BinTableHDU.from_columns([c_a])

# Remove output files
if os.path.exists(filename_out):
    os.remove(filename_out)

# Write to a file, then reload the file
hdu.writeto(filename_out)
hdul_b = pf.open(filename_out)
d_out = hdul_b[1].data["TEST"]

# Compare the two files
print d_in.shape, d_in.dtype
print d_out.shape, d_in.dtype

print d_in[0]
print d_out[0]

print np.allclose(d_in[0], d_out[0])

The print outputs are:

(7956, 1, 1, 1, 4, 3) >f4
(7956, 1, 1, 1, 4, 3) >f4
[[[[[ 12.4308672    0.56860745   3.99993873]
    [ 12.74043655   0.31398511   3.99993873]
    [  0.           0.           3.99993873]
    [  0.           0.           3.99993873]]]]]
[[[[[ -3.13737987e+13   7.20336838e+01              nan]
    [ -7.43180390e+12  -2.16078206e-14              nan]
    [  0.00000000e+00   0.00000000e+00              nan]
    [  0.00000000e+00   0.00000000e+00              nan]]]]]
False

It seems to be related with endianess, as if I change the line to say

d_in = a[0].data["DATA"].astype('<f4')

everything is dandy.

I would have thought that this kind of thing should be handled transparently?

PS: I'm running Python 2.7.6 on Mac yosemite, with astropy 1.0, numpy 1.9.2.

See astropy/astropy#3561 for why I was going crazy for several hours

embray · 2015-03-02T22:25:40Z

Interesting--normally that should work fine, though it's not the first time there have been bugs involving random groups data, as the feature is just not as well tested. If I repeat a similar use case, but using a normal binary table as the original file, this works as it should. So it seems somewhere along the way it's byte swapping when it shouldn't be.

embray · 2015-03-02T22:38:40Z

I think I see now. This actually doesn't have to do with it being random groups data for once, and has to do with the fact that it uses TDIMn, and that is going down a code path that's resulting in an unnecessary byte swap in this particular case.

telegraphic · 2015-03-02T23:41:48Z

Glad to hear you're narrowing it down! To add a related enhancement idea: shouldn't the column constructor be able to figure out the dim and format from the array data shape and dtype, if an array is passed?

I'm currently figuring these out with this code (which probably has a nicer equivalent in the bowels of pyfits):

def fits_format_code_lookup(numpy_dtype, numpy_shape):
    """ Return a FITS format code from a given numpy dtype

    Parameters
    ----------
    numpy_dtype: a numpy dtype object
        Numpy dtype to lookup
    numpy_shape: tuple
        shape of data array to be converted into a FITS format

    Returns
    -------
    fmt_code: string
        FITS format code, e.g. 8A for character string of length 8
    fits_dim: string or None
        Returns fits dimension for TDIM keyword

    Notes
    -----

    FITS format code         Description                     8-bit bytes
    --------------------------------------------------------------------
    L                        logical (Boolean)               1
    X                        bit                             *
    B                        Unsigned byte                   1
    I                        16-bit integer                  2
    J                        32-bit integer                  4
    K                        64-bit integer                  4
    A                        character                       1
    E                        single precision floating point 4
    D                        double precision floating point 8
    C                        single precision complex        8
    M                        double precision complex        16
    P                        array descriptor                8
    Q                        array descriptor                16
    """

    np_type = numpy_dtype.type


    # Get length of array last axis
    fmt_len = 1
    if len(numpy_shape) > 1:
        for ii in numpy_shape[1:]:
            fmt_len *=  ii
        #fmt_len = numpy_shape[-1]
    else:
        fmt_len = 1

    if np_type is np.string_ or "string" in str(np_type):
        fits_dim = None
        if numpy_dtype.itemsize == 1:
            fits_fmt = "A"
        else:
            fits_fmt = '%iA' % numpy_dtype.itemsize
        return fits_fmt, fits_dim

    else:
        #TODO: handle special case of uint16,32,64, requires BZERO / scale
        dtype_dict = {
            np.uint8:   'B',
            np.uint16:  'J',    # No native FITS equivalent
            np.uint32:  'K',
            np.uint64:  'K',
            np.int8:    'I',    # No native FITS equivalent
            np.int16:   'I',
            np.int32:   'J',
            np.int64:   'K',
            np.float16: 'E',    # No native FITS equivalent
            np.float32: 'E',
            np.float64: 'D',
            np.complex64: 'C',
            np.complex128: 'M',
            np.bool_:      'L'
        }

        fmt_code = dtype_dict.get(np_type)

        if fmt_len > 1:
            fmt_code = "%i%s" % (fmt_len, fmt_code)

        if len(numpy_shape) == 1:
            fits_dim = None
        else:
            #NOTE: Needs to be flipped to be FORTRAN order
            # From FITS spec: The value field of this indexed keyword shall contain a
            # character string describing how to interpret the contents of field n as
            # a multi-dimensional array with a format of “(l, m, n ...)” where l, m, n, ...
            # are the dimensions of the array. The data are ordered such that the array index
            # of the first dimension given (l) is the most rapidly varying and that of the last
            # dimension given is the least rapidly varying.
            fits_dim = str(numpy_shape[1:][::-1])

        return fmt_code, fits_dim

Probably not the most concise, but if you agree that being able to deduce dim and format from the array would be benefical, please feel to grab chunks of it. There are no doubt edge cases and things like BZERO and BSCALE to worry about though :/

embray · 2015-03-03T18:00:29Z

For the most part it does do this internally. You can, for example, eschew the Column interface entirely and initialize a table HDU from an existing structured Numpy array:

In [150]: a = np.array([(np.arange(12).reshape(1, 1, 1, 3, 4),)], dtype=[('test', '(1,1,1,3,4)f8')])

In [151]: table = fits.BinTableHDU.from_columns(a)

In [152]: table.header
Out[152]: 
XTENSION= 'BINTABLE'           / binary table extension                         
BITPIX  =                    8 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                   96 / length of dimension 1                          
NAXIS2  =                    1 / length of dimension 2                          
PCOUNT  =                    0 / number of group parameters                     
GCOUNT  =                    1 / number of groups                               
TFIELDS =                    1 / number of table fields                         
TTYPE1  = 'test    '                                                            
TFORM1  = '12D     '                                                            
TDIM1   = '(4,3,1,1,1)'                                                         

In [153]: table.data['TEST']
Out[153]: 
array([[[[[[  0.,   1.,   2.,   3.],
           [  4.,   5.,   6.,   7.],
           [  8.,   9.,  10.,  11.]]]]]])

In [154]: table.data['TEST'].shape
Out[154]: (1, 1, 1, 1, 3, 4)

The old Column interface is useful for users familiar with FITS, but IMO should be eventually de-emphasized (the plan is to deprecate it anyways). Another option is to directly use Astropy's Table interface and write it directly to a FITS file. Unfortunately this just doesn't support all kinds of FITS data yet and can't write, for example, random groups arrays. Anyways, thank you for the code sample, but this already exists. It just isn't used through the Column interface which assumes you're sort of designing a FITS format from scratch, rather than just trying to write some existing Numpy arrays to a FITS file.

embray · 2015-03-03T18:22:11Z

Relatedly, see #3552

telegraphic · 2015-03-03T19:30:18Z

Ah, thanks for pointing me toward the structured array route; I suspected that I was probably reinventing the wheel when I wrote it.

embray · 2015-03-03T21:47:24Z

It's a shame that isn't better documented--that much should be fixed,.

…special case for string types where the string lengths don't match up (i.e. we are cramming a 2x2 array of 2 character strings into a single 8 character string). It shouldn't apply to numeric types, and the 'view' here also just screws up byteorder

… summarize for the changelog =_=

embray · 2015-03-28T00:04:05Z

Fix for the underlying bug reported in this issue attached (the existing workaround should still be considered preferable to the original use case in the first place though).

Data mismatch when converting FITS random group data to binary table

… summarize for the changelog =_=

Data mismatch when converting FITS random group data to binary table

telegraphic added a commit to telegraphic/fits2hdf that referenced this pull request Mar 2, 2015

Improving random group round trip

c0d636e

See astropy/astropy#3561 for why I was going crazy for several hours

astrofrog added the io.fits label Mar 2, 2015

astrofrog assigned embray Mar 2, 2015

embray added Bug Affects-release labels Mar 2, 2015

embray mentioned this pull request Mar 3, 2015

Mention Table interface in "Creating a New Table File" in FITS docs #3552

Closed

embray added 3 commits March 27, 2015 19:53

Add a regression test demonstrating the problem

b769a0c

Add changelog entry for astropy#3561. Another bug that's difficult to…

2d7330c

… summarize for the changelog =_=

embray added this to the v1.0.2 milestone Mar 28, 2015

embray added a commit that referenced this pull request Mar 28, 2015

Merge pull request #3561 from embray/fits/issue-3561

09aab69

Data mismatch when converting FITS random group data to binary table

embray merged commit 09aab69 into astropy:master Mar 28, 2015

embray deleted the fits/issue-3561 branch March 28, 2015 21:03

dhomeier pushed a commit to dhomeier/astropy that referenced this pull request Mar 31, 2015

Add changelog entry for astropy#3561. Another bug that's difficult to…

696157f

… summarize for the changelog =_=

embray added a commit that referenced this pull request Apr 13, 2015

Merge pull request #3561 from embray/fits/issue-3561

057e85b

Data mismatch when converting FITS random group data to binary table

Navigation Menu

Repository navigation

Data mismatch when converting FITS random group data to binary table - #3561

Data mismatch when converting FITS random group data to binary table#3561
embray merged 3 commits into
astropy:masterfrom
embray:fits/issue-3561

embray commented Mar 28, 2015

Uh oh!

embray commented Mar 2, 2015

Uh oh!

embray commented Mar 2, 2015

Uh oh!

telegraphic commented Mar 2, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

telegraphic commented Mar 3, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

embray commented Mar 28, 2015

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

Conversation

embray commented Mar 28, 2015

Uh oh!

embray commented Mar 2, 2015

Uh oh!

embray commented Mar 2, 2015

Uh oh!

telegraphic commented Mar 2, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

telegraphic commented Mar 3, 2015

Uh oh!

embray commented Mar 3, 2015

Uh oh!

embray commented Mar 28, 2015

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants