python - 用astropy包修改python脚本
问题描述
我有 python 脚本,它采用合适的图像并绘制通量与像素的关系。它从对应于 0、30、60、90、120 和 150 度的选定中心点绘制 6 条线。然后它计算沿每条线的通量。我想使用 2 个图像并叠加它们,并分别为这两个图像绘制每一行以比较它们。即第一张图像的 0 度绘制在第二张图像的 0 度上。该脚本有一个用于天文学的 astropy 包,可以在线下载,试试 anaconda。
脚本如下:
#CODE CREATED BY A.AVISON 20-APR-2018#
import numpy as np
import scipy.ndimage
import matplotlib.pyplot as plt
from astropy.io import fits
#--- FITSNAME ---#
fitsName="fix_N7027.JVLA+eMer.eq.image.fits"
#---get fits image data---#
hdulist = fits.open(fitsName)
imData=hdulist[0].data
imData[np.isnan(imData)] = -0.0005 #== set NaNs to something which is a real number
roughSigma=np.median(imData[np.where(imData>0.0)])#=== takes the median non-zero values as a pseudoRMS
levs=np.arange(5.0,55.0,10.0)*roughSigma
#=== DEFINE CeNTRAL POSTION
y_peak=1043 #in pixels
x_peak=990 #in pixels
cutLen=500.0 #in pixels
#-- Plot...
fig, axes = plt.subplots(ncols=2, figsize=[12,6])
axes[0].imshow(imData,cmap='bone',origin='lower')#'gist_ncar'
axes[0].contour(imData,colors='w',linewidths=0.5,levels=levs)
#-- Extract data values along the line...
# Make a line with "num" points...
numCuts=6
angPerCuts=180.0/float(numCuts)
for ang in range(numCuts):
x1=x_peak-float(np.floor(cutLen*np.cos(np.radians(float(ang)*angPerCuts))))
y1=y_peak-float(np.floor(cutLen*np.sin(np.radians(float(ang)*angPerCuts))))
x0=x_peak+float(np.floor(cutLen*np.cos(np.radians(float(ang)*angPerCuts))))
y0=y_peak+float(np.floor(cutLen*np.sin(np.radians(float(ang)*angPerCuts))))
num = int(cutLen)
x, y = np.linspace(x0, x1, num), np.linspace(y0, y1, num)
# Extract the values along the line, using cubic interpolation
zi = 0
zi = scipy.ndimage.map_coordinates(imData, np.vstack((x,y)))
axes[0].plot(y,x, '-',linewidth=1.0,alpha=0.8,label=str(float(ang)*angPerCuts)+"degrees")
axes[0].plot(y1,x1, 'g.',linewidth=1.0,alpha=0.05)
axes[1].plot(zi,'-',alpha=0.8,linewidth=1.0,label=str(float(ang)*angPerCuts)+"degrees")
axes[1].set_xlim([0,cutLen])
#=== OVERLAY 1/r^x models ===#
axes[1].legend()
plt.show()
hdulist.close()
附上剧情。
有人建议将代码加倍(再次重复)并修改轴!!!。我怎样才能做到这一点?。
解决方案
推荐阅读
- javascript - Ember.js - classNameBindings 正在删除从外部添加的类
- c# - 如何针对 Active Directory 正确验证用户身份?
- serverless-framework - 使用无服务器框架通过部署调用 lambda 函数
- typescript - 如何通过在打字稿中分组将数组转换为矩阵
- php - 碳解析到 ISO8601
- excel - 启动画面导致 - 运行时错误 453:“在 user32 中找不到 DLL 入口点 FindWindowA”
- erlang - 列表:foreach 在 Erlang 中返回 ok
- c# - 使用捕获区域滑动 Unity 2D
- android - 安卓。我可以将资源 xml 读取为字符串吗?在运行时
- html - 带有旁白的 Edge 中的段落不可读