处理不同FITS文件时出现的WCS错误排查求助
问题:处理Orion FITS文件时WCS初始化失败
我在处理两个FITS文件,示例文件能正常运行,但Orion文件初始化WCS时报错。
核心处理代码
from astropy.io import fits from astropy.wcs import WCS hdul = fits.open(cubefits) hdul.info() hdul.readall() hd = hdul[0].header wcs = WCS(hd)
示例文件WCS正常输出
WCS Keywords Number of WCS axes: 3 CTYPE : 'RA---TAN' 'DEC--TAN' 'FREQ' CRVAL : 0.70395330810396 0.69721933449752 125.0 CRPIX : 507.90474039695 503.02270076207 1.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.0013888888888889 0.0013888888888889 1.0 NAXIS : 1009 1009 181
Orion文件的错误信息
In [94]: wcs = WCS(hd) WARNING: FITSFixedWarning: The WCS transformation has more axes (3) than the image it is associated with (0) [astropy.wcs.wcs] WARNING: FITSFixedWarning: 'celfix' made the change 'Ill-conditioned coordinate transformation parameter'. [astropy.wcs.wcs] --------------------------------------------------------------------------- InvalidTransformError Traceback (most recent call last) Input In [94], in <cell line: 1>() ----> 1 wcs = WCS(hd) InvalidTransformError: ERROR 7 in wcsset() at line 2660 of file cextern/wcslib/C/wcs.c: Ill-conditioned coordinate transformation parameter. ERROR 4 in celset() at line 441 of file cextern/wcslib/C/cel.c: Ill-conditioned coordinate transformation parameters No valid solution for latp for these values of phip, phi0, and theta0.
两个FITS文件表头信息
示例文件(EXAMPLE)
SIMPLE = T / conforms to FITS standard BITPIX = -64 / array data type NAXIS = 3 / number of array dimensions NAXIS1 = 1009 NAXIS2 = 1009 NAXIS3 = 181 WCSAXES = 3 / Number of coordinate axes CRPIX1 = 507.90474039695 / Pixel coordinate of reference point CRPIX2 = 503.02270076207 / Pixel coordinate of reference point CRPIX3 = 1.0 / Pixel coordinate of reference point CDELT1 = 0.0013888888888889 / [deg] Coordinate increment at reference point CDELT2 = 0.0013888888888889 / [deg] Coordinate increment at reference point CDELT3 = 1.0 / [Hz] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CUNIT3 = 'Hz' / Units of coordinate increment and value CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection CTYPE3 = 'FREQ' / Frequency (linear) CRVAL1 = 0.70395330810396 / [deg] Coordinate value at reference point CRVAL2 = 0.69721933449752 / [deg] Coordinate value at reference point CRVAL3 = 125.0 / [Hz] Coordinate value at reference point LONPOLE = 180.0 / [deg] Native longitude of celestial pole LATPOLE = 0.69721933449752 / [deg] Native latitude of celestial pole RADESYS = 'ICRS' / Equatorial coordinate system CUBE = DATAS = DATE = '2020-03-30 15:48:36.001328' / date of the creation
Orion文件
SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T / HIERARCH LASER_SHIFT = '-5.603303303303305' / [sample] Shift between kids and la CONTINUE '100_S49749_NA_NA' SCAN = '49749 ' OBJECT = 'OrionSim' / OBSTYPE = 'OTF ' / Observation Type OBSDATE = '2022-08-31T00:00:00.000' SITELAT = '-23.00575' / [deg] Site latitude SITELONG= '-67.7592222222222' SITEELEV= '5105.0 ' / [m] Site elevation PROJID = 'E-0109.A-0442A-2022' SCANNUM = '49749 ' OBS-ID = '49749 ' DATE = '2023-02-15T11:57:57.513287' EXPTIME = '395.13672846369445' / [s] Exposure time DATE-OBS= '2022-08-31T11:01:12.421' DATE-END= '2022-08-31T11:07:47.558' ORIGIN = 'None ' NINT = '736 ' / Total number of blocs NPTINT = '2048 ' / Number of points per bloc NSAMPLES= '1507328 ' LABEL = 'None ' INTTIME = '0.000262144' / [s] Integration Time per sample WCSAXES = 3 / Number of coordinate axes CRPIX1 = 381 / Pixel coordinate of reference point CRPIX2 = 649 / Pixel coordinate of reference point CRPIX3 = 0.00000000000 / Pixel coordinate of reference point CDELT1 = 0.0011111111111111 / [deg] Coordinate increment at reference point CDELT2 = 0.0011111111111111 / [deg] Coordinate increment at reference point CDELT3 = 600000000.000 / [Hz] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CUNIT3 = 'Hz' / Units of coordinate increment and value CTYPE1 = 'RA---SFL' / Right ascension, Sanson-Flamsteed projection CTYPE2 = 'DEC--SFL' / Declination, Sanson-Flamsteed projection CTYPE3 = 'FREQ' / Coordinate type code CRVAL1 = 83.8402938843 / [deg] Coordinate value at reference point CRVAL2 = -5.18931388855 / [deg] Coordinate value at reference point CRVAL3 = 0.600000000000 / [Hz] Coordinate value at reference point LONPOLE = 0.0 / [deg] Native longitude of celestial pole LATPOLE = 69.997710986349 / [deg] Native latitude of celestial pole MJDREF = 0.0 / [d] MJD of fiducial time RADESYS = 'ICRS' / Equatorial coordinate system CRDELT1 = 0 / CRDELT2 = 0 /
错误原因分析
- NAXIS与WCS轴不匹配:Orion文件
NAXIS=0(无数据轴),但WCS声明WCSAXES=3,导致WCS无法关联到实际数据,触发第一个警告。 - SFL投影参数异常:Sanson-Flamsteed(SFL)投影对
LATPOLE等参数有严格要求,当前LATPOLE、CRVAL1、CRVAL2的组合无法生成有效坐标转换矩阵,导致No valid solution for latp错误。 - 冗余参数干扰:表头中
CRDELT1=0和CRDELT2=0属于冗余无效参数,与已定义的CDELT冲突,造成转换矩阵病态。
解决方案
1. 修复数据轴信息
如果Orion文件应为三维数据立方体,补充NAXIS1、NAXIS2、NAXIS3并设置NAXIS=3;如果数据存储在扩展HDU中,尝试读取对应扩展的表头:
hdul = fits.open(cubefits) # 遍历HDU找包含有效数据的扩展 for hdu in hdul: if hdu.data is not None and hdu.header.get('NAXIS', 0) > 0: hd = hdu.header wcs = WCS(hd) break
2. 修正SFL投影参数
将LATPOLE改为与CRVAL2一致,或标准天极值(ICRS北天极为90°):
hd = hdul[0].header # 方案1:LATPOLE与CRVAL2保持一致 hd['LATPOLE'] = hd['CRVAL2'] # 方案2:设置为北天极 # hd['LATPOLE'] = 90.0 wcs = WCS(hd)
3. 移除冗余无效参数
删除冲突的CRDELT1和CRDELT2:
hd = hdul[0].header for key in ['CRDELT1', 'CRDELT2']: if key in hd: del hd[key] wcs = WCS(hd)
4. 禁用自动修复并手动指定轴
若自动修复引发问题,可禁用FITS固定功能并手动提取所需WCS轴:
# 禁用自动修复 wcs = WCS(hd, fix=False) # 仅保留天球坐标轴 wcs = wcs.sub(['celestial'])
内容的提问来源于stack exchange,提问作者Raizen
相关产品推荐
相关产品推荐

