Compression settings ignored when creating empty CompImageHDU

Description

When creating a HDU without specifying any data, the quantize_level and quantize_method arguments (at least) passed to the HDU class during initialisation are not honoured when data is later set via the HDU's data attribute. Initialising the HDU with even data=np.array([0.]) is enough to (seemingly) solve this issue.

Expected behavior

I would expect the compression settings passed as keyword parameters during initialisation to be remembered in the HDU object and then taken into account when the HDU's data is later updated by the user.

Actual behavior

It seems that if no data is passed during initialisation, the compression settings are assumed to be irrelevant and ignored. I think this should at least be documented (if not already).

Steps to Reproduce

import numpy as np
from astropy.io import fits

tmpfile = 'test_fits.fits'


def diff(hdu):
    # Generate the data to write to the FITS file
    np.random.seed(0)
    data = np.random.rand(50, 30).astype('float32')

    # Set the actual data
    hdu.data = data.copy()

    # Write compressed data to file
    hdul = fits.HDUList([fits.PrimaryHDU(), hdu])
    hdul.writeto(tmpfile, overwrite=True)

    # Print max abs. difference due to compression
    hdul = fits.open(tmpfile)
    print(np.abs(hdul[1].data - data).max())


diff(fits.CompImageHDU(quantize_level=-32))
diff(fits.CompImageHDU(quantize_level=-32, data=np.array([0.])))

kwargs = {'quantize_method': 2, 'dither_seed': 1}
diff(fits.CompImageHDU(**kwargs))
diff(fits.CompImageHDU(**kwargs, data=np.array([0.])))

Actual Values Returned

0.014079988
0.99980855
0.014079988
0.01334089

Expected Values Returned

i.e. initialising with data=None or data=np.array([0.]) should ultimately produce identical compressed values.

System Details

macOS-11.6-x86_64-i386-64bit
Python 3.9.6 | packaged by conda-forge | (default, Jul 11 2021, 03:36:15) 
[Clang 11.1.0 ]
Numpy 1.21.1
astropy 5.0.dev703+g155f428e0.d20210903
Scipy 1.7.1
Matplotlib 3.4.3

@embray

Welcome to Astropy 👋 and thank you for your first issue!

A project member will respond to you as soon as possible; in the meantime, please double-check the guidelines for submitting issues and make sure you've provided the requested details.

GitHub issues in the Astropy repository are used to track bug reports and feature requests; If your issue poses a question about how to use Astropy, please instead raise your question in the Astropy Discourse user forum and close this issue.

If you feel that this issue has not been responded to in a timely manner, please leave a comment mentioning our software support engineer @embray, or send a message directly to the development mailing list. If the issue is urgent or sensitive in nature (e.g., a security vulnerability) please send an e-mail directly to the private e-mail [email protected].

This issue may be related to the _update_header_data method not being called when the data attribute is updated. #12058 may be a partial solution to this, however it may also be necessary to save the keyword parameters passed to __init__ and also pass them to _update_header_data every time the data setter method is called (i.e. data is updated).

I'm trying to dig into this. I'll start with a summary of my initial findings

This issue may be related to the _update_header_data method not being called when the data attribute is updated.

I came to the same conclusion. To give a little bit more details:

_update_header_data is indeed called when self.data is initialised:

astropy/astropy/io/fits/hdu/compressed/compressed.py

Lines 321 to 339 in b299eb5

    
           self.data = data 
        
           # Update the table header (_header) to the compressed 
        
           # image format and to match the input data (if any); 
        
           # Create the image header (_image_header) from the input 
        
           # image header (if any) and ensure it matches the input 
        
           # data; Create the initially empty table data array to 
        
           # hold the compressed data. 
        
           self._update_header_data( 
        
               header, 
        
               name, 
        
               compression_type=compression_type, 
        
               tile_shape=tile_shape, 
        
               hcomp_scale=hcomp_scale, 
        
               hcomp_smooth=hcomp_smooth, 
        
               quantize_level=quantize_level, 
        
               quantize_method=quantize_method, 
        
               dither_seed=dither_seed, 
        
           )

Within that function the header gets initialised

astropy/astropy/io/fits/hdu/compressed/compressed.py

Lines 519 to 534 in b299eb5

    
           self._header, self.columns = _image_header_to_bintable_header_and_coldefs( 
        
               image_header, 
        
               self._image_header, 
        
               self._header, 
        
               name=name, 
        
               huge_hdu=huge_hdu, 
        
               compression_type=compression_type, 
        
               tile_shape=tile_shape, 
        
               hcomp_scale=hcomp_scale, 
        
               hcomp_smooth=hcomp_smooth, 
        
               quantize_level=quantize_level, 
        
               quantize_method=quantize_method, 
        
               dither_seed=dither_seed, 
        
               axes=self._axes, 
        
               generate_dither_seed=self._generate_dither_seed, 
        
           )

and it's apparent that this operation depends on the state of self.data:

astropy/astropy/io/fits/hdu/compressed/compressed.py

Lines 496 to 497 in b299eb5

    
           image_hdu = ImageHDU(data=self.data, header=self._header) 
        
           self._image_header = CompImageHeader(self._header, image_hdu.header)

Now, the important part is that self._image_header["BITPIX"] is set to 8 ("Character or unsigned binary integer") by default (i.e., if self.data is None), but it's set to -32 (resp. -64) if self.data is a numpy array with dtype float32 (resp. float64). It so happens that BITPIX < 0 (which I read as a proxy for "data is IEEE floating point") is a necessary condition in _image_header_to_bintable_header_and_coldefs for quantize_level to have any effect:

astropy/astropy/io/fits/hdu/compressed/header.py

Line 820 in b299eb5

if generated_image_header["BITPIX"] < 0: # floating point image

Indeed, ensuring that the generated header stays in sync with self.data by calling _update_header_data (combined with a way to memorise user-provided options) in the property setter seems like a potential solution, however it gets more complicated and this approach doesn't work out of the box because self.data = ... is actually called before self.header is in a usable state, so it results in an exception. For reproducibility, I'll put here the exact patch I used for this experiment, the test script, and the traceback:

Naive approach

diff --git a/astropy/io/fits/hdu/compressed/compressed.py b/astropy/io/fits/hdu/compressed/compressed.py
index a8db2a5d40..9d948176f3 100644
--- a/astropy/io/fits/hdu/compressed/compressed.py
+++ b/astropy/io/fits/hdu/compressed/compressed.py
@@ -593,6 +593,7 @@ class CompImageHDU(BinTableHDU):
                     type(data), data.dtype.fields
                 )
             )
+        self._update_header_data(self.header)
 
     @lazyproperty
     def compressed_data(self):

from astropy.io import fits
hdu = fits.CompImageHDU()

Traceback (most recent call last):
  File "/Users/clm/dev/astropy-project/coordinated/astropy/bugs/12216/g.py", line 2, in <module>
    hdu = fits.CompImageHDU()
          ^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/utils/decorators.py", line 604, in wrapper
    return function(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/hdu/compressed/compressed.py", line 321, in __init__
    self.data = data
    ^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/utils/decorators.py", line 848, in __set__
    ret = self.fset(obj, val)
          ^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/hdu/compressed/compressed.py", line 596, in data
    self._update_header_data(self.header)
                             ^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/utils/decorators.py", line 837, in __get__
    val = self.fget(obj)
          ^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/hdu/compressed/compressed.py", line 657, in header
    self._image_header = _bintable_header_to_image_header(self._header)
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/hdu/compressed/header.py", line 400, in _bintable_header_to_image_header
    "BITPIX", bintable_header["ZBITPIX"], hcomments["ZBITPIX"], before=1
              ~~~~~~~~~~~~~~~^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/header.py", line 170, in __getitem__
    card = self._cards[self._cardindex(key)]
                       ^^^^^^^^^^^^^^^^^^^^
  File "/Users/clm/dev/astropy-project/coordinated/astropy/astropy/io/fits/header.py", line 1726, in _cardindex
    raise KeyError(f"Keyword {keyword!r} not found.")
KeyError: "Keyword 'ZBITPIX' not found."

It's still possible that the proposed approach is viable and just requires more work, but I think there might be other avenues worth considering

decouple this class' initialisation process (easier said than done)
make CompressedImageHDU.data immutable (this wouldn't be as simple as removing the property setter). I think this would require a deprecation period, at the end of which data would become a required parameter of CompressedImageHDU.__init__.
add an optional dtype argument to CompressedImageHDU.__init__, that'd be required if and only if data isn't provided. I anticipate this would be the simplest implementation but it would increase the class complexity instead of decreasing it, which doesn't meet our goals for maintainability.

I think the next step is to actually try to follow through with the initial approach, just to get a proof-of-concept working example that we can use to judge wether it's actually viable or not. I'll give it a try.

Actually, taking a slightly different approach, I think the issue can be fixed by just memorizing compression options and re-injecting them on each header update. I'll open a PR and see if it gets traction.

TST: add regression tests for issue astropy#12216

Merge pull request #16019 from neutrinoceros/io/fits/bug/remember_com…

TST: add regression tests for issue astropy#12216

Navigation Menu

Repository navigation

Compression settings ignored when creating empty `CompImageHDU` #12216

Description

Expected behavior

Actual behavior

Steps to Reproduce

Actual Values Returned

Expected Values Returned

System Details

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

	self.data = data

	# Update the table header (_header) to the compressed
	# image format and to match the input data (if any);
	# Create the image header (_image_header) from the input
	# image header (if any) and ensure it matches the input
	# data; Create the initially empty table data array to
	# hold the compressed data.
	self._update_header_data(
	header,
	name,
	compression_type=compression_type,
	tile_shape=tile_shape,
	hcomp_scale=hcomp_scale,
	hcomp_smooth=hcomp_smooth,
	quantize_level=quantize_level,
	quantize_method=quantize_method,
	dither_seed=dither_seed,
	)

	self._header, self.columns = _image_header_to_bintable_header_and_coldefs(
	image_header,
	self._image_header,
	self._header,
	name=name,
	huge_hdu=huge_hdu,
	compression_type=compression_type,
	tile_shape=tile_shape,
	hcomp_scale=hcomp_scale,
	hcomp_smooth=hcomp_smooth,
	quantize_level=quantize_level,
	quantize_method=quantize_method,
	dither_seed=dither_seed,
	axes=self._axes,
	generate_dither_seed=self._generate_dither_seed,
	)

	image_hdu = ImageHDU(data=self.data, header=self._header)
	self._image_header = CompImageHeader(self._header, image_hdu.header)

Compression settings ignored when creating empty CompImageHDU #12216

Description

Description

Expected behavior

Actual behavior

Steps to Reproduce

Actual Values Returned

Expected Values Returned

System Details

Activity

github-actions commented on Sep 29, 2021

ConorMacBride commented on Sep 30, 2021

neutrinoceros commented on Feb 9, 2024

neutrinoceros commented on Feb 9, 2024

neutrinoceros commented on Feb 9, 2024

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

Compression settings ignored when creating empty `CompImageHDU` #12216