用python绘制SDSS图像

Question

我需要用 python 绘制一些 SDSS 图像。我从 SDSS DR7 下载了原始帧，并尝试使用 WCS 坐标投影显示（缩放）图像。 原始图像右侧为北，顶部为东。根据 DS9，对齐并不完美。 我想要最终的图像，顶部是北，左东。 我的尝试是这个：

import numpy as np
import matplotlib.pyplot as plt
import astropy.visualization as vis
import astropy.units as un
from astropy.io import fits
from astropy.wcs import WCS
from astropy.nddata.utils import Cutout2D

def open_image(imagename):
    hdu = fits.open(imagename)
    image = hdu[0].data
    header = hdu[0].header
    wcs = WCS(header)
    return image, header, wcs

def plot_image(imagename):
    image, header, wcs = open_image(imagename) #loading image

    print(wcs)

    #selecting a region centered on the target
    zoom = Cutout2D(image,(412,559),51*un.pixel, wcs = wcs, copy=True)

    #setting axis format
    fig,ax = plt.subplots(1,1, subplot_kw=dict(projection=zoom.wcs))
    ra = ax.coords[0]
    dec = ax.coords[1]
    ra.set_major_formatter('hh:mm:ss.s')
    dec.set_major_formatter('dd:mm:ss')

    #setting image scale
    interval = vis.PercentileInterval(99.9)
    vmin,vmax = interval.get_limits(zoom.data)
    norm = vis.ImageNormalize(vmin=vmin, vmax=vmax, stretch=vis.LogStretch(1000))
    ax.imshow(zoom.data, cmap =plt.cm.Reds, norm = norm, origin = 'lower')          
    ax.set_ylabel('Dec.')
    ax.set_xlabel('RA')
    plt.show()

if __name__ == '__main__':
    imagename = '../SDSS_g_orig.fits'
    plot_image(imagename)

您可以在附件中找到结果。坐标轴上不存在坐标。图像也保持其原始方向。我假设 WCS 投影会自动旋转图像以对齐它，但显然不是这样。

如何将图像与左上北对齐？

P.s.：这是 print(wcs) 行的输出：

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 195.77014030000001  16.47383778  
CRPIX : 1024.5  744.5  
CD1_1 CD1_2  : 6.1883041204266702e-06  0.000109834787851833  
CD2_1 CD2_2  : 0.000109820625  -6.2125672043002999e-06  
NAXIS : 2048  1489

output image

编辑：我尝试使用SkyView 模块，但结果并不令人满意。图像的强度轮廓（图像+轮廓是我的最终目标）是不同的，并且天空视图图像在目标中间显示了某种伪影。

Original with contours Skyview with contours

Answer 1

在天文学 Facebook 组中的 Python 用户的帮助下，我找到了答案。 它需要安装Montage 和montage-wrapper，但它允许旋转图像并继续使用astropy 和matplotlib。

我只需要以这种方式修改 open_image() 函数：

import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
import montage_wrapper as montage

def open_image(imagename):

    hdu = fits.open(imagename)
    hdu = montage.reproject_hdu(hdu[0], north_aligned=True) 
    image = hdu.data
    nans = np.isnan(image)
    image[nans] = 0
    header = hdu.header
    wcs = WCS(header)
    return image, header, wcs

生成的图像将以北向上的方式旋转。脚本的其余部分不需要修改。只有 Cutout2D 对象必须在目标上居中。

我正在附加带有轮廓的最终图像。

final+contours

Answer 2

更简单的方法可能是使用astroquery's SkyView module。例如：

import matplotlib.pyplot as plt
from astroquery.skyview import SkyView
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

# Query for SDSS g images centered on target name
hdu = SkyView.get_images("M13", survey='SDSSg')[0][0]

# Tell matplotlib how to plot WCS axes
wcs = WCS(hdu.header)
ax = plt.gca(projection=wcs)

# Plot the image
ax.imshow(hdu.data)
ax.set(xlabel="RA", ylabel="Dec")
plt.show()

生成这样的图像：

您可以将“M13”替换为您的目标名称或坐标——参见SkyView.get_images docs for all of the finer knobs you can tune。

然后，如果你想翻转任一轴的方向，你可以调用

ax.invert_xaxis()

或

ax.invert_yaxis()