Repository navigation
Pyfits: Incorrect writing of non contiguous data in compressed FITS #2150
Description
Activity
That makes sense--CFITSIO is going to assume the data is stored in C-order and operate on it as such. Since there's no way (I didn't think) to change how CFITSIO is going to work on the data it would have to be transposed before going through the compression routines. This of course introduces some extra overhead but is not insurmountable. Still, better to just use C-order in the first place if possible.
Well, the user should not have to worry about his data being strided or not when passing them to pyfits. The ImageHDU class works flawlessly in that regard. You might consider passing np.ascontiguousarray(the_input_image) to CFITSIO. It will make a copy only if the data is not C contiguous.
Oh I agree completely. I was just wondering if in the meantime any trouble caused by this could be avoided.
I wonder whether it would be better to make a copy or to just transpose the array in-place. It's a tradeoff...
Modifying the data in-place seems like an undesirable side-effect if the user needs the data in fortran-order, no? It seems that a copy would make sense, maybe with a warning that this is inefficient?
The issue is more general than that: a non C-contiguous array is not necessarily F-contiguous. Writing a[::2] would fail too I'm pretty sure. Pyfits is already quite verbose, I don't think a warning would be good. There are hundreds of places in scipy where the data is made contiguous in some ordering, without warnings. Pyfits/CFITSIO behaviour should be documented though. And ImageHDU doesn't issue a warning. Doesn't it make a copy for non C-contiguous arrays ?
No. In the example you gave they're really two very different arrays:
>>> a = np.arange(100, dtype='uint8').reshape((10, 10)) >>> a array([[ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9], [10, 11, 12, 13, 14, 15, 16, 17, 18, 19], [20, 21, 22, 23, 24, 25, 26, 27, 28, 29], [30, 31, 32, 33, 34, 35, 36, 37, 38, 39], [40, 41, 42, 43, 44, 45, 46, 47, 48, 49], [50, 51, 52, 53, 54, 55, 56, 57, 58, 59], [60, 61, 62, 63, 64, 65, 66, 67, 68, 69], [70, 71, 72, 73, 74, 75, 76, 77, 78, 79], [80, 81, 82, 83, 84, 85, 86, 87, 88, 89], [90, 91, 92, 93, 94, 95, 96, 97, 98, 99]], dtype=uint8) >>> a.tostring() '\x00\x01\x02\x03\x04\x05\x06\x07\x08\t\n\x0b\x0c\r\x0e\x0f\x10\x11\x12\x13\x14\x15\x16\x17\x18\x19\x1a\x1b\x1c\x1d\x1e\x1f !"#$%&\'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_`abc' >>> b = np.arange(100, dtype='uint8').reshape((10, 10), order='F') >>> b array([[ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90], [ 1, 11, 21, 31, 41, 51, 61, 71, 81, 91], [ 2, 12, 22, 32, 42, 52, 62, 72, 82, 92], [ 3, 13, 23, 33, 43, 53, 63, 73, 83, 93], [ 4, 14, 24, 34, 44, 54, 64, 74, 84, 94], [ 5, 15, 25, 35, 45, 55, 65, 75, 85, 95], [ 6, 16, 26, 36, 46, 56, 66, 76, 86, 96], [ 7, 17, 27, 37, 47, 57, 67, 77, 87, 97], [ 8, 18, 28, 38, 48, 58, 68, 78, 88, 98], [ 9, 19, 29, 39, 49, 59, 69, 79, 89, 99]], dtype=uint8) >>> b.tostring() '\x00\n\x14\x1e(2<FPZ\x01\x0b\x15\x1f)3=GQ[\x02\x0c\x16 *4>HR\\\x03\r\x17!+5?IS]\x04\x0e\x18",6@JT^\x05\x0f\x19#-7AKU_\x06\x10\x1a$.8BLV`\x07\x11\x1b%/9CMWa\x08\x12\x1c&0:DNXb\t\x13\x1d\'1;EOYc'At a certain level PyFITS doesn't care what the shape or index order of the array is. It's just a raw sequence of bytes. This is where things run into trouble with the compression code in particular, because once you've written to the file and read it back again it has no way of knowing that the array that was passed in was "Fortran order". That information is not in any way preserved.
@astrofrog In place transposition of the data would only happen before writing, and then restored after. We already do this with regard to byte order. Unfortunately transposing a large array is a lot more expensive that byte swapping (which is about as fast as it can be on modern architectures).
@embray not sure what you're answering no, but the fix you're proposing is not general enough. As I anticipated:
>>> orig = np.arange(10, dtype=float)[::2] >>> primary = fits.PrimaryHDU() >>> hdu = fits.CompImageHDU(orig) >>> hdulist = fits.HDUList([primary, hdu]) >>> hdulist.writeto('test_strided.fits', clobber=True) >>> print 'orig:', orig orig: [ 0. 2. 4. 6. 8.] >>> print 'actual:', fits.open('test_strided.fits')[1].data actual: [ 0. 1. 2. 3. 4.]What matters is whether or not the data is contiguous in memory or not.
Ah I see, good point. In that case it would be absolutely necessary to copy the data into a contiguous buffer. No way around it there :/ At least in the C-contiguous vs. F-contiguous case the issue could be fixed without copying by transposing the data in-place. But that wouldn't work when it comes to arbitrary strides (well, in principle they could be handled by shuffling data around in-place as well, but at that point a copy is much simpler).
4 remaining items
We ran into this issue with Gammapy, see here:
gammapy/gammapy#1346[updated test case and output]
Minimal test case to reproduce:
import numpy as np from astropy.io import fits data = np.arange(6).reshape((3, 2)).T hdu = fits.PrimaryHDU(data) filename = 'tmp.fits.gz' hdu.writeto(filename, overwrite=True) hdu2 = fits.open(filename)[0] data2 = hdu2.data print('data:\n', data) print('data2:\n', data2)gives the incorrect result (i.e. data is different after write / read):
data: [[0 2 4] [1 3 5]] data2: [[0 1 2] [3 4 5]]When removing the
.gz, i.e. writing just totmp.fits, the data comes back OK.This is very confusing for users, either an error, or making it work via a copy would be preferable.
I think I'm blind. I can't see any difference between those two. Could you clarify?
@MSeifert04 - I had copy & pasted the wrong test case and output. I've updated my last comment just now, please refresh.
Ah thanks!
That's not good ... I will look into this (although probably not this week) but if it's indeed effort-high (label) raising an exception might make sense. :)
Reading the discussion above,
np.ascontiguousarrayseems the easiest way to go for now, not the most efficient but at least it would fix this usecase.- added 2 commits that reference this issue
on Feb 19, 2020
Using astropy 0.3.dev4315, the following code prints False instead of True. If ImageHDU is substituted to CompImageHDU, it works correctly and prints True.