问题描述
我已经在python中加载并绘制了FITS文件。 在上一篇文章的帮助下,我设法将轴从像素转换为天体坐标。但是我无法正确地以毫秒为单位获取它们。 代码如下
import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
filename = get_pkg_data_filename('hallo.fits')
hdu = fits.open(filename)[0]
wcs = WCS(hdu.header).celestial
wcs.wcs.crval = [0,0]
plt.subplot(projection=wcs)
plt.imshow(hdu.data[0][0],origin='lower')
plt.xlim(200,800)
plt.ylim(200,800)
plt.xlabel('Relative R.A ()')
plt.ylabel('Relative Dec ()')
plt.colorbar()
我不知道y标签被切掉是什么原因。
正如另一篇文章中所示,可以使用
wcs.wcs.ctype = [ 'XOFFSET','YOFFSET' ]
将其切换为毫秒,我得到
但是比例不正确!。 例如,0deg00min00.02sec应该是20 mas,而不是0.000002! 我在这里想念什么吗?
解决方法
看起来像光谱索引图。真好! 我认为问题可能在于FITS隐式地将度数用于CDELT之类的值。并且应将其明确转换为mas用于该图。 最直接的方法是将CDELT值乘以3.6e6,以将度转换为mas。 但是,如果您想在某个时候转换为不同的单位,则可以使用一种更通用的方法:
import astropy.units as u
w.wcs.cdelt = (w.wcs.cdelt * u.deg).to(u.mas)
因此,它基本上首先说CDELT的单位是度,然后将其转换为mas。
整个工作流程是这样的:
def make_transform(f):
'''use already read-in FITS file object f to build pixel-to-mas transformation'''
print("Making a transformation out of a FITS header")
w = WCS(f[0].header)
w = w.celestial
w.wcs.crval = [0,0]
w.wcs.ctype = [ 'XOFFSET','YOFFSET' ]
w.wcs.cunit = ['mas','mas']
w.wcs.cdelt = (w.wcs.cdelt * u.deg).to(u.mas)
print(w.world_axis_units)
return w
def read_fits(file):
'''read fits file into object'''
try:
res = fits.open(file)
return res
except:
return None
def start_plot(i,df=None,w=None,xlim = [None,None],ylim=[None,None]):
'''starts a plot and returns fig,ax .
xlim,ylim - axes limits in mas
'''
# make a transformation
# Using a dataframe
if df is not None:
w = make_transform_df(df)
# using a header
if w is not None:
pass
# not making but using one from the arg list
else:
w = make_transform(i)
# print('In start_plot using the following transformation:\n {}'.format(w))
fig = plt.figure()
if w.naxis == 4:
ax = plt.subplot(projection = w,slices = ('x','y',0 ))
elif w.naxis == 2:
ax = plt.subplot(projection = w)
# convert xlim,ylim to coordinates of BLC and TRC,perform transformation,then return back to xlim,ylim in pixels
if any(xlim) and any(ylim):
xlim_pix,ylim_pix = limits_mas2pix(xlim,ylim,w)
ax.set_xlim(xlim_pix)
ax.set_ylim(ylim_pix)
fig.add_axes(ax) # note that the axes have to be explicitly added to the figure
return fig,ax
rm = read_fits(file)
wr = make_transform(rm)
fig,ax = start_plot(RM,w=wr,xlim = xlim,ylim = ylim)
然后仅用imshow或轮廓或其他内容绘制到轴ax上。 当然,可以减少这段代码来满足您的特定需求。