切片多维数据

WCSAxes可以绘制一维或二维数据。如果我们有一个数据集的维数比我们要绘制的图高，我们必须选择哪些维度用于绘图的x或x和y轴。本例将演示如何切片FITS数据立方体并从中绘制图像。

切片WCS对象

就像中介绍的例子一样 使用世界坐标初始化轴 ，我们将使用 astropy.io.fits 并解析WCS信息。原始的FITS文件可以从 here .

import matplotlib.pyplot as plt
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename
filename = get_pkg_data_filename('l1448/l1448_13co.fits')
hdu = fits.open(filename)[0]
wcs = WCS(hdu.header)
image_data = hdu.data

这是一个三维数据集，您可以通过以下方式检查标题信息：

>>> hdu.header  
...
NAXIS = 3 /number of axes
CTYPE1  = 'RA---SFL'           /
CTYPE2  = 'DEC--SFL'           /
CTYPE3  = 'VELO-LSR'           /
...

头关键字“NAXIS”给出数据集的维度数。关键字“CTYPE1”、“CTYPE2”和“CTYPE3”分别给出了这些维度的数据类型：赤经、赤纬和速度。

然后实例化 WCSAxes 使用 WCS 对象并选择要绘制的切片：

import matplotlib.pyplot as plt
ax = plt.subplot(projection=wcs, slices=(50, 'y', 'x'))

通过设置 slices=(50, 'y', 'x') ，我们选择在y轴上绘制第二维度，在x轴上绘制第三维度。即使我们没有绘制所有的尺寸，我们也必须指定要为未显示的尺寸选择哪些切片。在本例中，我们没有绘制第一个尺寸，所以我们选择了要显示的切片50。您可以通过更改所选切片并查看打印图像的变化来进行试验。

绘制图像

然后我们将轴添加到图像中并使用该方法绘制它 imshow() .

ax.coords[2].set_ticklabel(exclude_overlapping=True)
ax.imshow(image_data[:, :, 50].transpose())

(png, svg, pdf)

在这里， image_data 是一个 ndarray 对象。在Numpy中，轴的顺序颠倒，因此FITS文件中的第一个尺寸显示在最后，最后一个尺寸显示在第一个，依此类推。因此，索引传递给 imshow() 应该与传递给 slices 但顺序相反。我们也需要 transpose() image_data 因为我们已经颠倒了在切片的x轴和y轴上绘制的尺寸。

如果不想反转绘制的尺寸，可以简单地执行以下操作：

import matplotlib.pyplot as plt
ax = plt.subplot(projection=wcs, slices=(50, 'x', 'y'))
ax.imshow(image_data[:, :, 50])

(png, svg, pdf)

绘制一维数据

如果我们想绘制一个像素的光谱轴，我们可以将其切片到一维。

import matplotlib.pyplot as plt
ax = plt.subplot(projection=wcs, slices=(50, 50, 'x'))

这里我们选择了第一和第二维度的50个像素，并将第三个维度作为我们的x轴。

我们现在可以画出这个像素的光谱轴。注意，我们在调用 ax.plot ， WCSAxes 将为我们显示世界坐标。

ax.plot(image_data[:, 50, 50])

因为这仍然是一个 WCSAxes 绘图，我们可以设置x轴的显示单位

ra, dec, vel = ax.coords
vel.set_format_unit(u.km/u.s)

(png, svg, pdf)

如果我们想沿着一个空间维度绘制一个一维的图，即沿着图像中一行的强度， WCSAxes 默认为显示此绘图的两个世界坐标。我们可以定制颜色并为每个空间轴添加网格线。

import matplotlib.pyplot as plt
ax = plt.subplot(projection=wcs, slices=(50, 'x', 0))

ax.plot(image_data[0, :, 50])

ra, dec, wave = ax.coords
ra.set_ticks(color="red")
ra.set_ticklabel(color="red")
ra.grid(color="red")

dec.set_ticks(color="blue")
dec.set_ticklabel(color="blue")
dec.grid(color="blue")

(png, svg, pdf)