- c - 在位数组中找到第一个零
- linux - Unix 显示有关匹配两种模式之一的文件的信息
- 正则表达式替换多个文件
- linux - 隐藏来自 xtrace 的命令
我一直在尝试使用我在网上找到的修改后的脚本来绘制拟合图像的径向轮廓。我总是得到与预期完全不同的 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()
编辑:
我添加了注释行以从标题中检索中心。但是在选择使用 argmax
或标题信息来查找中心之前,您必须测试更多的拟合文件。
头信息的第一部分:
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
关于python - 来自 .fits 图像的径向轮廓,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/34965275/
我正在处理一组标记为 160 个组的 173k 点。我想通过合并最接近的(到 9 或 10 个组)来减少组/集群的数量。我搜索过 sklearn 或类似的库,但没有成功。 我猜它只是通过 knn 聚类
我有一个扁平数字列表,这些数字逻辑上以 3 为一组,其中每个三元组是 (number, __ignored, flag[0 or 1]),例如: [7,56,1, 8,0,0, 2,0,0, 6,1,
我正在使用 pipenv 来管理我的包。我想编写一个 python 脚本来调用另一个使用不同虚拟环境(VE)的 python 脚本。 如何运行使用 VE1 的 python 脚本 1 并调用另一个 p
假设我有一个文件 script.py 位于 path = "foo/bar/script.py"。我正在寻找一种在 Python 中通过函数 execute_script() 从我的主要 Python
这听起来像是谜语或笑话,但实际上我还没有找到这个问题的答案。 问题到底是什么? 我想运行 2 个脚本。在第一个脚本中,我调用另一个脚本,但我希望它们继续并行,而不是在两个单独的线程中。主要是我不希望第
我有一个带有 python 2.5.5 的软件。我想发送一个命令,该命令将在 python 2.7.5 中启动一个脚本,然后继续执行该脚本。 我试过用 #!python2.7.5 和http://re
我在 python 命令行(使用 python 2.7)中,并尝试运行 Python 脚本。我的操作系统是 Windows 7。我已将我的目录设置为包含我所有脚本的文件夹,使用: os.chdir("
剧透:部分解决(见最后)。 以下是使用 Python 嵌入的代码示例: #include int main(int argc, char** argv) { Py_SetPythonHome
假设我有以下列表,对应于及时的股票价格: prices = [1, 3, 7, 10, 9, 8, 5, 3, 6, 8, 12, 9, 6, 10, 13, 8, 4, 11] 我想确定以下总体上最
所以我试图在选择某个单选按钮时更改此框架的背景。 我的框架位于一个类中,并且单选按钮的功能位于该类之外。 (这样我就可以在所有其他框架上调用它们。) 问题是每当我选择单选按钮时都会出现以下错误: co
我正在尝试将字符串与 python 中的正则表达式进行比较,如下所示, #!/usr/bin/env python3 import re str1 = "Expecting property name
考虑以下原型(prototype) Boost.Python 模块,该模块从单独的 C++ 头文件中引入类“D”。 /* file: a/b.cpp */ BOOST_PYTHON_MODULE(c)
如何编写一个程序来“识别函数调用的行号?” python 检查模块提供了定位行号的选项,但是, def di(): return inspect.currentframe().f_back.f_l
我已经使用 macports 安装了 Python 2.7,并且由于我的 $PATH 变量,这就是我输入 $ python 时得到的变量。然而,virtualenv 默认使用 Python 2.6,除
我只想问如何加快 python 上的 re.search 速度。 我有一个很长的字符串行,长度为 176861(即带有一些符号的字母数字字符),我使用此函数测试了该行以进行研究: def getExe
list1= [u'%app%%General%%Council%', u'%people%', u'%people%%Regional%%Council%%Mandate%', u'%ppp%%Ge
这个问题在这里已经有了答案: Is it Pythonic to use list comprehensions for just side effects? (7 个答案) 关闭 4 个月前。 告
我想用 Python 将两个列表组合成一个列表,方法如下: a = [1,1,1,2,2,2,3,3,3,3] b= ["Sun", "is", "bright", "June","and" ,"Ju
我正在运行带有最新 Boost 发行版 (1.55.0) 的 Mac OS X 10.8.4 (Darwin 12.4.0)。我正在按照说明 here构建包含在我的发行版中的教程 Boost-Pyth
学习 Python,我正在尝试制作一个没有任何第 3 方库的网络抓取工具,这样过程对我来说并没有简化,而且我知道我在做什么。我浏览了一些在线资源,但所有这些都让我对某些事情感到困惑。 html 看起来
我是一名优秀的程序员,十分优秀!