Как получить 2D вырез изображения из позиции SkyCoord? - PullRequest
Новая
       59

Как получить 2D вырез изображения из позиции SkyCoord?

1 голос
/ 12 января 2020

Я следую примеру из документации Astropy для 2D Cutout.

Заголовок моего файла FITS: 

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. Хотя я не уверен, является ли это настоящей проблемой здесь. Кто-нибудь может помочь исправить эту ошибку?

1 Ответ

0 голосов
/ 13 января 2020

Так как ваш WCS 3D, но вы получаете 2d вырез, вам нужно отбросить 3-е измерение. Попробуйте cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1.celestial ), где .celestial - это удобный инструмент для удаления третьего измерения в WCS.

...