如何从 SkyCoord 位置获取图像的 2D Cutout？

Question

我正在按照 Astropy 文档中的 2D Cutout 示例进行操作。

我的 FITS 文件的header：

SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    3 / number of data axes                            
NAXIS1  =                  512 / length of data axis 1                          
NAXIS2  =                  512 / length of data axis 2                          
NAXIS3  =                    3 / length of data axis 3                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
SURVEY  = 'DECaLS  '                                                            
VERSION = 'DR8-south'                                                           
IMAGETYP= 'image   '                                                            
BANDS   = 'grz     '                                                            
BAND0   = 'g       '                                                            
BAND1   = 'r       '                                                            
BAND2   = 'z       '                                                            
CTYPE1  = 'RA---TAN'           / TANgent plane                                  
CTYPE2  = 'DEC--TAN'           / TANgent plane                                  
CRVAL1  =            186.11382 / Reference RA                                   
CRVAL2  =           0.15285422 / Reference Dec                                  
CRPIX1  =                256.5 / Reference x                                    
CRPIX2  =                256.5 / Reference y                                    
CD1_1   = -7.27777777777778E-05 / CD matrix                                     
CD1_2   =                   0. / CD matrix                                      
CD2_1   =                   0. / CD matrix                                      
CD2_2   = 7.27777777777778E-05 / CD matrix                                      
IMAGEW  =                 512. / Image width                                    
IMAGEH  =                 512. / Image height

到目前为止我已经尝试过：

from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

position = SkyCoord(hdu[0].header['CRVAL1']*u.deg,hdu[0].header['CRVAL2']*u.deg)
size = 200*u.pixel

wcs1 = WCS(hdu[0].header)

cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )

我运行最后一行出错了。错误：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-142-7cc21a13e941> in <module>
----> 1 cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )

/Applications/anaconda3/lib/python3.7/site-packages/astropy/nddata/utils.py in __init__(self, data, position, size, wcs, mode, fill_value, copy)
    714         if wcs is not None:
    715             self.wcs = deepcopy(wcs)
--> 716             self.wcs.wcs.crpix -= self._origin_original_true
    717             self.wcs.array_shape = self.data.shape
    718             if wcs.sip is not None:

ValueError: operands could not be broadcast together with shapes (3,) (2,) (3,)

我的猜测是因为我的文件中的 naxis =3 并且文档假定 naxis = 2。虽然我不确定这是否是这里的实际问题。有人可以帮助解决这个错误吗？

Answer 1

由于您的 WCS 是 3d，但您得到的是 2d 切口，因此您需要删除第 3 维。尝试 cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1.celestial )，其中 .celestial 是在 WCS 中删除第三维的便捷工具。

如何从 SkyCoord 位置获取图像的 2D Cutout？

How to obtain 2D Cutout of an image from a SkyCoord position?

fits

astropy