Concatenate FITS spectra together

fitscat

#!/usr/bin/env python
# fitscat -- concatenate several FITS spectra together

import sys
import os
import getopt
from astropy.io import fits as pyfits
from numpy import *
from scipy import interpolate

speed_of_light = 299792458. # m/s
VERSION = "0.1"

def usage():
    """
    Display usage

    """
    
    print """Usage: fitscat [option] fitsfile1 fitsfile2 [fitsfile3 [...]]
Contatenate several FITS cubes

Arguments:   
   fitsfile1              a FITS cube

Options:
   -h,--help              display this help
   -V,--version           display version
   -o,--output=filename   output file name (default cat.fits)
   -l,--vlsr=v            source velocity in the LSR frame in km/s (default 0)
   -d,--delete            delete original files
   -v,--verbose           verbose mode
"""

def version():
    """
    Display version number

    """

    print "This is fitscat, version %s" % VERSION
    print """Copyright (c) 2010-2015 Sebastien Maret

This is free software. You may redistribute copies of it under the terms
of the GNU General Public License. There is NO WARRANTY, to the extent
permitted by law."""

def readfitscube(fitsfile, vlsr = 0):
    """
    Read a FITS cube file

    This function reads a cube in FITS format. The 3rd axis units
    should be in velocity units. The velocity is converted into rest
    frequency, after taking into account the source VLSR.
    
    Arguments:
    fitsfile  -- fits file name
    vlsr      -- source velocity in the LSR frame in km/s (default 0)
    
    """

    vlsr *= 1e3 # m/s -> km/s

    fits = pyfits.open(fitsfile)

    header = fits[0].header# Primary HDU header

    # Check that we have a valid fits file

    if header['naxis'] != 3:
        raise ValueError, "%s does not have 3 axis" % fitsfile
    if header['ctype3'] != 'VELO-LSR':
        raise ValueError, "%s has an invalid type for axis 3" % fitsfile
    try:
        rest_freq = header['restfreq']
    except:
        raise ValueError, "%s has no RESTFREQ keyword" % fitsfile

    # Compute the frequency values

    naxis3 = header['naxis3']
    cdelt3 = header['cdelt3']
    crpix3 = header['crpix3']
    crval3 = header['crval3']

    freq = arange(naxis3, dtype = float32)
    for i in range(naxis3):
        veloc = (i + 1 - crpix3) * cdelt3 + crval3 + vlsr
        freq[i] = rest_freq * (1 - veloc / speed_of_light) # Hz -- check formulae !

    # Now read the brightness temperatures

    tb = fits[0].data

    # Reverse the frequency axis

    freq = freq[::-1]
    tb = tb[::-1,:,:]
    
    return (freq, tb, header)
    
def main ():
    """
    Main program for fitscat

    """

    # Parse command line options and arguments

    verbose = False
    fileout = "cat.fits"
    vlsr = 0.
    delete = False

    try:
	opts, args = getopt.getopt(sys.argv[1:], "hvVo:l:d",
				   ["help", "version", "verbose", "output", "vlsr", "delete"]) 
    except getopt.GetoptError:
	print >> sys.stderr, "fitscat: error: unknown option"
        print >> sys.stderr, "Type 'fitscat --help' for more information"
	sys.exit(1)

    for opt, arg in opts:
	if opt in ("-h", "--help") :
	    usage()
	    sys.exit(0)
	if opt in ("-v", "--version"):
	    version()
            sys.exit(0)
	if opt in ("-V", "--verbose"):
	    verbose = True
        if opt in ("-o", "--output"):
            fileout = arg
        if opt in ("-l", "--vlsr"):
            try:
                vlsr = float(arg)
            except ValueError:
                print >> sys.stderr, "fitscat: error: incorrect value for %s option" % opt
                print >> sys.stderr, "Type 'fitscat --help' for more information"
        if opt in ("-d", "--delete"):
            delete = True

    if len(args) < 2:
	print >> sys.stderr, "fitscat: error: incorrect number of arguments"
        print >> sys.stderr, "Type 'fitscat --help' for more information"
	sys.exit(1)

    fitsfile = args

    # We need to know the frequency range that is covered by the
    # cubes, so we need to read them all first
    
    input_cubes = []
    header = None
    min_freq, max_freq, min_freq_step = 1e99, 0., 1e99

    for f in fitsfile:

        if verbose:
            print "Reading %s..." % f
        try:
            freq, tb, hdr = readfitscube(f, vlsr = vlsr)
        except ValueError, err:
            print >> sys.stderr, "fitscat: error: %s" % err
            sys.exit(1)

        # Save the first header as well as the minimum frequency,

        if not(header): header = hdr
        if min(freq) < min_freq: min_freq = min(freq)
        if max(freq) > max_freq: max_freq = max(freq)
        if abs(freq[1] - freq[0]) < min_freq_step: min_freq_step = abs(freq[1] - freq[0])

        input_cubes.append([freq, tb])

    # Remove the continuum on each cubes
    cont = input_cubes[0][1][0] # assumes that all cubes have the same continuum
    for c in input_cubes:
        c[1] -= cont
    
    # Create an output cube covering the entire frequency range
    # and with the smallest frequency sampling

    freq_step = min_freq_step
    output_freq = arange(min_freq, max_freq, freq_step, dtype = float32)
    npix1, npix2 = len(cont[0]), len (cont[1])
    output_tb = zeros((len(output_freq), npix1, npix2), dtype = float32)

    # Resample each FITS cube on the frequency grid and sum them

    if verbose:
        print "Concatenating cubes..."

    for c in input_cubes:
        for i in range(npix1):
            for j in range(npix2):
                tb_interp = interpolate.interp1d(c[0], c[1][:,i,j], bounds_error = False, fill_value = 0.)
                output_tb[:,i,j] += tb_interp(output_freq)

    # Add the continuum
    output_tb += cont

    # Write the FITS

    central_freq = (output_freq[-1] + output_freq[0])/2.
    crpix3 = len(output_freq)/2. + 1.

    header['ctype3'] = "FREQ"
    header['crval3'] = 0.
    header['crpix3'] = crpix3
    header['cdelt3'] = output_freq[1] - output_freq[0]
    header['restfreq'] = central_freq

    if os.path.exists(fileout):
        os.unlink(fileout)

    hdu = pyfits.PrimaryHDU(output_tb, header = header)
    hdu.writeto(fileout)

    if verbose:
        print "%s created." % fileout

    if delete:
        for f in fitsfile:
            os.unlink(f)

if __name__ == "__main__":
    main()

I'm not sure if it works for cubes or only spectra. In the later case the program can be easily adapted to handle cubes.

The last version should work with cubes. I'm not 100% that the frequency scale is correct thought.