Making RGB images from FITS files with python/matplotlib.

Converting astronomical data taken in multiple filters into representative-color RGB images often provides one of the most visually appealing (and informative) views of a target. This can be done in ds9 very easily; however, if you want a little more control, then python/matplotlib can be used in a very similar fashion. I also note that APLpy has excellent capabilities for images with proper World Coordinate System headers already (see the APLpy RGB tutorial). APLpy is far easier to use, but if you don’t have that choice, then you can use matplotlib as I am going to show below.

First, lets assume we are starting with three images, say J.fits, H.fits, and K.fits and that these images all have the same plate scale, the same angle, and the same position on the sky but lacks WCS information (as is often the case for ground-based OIR observations). Load up the images into numpy arrays using the pyfits module.

import pyfits
import numpy as np
import pylab as py
import img_scale
j_img = pyfits.getdata('J.fits')
h_img = pyfits.getdata('H.fits')
k_img = pyfits.getdata('K.fits')

You can use pylab.imshow on a NxMx3 array; however, each layer of the image must first be mapped onto the luminosity scale. This is done with a module called that I borrowed from Min-Su Shin at U. of Michigan who very generously freely distributed his code (thanks!). You can get from this AstroBetter wiki page. Using this module you can linear, log, sqrt, or asinh stretch each color channel independently and set minimum or maximum flux values.

img = np.zeros((j_img.shape[0], j_img.shape[1], 3), dtype=float)
img[:,:,0] = img_scale.sqrt(k_img, scale_min=0, scale_max=10000)
img[:,:,1] = img_scale.sqrt(h_img, scale_min=0, scale_max=10000)
img[:,:,2] = img_scale.sqrt(j_img, scale_min=0, scale_max=10000)

Finally, we can display the image:

py.imshow(img, aspect='equal')
py.title('Blue = J, Green = H, Red = K')

If your input images don’t have the same rotation angles, plate scales, or positions, then IRAF’s handy imlintran can be used to rescale, rotate, and translate. This can still be done in python via pyraf. Here is an example of calling imlintran from within python for a J.fits image that I want to rotate, shift, and rescale to match H and K images. I will assume the H and K images have a position angle of zero.

from pyraf import iraf as ir
h_img = pyfits.getdata('H.fits')
k_img = pyfits.getdata('K.fits')
scaleK = 0.1   # arcsec per pixel
scaleJ = 0.15 # arcsec per pixel)
paJ = 30.0  # degrees
ir.imlintran.boundary = 'constant'
ir.imlintran.constant = 0
ir.imlintran.interpolant = 'spline3'
ir.imlintran.fluxconserve = 'yes'
# the input pixel position of a reference source = 0
ir.imlintran.yin = 0
# the output pixel position of a reference source
ir.imlintran.xout = 10
ir.imlintran.yout = 10
# Set the output size of the final image
ir.imlintran.ncols = k_img.shape[1]
ir.imlintran.nlines = k_img.shape[0]
ir.imlintran('j.fits', 'j_rot_scale.fits', angleJ, angleJ, scaleK/scaleJ, scaleK/scaleJ)
j_img = pyfits.getdata('j_rot_scale.fits')

Alternatively, you can use scipy.ndimage.interpolate to shift images.

from scipy.ndimage import interpolation as interp
old_j = pyfits.getdata('j.fits')
new_j = interp.shift(old_j, [shiftY, shiftX])

The image attached in this post was plotted using imshow as described above and is a JHK color composite from the UKIDSS galactic plane survey.

6 comments… add one
  • John Gizis Oct 22, 2010 @ 16:18

    This is really helpful! Thanks.

  • Jonathan Whitmore Jan 11, 2011 @ 19:23

    The APLpy documentation link is broken in the entry, you mean to link to:

  • Paul Anthony Wilson Feb 15, 2011 @ 17:45

    This is really cool. I would also urge people to try FITS Liberator which is provided by ESA/ESO/NASA together with Photoshop (for those people who can afford it, or get it through other means). I like the method above as it is open source.

  • Brandon Doyle Aug 12, 2014 @ 11:38

    Out of curiosity, how would one create an RGB if the .fits data is from different telescopes, of different resolutions, alignments etc.? I’m trying to use aplpy but it doesn’t seem to want to take my HST and VLA data…

  • nicolás Feb 20, 2015 @ 16:27

    When i put import img_scale in my script, the code doesn’t work, you hace an idea why?

  • Lana Oct 28, 2020 @ 12:04

    Thanks for this. Really helped a lot.

Leave a Reply

Your email address will not be published. Required fields are marked *