来自 .fits 图像的径向轮廓
Radial Profile from a .fits image
我一直在尝试使用我在网上找到的修改后的脚本来绘制拟合图像的径向轮廓。我总是得到与预期完全不同的 y 轴单位。我什至不确定 y 轴单位是什么。我附上了拟合文件和我不断获得的配置文件以及我使用另一个程序绘制的正确径向配置文件。
我是 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 轴单位错误的配置文件
使用 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()
编辑：
我添加了一条注释行，用于从 header 中检索中心。但是在选择使用 argmax 或 header 信息来找到中心之前，您必须测试更多的拟合文件。
header 信息的第一部分：
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
我一直在尝试使用我在网上找到的修改后的脚本来绘制拟合图像的径向轮廓。我总是得到与预期完全不同的 y 轴单位。我什至不确定 y 轴单位是什么。我附上了拟合文件和我不断获得的配置文件以及我使用另一个程序绘制的正确径向配置文件。
我是 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 轴单位错误的配置文件
使用 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()
编辑：
我添加了一条注释行，用于从 header 中检索中心。但是在选择使用 argmax 或 header 信息来找到中心之前，您必须测试更多的拟合文件。
header 信息的第一部分：
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