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.
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