使用 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

更简单的方法可能是使用 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()

Answer 2

在 Astronomy Facebook 群组中 Python 用户的帮助下，我找到了答案。它需要安装 Montage and montage-wrapper but it allows to rotate the image and to continue work with astropy and 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

使用 python 绘制 SDSS 图像

Plotting SDSS images with python

python

astronomy

fits