astropy WCS 修复：备选 WCS（key=‘A‘/‘B‘）畸变关键字按后缀独立读取的完整解析

• 科学计算
• 科研
• 数据分析

导读

本文围绕 astropy 仓库中的修复条目 19228.bugfix.rst 展开：当一张 FITS 图像同时包含主 WCS（无后缀）与备选 WCS（key='A'、'B'等）时，原先的畸变（distortion）关键字解析会错误地把主 WCS 的畸变信息加载到备选 WCS 上，甚至导致加载崩溃。修复后，备选 WCS 会严格依据自身关键字后缀（如CPDIS1A、DP1A）独立读取畸变查找表。读完本文，你将掌握 FITS 备选 WCS 畸变关键字的命名规则、astropy.wcs.WCS的底层解析逻辑，以及如何构造、校验和复现这一修复场景。

背景：FITS 备选 WCS 与关键字后缀约定

FITS 标准允许在同一图像头中存储多套世界坐标系统（WCS）。主 WCS 的关键字不带后缀（如CTYPE1、CRPIX1、CDELT1），而备选 WCS 通过在关键字名后追加单字母后缀（A–Z）来区分，例如CTYPE1A、CRPIX1A、CDELT1A对应key='A'的 WCS。

在 astropy 中，通过wcs.WCS(header, fobj, key='A')即可选择读取哪一套 WCS（key的取值是空格' '或'A'–'Z'，见 wcs.py 构造签名）。关键问题在于：畸变关键字同样遵循这一后缀规则，而旧版本实现没有为畸变解析正确拼接后缀，导致张冠李戴。

Bug 描述：备选 WCS 误载主畸变或直接崩溃

修复条目 19228.bugfix.rst 原文明确了问题与修复效果：

Fixed WCS distortion keyword handling to respect alternate WCS key suffixes. Secondary WCS (key='A', 'B', etc.) now correctly read their own distortion keywords (CPDIS1A, DP1A, etc.) instead of incorrectly loading primary WCS distortion or crashing.

即此前存在两类故障表现：

1. 错误加载主畸变：当备选 WCS（如key='A'）没有自己的畸变关键字（不存在CPDIS1A）时，解析器会错误地回退读取主 WCS 的CPDIS1/CPDIS2畸变查找表，从而把主 WCS 的畸变错误地施加到备选 WCS 坐标变换上，产生错误的天球坐标；
2. 解析崩溃：在部分关键字组合（例如存在DP*系列辅助关键字但缺少带正确后缀的入口）下，解析过程会抛出异常中断加载。

修复的核心，是让_read_distortion_kw在构造畸变关键字名时把key后缀（如'A'）拼入CPDIS/DP等关键字，使备选 WCS 只认自己的畸变关键字。

修复原理：源码级解析_read_distortion_kw的后缀拼接逻辑

修复落点在WCS.__init__中调用畸变解析的位置：构造时会以key=key传入_read_distortion_kw，见 wcs.py#L632-L635：

det2im = self._read_det2im_kw(header, fobj, err=minerr) cpdis = self._read_distortion_kw( header, fobj, dist="CPDIS", err=minerr, key=key ) self._fix_pre2012_scamp_tpv(header) sip = self._read_sip_kw(header, wcskey=key)

而 _read_distortion_kw 的实现 是整个修复的核心，其关键步骤如下：

key = key.strip().upper() key_suffix = key if dist == "CPDIS": d_kw = "DP" err_kw = "CPERR" else: d_kw = "DQ" err_kw = "CQERR" tables = {} for i in range(1, self.naxis + 1): axis_key_suffix = str(i) + key_suffix d_error_key = err_kw + axis_key_suffix ... distortion = dist + axis_key_suffix

逐一拆解这段逻辑：

• 后缀统一归一化：key先做strip().upper()，保证'a'与'A'等价；
• 主/备选切换：key=' '（主 WCS）时key_suffix为空，关键字保持CPDIS1、CPERR1、DP1.EXTVER形式；key='A'时生成CPDIS1A、CPERR1A、DP1A.EXTVER；key='B'时则对应CPDIS1B、DP1B.*，依此类推；
• 按轴遍历：对每个轴i分别构造err_kw + str(i) + key_suffix与dist + str(i) + key_suffix，即每个轴独立判断；
• 误差阈值过滤：err参数（构造器中的minerr，默认0.0）——若CPERRiX值小于该阈值，则该轴畸变被跳过（tables[i] = None）；
• 回退读取扩展：对DPiX.AXIS.i指示的轴序判断是否需要对WCSDVARR扩展数据做转置，再从扩展头读取CRPIX/CRVAL/CDELT构建DistortionLookupTable对象；
• 清除已消费关键字：读取后从 header 中删除对应的畸变关键字（del header[...]），避免它们被后续传给 WCSLIB 的核心 WCS 解析流程。

关键关键字速查表

以key='A'的备选 WCS 为例，修复后的正确关键字集合如下（主 WCS 去掉A后缀即可）：

_read_distortion_kw对dist参数还支持"CQDIS"（sequent 畸变，对应DQ/CQERR关键字族），这是同一套后缀逻辑的另一应用分支，见 wcs.py#L1207-L1212。

回归测试：两条测试用例如何锁定修复行为

修复的正确性由 test_wcs.py 中的两条测试 严格锁定，它们分别覆盖了 bug 的两种故障表现：

1.test_alternate_wcs_distortion_key_isolation（隔离性）

构造一个只有主 WCS 带畸变（CPDIS1/CPDIS2=LOOKUP，并附DP1.EXTVER、DP1.NAXES、DP1.AXIS.1等）、备选 WCS 没有任何畸变关键字（注释明确写着 "NOTE: No CPDIS1A, CPDIS2A keywords"）的头部。随后：

wcs_primary = wcs.WCS(header, hdulist, key=" ") assert wcs_primary.cpdis1 is not None assert wcs_primary.cpdis2 is not None wcs_secondary = wcs.WCS(header, hdulist, key="A") assert wcs_secondary.cpdis1 is None, \ "Secondary WCS should NOT load primary WCS distortion when CPDIS1A is missing" assert wcs_secondary.cpdis2 is None, \ "Secondary WCS should NOT load primary WCS distortion when CPDIS2A is missing"

这正是对"备选 WCS 不得回退加载主畸变"这一修复点位的直接验证。

2.test_alternate_wcs_distortion_independent（独立性）

构造主 WCS 与备选 WCS 各自拥有独立畸变的头部：主 WCS 用CPDIS1/DP1.EXTVER = 1指向WCSDVARRver 1；备选 WCS 用CPDIS1A/DP1A.EXTVER = 2指向WCSDVARRver 2，且两份查找表数据内容不同。测试断言：

wcs_secondary = wcs.WCS(header, hdulist, key="A") assert wcs_secondary.cpdis1 is not None assert wcs_secondary.cpdis2 is not None assert wcs_primary.cpdis1 is not wcs_secondary.cpdis1 assert not np.array_equal(wcs_primary.cpdis1.data, wcs_secondary.cpdis1.data)

即备选 WCS 必须从带A后缀的关键字独立加载自己的畸变，且数据与主 WCS 完全隔离。修复前的实现无法通过这两条测试。

实战：构造一个可运行的备选 WCS 畸变加载示例

以下示例可直接复制运行（需 astropy 已安装），完整复现上述两条测试的构造与验证思路：

import numpy as np from astropy.io import fits from astropy import wcs # 1. 主 HDU，写入主 WCS 与备选 WCS (key='A') 的核心关键字 primary = fits.PrimaryHDU() header = primary.header header["NAXIS"] = 2 header["NAXIS1"] = 1024 header["NAXIS2"] = 1024 # 主 WCS（无后缀） header["CTYPE1"] = "RA---TAN" header["CTYPE2"] = "DEC--TAN" header["CRPIX1"] = 512 header["CRPIX2"] = 512 header["CRVAL1"] = 150.0 header["CRVAL2"] = -35.0 header["CDELT1"] = -0.0002777778 header["CDELT2"] = 0.0002777778 # 备选 WCS（带 A 后缀） header["CTYPE1A"] = "HPLN-TAN" header["CTYPE2A"] = "HPLT-TAN" header["CRPIX1A"] = 512 header["CRPIX2A"] = 512 header["CRVAL1A"] = 0.0 header["CRVAL2A"] = 0.0 header["CDELT1A"] = 0.5 header["CDELT2A"] = 0.5 # 2. 主 WCS 的畸变关键字（无后缀，指向 WCSDVARR ver 1） header["CPDIS1"] = "LOOKUP" header["CPDIS2"] = "LOOKUP" header["DP1.EXTVER"] = 1 header["DP1.NAXES"] = 2 header["DP1.AXIS.1"] = 1 header["DP1.AXIS.2"] = 2 header["DP2.EXTVER"] = 1 header["DP2.NAXES"] = 2 header["DP2.AXIS.1"] = 1 header["DP2.AXIS.2"] = 2 # 3. 备选 WCS 的畸变关键字（带 A 后缀，指向 WCSDVARR ver 2） header["CPDIS1A"] = "LOOKUP" header["CPDIS2A"] = "LOOKUP" header["DP1A.EXTVER"] = 2 header["DP1A.NAXES"] = 2 header["DP1A.AXIS.1"] = 1 header["DP1A.AXIS.2"] = 2 header["DP2A.EXTVER"] = 2 header["DP2A.NAXES"] = 2 header["DP2A.AXIS.1"] = 1 header["DP2A.AXIS.2"] = 2 # 4. 两份不同的畸变查找表数据（129x129，值不同以区分） dist_data_primary = np.random.randn(129, 129).astype(np.float32) * 0.01 dist_data_secondary = np.random.randn(129, 129).astype(np.float32) * 0.02 def make_dist_hdu(data, ver): hdu = fits.ImageHDU(data, name="WCSDVARR") hdu.ver = ver hdu.header["CRPIX1"] = 65.0 hdu.header["CRPIX2"] = 65.0 hdu.header["CRVAL1"] = 513.0 hdu.header["CRVAL2"] = 1.0 hdu.header["CDELT1"] = 8.0 hdu.header["CDELT2"] = 8.0 return hdu hdulist = fits.HDUList([ primary, make_dist_hdu(dist_data_primary, 1), make_dist_hdu(dist_data_secondary, 2), ]) # 5. 分别加载并断言隔离性 wcs_primary = wcs.WCS(header, hdulist, key=" ") wcs_secondary = wcs.WCS(header, hdulist, key="A") assert wcs_primary.cpdis1 is not None assert wcs_secondary.cpdis1 is not None assert wcs_primary.cpdis1 is not wcs_secondary.cpdis1 assert not np.array_equal(wcs_primary.cpdis1.data, wcs_secondary.cpdis1.data) print("主 WCS 与备选 WCS 畸变独立加载成功")

关键点：fobj必须传入fits.HDUList，否则LOOKUP类型畸变会抛出ValueError（"an astropy.io.fits.HDUList is required for Lookup table distortion"，见 wcs.py#L1236-L1240）。

相关机制：SIP 畸变同样遵循后缀规则

本次修复还联动保证了一条既有约定：SIP 畸变关键字的读取也支持备选后缀。_read_sip_kw接收wcskey=key参数（见 wcs.py#L638 与 wcs.py#L1416-L1480），并按CTYPE{nax}{wcskey}检查-SIP后缀一致性；若检测到 SIP 系数但CTYPE缺少-SIP后缀，会给出坐标可能不正确的提示日志。相关回归测试包括test_sip_with_altkey（验证key='A'时读取CTYPE1A = "RA---SIN-SIP"）与test_keyedsip（key="A"时正确读取A_ORDER/A_i_j系数），见 test_wcs.py#L1449-L1501。

底层执行链：畸变查找表如何进入坐标变换管线

从 C 层看，cpdis1/cpdis2是DistortionLookupTable对象，最终以distortion_lookup_t*指针数组的形式挂入 WCS 管线（pipeline.c#L22-L40）。坐标变换时：

• pix2foc阶段调用p4_pix2foc将像素坐标先经 prior 畸变（cpdis）映射到"焦点平面"（astropy_wcs.c#L373）；
• 管线中has_p4判断cpdis[0]/cpdis[1]是否非空，决定是否启用该阶段（pipeline.c#L92）。

因此cpdis1/cpdis2若被错误赋值（旧 bug 的根源），备选 WCS 的all_pix2world/all_world2pix坐标结果就会混入主 WCS 的畸变校正。

使用注意事项与限制

1. 仅支持 2 维：FITS WCS distortion paper 查找表与 SIP 畸变只支持 2 维；若核心 WCS 维度大于 2，构造会抛出ValueError，提示通过naxis参数降维后使用（wcs.py#L681-L692）；
2. to_header不写出畸变：to_header明确不写 distortion paper 信息（需多 HDU 表示），完整写出应使用to_fits（经由_write_distortion_kw写出CPDISi、DPi.EXTVER、DPi.NAXES、DPi.AXIS.i及WCSDVARR扩展，见 wcs.py#L1285-L1341）；
3. 多项式畸变未实现：若CPDISiX的值不是LOOKUP（如TPD等多项式类型），当前实现仅发出AstropyUserWarning（"Polynomial distortion is not implemented"），不加载畸变（wcs.py#L1272-L1276）；
4. minerr阈值：构造函数minerr（默认0.0）可跳过最大误差小于阈值的畸变，可用于条件性忽略低精度畸变表。

总结

19228 号修复从根本上一并解决了备选 WCS 畸变关键字的命名解析与数据隔离两个问题：_read_distortion_kw现在严格按照key后缀拼接CPDIS/CPERR/DP等关键字族，使key='A'/'B'等备选 WCS 只能读取自己的畸变，绝不回退到主 WCS；配套的隔离性与独立性两条回归测试保证了该行为长期稳定。对于处理多 WCS 天文图像（如 HST/地面望远镜多 CCD 数据）的开发者而言，这套后缀规则与验证方法是构造合规 FITS 头、排查坐标偏差的重要依据。

• 科学计算
• 科研
• 数据分析