.fits图像的径向轮廓

Question

我一直试图用我在网上找到的修改过的脚本来绘制一个适合的图像的径向轮廓。我总是得到y轴单位，它与预期完全不同。我甚至不知道y轴单位是什么。我已经附加了fits文件和我不断得到的配置文件，以及我用另一个程序绘制的正确的径向轮廓。

我对python非常陌生，所以我不知道为什么会这样。对此的任何帮助都将受到极大的感谢。

这是我一直在使用的代码：

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

带错y轴单元的轮廓

​

​

轮廓与预期正确的单位创建使用凯尔特孔径测光工具。

​

​

Answer 1

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

​

​

编辑：

我添加了一个注释行，用于从标头中检索中心。但是，在选择使用argmax或头信息找到中心之前，您必须测试更多的fits文件。

标题信息的第一部分：

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