我一直在寻找并寻找解决此问题的方法,但一无所获。
我正在通过 matplotlib 生成矩形 FITS 图像,然后使用 AstroPy(或 PyFITS)将 WCS 坐标应用于它们。我的图像是银河经纬度,所以适合我的地图的标题关键字应该是GLON-CAR
和GLAT-CAR
(对于笛卡尔投影)。我查看了在SAO DS9中使用相同地图投影的其他地图,并且坐标效果很好......网格完全正交,因为它应该是。可以在此处找到 FITS 标准投影。
但是当我生成我的地图时,坐标根本不是笛卡尔坐标。这是我的地图(左)和另一张大致相同区域的参考地图(右)的并排比较。两者都在 FITS 标头中列出GLON-CAR
,GLAT-CAR
但在 SAO DS9 中查看时我的很奇怪(请注意,坐标网格是 SAO DS9 根据 FITS 标头中的数据生成的,或者至少存储在 FITS 文件中的某个位置):
这是有问题的,因为如果投影错误,坐标分配算法将为每个像素分配不正确的坐标。
有没有人遇到过这个问题,或者知道可能是什么问题?
我已经尝试应用其他投影(只是为了看看它们在 SAO DS9 中的表现如何)并且结果很好……但是我的笛卡尔和墨卡托投影并没有像他们应该的那样使用正交网格。
我不敢相信这会是 AstroPy 中的错误,但我找不到任何其他原因......除非我在标题中的参数格式不正确,但我仍然不明白这会如何导致我的问题我正在经历。或者你会推荐使用别的东西吗?(我查看了 matplotlib 底图,但在我的计算机上运行时遇到了一些麻烦)。
我的标题代码如下:
from __future__ import division
import numpy as np
from astropy.io import fits as pyfits # or use 'import pyfits, same thing'
#(lots of code in between: defining variables and simple calculations...
#probably not relevant)
header['BSCALE'] = (1.00000, 'REAL = TAPE*BSCALE + BZERO')
header['BZERO'] = (0.0)
header['BUNIT'] = ('mag ', 'UNIT OF INTENSITY')
header['BLANK'] = (-100.00, 'BLANK VALUE')
header['CRVAL1'] = (glon_center, 'REF VALUE POINT DEGR') #FIRST COORDINATE OF THE CENTER
header['CRPIX1'] = (center_x+0.5, 'REF POINT PIXEL LOCATION') ## REFERENCE X PIXEL
header['CTYPE1'] = ('GLON-CAR', 'COORD TYPE : VALUE IS DEGR')
header['CDELT1'] = (-glon_length/x_length, 'COORD VALUE INCREMENT WITH COUNT DGR') ### degrees per pixel
header['CROTA1'] = (0, 'CCW ROTATION in DGR')
header['CRVAL2'] = (glat_center, 'REF VALUE POINT DEGR') #Y COORDINATE OF THE CENTER
header['CRPIX2'] = (center_y+0.5, 'REF POINT PIXEL LOCATION') #Y REFERENCE PIXEL
header['CTYPE2'] = ('GLAT-CAR', 'COORD TYPE: VALUE IS DEGR') # WAS CAR OR TAN
header['CDELT2'] = (glat_length/y_length, 'COORD VALUE INCREMENT WITH COUNT DGR') #degrees per pixel
header['CROTA2'] = (rotation, 'CCW ROTATION IN DEGR') #NEGATIVE ROTATES CCW around origin (bottom left).
header['DATAMIN'] = (data_min, 'Minimum data value in the file')
header['DATAMAX'] = (data_max, 'Maximum data value in the file')
header['TELESCOP'] = ("Produced from 2MASS")
pyfits.update(filename, map_data, header)
感谢您的任何帮助,您可以提供。