2.5.5 3D measurements and spectra ¶

In the context of optical IFUs or radio IFUs in astronomy, a “Spectrum” is defined as separate measurements on each 2D slice of the 3D cube. Each 2D slice is defined by the first two FITS dimensions: the first FITS dimension is the horizontal axis and the second is the vertical axis. As with the tutorial on 2D image analysis (in Segmentation and making a catalog), let’s run Segment to see how it works in 3D. Like NoiseChisel above, to simplify the commands, let’s make an alias (3D detection with NoiseChisel):

$ alias astsegment-3d="astsegment \
           --config=/usr/local/etc/astsegment-3d.conf"

$ astsegment-3d det.fits --output=seg.fits

$ astscript-fits-view seg.fits

You see that we now have 3D clumps and 3D objects. So we can go ahead to do measurements. MakeCatalog can do single-valued measurements (as in 2D) on 3D datasets also. For example, with the command below, let’s get the flux-weighted center (in the three dimensions) and sum of pixel values. There isn’t usually a standard name for the third WCS dimension (unlike Ra/Dec). So in Gnuastro, we just call it --w3. With the second command, we are having a look at the first 5 rows. Note that we are not using -Y with asttable anymore because it the wavelength column would only be shown as zero (since it is in meters!).

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum --output=cat.fits

$ asttable cat.fits -h1 -O --txtf64p=5 --head=5
# Column 1: OBJ_ID [counter    ,i32,] Object identifier.
# Column 2: RA     [deg        ,f64,] Flux weighted center (WCS axis 1).
# Column 3: DEC    [deg        ,f64,] Flux weighted center (WCS axis 2).
# Column 4: AWAV   [m          ,f64,] Flux weighted center (WCS axis 3).
# Column 5: SUM    [input-units,f32,] Sum of sky subtracted values.
1  3.99677e+01   -1.58660e+00   4.82994e-07   7.311189e+02
2  3.99660e+01   -1.58927e+00   4.86411e-07   7.872681e+03
3  3.99682e+01   -1.59141e+00   4.90609e-07   1.314548e+03
4  3.99677e+01   -1.58666e+00   4.90816e-07   7.798024e+02
5  3.99659e+01   -1.58930e+00   4.93657e-07   3.255210e+03

Besides the single-valued measurements above (that are shared with 2D inputs), on 3D cubes, MakeCatalog can also do per-slice measurements. The options for these measurements are formatted as --*in-slice. With the command below, you can check their list:

$ astmkcatalog --help | grep in-slice
 --area-in-slice        [3D input] Number of labeled in each slice.
 --area-other-in-slice  [3D input] Area of other lab. in projected area.
 --area-proj-in-slice   [3D input] Num. voxels in '--sum-proj-in-slice'.
 --sum-err-in-slice     [3D input] Error in '--sum-in-slice'.
 --sum-in-slice         [3D input] Sum of values in each slice.
 --sum-other-err-in-slice   [3D input] Area in '--sum-other-in-slice'.
 --sum-other-in-slice   [3D input] Sum of other lab. in projected area.
 --sum-proj-err-in-slice   [3D input] Error of '--sum-proj-in-slice'.
 --sum-proj-in-slice    [3D input] Sum of projected area in each slice.

For every label and measurement, these options will give many values in a vector column (see Vector columns). Let’s have a look by asking for the sum of values and area of each label in each slice associated to each label with the command below. There is just one important point: in 3D detection with NoiseChisel, we ran NoiseChisel on the signal-to-noise image, not the continuum-subtracted image! So the values to use for the measurement of each label should come from the no-continuum.fits file (not seg.fits).

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum  \
               --area-in-slice --sum-in-slice --output=cat.fits \
               --valuesfile=no-continuum.fits

$ asttable -i cat.fits
--------
seg_cat.fits (hdu: 1)
-------          -----       ----          -------
No.Name          Units       Type          Comment
-------          -----       ----          -------
1  OBJ_ID        counter     int32         Object identifier.
2  RA            deg         float64       Flux wht center (WCS 1).
3  DEC           deg         float64       Flux wht center (WCS 2).
4  AWAV          m           float64       Flux wht center (WCS 3).
5  SUM           input-units float32       Sum of sky-subed values.
6  AREA-IN-SLICE counter     int32(3681)   Number of pix. in each slice.
7  SUM-IN-SLICE  input-units float32(3681) Sum of values in each slice.
--------
Number of rows: 211
--------

You can see that the new AREA-IN-SLICE and SUM-IN-SLICE columns have a (3681) in their types. This shows that unlike the single-valued columns before them, in these columns, each row has 3681 values (a “vector” column). If you are not already familiar with vector columns, please take a few minutes to read Vector columns. Since a MUSE data cube has 3681 slices, this is effectively the spectrum of each object.

Let’s find the object that corresponds to the H-alpha emission of the brightest galaxy (that we found in Viewing spectra and redshifted lines). That emission line was around 8565.93 Angstroms, so let’s look for the objects within \(\pm5\) Angstroms of that value (between 8560 to 8570 Angstroms):

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -cobj_id,ra,dec -Y
198    39.965897   -1.589279

From the command above, we see that at this wavelength, there was only one object. Let’s extract its spectrum by asking for the sum-in-slice column:

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 \
           -carea-in-slice,sum-in-slice

If you look into the outputs, you will see that it is a single line! It contains a long list of 0 values at the start and nan values in the end. If you scroll slowly, in the middle of each you will see some non-zero and non-NaN numbers. To help interpret this more easily, let’s transpose these vector columns (so each value of the vector column becomes a row in the output). We will use the --transpose option of Table for this (just note that since transposition changes the number of rows, it can only be used when your table only has vector columns and they all have the same number of elements (as in this case, for more):

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 \
           -carea-in-slice,sum-in-slice --transpose

We now see the measurements on each slice printed in a separate line (making it much more easier to visually read). However, without a counter, it is very hard to interpret them. Let’s pipe the output to a new Table command and use column arithmetic’s counter operator for displaying the slice number (see Size and position operators). Note that since we are piping the output, we also added -O so the column metadata are also passed to the new instance of Table:

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -O \
           -carea-in-slice,sum-in-slice --transpose \
           | asttable -c'arith $1 counter swap',2
...[[truncated]]...
3040   0       nan
3041   0       nan
3042   0       nan
3043   0       nan
3044   1       4.311140e-01
3045   18      3.936019e+00
3046   161    -5.800080e+00
3047   360     2.967184e+02
3048   625     1.912855e+03
3049   823     5.140487e+03
3050   945     7.174101e+03
3051   999     6.967604e+03
3052   1046    6.468591e+03
3053   1025    6.457354e+03
3054   996     6.599119e+03
3055   966     6.762280e+03
3056   873     5.014052e+03
3057   649     2.003334e+03
3058   335     3.167579e+02
3059   131     1.670975e+01
3060   25     -2.953789e+00
3061   0       nan
3062   0       nan
3063   0       nan
3064   0       nan
...[[truncated]]...

$ astscript-fits-view seg.fits

After DS9 opens with the last command above, go to slice 3044 (which is the first non-NaN slice in the spectrum above). In the OBJECTS extension of this slice, you see several non-zero pixels. The few non-zero pixels on the bottom have a label of 197 and the single non-zero pixel at a higher Y axis position has a label of 198 (which as we saw above, was the label of the H-alpha emission of this galaxy). The few 197 labeled pixels in this slice are the last voxels of the NII emission that is just blue-ward of H-alpha.

The single pixel you see in slice 3044 is why you see a value of 1 in the AREA-IN-SLICE column. As you go to the next slices, if you count the pixels, you will see they add up to the same number you see in that column. The values in the SUM-IN-SLICE are the sum of values in the continuum-subtracted cube for those same voxels. You should now be able to understand why the --sum-in-slice column has NaN values in all other slices: because this label doesn’t exist in any other slice! Also, within slices that contain label 198, this column only uses the voxels that have the label. So as you see in the second column above, the area that is used in each changes.

Therefore --sum-in-slice or area-in-slice are the raw 3D spectrum of each 3D emission-line. This is a different concept from the traditional “spectrum” where the same area is used over all the slices. To get that you should use the --sumprojinslice column of MakeCatalog. All the --*in-slice options that contain a proj in their name are measurements over the fixed “projection” of the 3D volume on the 2D surface of each slice. To see the effect, let’s also ask MakeCatalog to measure this projected sum column:

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum  \
               --area-in-slice --sum-in-slice --sum-proj-in-slice \
               --output=cat.fits --valuesfile=no-continuum.fits
$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -O \
           -carea-in-slice,sum-in-slice,sum-proj-in-slice \
           --transpose \
           | asttable -c'arith $1 counter swap',2,3
...[[truncated]]...
3040   0       nan            8.686357e+02
3041   0       nan            4.384907e+02
3042   0       nan            4.994813e+00
3043   0       nan           -1.595918e+02
3044   1       4.311140e-01  -2.793141e+02
3045   18      3.936019e+00  -3.251023e+02
3046   161    -5.800080e+00  -2.709914e+02
3047   360     2.967184e+02   1.049625e+02
3048   625     1.912855e+03   1.841315e+03
3049   823     5.140487e+03   5.108451e+03
3050   945     7.174101e+03   7.149740e+03
3051   999     6.967604e+03   6.913166e+03
3052   1046    6.468591e+03   6.442184e+03
3053   1025    6.457354e+03   6.393185e+03
3054   996     6.599119e+03   6.572642e+03
3055   966     6.762280e+03   6.716916e+03
3056   873     5.014052e+03   4.974084e+03
3057   649     2.003334e+03   1.870787e+03
3058   335     3.167579e+02   1.057906e+02
3059   131     1.670975e+01  -2.415764e+02
3060   25     -2.953789e+00  -3.534623e+02
3061   0       nan           -3.745465e+02
3062   0       nan           -2.532008e+02
3063   0       nan           -2.372232e+02
3064   0       nan           -2.153670e+02
...[[truncated]]...

As you see, in the new SUM-PROJ-IN-SLICE column, we have a measurement in each slice: including slices that do not have the label of 198 at all. Also, the area used to measure this sum is the same in all slices (similar to a classical spectrometer’s output).

However, there is a big problem: have a look at the sums in slices 3040 and 3041: the values are increasing! This is because of the emission in the NII line that also falls over the projected area of H-alpha. This shows the power of IFUs as opposed to classical spectrometers: we can distinguish between individual lines based on spatial position and do measurements in 3D!

Finally, in case you want the spectrum with the continuum, you just have to change the file given to --valuesfile:

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum  \
               --area-in-slice --sum-in-slice --sum-proj-in-slice \
               --valuesfile=a370-crop.fits \
               --output=cat-with-continuum.fits