Converting an RGB image to fits (astropy)

I am trying to save a NumPy RGB array (a 3d array with dimension 512x512x3), to a .fits file using astropy.io. Saving this array to a .png file works fine, but when I'm trying to save it to a .fits file using HDUList.writeto, it produces just garbage colors. Also, I have to reshape the image like so:

img = img.reshape(img.shape[2], img.shape[0], img.shape[1])


or it will try to store just 3 column of the image instead of 512.

Each element in the 3d array contains a value between 0 and 255.

• What does "it produces just garbage colors" mean? FITS files aren't meant for storing RGB "colors". – Peter Erwin Jun 10 '19 at 13:19

I haven't done much astronomical image processing before, but as this question is unanswered I'll give it a shot - hopefully to some avail. If the problem is more specific, a code sample/image sample would probably be useful for further diagnosis, but otherwise this example may help. It discusses the process of writing a 3-channel image to separate FITS images. I would try this first, and check if this output is also garbage - it might just be an issue with how you're using the HDUList.writeto() function. The relevant code is pasted below (with a couple of edits, as the example uses the Pillow Image class at first):

Split the three channels (RGB):

r, g, b = img[:, :, 0], img[:, :, 1], img[:, :, 2]


Write out the channels as separate FITS images

red = fits.PrimaryHDU(data=r)
red.writeto('red.fits')

green = fits.PrimaryHDU(data=g)
green.writeto('green.fits')

blue = fits.PrimaryHDU(data=b)
blue.writeto('blue.fits')


Also out of curiosity - why is your "numpy RGB array" storing values between 0 and 512? The standard for images is usually integers between 0-255 or floats between 0 and 1; but my guess is that astronomy requires a larger bit depth? In any case, you should also make sure the dtype of your numpy array is large enough for the desired bit depth - often images are stored with dtype="uint8" which ranges from 0 to 255.

• Hello, thanks for your reply. I tried your splitting the r,g,b channels and saving a .fits file for each channel like you described. It works, but the observed intensities are clearly quite different than the .png file. It makes me think I need to do some kind of transformation to the values before storing them to give a good result, but I don't really know what that is. Also, you're right about the 0-255 range, that was my mistake. – user3042966 Jun 10 '19 at 16:58
• I uploaded a sample png with a .fits for the r,g,b channels here: link – user3042966 Jun 10 '19 at 17:09
• @user3042966 What are you observing the intensities with? Astronomical software such as ds9 may just be setting the colorscale to log rather than linear values, giving the impression that the intensities are off. – usernumber Jan 7 at 9:30