WCS comme projection matplotlib pour une tranche de cube de données chargée à l'aide de l'astrophotographie?

Question

J'ai un fichier FITS nommé 'my_cube.fits' avec WCS. Le fichier contient des informations spatiales sur les axes 1 et 2 (X et Y) et des informations spectrales sur l'axe 3 (Z). Quand je charge à l'aide de astropy.io.fits, l'axe spectral est égal à 0 et sont des axes spatiaux 1 et 2. Le fichier est chargé comme suit:

import astropy.io.fits as pyfits 

    filename = 'my_cube.fits' 
    my_data = pyfits.getdata(filename) 
    my_header = pyfits.getheader(filename) 

J'utilise matplotlib pour afficher des données et I Je veux savoir comment afficher une seule trame spectrale de mon cube de données en utilisant son WCS. Disons que:

from astropy.wcs import WCS 
    from matplotlib import pyplot as plt 

    my_wcs = WCS(my_header) 
    fig = plt.figure() 
    ax = fig.add_subplot(111, projection=my_wcs) 
    ax.imshow(my_data[5, :, :]) 

    plt.show() 

Si je fais juste ça, j'ai:

... 
    File "/usr/local/lib/python3.4/dist-packages/matplotlib/figure.py", line 1005, in add_subplot 
    a = subplot_class_factory(projection_class)(self, *args, **kwargs) 
    File "/usr/local/lib/python3.4/dist-packages/matplotlib/axes/_subplots.py", line 73, in __init__ 
    self._axes_class.__init__(self, fig, self.figbox, **kwargs) 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/core.py", line 49, in __init__ 
    self.patch = self.coords.frame.patch 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/frame.py", line 129, in patch 
    self._update_patch_path() 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/frame.py", line 115, in _update_patch_path 
    self.update_spines() 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/frame.py", line 192, in update_spines 
    self['b'].data = np.array(([xmin, ymin], [xmax, ymin])) 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/frame.py", line 40, in data 
    self._world = self.transform.transform(self._data) 
    File "/usr/local/lib/python3.4/dist-packages/wcsaxes/transforms.py", line 166, in transform 
    world = self.wcs.wcs_pix2world(pixel_full, 1) 
    File "/usr/local/lib/python3.4/dist-packages/astropy/wcs/wcs.py", line 1329, in wcs_pix2world 
'output', *args, **kwargs) 
    File "/usr/local/lib/python3.4/dist-packages/astropy/wcs/wcs.py", line 1231, in _array_converter 
    return _return_single_array(xy, origin) 
    File "/usr/local/lib/python3.4/dist-packages/astropy/wcs/wcs.py", line 1212, in _return_single_array 
"of shape (N, {0})".format(self.naxis)) 
    ValueError: When providing two arguments, the array must be of shape (N, 3) 

Answer 1

Votre objet WCS contient les deux axes spatiaux et l'axe spectral, puisque vous avez initialisé avec l'en-tête complet (Je suppose que votre en-tête contient également une solution spectrale WCS). Ainsi, créer un tracé 2D avec seulement l'espace ne fonctionnera pas, puisque votre projection donnée à la sous-placette est en 3D.

Je ne l'ai pas fait moi-même, je n'ai pas de fichier d'exemple, mais la documentation WCS mentionne que vous pouvez utiliser the sub function, or even the celestial functions pour obtenir les axes individuels ou les axes célestes (spatiaux); ces fonctions retourneront un objet WCS avec juste ces axes.

Ainsi, vous pouvez probablement utiliser les éléments suivants:

my_wcs = WCS(my_header).celestial 
fig = plt.figure() 
ax = fig.add_subplot(111, projection=my_wcs) 

Answer 2

Ceci est un bon cas d'utilisation pour spectral-cube, qui enveloppe efficacement astropy.io.fits pour des utilisations du cube.

La solution d'Evert, une fois corrigée, fonctionnera, en supposant que vous avez wcsaxes installé.

Pour utiliser spectral_cube pour cela, un exemple simple est:

cube = SpectralCube.read(filename) 
cube[5,:,:].quicklook() # uses aplpy 

# using wcsaxes 
fig = plt.figure() 
ax = fig.add_axes([0.1,0.1,0.8,0.8], projection=cube[5,:,:].wcs) 
ax.imshow(cube[5,:,:]) # you may need cube[5,:,:].value depending on mpl version