Making RGB images from FITS files with python/matplotlib.

by Jessica Lu on October 22, 2010

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.

{ 5 comments… read them below or add one }

1 John Gizis October 22, 2010 at 4:18 pm

This is really helpful! Thanks.

Reply

2 Jonathan Whitmore January 11, 2011 at 7:23 pm

The APLpy documentation link is broken in the entry, you mean to link to: http://aplpy.github.com/documentation/index.html

Reply

3 Paul Anthony Wilson February 15, 2011 at 5:45 pm

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.

Reply

4 Brandon Doyle August 12, 2014 at 11:38 am

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…

Reply

5 nicolás February 20, 2015 at 4:27 pm

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

Reply

Leave a Comment

Previous post:

Next post: