# 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.44937493,01:20:00.43755249) 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/cm^{2}/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:

(see Spectrum, Image and Cube format for more information).

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 0x7fe703df4e50>
```

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)
[INFO] loop_spe_multiprocessing (poly_spec): 200 tasks
```

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 0x7fe7014cb880>
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 0x7fe70174f250>
```

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 0x7fe7015b5c10>
```

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.44937493,01:20:00.43755249) 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.
In [52]: sp1.mask_region(9050, 9125)
# 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 [54]: sp1.unmask()
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
....: sp.mask_region(9050, 9125)
....: cont = sp.poly_spec(1)
....: sp.unmask()
....: 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 [74]: ha_ew.mask_selection(np.where((ima1.data)<4000))
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 0x7fe7011dda90>
```

## 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)
[INFO] loop_ima_multiprocessing (gaussian_filter): 10 tasks
```

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 0x7fe70132d160>
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 0x7fe70102e250>
```

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)
[INFO] loop_ima_multiprocessing (<function remove_background_gradient at 0x7fe701075b80>): 10 tasks
```

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 0x7fe700fd2df0>
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 0x7fe700f26730>
```

## 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 0x7fe7009ba940>
```

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.