# Cube object¶

MPDAF Cube objects contain a series of images at neighboring wavelengths. These are stored as a 3D array of pixel values. The wavelengths of the images increase along the first dimension of this array, while the Y and X axes of the images lie along the second and third dimensions. In addition to the data array, Cubes contain a WaveCoord object that describes the wavelength scale of the spectral axis of the cube, and a WCS object that describes the spatial coordinates of the spatial axes. Optionally a 3D array of variances can be provided to give the statistical uncertainties of the fluxes. Finally a mask array is provided for indicating bad pixels.

The fluxes and their variances are stored in numpy masked arrays, so virtually all numpy and scipy functions can be applied to them.

Preliminary imports:

In [1]: import numpy as np

In [2]: import matplotlib.pyplot as plt

In [3]: import astropy.units as u

In [4]: from mpdaf.obj import Cube, WCS, WaveCoord


## Cube creation¶

There are two common ways to obtain a Cube object:

• A cube can be created from a user-provided 3D array of pixel values, or from both a 3D array of pixel values, and a corresponding 3D array of variances. These can be simple numpy arrays, or they can be numpy masked arrays in which bad pixels have been masked. For example:

In [5]: wcs1 = WCS(crval=0, cdelt=0.2)

In [6]: wave1 = WaveCoord(cdelt=1.25, crval=4000.0, cunit=u.angstrom)

# numpy data array
In [7]: MyData = np.ones((400, 30, 30))

# cube filled with MyData
In [8]: cube = Cube(data=MyData, wcs=wcs1, wave=wave1) # cube 400X30x30 filled with data

In [9]: cube.info()
[INFO] 400 x 30 x 30 Cube (no name)
[INFO] .data(400 x 30 x 30) (no unit), no noise
[INFO] spatial coord (pix): min:(0.0,0.0) max:(5.8,5.8) step:(0.2,0.2) rot:-0.0 deg
[INFO] wavelength: min:4000.00 max:4498.75 step:1.25 Angstrom

• Alternatively, a cube can be read from a FITS file. In this case the flux and variance values are read from specific extensions:

# data and variance arrays read from the file (extension DATA and STAT)
In [10]: obj1 = Cube('obj/CUBE.fits')

In [11]: obj1.info()
[INFO] 1595 x 10 x 20 Cube (obj/CUBE.fits)
[INFO] .data(1595 x 10 x 20) (1e-40 erg / (Angstrom cm2 s)), .var(1595 x 10 x 20)
[INFO] center:(-30:00:00.4494,01:20:00.4376) size:(2.000",4.000") step:(0.200",0.200") rot:-0.0 deg frame:FK5
[INFO] wavelength: min:7300.00 max:9292.50 step:1.25 Angstrom


By default, if a FITS file has more than one extension, then it is expected to have a ‘DATA’ extension that contains the pixel data, and possibly a ‘STAT’ extension that contains the corresponding variances. If the file doesn’t contain extensions of these names, the “ext=” keyword can be used to indicate the appropriate extension or extensions.

The spatial and spectral world-coordinates of a Cube are recorded within WCS and WaveCoord objects, respectively. When a cube is read from a FITS file, these are automatically generated, based on FITS header keywords. Alternatively, when an cube is extracted from another cube, they are derived from the world coordinate information of the original cube.

In the above example, information about a cube is printed using the info method. In this case the information is as follows:

• The cube has dimensions, 1595 x 10 x 20, which means that there are 1595 images of different wavelengths, each having 10 x 20 spatial pixels.

• In addition to the pixel array (.data(1595 x 10 x 20)), there is also a variance array of the same dimensions (.var(1595 x 10 x 20)).

• The flux units of the pixels are 10-20 erg/s/cm2/Angstrom.

• The center of the field of view is at DEC: -30° 0’ 0.45” and RA: 1°20’0.437”.

• The size of the field of view is 2 arcsec x 4 arcsec. The pixel dimensions are 0.2 arcsec x 0.2 arcsec.

• The rotation angle of the field is 0°, which means that North is along the positive Y axis.

• The wavelength range is 7300-9292.5 Angstrom with a step of 1.25 Angstrom

The indexing of the cube arrays follows the Python conventions for indexing a 3D array. The pixel in the bottom-lower-left corner is referenced as [0,0,0] and the pixel [k,p,q] refers to the horizontal position q, the vertical position p, and the spectral position k, as follows:

The following example computes a reconstructed white-light image and displays it. The white-light image is obtained by summing each spatial pixel of the cube along the wavelength axis. This converts the 3D cube into a 2D image. The cube in this examples contains an observation of a single galaxy.

In [12]: ima1 = obj1.sum(axis=0)

In [13]: plt.figure(figsize=(8, 4))
Out[13]: <Figure size 800x400 with 0 Axes>

In [14]: ima1.plot(scale='arcsinh', colorbar='v')
Out[14]: <matplotlib.image.AxesImage at 0x7fbe8a02b9e8>


The next example computes the overall spectrum of the cube by taking the cube and summing along the X and Y axes of the image plane. This yields the total flux per spectral pixel.

In [15]: sp1 = obj1.sum(axis=(1,2))

In [16]: plt.figure()
Out[16]: <Figure size 640x480 with 0 Axes>

In [17]: sp1.plot()


## Loops over all spectra¶

The examples in this section will demonstrate how a procedure can be applied iteratively to the spectra of every image pixel of a cube. The goal of the examples will be to create a version of the above data-cube that has had the continuum background subtracted. For each image pixel, a low-order polynomial will be fitted to the spectrum of that pixel. This results in a polynomial curve that approximates the continuum spectrum of the pixel. This polynomial is then subtracted from the spectrum of that pixel, and the difference spectrum is recorded in a new output cube.

To illustrate the procedure, we start by fitting the continuum to the overall spectrum that was obtained in the previous example:

In [18]: plt.figure()
Out[18]: <Figure size 640x480 with 0 Axes>

In [19]: cont1 = sp1.poly_spec(5)

In [20]: sp1.plot()

In [21]: cont1.plot(color='r')


Next we do the same to a single pixel at the edge of the galaxy:

In [22]: plt.figure()
Out[22]: <Figure size 640x480 with 0 Axes>

In [23]: sp1 = obj1[:,5,2]

In [24]: sp1.plot()

In [25]: sp1.poly_spec(5).plot(color='r')


In principle, the above procedure could be performed to each pixel by writing a nested loop over the X and Y axes of the cube. However, instead of using two loops, one can use the spectrum iterator method, iter_spe of the Cube object. In the following example this is used to iteratively extract the six spectra of a small 2 x 3 pixel sub-cube, and determine their peak values:

In [26]: from mpdaf.obj import iter_spe

In [27]: small = obj1[:,0:2,0:3]

In [28]: small.shape
Out[28]: (1595, 2, 3)

In [29]: for sp in iter_spe(small):
....:     print(sp.data.max())
....:


Now let’s use the same approach to do the continuum subtraction procedure. We start by creating an empty datacube with the same dimensions as the original cube, but without variance information (using the clone method). Using two spectrum iterators we iteratively extract the spectra of each image pixel of the input cube and the empty output cube. At each iteration we then fit a polynomial spectrum to the input spectrum and record it in the output spectrum.

In [30]: cont1 = obj1.clone(data_init=np.empty, var_init=np.zeros)

In [31]: for sp, co in zip(iter_spe(obj1), iter_spe(cont1)):
....:     co[:] = sp.poly_spec(5)
....:


The result is a continuum datacube. Note that we have used the co[:] = sp.poly_spec(5) assignment rather than the more intuitive co = sp.poly_spec(5) assignment. The difference is that in python co=value changes the object that the co variable refers to, whereas co[:] changes the contents of the object that it currently points to. We want to change the contents of the spectrum in the output cube, so the latter is needed.

There is another way to compute the continuum datacube that can be much faster when used on a computer with multiple processors. This is to use the loop_spe_multiprocessing function. This uses multiple processors to apply a specified function to each spectrum of a cube and return a new cube that contains the resulting spectra:

In [32]: from mpdaf.obj import Spectrum

In [33]: cont2 = obj1.loop_spe_multiprocessing(f=Spectrum.poly_spec, deg=5)

[K
Completed 180 of 200 tasks (90%) in 5.0s
[K


To compare the results of the two methods, the following example sums the images of the two continuum cubes over the wavelength axis and displays the resulting white-light images of the continuum:

In [34]: rec1 = cont1.sum(axis=0)

In [35]: plt.figure()
Out[35]: <Figure size 640x480 with 0 Axes>

In [36]: rec1.plot(scale='arcsinh', colorbar='v', title='method 1')
Out[36]: <matplotlib.image.AxesImage at 0x7fbe8b02e3c8>

In [37]: rec2 = cont2.sum(axis=0)

In [38]: plt.figure()
Out[38]: <Figure size 640x480 with 0 Axes>

In [39]: rec2.plot(scale='arcsinh', colorbar='v', title='method2')
Out[39]: <matplotlib.image.AxesImage at 0x7fbe896f8080>


Next we subtract the continuum cube from the original cube to obtain a cube of the line emission of the galaxy. For display purposes this is then summed along the wavelength axis to yield an image of the sum of all of the emission lines in the cube:

In [40]: line1 = obj1 - cont1

In [41]: plt.figure()
Out[41]: <Figure size 640x480 with 0 Axes>

In [42]: line1.sum(axis=0).plot(scale='arcsinh', colorbar='v')
Out[42]: <matplotlib.image.AxesImage at 0x7fbe89693860>


Next we compute the equivalent width of the Hα emission in the galaxy. First we isolate the emission line by truncating the object datacube in wavelength:

In [43]: plt.figure()
Out[43]: <Figure size 640x480 with 0 Axes>

# Obtain the overall spectrum of the cube.
In [44]: sp1 = obj1.sum(axis=(1,2))

In [45]: sp1.plot()

# Obtain the spectral pixel indexes of wavelengths 9000 and 9200
In [46]: k1,k2 = sp1.wave.pixel([9000,9200], nearest=True)

# Extract a sub-cube restricted to the above range of wavelengths.
In [47]: emi1 = obj1[k1:k2+1,:,:]

In [48]: emi1.info()
[INFO] 161 x 10 x 20 Cube (obj/CUBE.fits)
[INFO] .data(161 x 10 x 20) (1e-40 erg / (Angstrom cm2 s)), .var(161 x 10 x 20)
[INFO] center:(-30:00:00.4494,01:20:00.4376) size:(2.000",4.000") step:(0.200",0.200") rot:-0.0 deg frame:FK5
[INFO] wavelength: min:9000.00 max:9200.00 step:1.25 Angstrom

# Obtain the overall spectrum of the above sub-cube.
In [49]: sp1 = emi1.sum(axis=(1,2))

# Plot the sub-spectrum in red over the original spectrum.
In [50]: sp1.plot(color='r')


Next we fit and subtract the continuum. Before doing the polynomial fit we mask the region of the emission lines (sp1.mask), so that the lines don’t affect the fit, and then we perform a linear fit between the continuum on either side of the masked region. Then the spectrum is unmasked and the continuum subtracted:

In [51]: plt.figure()
Out[51]: <Figure size 640x480 with 0 Axes>

# Mask the region containing the line emission.

# Fit a line to the continuum on either side of the masked region.
In [53]: cont1 = sp1.poly_spec(1)

# Unmask the region containing the line emission.

In [55]: plt.figure()
Out[55]: <Figure size 640x480 with 0 Axes>

In [56]: sp1.plot()

In [57]: cont1.plot(title="Fitted continuum")

In [58]: plt.figure()
Out[58]: <Figure size 640x480 with 0 Axes>

# Subtract the continuum from the spectrum to leave the line emission.
In [59]: line1 = sp1 - cont1

In [60]: line1.plot(title="Continuum subtracted")


Next we compute the total Hα line flux by simple integration (taking into account the pixel size in Angstrom) over the wavelength range centered around the Hα line and the continuum mean flux at the same location:

In [61]: plt.figure()
Out[61]: <Figure size 640x480 with 0 Axes>

# Find the spectral pixel index of the peak flux.
In [62]: k = line1.data.argmax()

In [63]: line1[55-10:55+11].plot(color='r')

# Integrate by summing pixels, multiplied by the pixel width.
In [64]: fline = (line1[55-10:55+11].sum()*line1.unit) * (line1.get_step(unit=line1.wave.unit)*line1.wave.unit)

# Obtain the mean continuum flux.
In [65]: cline = cont1[55-10:55+11].mean()*cont1.unit

# Compute the equivalent width of the line.
In [66]: ew = fline/cline

In [67]: print(fline, cline, ew)
[9224.61872579  507.99642777] 1e-40 erg / (cm2 s) [1936.55733796           inf] 1e-40 erg / (Angstrom cm2 s) [4.7634111 0.       ] Angstrom

In [68]: plt.close("all")


Finally we repeat this for all datacube spectra, and we save the Hα flux and equivalent width in two images. We start by creating two images with identical shapes and world-coordinates for the reconstructed image and then use the spectrum iterator iter_spe:

In [69]: ha_flux = ima1.clone(data_init=np.empty)

In [70]: cont_flux = ima1.clone(data_init=np.empty)

In [71]: ha_ew = ima1.clone(data_init=np.empty)

In [72]: for sp,pos in iter_spe(emi1, index=True):
....:     p,q = pos
....:     cont = sp.poly_spec(1)
....:     line = sp - cont
....:     fline = line[55-10:55+11].sum()[0] * line.get_step(unit=line.wave.unit)
....:     cline = cont[55-10:55+11].mean()[0]
....:     ew = fline/cline
....:     cont_flux[p,q] = cline
....:     ha_flux[p,q] = fline
....:     ha_ew[p,q] = ew
....:

In [73]: ha_flux.unit = sp.unit * sp.wave.unit

In [75]: ha_ew.unit = ha_flux.unit / cont_flux.unit

In [76]: fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16, 4))

In [77]: cont_flux.plot(ax=ax1, title="continuum mean flux (%s)"%cont_flux.unit, colorbar='v')
....: ha_flux.plot(ax=ax2, title="H$\\alpha$ line total flux (%s)"%ha_flux.unit, colorbar='v')
....: ha_ew.plot(ax=ax3, title="H$\\alpha$ line ew (%s)"%ha_ew.unit, colorbar='v')
....:
Out[77]: <matplotlib.image.AxesImage at 0x7fbe895d6eb8>


## Loops over all images¶

In this section, we demonstrate ways to perform of procedure iteratively to each image of a cube. We consider the datacube as a collection of monochromatic images and we process each of them. For each monochromatic image we apply a convolution by a gaussian kernel.

# Data and variance arrays read from the file (extension DATA and STAT)
In [78]: cube = Cube('obj/Central_Datacube_bkg.fits')


First, we use the image iterator iter_ima, which operates similarly to the spectrum iterator described earlier on this page, except that it returns successive images instead of successive spectra:

In [79]: from mpdaf.obj import iter_ima

In [80]: cube2 = cube.clone(data_init=np.empty, var_init=np.empty)

In [81]: for ima,k in iter_ima(cube, index=True):
....:     cube2[k,:,:] = ima.gaussian_filter(sigma=3)
....:


We can also make use of multiple processors to do the same thing faster, by using the loop_ima_multiprocessing method. This applies a specified procedure to all images within a cube and returns a new cube of the processed images:

In [82]: from mpdaf.obj import Image

In [83]: cube2 = cube.loop_ima_multiprocessing(f=Image.gaussian_filter, sigma=3)


We then plot the results:

In [84]: plt.figure()
Out[84]: <Figure size 640x480 with 0 Axes>

In [85]: cube.sum(axis=0).plot(title='before Gaussian filter', scale='log')
Out[85]: <matplotlib.image.AxesImage at 0x7fbe894f4208>

In [86]: plt.figure()
Out[86]: <Figure size 640x480 with 0 Axes>

In [87]: cube2.sum(axis=0).plot(title='after Gaussian filter', scale='log')
Out[87]: <matplotlib.image.AxesImage at 0x7fbe894c0278>


Next we will use the loop_ima_multiprocessing method to fit and remove a background gradient from a simulated datacube. For each image of the cube, we fit a 2nd order polynomial to the background values (selected here by simply applying a flux threshold to mask all bright objects). We do so by doing a chi^2 minimization over the polynomial coefficients using the numpy recipe np.linalg.lstsq(). For this, we define a function that takes an image as its sole parameter and returns a background-subtracted image:

In [88]: def remove_background_gradient(ima):
....:     ksel = np.where(ima.data.data<2.5)
....:     pval = ksel[0]
....:     qval = ksel[1]
....:     zval = ima.data.data[ksel]
....:     degree = 2
....:     Ap = np.vander(pval,degree)
....:     Aq = np.vander(qval,degree)
....:     A = np.hstack((Ap,Aq))
....:     (coeffs,residuals,rank,sing_vals) = np.linalg.lstsq(A,zval)
....:     fp = np.poly1d(coeffs[0:degree])
....:     fq = np.poly1d(coeffs[degree:2*degree])
....:     X,Y = np.meshgrid(list(range(ima.shape[0])), list(range(ima.shape[1])))
....:     ima2 = ima - np.array(list(map(lambda q,p: fp(p)+fq(q),Y,X)))
....:     return ima2
....:


We can then create the background-subtracted cube:

In [89]: cube2 = cube.loop_ima_multiprocessing(f=remove_background_gradient)


Finally, we compare the results for one of the slices:

In [90]: plt.figure()
Out[90]: <Figure size 640x480 with 0 Axes>

In [91]: cube[5,:,:].plot(vmin=-1, vmax=4)
Out[91]: <matplotlib.image.AxesImage at 0x7fbe8944ca90>

In [92]: plt.figure()
Out[92]: <Figure size 640x480 with 0 Axes>

In [93]: cube2[5,:,:].plot(vmin=-1, vmax=4)
Out[93]: <matplotlib.image.AxesImage at 0x7fbe8940eb38>


## Creation of images¶

The Cube class contains several methods to create images from a specific part of the cube. get_image extracts an image as the mean or sum of a given wavelenth band. bandpass_image sums the images of a cube after multiplying the cube by a given spectral bandpass curve.

get_band_image can be used to create an image with the filters included in MPDAF:

In [94]: cube = Cube('sdetect/minicube.fits')

In [95]: imV = cube.get_band_image('Johnson_V')

In [96]: imR = cube.get_band_image('Cousins_R')

In [97]: imI = cube.get_band_image('Cousins_I')

In [98]: fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))

In [99]: imV.plot(ax=ax1, title="Johnson_V")
....: imR.plot(ax=ax2, title="Cousins_R")
....: imI.plot(ax=ax3, title="Cousins_I")
....:
Out[99]: <matplotlib.image.AxesImage at 0x7fbe88f83198>


The filters are taken from the filter_list.fits file from the MUSE DRS:

In [100]: plt.figure()
Out[100]: <Figure size 640x480 with 0 Axes>

In [101]: from astropy.io import fits

In [102]: with fits.open('../obj/filters/filter_list.fits') as hdul:
.....:     fig, axes = plt.subplots(2, 2, figsize=(12, 12))
.....:     axes = axes.ravel()
.....:     for hdu in hdul[1:]:
.....:         if hdu.name.startswith('SDSS'):
.....:             ax = axes[1]
.....:         elif hdu.name.startswith('ACS'):
.....:             ax = axes[2]
.....:         elif hdu.name.startswith('WF'):
.....:             ax = axes[3]
.....:         else:
.....:             ax = axes[0]
.....:         ax.plot(hdu.data['lambda'], hdu.data['throughput'], label=hdu.name)
.....:     for ax in axes:
.....:         ax.legend()
.....:


## Sub-cube extraction¶

Warning

To be written.

select_lambda returns the sub-cube corresponding to a wavelength range.

subcube extracts a sub-cube around a position.

aperture

subcube_circle_aperture extracts a sub-cube from an circle aperture of fixed radius.