Masking a 3-D cube based on 5 sigma detection of the peak¶
Although NDF cubes should already have a variance array present, allowing you to easily mask the cube based on the SNR in each channel, its sometimes often useful to mask a cube based on the detection of the peak pixel. This quick how-to shows you how to do this in Starlink.
First of all, as always, you need to initialise Starlink and set up the individual package we will use. In this case, that is KAPPA. (Please note that if you use the tcsh shell instead of the bash shell these commands will be slightly different.) You should replace /star in the following commands with the location of your Starlink package, e.g. ~/star-2018A.
mkdir CubeAnalysis
cd CubeAnalysis
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
kappa
As we're using Starlink within a jupyter notebook, we don't want it to prompt us for input, so we will set the environmental variable ADAM_NOPROMPT to any value. If you're running this in an interactive bash shell you don't need to set this, but if you write shell scripts you don't want to interact with please set this variable at the top of the script.
export ADAM_NOPROMPT=1
For the purposes of this how-to, we are using a reduced HARP spectral cube named g34-3-12co_trim.sdf located in the ../Data directory (relative to our current directory). First of all, we are going to collapse the cube across the velocity axis, leaving the peak emission from each spectra. Currently, we don't care which velocity the peak occured at. We use the command collapse with the estimator=max option.
collapse in=../Data/g34-3-12co_small_trim.sdf axis=VRAD estimator=max out=g34-3-12co_peakcoll.sdf WLIM=0.0
Now lets see what this map looks like, using the KAPPA command display, with the default options. We can bring it up in a separate XWindows window:
display g34-3-12co_peakcoll.sdf dev=xw mode=faint style=!
lutwarm device=PNG
palentry 0 White device=PNG
palentry 1 Black device=PNG
And lets also have a black border and numeric labels, as well as a light grey tick marks. We'll alter these by creating a style file:
echo 'border=1' > style.dat
echo 'color(border) = black' >>style.dat
echo 'color(numlab) = black' >>style.dat
echo 'color(ticks) = grey' >>style.dat
We'll use this to plot our png like so:
display g34-3-12co_peakcoll.sdf dev=peakcoll.png/PNG style=^style.dat mode=faint
errclip in=g34-3-12co_peakcoll.sdf limit=5 mode=SNR out=g34-3-12co_5sigpeak.sdf
display g34-3-12co_5sigpeak.sdf dev=peakcoll_clip.png/PNG style=^style.dat mode=faint
ndftrace ../Data/g34-3-12co_small_trim.sdf
The part we're interested in is the pixel bounds:
ndftrace ../Data/g34-3-12co_small_trim.sdf|grep 'Pixel bounds'
If you were doing this in a script, you could get the values programatically using parget) as BASH arrays:
LOWER_BOUND=(`parget lbound ndftrace`)
UPPER_BOUND=(`parget ubound ndftrace`)
And you can then get the third element of the array like so:
LBOUND=${LOWER_BOUND[2]}
UBOUND=${UPPER_BOUND[2]}
This would let you grow your cube using manic into a cube pixel-aligned with the original data and having the same dimension.
Manic will do this by stacking identical copies of our 2 dimensional image repeatedly in the 3rd dimension. This means it will have the correct bad-pixel mask to blank out all pixels in our input cube that do not have 5$\sigma$ detection of the peak pixel.
To call manic, we specify that we want the 1st and 2nd axis to remain the same, and the third one to be a number of copies of the same image via axes=[1,2,0]. We then set the upper bound and lower bound parameters to be the same pixel bounds as the original data using ubound and lbound.
manic g34-3-12co_5sigpeak.sdf axes=[1,2,0] lbound=$LBOUND ubound=$UBOUND out=5sigpeakmask_3d.sdf
Finally, we use the command copybad to copy the bad pixel mask from our mask to our original data:
copybad in=../Data/g34-3-12co_small_trim.sdf ref=5sigpeakmask_3d.sdf out=g34-3-12co_5sigpeakmasked.sdf
You could now examine this cube in Gaia if you want to see the output cube.