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 img_scale.py that I borrowed from Min-Su Shin at U. of Michigan who very generously freely distributed his code (thanks!). You can get img_scale.py 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.clf()
py.imshow(img, aspect='equal')
py.title('Blue = J, Green = H, Red = K')
py.savefig('my_rgb_image.png')

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.unlearn('imlintran')
 
ir.imlintran.boundary = 'constant'
ir.imlintran.constant = 0
ir.imlintran.interpolant = 'spline3'
ir.imlintran.fluxconserve = 'yes'
 
# the input pixel position of a reference source
ir.imlintran.xin = 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.