Spectrum, Image and Cube format


Spectrum, Image and Cube objects all have the following items, where O denotes the name of the object:

Component Description
O.filename A FITS filename if the data were loaded from a FITS file
O.primary_header A FITS primary header instance
O.wcs World coordinate spatial information (WCS object)
O.wave World coordinate spectral information (WaveCoord object)
O.shape An array of the dimensions of the cube
O.data A numpy masked array of pixel values
O.data_header A FITS data header instance
O.unit The physical units of the data values
O.dtype The data-type of the data array (int, float)
O.var An optional numpy masked array of pixel variances
O.mask An array of the masked state of each pixel
O.ndim The number of dimensions in the data, variance and mask arrays

Masked arrays are arrays that can have missing or invalid entries. The numpy.ma module provides a nearly work-alike replacement for numpy that supports data arrays with masks. See the DataArray documentation for more details.

When an object is constructed from a MUSE FITS file, O.data will contain the DATA extension, O.var will contain the STAT extension and O.mask will contain the DQ extension if it exists. A DQ extension contains the pixel data quality. By default all bad pixels are masked. But it is possible for the user to create his mask by using the Python interface for euro3D convention.


The format of each numpy array follows the indexing used by Python to handle 2D and 3D arrays. For an MPDAF cube C, the pixel in the bottom-lower-left corner is C[0,0,0] and pixel C[k,p,q] refers to the horizontal position q, the vertical position p, and the spectral position k, as follows:


In total, this cube C contains nq pixels in the horizontal direction, np pixels in the vertical direction and nk channels in the spectral direction. In the cube, each numpy masked array has 3 dimensions, array[k,p,q], where k is the spectral axis, and p and q are the spatial axes:

In [1]: from mpdaf.obj import Cube

# data array is read from the file (extension number 0)
In [2]: cube = Cube(filename='sdetect/minicube.fits')

# The 3 dimensions of the cube:
In [3]: cube.shape
Out[3]: (3681, 40, 40)

In [4]: cube.data.shape
Out[4]: (3681, 40, 40)

In [5]: cube.var.shape
Out[5]: (3681, 40, 40)

In [6]: cube.mask.shape
Out[6]: (3681, 40, 40)

Cube[k,p,q] returns the value of Cube.data[k,p,q]:

In [7]: cube[3659, 8, 28]
Out[7]: 493.9726257324219

Similarly Cube[k,p,q] = value sets the value of Cube.data[k,p,q], and Cube[k1:k2,p1:p2,q1:q2] = array sets the corresponding subset of Cube.data. Finally Cube[k1:k2,p1:p2,q1:q2] returns a sub-cube, as demonstrated in the following example:

In [8]: cube.info()
[INFO] 3681 x 40 x 40 Cube (sdetect/minicube.fits)
[INFO] .data(3681 x 40 x 40) (1e-20 erg / (Angstrom cm2 s)), .var(3681 x 40 x 40)
[INFO] center:(10:27:56.3962,04:13:25.3588) size:(8.000",8.000") step:(0.200",0.200") rot:-0.0 deg frame:FK5
[INFO] wavelength: min:4749.89 max:9349.89 step:1.25 Angstrom

In [9]: cube[3000:,10:20,25:40].info()
[INFO] 681 x 10 x 15 Cube (sdetect/minicube.fits)
[INFO] .data(681 x 10 x 15) (1e-20 erg / (Angstrom cm2 s)), .var(681 x 10 x 15)
[INFO] center:(10:27:55.3962,04:13:25.1893) size:(2.000",3.000") step:(0.200",0.200") rot:-0.0 deg frame:FK5
[INFO] wavelength: min:8499.89 max:9349.89 step:1.25 Angstrom

Likewise, Cube[k,:,:] returns an Image, as demonstrated below:

In [10]: ima1 = cube[1000, :, :]

In [11]: plt.figure()
Out[11]: <matplotlib.figure.Figure at 0x7f5811cdd7d0>

In [12]: ima1.plot(colorbar='v', title = '$\lambda$ = %.1f (%s)' %(cube.wave.coord(1000), cube.wave.unit))
Out[12]: <matplotlib.image.AxesImage at 0x7f5811e79090>

In [13]: ima2 = cube[3000, :, :]

In [14]: plt.figure()
Out[14]: <matplotlib.figure.Figure at 0x7f5811cdd8d0>

In [15]: ima2.plot(colorbar='v', title = '$\lambda$ = %.1f (%s)' %(cube.wave.coord(3000), cube.wave.unit))
Out[15]: <matplotlib.image.AxesImage at 0x7f581d7f2310>

In [16]: plt.figure()
Out[16]: <matplotlib.figure.Figure at 0x7f5818dc0950>

In [17]: ima2[5:25, 15:35].plot(colorbar='v',title = 'Zoom $\lambda$ = %.1f (%s)' %(cube.wave.coord(3000), cube.wave.unit))
Out[17]: <matplotlib.image.AxesImage at 0x7f5819885990>
_images/ObjFormatIma1.png _images/ObjFormatIma2.png _images/ObjFormatZommIma2.png

In the Image objects extracted from the cube, Image[p1:p2,q1:q2] returns a sub-image, Image[p,q] returns the value of pixel (p,q), Image[p,q] = value sets value in Image.data[p,q], and Image[p1:p2,q1:q2] = array sets the corresponding part of Image.data.

Finally, Cube[:,p,q] returns a Spectrum:

In [18]: spe = cube[:, 8, 28]

In [19]: import astropy.units as u

In [20]: from mpdaf.obj import deg2sexa

In [21]: coord_sky = cube.wcs.pix2sky([8, 28], unit=u.deg)

In [22]: dec, ra = deg2sexa(coord_sky)[0]

In [23]: plt.figure()
Out[23]: <matplotlib.figure.Figure at 0x7f5818dc0a50>

In [24]: spe.plot(title = 'Spectrum ra=%s dec=%s' %(ra, dec))

In [25]: plt.figure()
Out[25]: <matplotlib.figure.Figure at 0x7f5818d7fed0>

In [26]: spe[1640:2440].plot(title = 'Zoom Spectrum ra=%s dec=%s' %(ra, dec))
_images/ObjFormatSpe.png _images/ObjFormatZoomSpe.png

Getters and setters

Cube.get_step, Image.get_step and Spectrum.get_step return the world-coordinate separations between pixels along each axis of a cube, image, or spectrum, respectively:

In [27]: cube.get_step(unit_wave=u.nm, unit_wcs=u.deg)
Out[27]: array([1.25000000e-01, 5.55555556e-05, 5.55555556e-05])

In [28]: ima1.get_step(unit=u.deg)
Out[28]: array([5.55555556e-05, 5.55555556e-05])

In [29]: spe.get_step(unit=u.angstrom)
Out[29]: 1.25

Cube.get_range returns the range of wavelengths, declinations and right ascensions in a cube. Similarly, Image.get_range returns the range of declinations and right ascensions in an image, and Spectrum.get_range returns the range of wavelengths in a spectrum, as demonstrated below:

In [30]: cube.get_range(unit_wave=u.nm, unit_wcs=u.deg)
array([474.9890625 ,  10.46458228,  63.35455983, 934.9890625 ,
        10.46674895,  63.35676316])

In [31]: ima1.get_range(unit=u.deg)
Out[31]: array([10.46458228, 63.35455983, 10.46674895, 63.35676316])

In [32]: spe.get_range(unit=u.angstrom)
Out[32]: array([4749.890625, 9349.890625])

The get_start and get_end methods of cube, image and spectrum objects, return the world-coordinate values of the first and last pixels of each axis:

In [33]: print(cube.get_start(unit_wave=u.nm, unit_wcs=u.deg), cube.get_end(unit_wave=u.nm, unit_wcs=u.deg))
[474.9890625   10.46458228  63.35676315] [934.9890625   10.46674895  63.35455983]

In [34]: print(ima1.get_start(unit=u.deg), ima2.get_end(unit=u.deg))
[10.46458228 63.35676315] [10.46674895 63.35455983]

In [35]: print(spe.get_start(unit=u.angstrom), spe.get_end(unit=u.angstrom))
4749.890625 9349.890625

Note that when the declination axis is rotated away from the vertical axis of the image, the coordinates returned by get_start and get_end are not the minimum and maximum coordinate values within the image, so they differ from the values returned by get_range.

Cube.get_rot and Image.get_rot return the rotation angle of the declination axis to the vertical axis of the images within these objects:

In [36]: cube.get_rot(unit=u.deg)
Out[36]: -0.0

In [37]: ima1.get_rot(unit=u.rad)
Out[37]: -0.0

Updated flux and variance values can be assigned directly to the O.data and O.var attributes, respectively. Similarly, elements of the data can be masked or unmasked by assigning True or False values to the corresponding elements of the O.mask attribute. Changes to the spatial world coordinates must be performed using the set_wcs method:

In [38]: ima2.data[0:10,0:10] = 0

In [39]: ima2.mask[0:10,0:10] = True

In [40]: plt.figure()
Out[40]: <matplotlib.figure.Figure at 0x7f5811dffa10>

In [41]: ima2.plot()
Out[41]: <matplotlib.image.AxesImage at 0x7f58191e9290>