Repository navigation
io.fits writes incorrect data values in some cases if original image was uint16 #4974
Description
The following sequence of actions results in a FITS file in which the data values are arguably incorrect, and in which reading the file with io.fits definitely produces incorrect values.
I didn't set a milestone but would argue this should be fixed for 1.2.
I believe it is due to the changes in #3916 (see the comments for some discussion that happened at the time).
Code to reproduce is below, but conceptually:
- Read in a FITS file that has unsigned int16 data in it. After the read the dtype of the data is
uint16,BITPIX=16andBZERO=32768. - Do something that converts the data to a float (e.g. multiply by
1.0). After this operation,BITPIX=-64but alsoBZERO=32768. This is the key error, I think. - Write the data to a new fits file. The correct value is written for each pixel, and the value of
BZEROis written to the header of the fits file. - That new fits file, when reopened in astropy or in the couple of other programs I tried, has data value 32768 bigger than the original data because it (sensibly) follows the instructions in the FITS header, adding BZERO to each value.
To reproduce (ipython session, plain code is the block following):
In [1]: from astropy.io import fits
In [2]: import numpy as np
In [3]: # make a FITS file and save it, all data values are 1
In [4]: fake_hdu = fits.PrimaryHDU(data=np.ones((100, 100), dtype='uint16'))
In [5]: fake_hdu.writeto('fake_uint16.fits')
In [6]: f = fits.open('fake_uint16.fits') # Open the fits file we made
In [7]: print(f[0].header['bitpix'], f[0].header['bzero'])
16 32768
In [8]: assert (f[0].data == 1).all() # Is the data still 1?
In [9]: f[0].data = f[0].data * 1.0 # Make the data floating point
In [10]: assert (f[0].data == 1).all() # Is the data still 1?
In [11]: print(f[0].header['bitpix'], f[0].header['bzero'])
-64 32768
In [12]: # ^^^^^^^^^^ bad! ^^^^^^^^^^^^^
In [13]: f.writeto('fake_bad.fits') # Write data out...maybe it will round trip.
In [14]: bad = fits.open('fake_bad.fits')
In [15]: assert (bad[0].data == 1).all() # Is the data still 1?
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-15-14ad38116ad0> in <module>()
----> 1 assert (bad[0].data == 1).all() # Is the data still 1?
AssertionError:
In [16]: # No, it is not...
In [17]: # Check header...
In [18]: print(bad[0].header['bitpix'], 'bzero' in bad[0].header)
-64 False
In [19]: assert(bad[0].data == 32768 + 1).all() # Data is BZERO + 1Plain code for reproducing:
from astropy.io import fits
import numpy as np
# make a FITS file and save it, all data values are 1
fake_hdu = fits.PrimaryHDU(data=np.ones((100, 100), dtype='uint16'))
fake_hdu.writeto('fake_uint16.fits')
f = fits.open('fake_uint16.fits') # Open the fits file we made
print(f[0].header['bitpix'], f[0].header['bzero'])
assert (f[0].data == 1).all() # Is the data still 1?
f[0].data = f[0].data * 1.0 # Make the data floating point
assert (f[0].data == 1).all() # Is the data still 1?
print(f[0].header['bitpix'], f[0].header['bzero'])
# ^^^^^^^^^^ bad! ^^^^^^^^^^^^^
f.writeto('fake_bad.fits') # Write data out...maybe it will round trip.
bad = fits.open('fake_bad.fits')
assert (bad[0].data == 1).all() # Is the data still 1?
# No, it is not...
# Check header...
print(bad[0].header['bitpix'], 'bzero' in bad[0].header)
assert(bad[0].data == 32768 + 1).all() # Data is BZERO + 1
Activity
Delete scaling keywords from header when unsigned int data converted …