fits.PrimaryHDU doesn't properly initialize from scaled data - #2305

I'm not sure I fully follow your example, as in what you were trying to do originally (what the inputs and expected outputs were) and what you got instead. All of the behavior shown in the example is by design, though admittedly potentially confusing.

Also worth noting that what you're really trying to do here is something involving unsigned ints--at this point by fair the most used usecase for BZERO. I need to work more on making unsigned int support just more transparent because it's really absurd that it isn't.

By the way, perhaps what you're missing is the uint=True option needed to enable consistent support of pseudo unsigned ints. The main change I need to make is just to make this the default.

ETA: Ah, okay. One other issue here--if creating a new HDU object with unsigned integer data it should probably go ahead and put the appropriate BZERO value in the header upon initialization. Is that what you're getting at here?

if creating a new HDU object with unsigned integer data it should probably go ahead and put the appropriate BZERO value in the header upon initialization. Is that what you're getting at here?

Yes, sorry that wasn't clear. That is the only problem I was trying to report. More detail:

Given a data array and a header that contains header['BZERO']=32768 I would expect that

new_hdu = fits.PrimaryHDU(data, header=header)

should return new_hdu with scaled data and no BZERO in the header.

If I ask PrimaryHDU not to scale the data with

new_hdu = fits.PrimaryHDU(data, header=header, do_not_scale_image_data=True)

I would expect that new_hdu.data == data and that new_hdu.header['BZERO'] is 32768.

I should never get unscaled data with no 'BZERO'.

I'm not sure I fully follow your example, as in what you were trying to do originally

I was trying to do something that, in retrospect, was nonsensical :-/ : allowing the user to ask that data not be scaled when initializing a ccdproc.CCDData object from a FITS file by passing do_not_scale_data to fits.open.

In my defense, the immediate goal was allowing keywords to pass through the CCDData reader to fits.open...then I checked to see how do_not_scale_image_data was handled and found this. Now the CCDData reader raises an error if the user asks to open without scaling.

I think that all makes sense now. Agreed this specific case should be handled better.

Still need to give this some thought, and some more explanation of what you're doing might be useful.

The problem with the current code is this: When you instantiate a PrimaryHDU (or really just about any other HDU type, but we'll stick with the one example) with some data array and some header, pyfits has no idea where that data and where that header came from. Even if they were already associated with each other in some existing HDU, this form of instantiation has no knowledge of that association, or if the header is even "valid" for that particular data.

Because of this, the header argument to the initializer is really meant only for "metadata" keywords that are not required to describe the physical format of the data itself. In fact, the header you pass in has all such "structural" keywords stripped out of it, and then is rebuilt into a guaranteed valid (for a Primary HDU) header, including the correct NAXIS/BITPIX keywords for the data array that was passed in, and then followed by any non-"structural" keywords that were in the original header passed in.

Unfortunately, what I'm referring to as "structural" keywords also includes BSCALE and BZERO. If given an array and a header with BSCALE/BZERO, where their association is unknown, it's often completely ambiguous if and how those scale keywords should be applied to the input data array. Maybe it's already been scaled according to those parameters. Maybe it hasn't. Maybe the input header came from a completely different HDU and the user just wanted to copy the metadata, but BSCALE/BZERO just came along for the ride unintentionally. In other words, the header argument should not have any effect on how the data argument is interpreted. I think that was the original intent in the current design.

Now there is one small set of cases where it's less ambiguous, and that's the unsigned int case. But say we're given an array of uint16--if all the values in the array are < 32768, then the array can just as easily be stored in FITS as a normal int16 array. It's only when given an array of uint16 containing values >= 32768 that the pseudo-unsigned int workaround is necessary. And I believe that case should just work transparently. Users should be able to store/retrieved unsigned integer arrays from a FITS file without having to fuss around with scaling or headers or anything like that. I'm not sure if that would help with your problem though...

Okay, here's an example that I think illustrates what can and should be fixed for now (before I go about completely overhauling uint support). I think this is almost identical to the example you gave, but going through it myself helped understand:

In [1]: from astropy.io import fits

In [2]: import numpy as np

In [3]: uints = np.uint16(np.random.random_integers(0, 65535, size=[100, 100]))

In [4]: uints_hdu = fits.PrimaryHDU(data=uints)

# Initializing an HDU with uint data happily keeps it in that format
In [5]: uints_hdu.data
Out[5]: 
array([[53308, 47610, 52202, ..., 25902,  2372,  6734],
       [ 2123,  9195, 50909, ..., 63669, 26111, 28117],
       [29503, 30208, 28154, ..., 59558,  2159, 36058],
       ..., 
       [40962,  5128,  6651, ..., 55110, 31358, 65326],
       [ 9323,  3586, 10685, ...,  5374, 22279,  8573],
       [ 3428,  9785, 51038, ..., 61065, 30025,  3280]], dtype=uint16)

# But it does not put BZERO in the header
In [6]: uints_hdu.header
Out[6]: 
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   16 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  100                                                  
NAXIS2  =                  100                                                  
EXTEND  =                    T                                                  

In [7]: uints_hdu.writeto('uints.fits')

# Writing the HDU out to a file doesn't change the header either
In [8]: uints_hdu.header
Out[8]: 
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   16 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  100                                                  
NAXIS2  =                  100                                                  
EXTEND  =                    T                                                  

In [9]: h = fits.open('uints.fits', uint=True)

# But the written file *does* preserve the uint scaling (uint=True is required to return integers, otherwise it converts to floats)
In [10]: h[0].data
Out[10]: 
array([[53308, 47610, 52202, ..., 25902,  2372,  6734],
       [ 2123,  9195, 50909, ..., 63669, 26111, 28117],
       [29503, 30208, 28154, ..., 59558,  2159, 36058],
       ..., 
       [40962,  5128,  6651, ..., 55110, 31358, 65326],
       [ 9323,  3586, 10685, ...,  5374, 22279,  8573],
       [ 3428,  9785, 51038, ..., 61065, 30025,  3280]], dtype=uint16)

# But the header for this HDU *does* have BZERO in it
In [11]: h[0].header
Out[11]: 
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   16 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  100                                                  
NAXIS2  =                  100                                                  
EXTEND  =                    T                                                  
BSCALE  =                    1                                                  
BZERO   =                32768 

So in short, the initial uints_hdu I created should have also had BZERO in its header.

I agree that your examples illustrates the essential problem: the header can end up in a state that, applied to the data, results in incorrect data values.

Thanks for the summary of how the header initialization is handled -- useful to know for ccdproc.

@mwcraig If you'd like, would you please try out the attached change and see how that works for your needs? It's a one-line fix, hopefully.

@embray -- I don't think this fixes...and it appears to break in some cases the current (correct) behavior for unsigned int images; I'm afraid I may be confusing the issue!

I'll try to clarify the issue below, but want to emphasize first that the original use case I had for this issue is gone, so there is no urgency to getting a fix on my end.

In case this is tl;dr -- would it be helpful if I wrote a couple tests that demonstrate what I see as the problems?

Looking over your example again, I think it demonstrates a separate issue than the one I originally reported.

Assume that you have already created uint.fits as in your example. Then suppose we open that file without scaling data so we can see what the stored data/header looks like (there is no error in what is immediately below):

In [15]: h = fits.open('uints.fits', do_not_scale_image_data=True)

In [16]: h[0].header  # Should contain BZERO bc we have not yet scaled the data
Out[16]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T
BSCALE  =                    1
BZERO   =                32768

In [17]: h[0].data  # Should be signed int16 because the data has not been scaled
Out[17]:
array([[ 30716,  26561, -14519, ...,   3678, -19177,    863],
       [ 20054,  22152,   2374, ...,  22282,   5299, -18784],
       [ 19872,   8897,  -8137, ...,   3089,  29558,  29516],
       ...,
       [  5971, -23966,  22772, ...,  28757,  32167, -27685],
       [-32188, -29295, -21282, ...,  20658,  18497,  14584],
       [-29243,  -8268,  17596, ...,   5903,  -9962,   4725]], dtype=int16)

In [18]: h[0].header # BZERO still here because I explicitly asked that the data not be scaled
Out[18]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T
BSCALE  =                    1
BZERO   =                32768

So far so good -- the header and data are consistent in the sense that to get the physical data values from the stored data values I use BZERO/BSCALE from the header, physical = BSCALE * stored + BZERO.

If you open the same file but allow the data to be scaled this happens (there is no error in what is immediately below):

In [18]: h = fits.open('uints.fits')

In [19]: h[0].header  # BZERO will be in here because scaling happens only when data is accessed
Out[19]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T
BSCALE  =                    1
BZERO   =                32768

In [20]: h[0].data  # accessing the data automatically uses the header BZERO/BSCALE to scale the data correctly
Out[20]:
array([[ 63484.,  59329.,  18249., ...,  36446.,  13591.,  33631.],
       [ 52822.,  54920.,  35142., ...,  55050.,  38067.,  13984.],
       [ 52640.,  41665.,  24631., ...,  35857.,  62326.,  62284.],
       ...,
       [ 38739.,   8802.,  55540., ...,  61525.,  64935.,   5083.],
       [   580.,   3473.,  11486., ...,  53426.,  51265.,  47352.],
       [  3525.,  24500.,  50364., ...,  38671.,  22806.,  37493.]], dtype=float32)

In [21]: h[0].header  #  BZERO/BSCALE are correctly automatically removed from the header because they have been used to scale the data
Out[21]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -32 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T

Now to the case in your example, opening the data file with the uint keyword:

In [30]: h = fits.open('uints.fits', uint=True)

In [31]: h[0].header # BZERO is in the header, which is expected because the data hasn't been accessed
Out[31]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T
BSCALE  =                    1
BZERO   =                32768

In [32]: h[0].data  # Data is scaled to correct (unsigned) values
Out[32]:
array([[63484, 59329, 18249, ..., 36446, 13591, 33631],
       [52822, 54920, 35142, ..., 55050, 38067, 13984],
       [52640, 41665, 24631, ..., 35857, 62326, 62284],
       ...,
       [38739,  8802, 55540, ..., 61525, 64935,  5083],
       [  580,  3473, 11486, ..., 53426, 51265, 47352],
       [ 3525, 24500, 50364, ..., 38671, 22806, 37493]], dtype=uint16)

In [33]: h[0].header #  ERROR HERE -- BZERO is still in the header even though the stored data values ARE the physical values
Out[33]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100
EXTEND  =                    T
BSCALE  =                    1
BZERO   =                32768

I think this is an error because it is no longer the case that physical = BSCALE * stored + BZERO.

So in short, the initial uints_hdu I created should have also had BZERO in its header.

I think the answer here should have been no...my bad. Doing so leads to the same error as opening with uint=True.

I'm not sure it is the same error that I reported initially, though. In that case my expectation was that the handling of data scaling when creating an HDU with PrimaryHDU would mirror that of fits.open. In the block below, I read uint.fits without scaling to get a header and stored values that I can use as arguments to PrimaryHDU.

In [35]: h = fits.open('uints.fits', do_not_scale_image_data=True)

In [36]: new_hdu = fits.PrimaryHDU(data=h[0].data, header=h[0].header)

In [37]: new_hdu.data  # Expect this to have been scaled because that is what fits.open would do with the same header/data
Out[37]:
array([[ 30716,  26561, -14519, ...,   3678, -19177,    863],
       [ 20054,  22152,   2374, ...,  22282,   5299, -18784],
       [ 19872,   8897,  -8137, ...,   3089,  29558,  29516],
       ...,
       [  5971, -23966,  22772, ...,  28757,  32167, -27685],
       [-32188, -29295, -21282, ...,  20658,  18497,  14584],
       [-29243,  -8268,  17596, ...,   5903,  -9962,   4725]], dtype=int16)

In [38]: new_hdu.header  # Given output immediately above, expect that BZERO will still be in the header so that I can reconstruct the physical values
Out[38]:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  100

Here the error is that it is no longer possible to reconstruct the physical values because BZERO has been dropped.

Ah, the third case you mentioned I agree is confusingly inconsistent with the others. But I think it's somewhat intentional that it's treated differently with unsigned ints. The reason being if you read the file with uint=True you're going to want to also write the file while keeping the data as unsigned ints. So when the file is written out it's going to add BZERO to the header anyway.

I guess maybe a more consistent approach would be to add BZERO to the header before writing and then remove it again afterward. But having in the header reflects what will be written when the file is written to a file. On the other hand with arbitrary BSCALE/BZERO, those are left out of the header because they will not be used when writing the file (unless scale_back=True was given.

I have some ideas to improve the situation somewhat, but no matter what it has to be a little confusing. The area that can be improved most (and honestly that's most common anyway) is with unsigned ints. I think as it is this already works pretty okay. The only major change is that uint=True needs to be the default.

And yes, do_not_scale_image_data and uint arguments only really apply when reading data from a file. They don't have as much (or any) influence when creating a new HDU from scratch. It might be better if they were handled somehow as entirely internal details.

Oh, and as for your last example, this is exactly as expected. When you call

new_hdu = fits.PrimaryHDU(data=h[0].data, header=h[0].header)

it has no way of knowing that that data was originally associated with that header. The BSCALE/BZERO are stripped from the header and are no longer relevant. I agree this is a problem though.

This would be a good place to mention what I'm planning to do to address this, and some other issues. I think I've mentioned elsewhere too that what is really needed is a FITSData class that keeps as internal attributes all the FITSisms such as scaling and other minor details. That information should all be tied to the data array itself, and not any specific HDU that it belongs to, much less any header the data parameters were originally gleaned from.

Would it make sense, then, to drop the arguments do_not_scale_image_data, etc., from PrimaryHDU?

Possibly, eventually, but only as part of a more significant refactoring.

In the meantime I'm not entirely clear where to go with this. I agree the behavior is confusing, but other than the one tweak I made in this PR it's all by design, and somewhat by necessity. There's no real consistent way to deal with these issues of FITS. In particular, the uint case needs to be treated as distinct from arbitrary "scaled" data. Really BITPIX=16 plus BZERO=32768 should be treated like a single keyword indicating unsigned 16-bit int. And of course the data can't actually be stored with the correct values so it has to be transformed. But it serves a completely different purpose than arbitrarily scaled data.

I would leave it as it is (even without the tweak) -- ultimately, my confusion came from the existence of the scaling arguments in PrimaryHDU. From what you have said it makes sense to treat the most common case -- unsigned int data -- separately from other scaled data, and there are changes to FITS i/o coming in the future, so it may as well wait.

I do regret not suggesting for the A&C FITS limitations paper that FITS cannot even handle simple data types like unsigned int....

Hah, that's a good point. It should have been in there. In fairness, it would be an easy enough thing to fix in FITS. It mystifies me that they haven't (we could add 64-bit ints, but not unsigned ints?) I think I am going to include this change though because it is more consistent with the intended behavior, regardless of how confusing the intended behavior is.

With this passing I'm going to go ahead and merge as is. It's no less confusing but it is more in the direction of where I want to go with this.