世界坐标系的共享Python接口#
背景#
这个 WCS 类实现了在FITS文件中表示世界坐标系的最常见的“标准”，但它不能表示任意复杂的转换，并且对于如何使用标准beyond FITS文件没有达成一致意见。因此，存在其他世界坐标系转换方法，例如 gwcs 正在为詹姆斯·韦伯太空望远镜开发的软件包（也适用于其他数据）。
由于Astropy项目的目标之一是改善包之间的互操作性，我们已经为Python中使用的世界坐标系对象协作定义了一个标准化的应用程序编程接口（API）。该API在Astropy Enhancements（APE）14: A shared Python interface for World Coordinate Systems .
核心astropy包提供了基类，这些基类定义了在 astropy.wcs.wcsapi 模块，这些都列在 参考/API 下面部分。
概述#
虽然API的全部细节和动机在APE14中有详细说明，但本文档总结了直接在astropy核心包中实现的元素。高级界面可能是普通用户最感兴趣的。尤其是，最重要的方法是 pixel_to_world() 和 world_to_pixel() 方法。它们提供了WCS的基本元素：映射到世界坐标系和从世界坐标系映射。其余部分通常提供有关 kind 关于WCS结构的世界坐标或类似信息。
更详细地说，这里实现的关键类是提供主用户界面的高级类 (BaseHighLevelWCS 和子类），以及较低级别的接口 (BaseLowLevelWCS 和子类）。这些可以是不同的对象 or 同一个。对于FITS-WCS WCS 用于FITS-WCS的对象遵循这两个接口，允许对已经包含FITS-WCS的文件立即使用此API。更具体的例子概述如下。
基本用法#
让我们先看一下WCS的共享Python接口，使用一个带有两个天轴（赤经和赤纬）的简单图像：
>>> from astropy.wcs import WCS
>>> from astropy.utils.data import get_pkg_data_filename
>>> from astropy.io import fits
>>> filename = get_pkg_data_filename('galactic_center/gc_2mass_k.fits')
>>> hdulist = fits.open(filename)
>>> hdu = hdulist[0]
>>> wcs = WCS(hdu.header)
>>> wcs
WCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 266.4 -28.93333
CRPIX : 361.0 360.5
NAXIS : 721 720
我们可以检查变换中有多少像素和世界坐标轴，以及WCS应用于的数据的形状：
>>> wcs.pixel_n_dim
2
>>> wcs.world_n_dim
2
>>> wcs.array_shape
(720, 721)
请注意，数组形状应与数据的形状匹配：
>>> hdu.data.shape
(720, 721)
如上所述 像素约定和定义 ，通常认为图像的y轴是第一维度，而图像的x轴是第二维度。因此 array_shape 返回中的形状 相反的 排序到FITS标题中的NAXIS关键字（在FITS-WCS的情况下）。如果您对相反顺序的数据形状感兴趣（在FITS-WCS的情况下，这将与NAXIS顺序相匹配），那么您可以使用 pixel_shape ：：
>>> wcs.pixel_shape
(721, 720)
现在让我们检查每个轴的物理类型：
>>> wcs.world_axis_physical_types
['pos.eq.ra', 'pos.eq.dec']
这确实是一个有两个天轴的图像。
新接口的主要部分定义了转换坐标的标准方法。最方便的方法是使用高级方法 pixel_to_world() 和 world_to_pixel() ，可以直接转换为天体：
>>> coord = wcs.pixel_to_world([1, 2], [4, 3])
>>> coord
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
[(266.97242993, -29.42584415), (266.97084321, -29.42723968)]>
同样，我们可以把天体物体转换回来——我们可以通过创建星系坐标来测试这一点，这些坐标将自动转换：
>>> from astropy.coordinates import SkyCoord
>>> coord = SkyCoord('00h00m00s +00d00m00s', frame='galactic')
>>> pixels = wcs.world_to_pixel(coord)
>>> pixels
(array(356.85179997), array(357.45340331))
如果您希望使用这些像素坐标索引原始数据，请确保使用 world_to_array_index() 它以正确的顺序返回坐标以索引Numpy数组，并舍入到最接近的整数值：
>>> index = wcs.world_to_array_index(coord)
>>> index
(357, 357)
>>> hdu.data[index]
563.7532
>>> hdulist.close()
高级用法#
现在让我们来看看光谱立方体的WCS（两个天轴和一个光谱轴）：
>>> filename = get_pkg_data_filename('l1448/l1448_13co.fits')
>>> hdulist = fits.open(filename)
>>> hdu = hdulist[0]
>>> wcs = WCS(hdu.header)
>>> wcs
WCS Keywords
Number of WCS axes: 3
CTYPE : 'RA---SFL' 'DEC--SFL' 'VOPT'
CRVAL : 57.6599999999 0.0 -9959.44378305
CRPIX : -799.0 -4741.913 -187.0
PC1_1 PC1_2 PC1_3 : 1.0 0.0 0.0
PC2_1 PC2_2 PC2_3 : 0.0 1.0 0.0
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0
CDELT : -0.006388889 0.006388889 66.42361
NAXIS : 105 105 53
如前所述，我们可以检查变换中有多少像素轴和世界轴，以及WCS应用于的数据的形状，以及每个轴的物理类型：
>>> wcs.pixel_n_dim
3
>>> wcs.world_n_dim
3
>>> wcs.array_shape
(53, 105, 105)
>>> wcs.world_axis_physical_types
['pos.eq.ra', 'pos.eq.dec', 'spect.dopplerVeloc.opt']
这确实是一个光谱立方体，有RA/Dec和一个速度轴。
如前所述，我们可以在像素和高层天体之间进行转换：
>>> celestial, spectral = wcs.pixel_to_world([1, 2], [4, 3], [2, 3])
>>> celestial
<SkyCoord (ICRS): (ra, dec) in deg
[(51.73115731, 30.32750025), (51.72414268, 30.32111136)]>
>>> spectral
<SpectralCoord
(target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
(57.66, 0., 1000.)
(pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
(0., 0., 0.)>)
[2661.04211695, 2727.46572695] m / s>
后面是：
>>> from astropy import units as u
>>> coord = SkyCoord('03h26m36.4901s +30d45m22.2012s')
>>> pixels = wcs.world_to_pixel(coord, 3000 * u.m / u.s)
>>> pixels
(array(8.11341207), array(71.0956641), array(7.10297292))
和之前一样，我们可以使用以下方法索引数组值：
>>> index = wcs.world_to_array_index(coord, 3000 * u.m / u.s)
>>> index
(7, 71, 8)
>>> hdu.data[index]
0.22262384
>>> hdulist.close()
如果您有兴趣在不使用高级astropy对象的情况下将世界值转换为简单的Python标量或Numpy数组，可以使用以下方法 pixel_to_world_values() 做这个-看 参考/API 有关详细信息，请参阅。
扩展FITS-WCS中的物理类型#
如上图所示， world_axis_physical_types 属性返回每个轴的物理类型列表。对于FITS-WCS，这是由标题中的CTYPE值确定的。在物理类型未知的情况下， None 返回。但是，可以使用 custom_ctype_to_ucd_mapping 上下文管理器。考虑具有以下CTYPE的WCS:：
>>> from astropy.wcs import WCS
>>> wcs = WCS(naxis=1)
>>> wcs.wcs.ctype = ['SPAM']
>>> wcs.world_axis_physical_types
[None]
我们可以指定对于这个CTYPE，物理类型应该是 'food.spam' ：：
>>> from astropy.wcs.wcsapi.fitswcs import custom_ctype_to_ucd_mapping
>>> with custom_ctype_to_ucd_mapping({'SPAM': 'food.spam'}):
... wcs.world_axis_physical_types
['food.spam']
WCS对象的切片#
在处理带有WCS信息的数据时，一个常见的操作是对WCS进行切片-这可以是为数据的一个子区域提取WCS，保留总的维度数（例如，从图像中剪切出来的），也可以降低数据和相关WCS的维数（例如从光谱中提取切片）立方体）。
这个 SlicedLowLevelWCS 类可用于切片任何符合 BaseLowLevelWCS 应用程序编程接口。为了证明这一点，让我们从一个光谱立方体文件开始：
>>> filename = get_pkg_data_filename('l1448/l1448_13co.fits')
>>> wcs = WCS(fits.getheader(filename, ext=0))
这个 wcs 对象是的实例 WCS 符合 BaseLowLevelWCS 应用程序编程接口。然后我们可以使用 SlicedLowLevelWCS 类来分割多维数据集：
>>> from astropy.wcs.wcsapi import SlicedLowLevelWCS
>>> slices = [10, slice(30, 100), slice(30, 100)]
>>> subwcs = SlicedLowLevelWCS(wcs, slices=slices)
这个 slices 参数接受切片、整数值和省略号的任意组合，这些组合通常会对Numpy数组进行切片。在上面的例子中，我们提取了一个光谱切片，在这个切片中我们提取了天空中的一个子区域。
如果您正在实现自己的WCS类，那么可以选择实现 __getitem__ 内部使用 SlicedLowLevelWCS . 事实上 WCS 类执行此操作-上面的示例可以更简洁地编写为：
>>> wcs[10, 30:100, 30:100]
<...>
SlicedFITSWCS Transformation
This transformation has 2 pixel and 2 world dimensions
Array shape (Numpy order): (70, 70)
Pixel Dim Axis Name Data size Bounds
0 None 70 None
1 None 70 None
World Dim Axis Name Physical Type Units
0 None pos.eq.ra deg
1 None pos.eq.dec deg
Correlation between pixel and world axes:
Pixel Dim
World Dim 0 1
0 yes yes
1 yes yes
这个切片基础设施能够处理具有相关轴的WCS对象的切片-在这种情况下，您可能最终得到具有不同数量像素和世界坐标的WCS。例如，如果我们对光谱立方体进行切片以提取与光谱切片图像平面中的一行相对应的1D数据集，则最终的WCS将具有一个像素维度和两个世界维度（因为RA/Dec随提取的1D切片而变化）：
>>> wcs[10, 40, :]
<...>
SlicedFITSWCS Transformation
This transformation has 1 pixel and 2 world dimensions
Array shape (Numpy order): (105,)
Pixel Dim Axis Name Data size Bounds
0 None 105 None
World Dim Axis Name Physical Type Units
0 None pos.eq.ra deg
1 None pos.eq.dec deg
Correlation between pixel and world axes:
Pixel Dim
World Dim 0
0 yes
1 yes