pythonunits-of-measurementastronomyastropyfits

Radial Profile from a .fits image


I have been trying to plot a radial profile of a fits image using a modified script I found on-line. I always get y axis units which are completely different to what's expected. I'm not even sure what the y axis units are. I have attached the fits file and a profile I keep getting and the correct radial profile I plotted using another program.

I am very new to python so I have no idea why this keeps happening. Any help to fix this will be so greatly appreciated.

This is the code I've been using:

import numpy as np
import pyfits
import matplotlib.pyplot as plt
import matplotlib.image as mpimg

def azimuthalAverage(image, center=None):
    """
    Calculate the azimuthally averaged radial profile.

    image - The 2D image
    center - The [x,y] pixel coordinates used as the center. The default is 
             None, which then uses the center of the image (including 
             fracitonal pixels).

    """
    # Calculate the indices from the image
    y, x = np.indices(image.shape)


    if not center:
        center = np.array([(x.max()-x.min())/2.0, (y.max()-y.min())/2.0])

    r = np.hypot(x - center[0], y - center[1])

    # Get sorted radii
    ind = np.argsort(r.flat)
    r_sorted = r.flat[ind]
    i_sorted = image.flat[ind]

    # Get the integer part of the radii (bin size = 1)
    r_int = r_sorted.astype(int)

    # Find all pixels that fall within each radial bin.
    deltar = r_int[1:] - r_int[:-1]  # Assumes all radii represented
    rind = np.where(deltar)[1]       # location of changed radius
    nr = rind[1:] - rind[:-1]        # number of radius bin

    # Cumulative sum to figure out sums for each radius bin
    csim = np.cumsum(i_sorted, dtype=float)
    tbin = csim[rind[1:]] - csim[rind[:-1]]

    radial_prof = tbin / nr
    print center
    print i_sorted
    print radial_prof
    return radial_prof

#read in image
hdulist = pyfits.open('cit6ndf2fitsexample.fits')
scidata = np.array(hdulist[0].data)[0,:,:]
center = None
radi = 10
rad = azimuthalAverage(scidata, center)

plt.xlabel('radius(pixels?)', fontsize=12)
plt.ylabel('image intensity', fontsize=12)
plt.xlim(0,10)
plt.ylim(0, 3.2)
plt.plot(rad[radi:])
plt.savefig('testfig1.png')
plt.show()

Profile with wrong y axis units

enter image description here

Profile with expected correct units created using Celtech Aperture Photometry Tool.

enter image description here


Solution

  • from astropy.io import fits
    import numpy as np
    import matplotlib.pyplot as plt
    from matplotlib.ticker import AutoMinorLocator
    
    minorLocator = AutoMinorLocator()
    
    
    def radial_profile(data, center):
        x, y = np.indices((data.shape))
        r = np.sqrt((x - center[0])**2 + (y - center[1])**2)
        r = r.astype(np.int)
    
        tbin = np.bincount(r.ravel(), data.ravel())
        nr = np.bincount(r.ravel())
        radialprofile = tbin / nr
        return radialprofile 
    
    
    fitsFile = fits.open('testfig.fits')
    img = fitsFile[0].data[0]
    img[np.isnan(img)] = 0
    
    #center = np.unravel_index(img.argmax(), img.shape)
    center = (-fitsFile[0].header['LBOUND2']+1, -fitsFile[0].header['LBOUND1']+1)
    rad_profile = radial_profile(img, center)
    
    fig, ax = plt.subplots()
    plt.plot(rad_profile[0:22], 'x-')
    
    ax.xaxis.set_minor_locator(minorLocator)
    
    plt.tick_params(which='both', width=2)
    plt.tick_params(which='major', length=7)
    plt.tick_params(which='minor', length=4, color='r')
    plt.grid()
    ax.set_ylabel(fitsFile[0].header['Label'] + " (" + fitsFile[0].header['BUNIT'] + ")")
    ax.set_xlabel("Pixels")
    plt.grid(which="minor")
    plt.show()
    

    enter image description here

    EDIT:

    I added a commented line for retrieving the center from the headers. But you would have to test more fits files before choosing to use argmax or the header info to find the center.

    First part of the header info:

    SIMPLE  =                    T / file does conform to FITS standard             
    BITPIX  =                  -64 / number of bits per data pixel                  
    NAXIS   =                    3 / number of data axes                            
    NAXIS1  =                  259 / length of data axis 1                          
    NAXIS2  =                  261 / length of data axis 2                          
    NAXIS3  =                    1 / length of data axis 3                          
    EXTEND  =                    T / FITS dataset may contain extensions            
    COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
    COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
    LBOUND1 =                 -133 / Pixel origin along axis 1                      
    LBOUND2 =                 -128 / Pixel origin along axis 2                      
    LBOUND3 =                    1 / Pixel origin along axis 3                      
    OBJECT  = 'CIT 6   '           / Title of the dataset                           
    LABEL   = 'Flux Density'       / Label of the primary array                     
    BUNIT   = 'mJy/arcsec**2'      / Units of the primary array                     
    DATE    = '2015-12-18T06:45:40' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
    ORIGIN  = 'East Asian Observatory' / Origin of file                             
    BSCALE  =                  1.0 / True_value = BSCALE * FITS_value + BZERO       
    BZERO   =                  0.0 / True_value = BSCALE * FITS_value + BZERO       
    HDUCLAS1= 'NDF     '           / Starlink NDF (hierarchical n-dim format)       
    HDUCLAS2= 'DATA    '           / Array component subclass                       
    HDSTYPE = 'NDF     '           / HDS data type of the component                 
    TELESCOP= 'JCMT    '           / Name of Telescope