You need to enable JavaScript to run this app.
优惠活动
大模型
产品
解决方案
定价
更多

处理不同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 /

错误原因分析

  1. NAXIS与WCS轴不匹配:Orion文件NAXIS=0(无数据轴),但WCS声明WCSAXES=3,导致WCS无法关联到实际数据,触发第一个警告。
  2. SFL投影参数异常:Sanson-Flamsteed(SFL)投影对LATPOLE等参数有严格要求,当前LATPOLE、CRVAL1、CRVAL2的组合无法生成有效坐标转换矩阵,导致No valid solution for latp错误。
  3. 冗余参数干扰:表头中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

相关产品推荐
方舟 Agent Plan

超全模态模型 × Harness 升级,最新支持 Deepseek-V4.1-Flash、GLM-5.3 系列、Doubao-Seedream-5.0-pro、Kimi-K3 (部分), 限时 9.9 元起

最近更新时间:2026.07.22 13:32:02