FITS images

1. What are FITS images?

FITS stands for Flexible Image Transport System. The complete documentation can be found on this webpage hosted by NASA.

In brief, FITS is a binary format to store astronomical data (either image(s) or certain types of table(s)). The reason why we need it is that our files are usually huge, making the ASCII format impractical. In addition, there usually needs some meta-data to explain how the data were taken and were processed (e.g., the angular resolution or spectral resolution of the image(s) stored in a specific FITS file, the date of the observations, etc), such that the data file can be self-contained without depending on external documentation.

In light of this, a FITS file usually contains two parts:

one or more binary data table
one or more header files that contain the necessary information to understand the data.

Since the binary data are not human readable, we cannot handle the FITS files with a text editor (e.g., vim) or other commercial or open resource software like MS Excel or Open Office. Here I provide some simple tips for manipulating FITS images that can be useful to radio astronomers. If you are a graduate student or researcher, I recommend at least installing the ds9, CASA, Miriad software and the Astropy and APLpy Python packages. If you have access to IDL, you should also install the IDL Astronomy User’s Library.

If you know something that is very useful, please recommend it to me.

2. (Pre)Viewers

The following software may be useful for a quick preview of your FITS images or tables, which are particularly handy during a casual discussion with other colleagues. The links to them are provided. They are all reasonably easy to install.

2.1 ds9

For any astronomer, ds9 would be one of the most handy FITS image viewers. It is trivial to install no matter whether you are using Linux, Windows, or Mac OSX. There is no reason not using it. After installation, you can launch it by double-clicking its icon, otherwise, by typing

> ds9

in a Linux or OSX comand line.

2.2 fv

fv is a simple graphical interface for you to preview FITS tables. It can also export data to ASCII format. It is also very easy to install.

2.3 CASA-Viewer and CARTA

If you are a radio astronomer, it is most likely that you have installed the CASA Software Package. The viewer function of CASA is handy when you are previewing the radio interferometric images, although it crashes once in a while. It may be replaced by CARTA at some point although I have not adapted myself to it. If you are new to astronomy, I recommend beginning with CARTA. When your image file is large (e.g., >100 GB), presently, only CARTA allows you to preview the images reasonably smoothly.

2.4 Karma kvis

We used the karma package a lot in the old days. The kvis function is very useful when you are interactively examining the gas kinematics of an astronomical object with a spectral line image cube. I myself still use it once in a while although not everybody uses it nowadays.

3. Converting images in various formats

3.1 Miriad

After sourcing the Miriad software package, which is often used in the imaging of the SMA data, you can use the following commands to convert between the Miriad and FITS format images:

# The file names in in='' and out='' (string variables) 
# need to be replaced with your filenames

# converting from Miriad format to FITS format
fits in='mytarget.image.miriad' op=xyout out='mytarget.image.fits'

# converting from FITS to Miriad format
fits in='mytarget.image.fits' op=xyin out='mytarget.image.miriad'

They can be executed either by directly entering in a Linux/OSX command line or by executing as a BASH or CSH script. We can then program other software languages to manipulate or analyze the FITS format images, such as C/Fortran, IDL, or Python. For online documentation of the Miriad tasks please check the Miriad task index.

You can regrid two Miriad images or image cubes to have identical spatial or spectral coordinates using the following command:

# The '\' symbols denotes line-breaking
regrid in='image1.image.miriad'  \
       out='image1.regrided.image.miriad' \
       axes=1,2 \
       tin='image2.image.miriad'

This is useful when you want to align the pixels of two images before any further analyses. In this example, axes=1,2 means that we are regridding the two spatial axes. When you want to also regrid the spectral axis, you can use axes=1,2,3. It loads ‘image1.image.miriad’ and regrids it according to the header information in ‘image2.image.miriad’, and then outout the regridded image1 as a new file ‘image1.image.miriad’.

3.2 CASA

After launching CASA, you can type viewer in the CASA command line environment to launch the GUI image previewer.

You can use the following CASA Software Package commands to convert between the CASA and FITS format images:

# The file names (string variables) in imagename and fitsimage need 
# to be replaced with your filenames.

# converting from CASA format to FITS format
exportfits(
           imagename = 'mytarget.image',
           fitsimage = 'mytarget.image.fits',
           overwrite = True
          )
          
# converting from CASA format to FITS format
importfits(
           imagename = 'mytarget.image',
           fitsimage = 'mytarget.image.fits',
           overwrite = True
          )
          
# If your CASA image have redundant spectral or polarization dimension,
# and if you want to drop those redundant dimension, you can include the
# relevant options (dropstokes or dropdeg) like:
exportfits(
           imagename = 'mytarget.image',
           fitsimage = 'mytarget.image.fits',
           overwrite = True, dropdeg = True, dropstokes = True
          )
# This can be very useful since not all Python package can recognize 
# high-dimensional FITS images.

To crop a sub-image in CASA, you can use the CASA imsubimage task like:

# # box[[x1, y1], [x2, y2]] region
subimregion = "box[[1500pix, 1500pix], [4499pix, 4499pix]]"

img_name = 'mytarget'
imsubimage(imagename = img_name + '.image',
           outfile   = img_name + '.subim.image',
           region    = subimregion
          )

We also often use the CASA imregrid and imsmooth tasks to regrid a set of CASA format images to the same coordiante system and then smooth same synthesized beam size before exporting to FITS files for further analyses. Example of these commands are:

imregrid(
         # name of the image to be regrided
         imagename = imagename,
         # an existing image to serve as a template of coordiante grids
         template  = template_imagename,
         # output filename
         output    = imagename + '.regrid',
         # options in the regridding
         axes      = [0, 1], interpolation = 'linear',
         # whether or not overwritting the output file
         overwrite = True
        )
          
imsmooth(
         # name of the image to be smoothed
         imagename = image,
         # Shape of the smoothing kernel or final synthesized beam.
         # This example specified that the final synthesized beam should
         # be a gaussian beam with 0.53 arcsecond FWHM in the major and minor axes
         kernel    = 'gauss',
             major = '0.53arcsec', minor = '0.53arcsec', pa = '0deg',
             targetres = True,
         # output file name
         outfile   = image + '.smooth', 
         # whether or not overwritting the output file
         overwrite = True
        )

These CASA commands can be executed either by directly entering in a Linux/OSX command line or by executing as a Python3 script.

To see a full list of available CASA tasks, please check this webpage. For online documentation of the CASA tasks please check the CASA docs. In the CASA command line, you can also check the documentation and input variables by entering help taskname and inp, e.g., in this case:

# check the documentation for importfits/exportfits
CASA> help importfits
CASA> help exportfits

# load the task importfits
CASA> tget importfits
# check the input variable of the loaded task
CASA> inp

After obtaining the FITS images, we can then program other software languages to manipulate or analyze the FITS format images, such as C/Fortran, IDL, or Python.

4. Working with FITS using software languages

If you are programming C or Fortran, then CFITSIO is the package you need. I do not think this is the way you will be going but am providing a reference to it just in case (e.g., sometimes you may need to check other people’s programs which utilized this package).

4.1 Python

If you are analyzing images, you might build a conda environment by the following Linux/OSX command line commands:

conda create --name astroimgAna python=3.9
pip install --upgrade pip
pip install jupyter
pip install astropy scipy
pip install jupyter matplotlib h5py aplpy pyregion PyAVM healpy
pip install pandas

In 2023, after updating astropy to version >5.0, it may become impossible to install aplpy successfully. If that is what happens to you, you can try the following command:

conda create --name test python=3.7
pip install --upgrade pip
pip install numpy
pip install astropy==4.3.1
pip install aplpy==2.0.3
pip install jupyter
pip install scipy
pip install matplotlib h5py
conda install -c conda-forge pyregion
pip install PyAVM healpy
pip install pandas

I often make most of the following imports:

import os
import numpy as np

from astropy.io.fits import getdata
from astropy import wcs
from astropy.io import fits
from astropy import units as u
from astropy import constants as con
from astropy.coordinates import SkyCoord

import matplotlib
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib

import aplpy

matplotlib.use('PDF')

The following code is an example of importing FITS image as a Header Data Unit (HDU) and extracting certain header information in Python:

intensitymap = 'mytarget_StokesI.fits'
intensity_scale = 1000.0
if_success = False
try:
  
  # importing FITS image to a HDU
  Ihdu   = fits.open(intensitymap)
  
  # editing the FITS image by multiplying a scaling factor
  Ihdu[0].data = Ihdu[0].data * intensity_scale
  
  if_success = True

except:
  print('Unable to read the intensity FITS image. Please double check the image file.')
  print(intensitymap)

if ( if_success == True ):
  # Reading FITS header
  try:
    naxis1 = Ihdu[0].header['naxis1']
    naxis2 = Ihdu[0].header['naxis2']
    crval1 = Ihdu[0].header['crval1']
    crpix1 = Ihdu[0].header['crpix1']
    cdelt1 = Ihdu[0].header['cdelt1']
    crval2 = Ihdu[0].header['crval2']
    crpix2 = Ihdu[0].header['crpix2']
    cdelt2 = Ihdu[0].header['cdelt2']
    hduwcs = wcs.WCS( Ihdu[0].header)
  except:
    print( 'Warning. No coordinate headers' )

  try:
    bmaj = Ihdu[0].header['bmaj']
    bmin = Ihdu[0].header['bmin']
    bpa  = Ihdu[0].header['bpa']
  except:
    print('Warnning. No header for synthesized beam size')

After having these header information, we can make conversion between the world coordiante system (i.e., astronomical coordinate systems) and the pixel coordinate system using the some methods of astropy.wcs. The following example convert the central pixel in the image to the R.A. and Decl. values and then convert them back ot the x and y pixel numbers:

xpix_temp = int(round( naxis1/2.0 ))
ypix_temp = int(round( naxis2/2.0 ))

world = hduwcs.wcs_pix2world([ [xpix_temp, ypix_temp] ], 0)
ra_center  = world[0][0]
dec_center = world[0][1]

pixcrd = hduwcs.wcs_world2pix([ [ra_center, dec_center] ], 0)
ra_center_pix  = pixcrd[0][0]
dec_center_pix = pixcrd[0][1]

The value of a certain pixel in a (two dimensional image) is simply:

value = Ihdu[0].data[ypix][xpix]

This allows you to analyze the image, for example, measuring aperature photometry.

The (two-dimensional) image can be output as a plot using the APLpy commands like:

figsize = [9.0, 9.0]
cmap           = 'viridis'
plot_ticks     = True
tick_font      = 25
tick_ypad      = 0
plot_colorbar  = True
colorbar_location = 'top'
colorbar_width = 0.3
colorbar_pad = 0.3
colorbar_font = 15
colorbar_label = 'Colorbar'
colorbar_labelpad = 0.3
colorbar_labelfont = 20
width = 1.0
height = 1.0
plot_scalebar = True
distance      = 140.0
scalebar_size = 100.0
scalebar_text = '100 au'
scalebar_color = (0,0,0)
scalebar_font = 25.0
scalebar_linewidth = 3.0
plot_beam = True
beam_color = 'black'

if_plot = False
try:
    self.fig = aplpy.FITSFigure(Ihdu, figsize=(figsize[0], figsize[1]))
    if_plot = True

except:
    print('Unable to load or plot hdu. Please double check the image file.')

if ( if_plot == True ):
    vmax = np.nanmax(Ihdu[0].data)*1.2
    vmin = 0.0

    fig.show_colorscale(
                        vmax = vmax, vmin = vmin,
                        interpolation = 'bicubic',
                        cmap = cmap
                       )

# ticks
fig.tick_labels.set_font(size=tick_font)
fig.axis_labels.set_font(size=tick_font)
fig.axis_labels.set_ypad(tick_ypad)
if (plot_ticks != True):
    fig.axis_labels.hide_x()
    fig.axis_labels.hide_y()
    fig.tick_labels.hide_x()
    fig.tick_labels.hide_y()

# recentering
if ( (ra_center != 0) or (dec_center !=0) ):
    print('recentering')
    fig.recenter(ra_center, dec_center, width=width, height=height)

# plot color bar
if (plot_colorbar == True):
    fig.add_colorbar()
    fig.colorbar.show()
    fig.colorbar.set_location(colorbar_location)
    fig.colorbar.set_width(colorbar_width)
    fig.colorbar.set_pad(colorbar_pad)
    fig.colorbar.set_font(size=colorbar_font, weight='medium', \
                               stretch='normal', family='sans-serif', \
                               style='normal', variant='normal')

    fig.colorbar.set_axis_label_text(colorbar_label)
    fig.colorbar.set_axis_label_font(size=colorbar_labelfont)
    fig.colorbar.set_axis_label_pad(colorbar_labelpad)
    
# plot scalebar
if (plot_scalebar == True):

    fig.add_scalebar( scalebar_size * (1.0/distance) / 3600.0)
    fig.scalebar.set_label(scalebar_text)
    fig.scalebar.set_color(scalebar_color)
    fig.scalebar.set_font(size=scalebar_font)
    fig.scalebar.set_linewidth(scalebar_linewidth)

If you would like to overplot another image as contours, try loading as another HDU (e.g., chdu in the following example) and complete the commands like:

fig.show_contour(chdu, colors=ccolor, linestyles=clinestyle,
                       levels=clevels,
                       linewidths=clinewidth)

Finally, to save an output figure, try ty complete the command like:

fig.save(outfigname, transparent=True)

You will make plots routinely. It is recommended to integrate the above commands in a script or a wrapper in a way that is convenient to you. If anything in this part is not clear to you, you can check the documentation for APLpy. You might need some patience and may need to try some things since not everything has been very clearly documented.

We can also use APLpy to create 3-color (RGB) figures, for example, see this page. However, my experience of trying this function was not very smooth. Following are some tips. First, we need three independent FITS images (e.g., the observations on the same area in the sky, at different wavelengths). Before using APLpy to generate a 3-color figure from them, we need to first regrid them to the same dimension (using either the CASA imregrid task or the Miriad regrid task, see above). Then we have to use the following CASA commands to produce the subimages (the area to be presented) and to remove the frequency/velocity and/or Stokes axes:

import os

path = './'
Rmap = path = '5ch3oh0_2_mnt0.img'  # my red image
Gmap = path = '5ch3oh2_2_mnt0.img'  # my green image
Bmap = path = '5ch3oh3_2_mnt0.img'  # my blue image

for imname in [Rmap, Gmap, Bmap]:
    importfits(
                fitsimage = imname + '.FITS',
                imagename = imname + '.image',
                overwrite = True
                )
                
    imsubimage(
                imagename = imname + 'image',
                outfile = imname + '.subim.image',
                box = '174,174,324,324', # the box sub-region to be plotted (blcx, blcy, trcx, trcy),
                dropdeg = True, overwrite = True
                )
                
    exportfits(
                imagename = imname + 'subim.image',
                fitsimage = imname + '.dropdeg.fits',
                overwrite = True, dropdeg = True, dropstokes = True
                )

Then the following Python code and produce the PDF figure:

# importing packages
import os
import numpy as np

from astropy.io.fits import getdata
from astropy import wcs
from astropy.io import fits
from astropy import units as u
from astropy import constants as con
from astropy.coordinates import SkyCoord

import matplotlib
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib

import aplpy
from astropy.io import fits

matplotlib.use('PDF')

# defining input image names
fits_r = '5CH3OH0_2.MNT0.IMG.dropdeg.fits'
fits_g = '5CH3OH2_2.MNT0.IMG.dropdeg.fits'
fits_b = '5CH3OH3_2.MNT0.IMG.dropdeg.fits'


# Combing the 3 independent FITS images to a single FITS image cube
os.system('rm -rf g10_rgb_cube.fits')
aplpy.make_rgb_cube([fits_r, fits_g, fits_b], 'g10_rgb_cube.fits')

# Make an RGB image (PNG figure)
aplpy.make_rgb_image('g10_rgb_cube.fits', 'g10_rgb.png')

# Insert a header item CTYPE3 such that APLpy will know what we are trying to do
fits.setval('g10_rgb_cube.fits','CTYPE3',value='RGB')

# Plot the RGB image using the 2d image to indicate the projection
fig = plt.figure(figsize=(7, 7))
f = aplpy.FITSFigure('g10_rgb_cube.fits', dimensions=[0,1], slices=[2],
                    figure=fig, subplot=[0.1, 0.1, 0.8, 0.8]
                    )
f.show_rgb('g10_rgb.png')
f.savefig('g10_rgb.pdf')

The commands to change/include the fontsize, ticks, labeling, beam, etc are identical to those for the single-color image, which have been introduced above.

4.2 IDL

(I used to use IDL to make plots. I just found that most of my old IDL codes can not anymore be compiled. Let’s forget about this part unless you are interested in archaeology or anthropology.)