WCS als Matplotlib-Projektion für eine Datenwürfelscheibe, die mit Hilfe von Astropie geladen wurde?

Question

Ich habe eine FITS-Datei namens 'my_cube.fits' mit WCS. Die Datei enthält räumliche Informationen zu Achse 1 und 2 (X und Y) und spektrale Informationen zu Achse 3 (Z). Als ich es laden astropy.io.fits Verwendung ist die spektrale Achse 0 und Raumachsen sind 1 und 2. Die Datei wie folgt geladen:

import astropy.io.fits as pyfits 

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

Ich habe matplotlib wurde unter Verwendung von Daten anzuzeigen und I möchte wissen, wie man einen einzelnen Spektralrahmen meines Datenwürfels mit seinem WCS anzeigt. Lassen Sie uns sagen:

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() 

Wenn ich das nur tun, ich habe:

... 
    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

Ihr WCS Objekt, um die zwei Raumachsen und die spektrale Achse enthält, da man es mit dem vollständigen Header initialisiert (Ich nehme an, dass Ihr Header auch eine spektrale WCS-Lösung enthält). Daher wird das Erstellen eines 2D-Diagramms mit nur dem räumlichen Format nicht funktionieren, da Ihre Projektion für das Unterdiagramm 3D ist.

Ich habe das selbst nicht gemacht, noch habe ich eine Beispieldatei, aber die WCS-Dokumentation erwähnt, dass Sie the sub function, or even the celestial functions verwenden können, um die einzelnen Achsen oder die Himmels- (räumlichen) Achsen zu erhalten; Diese Funktionen geben ein WCS-Objekt mit genau diesen Achsen zurück.

So können Sie wahrscheinlich die folgenden Befehle verwenden:

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

Answer 2

Dies ist ein guter Anwendungsfall für spectral-cube, die effektiv astropy.io.fits für Cube verwendet wickelt.

Die einmal korrigierte Evert-Lösung funktioniert, vorausgesetzt, Sie haben wcsaxes installiert.

Um spectral_cube für diese, ist ein einfaches Beispiel:

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