And in my own experiments it seems if a single tile contains a NaN anywhere but the first element of the tile (curiously) the entire tile will contain NaNs when decompressed.
Furthermore, there is a separate but related bug where if a NaN occurs anywhere in the tile (except the very first element) the ZSCALE for the tile is also set to NaN, causing every pixel in the tile to be "scaled" to NaN when decompressing.
The way NaN handling works in CFITSIO is a little strange to begin with, and there are probably some bugs on the Astropy side of not pushing the right buttons to get it to do this properly.
>>> ab = np.arange(10000, dtype=np.double).reshape(100, 100)
>>> ab[0, 4] = np.nan
>>> hdu = fits.CompImageHDU(data=ab)
>>> hdu.writeto('test.fits', overwrite=True)
>>> hdul = fits.open('test.fits')
>>> hdul[1].data[0]
array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan, nan,
nan, nan, nan, nan, nan, nan, nan, nan, nan])
>>> hdul[1].data[1]
array([100., 101., 102., 103., 104., 105., 106., 107., 108., 109., 110.,
111., 112., 113., 114., 115., 116., 117., 118., 119., 120., 121.,
122., 123., 124., 125., 126., 127., 128., 129., 130., 131., 132.,
133., 134., 135., 136., 137., 138., 139., 140., 141., 142., 143.,
144., 145., 146., 147., 148., 149., 150., 151., 152., 153., 154.,
155., 156., 157., 158., 159., 160., 161., 162., 163., 164., 165.,
166., 167., 168., 169., 170., 171., 172., 173., 174., 175., 176.,
177., 178., 179., 180., 181., 182., 183., 184., 185., 186., 187.,
188., 189., 190., 191., 192., 193., 194., 195., 196., 197., 198.,
199.])
>>> hdul[1].compressed_data[0]
(array([128, 0, 0, 0, 0, 0, 0], dtype=uint8), array([], dtype=uint8), nan, 49.5)
>>> hdul[1].compressed_data.dtype
dtype((numpy.record, [('COMPRESSED_DATA', '>i4', (2,)), ('GZIP_COMPRESSED_DATA', '>i4', (2,)), ('ZSCALE', '>f8'), ('ZZERO', '>f8')]))
>>> hdul[1].header
XTENSION= 'IMAGE ' / Image extension
BITPIX = -64 / data type of original image
NAXIS = 2 / dimension of original image
NAXIS1 = 100 / length of original image axis
NAXIS2 = 100 / length of original image axis
PCOUNT = 0 / number of parameters
GCOUNT = 1 / number of groups
>>> hdul[1]._header
XTENSION= 'BINTABLE' / binary table extension
BITPIX = 8 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 32 / width of table in bytes
NAXIS2 = 100 / number of rows in table
PCOUNT = 21317 / number of group parameters
GCOUNT = 1 / number of groups
TFIELDS = 4 / number of fields in each row
TTYPE1 = 'COMPRESSED_DATA' / label for field 1
TFORM1 = '1PB(7) ' / data format of field: variable length array
TTYPE2 = 'GZIP_COMPRESSED_DATA' / label for field 2
TFORM2 = '1PB(264)' / data format of field: variable length array
TTYPE3 = 'ZSCALE ' / label for field 3
TFORM3 = '1D ' / data format of field: 8-byte DOUBLE
TTYPE4 = 'ZZERO ' / label for field 4
TFORM4 = '1D ' / data format of field: 8-byte DOUBLE
ZIMAGE = T / extension contains compressed image
ZTENSION= 'IMAGE ' / Image extension
ZBITPIX = -64 / data type of original image
ZNAXIS = 2 / dimension of original image
ZNAXIS1 = 100 / length of original image axis
ZNAXIS2 = 100 / length of original image axis
ZPCOUNT = 0 / number of parameters
ZGCOUNT = 1 / number of groups
ZTILE1 = 100 / size of tiles to be compressed
ZTILE2 = 1 / size of tiles to be compressed
ZCMPTYPE= 'RICE_1 ' / compression algorithm
ZNAME1 = 'BLOCKSIZE' / compression block size
ZVAL1 = 32 / pixels per block
ZNAME2 = 'BYTEPIX ' / bytes per pixel (1, 2, 4, or 8)
ZVAL2 = 4 / bytes per pixel (1, 2, 4, or 8)
ZNAME3 = 'NOISEBIT' / floating point quantization level
ZVAL3 = 16.0 / floating point quantization level
ZQUANTIZ= 'NO_DITHER' / No dithering during quantization
EXTNAME = 'COMPRESSED_IMAGE' / name of this binary table extension
Description
Originally raised here: https://stackoverflow.com/questions/65571367/how-to-get-nans-to-work-in-rice-compressed-fits-files-with-astropy
And in my own experiments it seems if a single tile contains a NaN anywhere but the first element of the tile (curiously) the entire tile will contain NaNs when decompressed.
Expected behavior
The FITS tile compression standard says NaNs should be handled using the
ZBLANKcolumn/header keyword:Actual behavior
However, when writing an array that contains NaNs ZBLANK is not used for them at all.
Furthermore, there is a separate but related bug where if a NaN occurs anywhere in the tile (except the very first element) the
ZSCALEfor the tile is also set to NaN, causing every pixel in the tile to be "scaled" to NaN when decompressing.The way NaN handling works in CFITSIO is a little strange to begin with, and there are probably some bugs on the Astropy side of not pushing the right buttons to get it to do this properly.
Steps to Reproduce
Here's a simple example: