Radial Profile from a .fits image

Question

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

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

You say you attached the .fits file but I can't see it anywhere. — Gabriel, Jan 23 '16 at 15:32
Also, if you found the script on-line, add a link to the script. Is this the one you are using? http://image-tools.readthedocs.org/en/latest/api/image_tools.radialprofile.azimuthalAverage.html?highlight=azimuthalaverage#azimuthalaverage — Gabriel, Jan 23 '16 at 15:38
And one more thing: what program are you using to get the "correct" radial density profile? — Gabriel, Jan 23 '16 at 15:46
Possible duplicate of [Most efficient way to calculate radial profile](http://stackoverflow.com/questions/21242011/most-efficient-way-to-calculate-radial-profile) — M4rtini, Jan 23 '16 at 16:20
Where does the "correct" profile comes from? Why do you expect those values? Keep in mind that FITS has some scaling keywords that could be causing this. And why don't you plot the first 10 values (`plot(rad[radi:])`, with `radi = 10`)? — , Jan 23 '16 at 22:22
You should show the header of the SCI extension of the file you're using. — Iguananaut, Jan 24 '16 at 14:15
@Gabriel. Sorry I have left out the link to the fits file. Here it is: https://dl.dropboxusercontent.com/u/96130369/testfig.fits. And this is the link to the script [file: radialProfile.py]: http://www.astrobetter.com/wiki/python_radial_profiles. I have used the Caltech Aperture Photometry tool for it. Thank You! — BubbleGum, Jan 24 '16 at 16:12
@M4rtini. I did really try to solve my problem using that question but it really did not seem to help. That question was more of optimising a radial profile. This is more of getting a python code like that to run for a .fits file which I think is what's causing this problem. Thank You! — BubbleGum, Jan 24 '16 at 16:13
@Iguananaut. Sorry I forgot to include the .fits file I used to create this profile. The fits file includes the header with all the image details. Thank You! — BubbleGum, Jan 24 '16 at 16:15
@Evert. The 'correct' profile is something I generated from Aperture Photometry Tool. And that is exactly what I have done in my code (The third and ninth lines from the bottom of my code). Thank You! — BubbleGum, Jan 24 '16 at 16:21
@BubbleGum did you check the answer by M4rtini below? If it helped you solve your issue, don't forget to mark it as "accepted". — Gabriel, Jan 25 '16 at 21:36

M4rtini · Accepted Answer · 2016-01-26T06:56:43.693

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()

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

Sorry about my comment before. It does run for that one fits file, but for every other fits file there is still an error in the y axis values. I really can't figure out where this comes from. Here are the other test fits files: [https://www.dropbox.com/sh/ddg9roz9pio2s5x/AAC5fmozzhGbOH0ZlytUYu1ea?dl=0]. These are the profiles they give using the above code: [https://www.dropbox.com/sh/br93ss6r4xqqbqe/AAChqUGCbLoMBKlI5GlnVzpMa?dl=0] and these are the profiles I created using the Aperture Photometry Tool: [https://www.dropbox.com/sh/j2wodhes1cum37y/AAB1xPHqGFA7MYwX2y99Ka4Na?dl=0]. Thank You! — BubbleGum, Jan 26 '16 at 06:34
Hi, Here is a link to three of the other files I tried: https://www.dropbox.com/sh/ddg9roz9pio2s5x/AAC5fmozzhGbOH0ZlytUYu1ea?dl=0 It turns out when I run it using the center from the header as you suggested it does work fine. So maybe something goes wrong in the array calculation for the max value? Thanks! — BubbleGum, Jan 26 '16 at 06:49
Use the second option for calculating the `center` variable. i just tested `test2.fits` and it seems right to me. I'll comment out the first, and make taking it from the headers the main one. — M4rtini, Jan 26 '16 at 06:55
Yes that will work. I think the error for the other comes up because the arrays have all the noise values as well so if there is a noise spike which is the max value the code considers this to be the max array value. Thank You for all the help!! — BubbleGum, Jan 26 '16 at 07:23

score 0 · Answer 2 · answered Jun 20 '23 at 22:54

Instead of hard coding, you can use AstroPy's Photutils package. It has an in-built functionality to give radial profiles and also supports units. You can use it something like this:

from astropy.io import fits
import numpy as np
import matplotlib.pyplot as plt
from photutils.centroids import centroid_quadratic
from photutils.profiles import RadialProfile


# opening the fits   
with fits.open('test.fits') as hdu:
    data = hdu[0].data

# finding the centroid
xycen = centroid_quadratic(data, xpeak=48, ypeak=52)
print(xycen)

# getting the radial profile
edge_radii = np.arange(15)
rp = RadialProfile(data, xycen, edge_radii)
rp.plot(label='Radial Profile')

You can also perform a Gaussian fit of your profile too.

plt.plot(rp.radius,rp.gaussian_fit,label='Gaussian Fit')
plt.axvline(rp.gaussian_fwhm,label=f'Gaussian FWHM {rp.gaussian_fwhm} pix')

Radial Profile from a .fits image

2 Answers2