PixTable object¶
In the reduction approach of the MUSE pipeline, data need to be kept un-resampled until the very last step. The pixel tables used for this purpose can be saved at each intermediate reduction step and hence contain lists of pixels together with output coordinates and values. The pixel tables values and units change according to the reduction step. Please consult the data reduction user manual for further informations.
The PixTable
python object is used to handle the MUSE pixel tables created by
the data reduction system. The PixTable object can be read and write to disk
and a few functions can be performed on the object. Note that pixel tables are
saved as FITS binary tables or as multi-extension FITS images (data reduction
software version 0.08 or above). PixTable python object detects the file format
and acts accordingly. But by default PixTable object writes to disk using the
multi-extension FITS images option because it’s faster. A 8.2 GB pixel table
can be saved 3x faster as FITS images and loaded 5x faster compared to saving
and loading as FITS table.
Also pixtable can be very large and then used a lot of RAM. To use efficiently the memory, the file is open in memory mapping mode when a PixTable object is created from an input FITS file: i.e. the arrays are not in memory unless they are used by the script.
PixTable format¶
Array | Type | Description | Units |
---|---|---|---|
xpos | float | x position of the pixel within the field of view | pixel, rad, deg |
ypos | float | y position of the pixel within the field of view | pixel, red, deg |
lambda | float | wavelength assigned to the pixel | Angstrom |
data | float | data value | count, 10**-20 erg/s/cm**2/A |
dq | int | 32bit bad pixel status (in Euro3D convention) | |
stat | float | data variance estimation | count**2, (10**-20 erg/s/cm**2/A)**2 |
origin | int | pixel location on detector, slice and channel number |
The origin column is composed of the IFU and slice numbers and the x and
y coordinates on the originating CCD. Using bit shifting and information in the
FITS headers these four numbers are encoded in a single 32bit integer. Note
that the MUSE
package provides a Slicer
class
to convert the slicer number between various numbering schemes.
Read a pixtable, display informations and extract a smaller pixtable centered around an object¶
Preliminary imports:
In [1]: import numpy as np
In [2]: import matplotlib.pyplot as plt
In [3]: from mpdaf.drs import PixTable
We read the PixTable
from the disk and check its basic informations (info
) and FITS header content:
In [4]: pix = PixTable('PIXTABLE-MUSE.2014-07-26T04:37:08.541.fits')
In [5]: pix.info()
[INFO] 24 merged IFUs went into this pixel table
[INFO] This pixel table was flux-calibrated
[INFO] projected (intermediate) (Gnomonic proje)
Filename: PIXTABLE-MUSE.2014-07-26T04:37:08.541.fits
No. Name Type Cards Dimensions Format
0 PRIMARY PrimaryHDU 2232 ()
1 xpos ImageHDU 9 (1, 323885723) float32
2 ypos ImageHDU 9 (1, 323885723) float32
3 lambda ImageHDU 9 (1, 323885723) float32
4 data ImageHDU 9 (1, 323885723) float32
5 dq ImageHDU 8 (1, 323885723) int32
6 stat ImageHDU 9 (1, 323885723) float32
7 origin ImageHDU 8 (1, 323885723) int32
[INFO] None
In [6]: print(pix.nrows)
323885723
This is a pixtable containing a MUSE exposure of HDFS. Note that the current table has 323885723 pixels. It corresponds to the full MUSE field for a single exposure. Let’s look to the corresponding reconstructed image from the associated datacube:
In [7]: from mpdaf.obj import Image
In [8]: rec = Image('IMAGE-HDFS-v1.11.fits')
In [9]: rec.info()
[INFO] 331 x 326 Image (/muse/HDFS/public/dataproducts/HDFS-DataProducts-v1.11/IMAGE-HDFS-v1.11.fits)
[INFO] .data(331 x 326) (1e-20 ct), .var(331 x 326)
[INFO] center:(-60:33:49.0427,22:32:55.53) size in arcsec:(66.139,65.126) step in arcsec:(0.200,0.200) rot:0.1 deg
In [10]: plt.figure()
Out[10]: <matplotlib.figure.Figure at 0x7ff6d7b3a550>
In [11]: rec.plot(scale='arcsinh')
Out[11]: <matplotlib.image.AxesImage at 0x7ff6d7935910>
We are interested in the brightest object. Let’s find its position by using peak
:
In [12]: peak = rec.peak()
Now we will use extract
to extract a circular region of the pixtable centered around the object and we will restrict
the wavelength to the 6000:6100 Angstrom range:
In [13]: objpix = pix.extract(filename='Star_pixtable.fits',sky=(peak['y'], peak['x'], 2., 'C'), lbda=(6000,6100))
In [14]: objpix.nrows
Out[14]: 24441
Note that we have extracted a circular (‘C’) region of 2 arcsec around the object. The new pixtable (objpix) is much smaller, only 24441 pixels. The pixtable has been saved as a FITS file (Star_pixtable.fits).
The method extract
can extract a subset of a pixtable using the following criteria:
- aperture on the sky (center, size and shape),
- wavelength range,
- IFU numbers,
- slice numbers,
- detector pixels,
- exposure numbers,
- stack numbers.
extract
creates a mask columns for all criteria, merges the masks and returns a new pixtable extracted with the final mask.
These methods are also available to do the extraction step par step:
select_lambda
returns a mask corresponding to the given wavelength range,select_stacks
returns a mask corresponding to given stacks,select_slices
returns a mask corresponding to given slices,select_ifus
returns a mask corresponding to given ifus,select_exp
returns a mask corresponding to given exposures,select_xpix
andselect_ypix
return a mask corresponding to detector pixels,select_sky
returns a mask corresponding to the given aperture on the sky,extract_from_mask
returns a new pixtable extracted with the given mask.
Let’s investigate this pixtable.
get_xpos
/get_ypos
return the relative x/y position of the pixel to the center of the field of view.
We start by plotting the sky positions:
In [15]: x = objpix.get_xpos()
In [16]: y = objpix.get_ypos()
In [17]: plt.figure()
Out[17]: <matplotlib.figure.Figure at 0x7ff495869fd0>
In [18]: plt.plot(y, x, '.')
Out[18]: [<matplotlib.lines.Line2D at 0x7ff495a8bfd0>]
Ok, we have a circular location of pixels as expected. Note that the plotted points seems to be ‘thick’. We can check this by zooming. For example if we zoom to the two points on the left side, this what we obtain.
This is typical of the pixel table. Because of distortion each pixel on the
detector has not exactly the same location on the sky for the various
wavelength.
Let’s see if we have some bad pixel identified.
get_dq
gets the dq column.:
In [19]: dq = objpix.get_dq()
In [20]: k = np.where(dq > 0)
In [21]: k
Out[21]:
(array([ 9363, 13049, 14485, 14611, 14738, 15074, 15158, 15199, 15704,
15830, 21261, 21279]),)
In [22]: plt.plot(y[k], x[k], 'r.')
Out[22]: [<matplotlib.lines.Line2D at 0x7ff495ab68d0>]
Indeed there are 12 bad pixels localised in 6 areas of the detectors. We can see their location as the red points in the plot.
Let’s now investigate how this object is mapped on the detector. We start to get the origin array with get_origin
:
In [23]: origin = objpix.get_origin()
Several methods exists to decode it:
origin2ifu
returns the ifu number of each pixel,origin2slice
returns the slice number of each pixel,origin2xpix
returns the x coordinates of the pixels on the detector,origin2ypix
returns the y coordinates of the pixels on the detector,origin2coords
returns (ifu, slice, ypix, xpix).
For example we decode the origin array to get the IFU number:
In [24]: ifu = objpix.origin2ifu(origin)
In [25]: np.unique(ifu)
Out[25]: array([5, 6, 7], dtype=uint8)
In [26]: k = np.where(ifu == 5)
In [27]: plt.plot(y[k],x[k],'ob')
Out[27]: [<matplotlib.lines.Line2D at 0x7ff4957c0c50>]
In [28]: k = np.where(ifu == 6)
In [29]: plt.plot(y[k],x[k],'or')
Out[29]: [<matplotlib.lines.Line2D at 0x7ff495854850>]
In [30]: k = np.where(ifu == 7)
In [31]: plt.plot(y[k],x[k],'oc')
Out[31]: [<matplotlib.lines.Line2D at 0x7ff495854f10>]
We can see that the star is split into three IFUs (5, 6 and 7). We plot the sky location according to the IFU number.
Now we are going to display the data as located on the original exposure.
Firs we have to compute separately the corresponding pixtable for each IFU (extract
)
and then we use the sub-pixtable to reconstruct the originating CCD image (reconstruct_det_image
):
In [32]: objpix5 = pix.extract(filename='Star_pixtable.fits',sky=(peak['y'], peak['x'], 2., 'C'), lbda=(6000,6100), ifu=5)
In [33]: ima5 = objpix5.reconstruct_det_image()
In [34]: plt.figure()
Out[34]: <matplotlib.figure.Figure at 0x7ff495a670d0>
In [35]: ima5.plot(vmin=0, vmax=200)
Out[35]: <matplotlib.image.AxesImage at 0x7ff495712410>
In [36]: objpix6 = pix.extract(filename='Star_pixtable.fits',sky=(peak['y'], peak['x'], 2., 'C'), lbda=(6000,6100), ifu=6)
In [37]: ima6 = objpix6.reconstruct_det_image()
In [38]: plt.figure()
Out[38]: <matplotlib.figure.Figure at 0x7ff495854e90>
In [39]: ima6.plot(vmin=0, vmax=200)
Out[39]: <matplotlib.image.AxesImage at 0x7ff6d72dd090>
In [40]: objpix7 = pix.extract(filename='Star_pixtable.fits',sky=(peak['y'], peak['x'], 2., 'C'), lbda=(6000,6100), ifu=7)
In [41]: ima7 = objpix7.reconstruct_det_image()
In [42]: plt.figure()
Out[42]: <matplotlib.figure.Figure at 0x7ff49576b0d0>
In [43]: ima7.plot(vmin=0, vmax=200)
Out[43]: <matplotlib.image.AxesImage at 0x7ff4953c2cd0>
This give a good view of the pixels that comes into the object for the wavelength 6000:6100 Angstrom.
Note that we restricted the wavelength range in the extract
method.
It would be also possible to used reconstruct_det_waveimage
that reconstructs the image of wavelength values on the detector from the pixtable.
Use the pixtable’s data¶
We will see how to use the pixel table to fit a 2D gaussian for a restricted wavelength range.
We start to define a function that fit a 2D gaussian to a set of points (x, y, data):
In [44]: from scipy.optimize import leastsq
In [45]: def fitgauss(x, y, data, peak, center, fwhm):
....: p0 = np.array([peak, center[0], center[1], fwhm/2.355])
....: res = leastsq(gauss2D, p0, args=[x, y, data])
....: return res
....:
In [46]: def gauss2D(p, arglist):
....: x, y, data = arglist
....: peak, x0, y0, sigma = p
....: g = peak*np.exp(-((x-x0)**2 + (y-y0)**2)/(2*sigma**2))
....: residual = data - g
....: return residual
....:
Let’s check if it works:
In [47]: y, x = np.meshgrid(np.arange(10), np.arange(10))
In [48]: g = 2.0*np.exp(-((x-5)**2+(y-5)**2)/(2*1.7**2))
In [49]: gn = np.random.normal(g, 0.1*np.sqrt(g))
In [50]: xp = x.ravel()
In [51]: yp = y.ravel()
In [52]: gnp = gn.ravel()
In [53]: fitgauss(xp, yp, gnp, 1.0, (4.9,5.1), 2*2.355)
Out[53]: (array([ 1.97698318, 4.98185621, 4.94862697, 1.70022846]), 1)
OK, so now we can test it on our object pixtable:
In [54]: x = objpix.get_xpos()
In [55]: y = objpix.get_ypos()
In [56]: data = objpix.get_data()
In [57]: center = (0,0)
In [58]: res = fitgauss(y, x, data, data.max(), center, 0.7/3600.)
In [59]: print('Peak:',res[0][0], 'Center:',res[0][1:3], 'Fwhm:',res[0][3]*2.355*3600)
Peak: 1465.94006348 Center: [ 0. 0.] Fwhm: 0.7
We have used get_data
to have the data column.
It exists a getter and a setter for each column of the pixtable.
We recommend that you use these setters to update a pixtable because they preserve the consistency of the file by updating the FITS header.
In place of the relative coordinates, we can use the absolute position on the sky given by get_pos_sky
:
In [60]: y, x = objpix.get_pos_sky()
In [61]: center = (peak['y'], peak['x'])
In [62]: res = fitgauss(y, x, data, data.max(), center, 0.7/3600.)
In [63]: print('Peak:',res[0][0], 'Center:',res[0][1:3], 'Fwhm:',res[0][3]*2.355*3600)
Peak: 1465.94006348 Center: [ -60.56826963 338.23752675] Fwhm: 0.7
Self-calibration method for empty field¶
In this section, we will apply a self-calibration method on a single pixel table to bring all slices to the same median value. This will work on fields with small object, e.g. objects smaller than a slice length (15 arcsec).
The pixtable previously loaded contains a MUSE exposure of HDFS. This is a reduced pixtable produced by scipost, without sky subtraction.
We will mask out all bright continuum objects present in the FoV. We use a mask which has been produced by SExtractor on the corresponding white light image of this exposure.
mask_column
method returns a PixTableMask
object containing the mask as a new column.
We save this mask column as a FITS table:
In [64]: mask = pix.mask_column('Mask-HDF-110814.fits')
In [65]: mask.write('maskcol.fits')
Then, we estimate a reference sky spectrum from the masked pixel table:
In [66]: skyref = pix.sky_ref(pixmask=mask)
In [67]: skyref.write('skyref.fits')
sky_ref
returns a Spectrum
. Let’s look to it:
In [68]: skyref.plot()
This reference spectrum is used by the auto calibration method to normalise data values in each MUSE slice. In this example, we choose to use the additive correction:
In [69]: autocalib = pix.subtract_slice_median(skyref, pixmask=mask)
[INFO] pixtable PIXTABLE-MUSE.2014-07-26T04:37:08.541.fits updated
subtract_slice_median
is a python pixtable method but it has been coded in C for efficiency.
In the same way divide_slice_median
computes the median value for all slices and divides each slice by this correction to bring all slices to the same median value.
These methods return an PixTableAutoCalib
object that lets the user to save calibration information in a fits file:
In [70]: autocalib.write('autocalib-PIXTABLE-MUSE.2014-07-26T04:37:08.541.fits')
Finally, we save this corrected pixel table (write
):
In [71]: pix.write('corr-PIXTABLE-MUSE.2014-07-26T04:37:08.541.fits')
This non sky subtracted corrected pixtable can then be used to create a datacube with the appropriate pipeline recipe. Sky subtraction can then be performed with the zap software.
Warning
The use of OpenMP by MPDAF could hand an other process using OpenMP during the same python program. This problem only affects GCC; ICC is not affected. There is currently no workaround; the solution is to stop MPDAF before starting on other process using OpenMP.