# 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.

In [1]:
mkdir CubeAnalysis
cd CubeAnalysis
source $STARLINK_DIR/etc/profile  In [2]: kappa   KAPPA commands are now available -- (Version 2.5-8) Type kaphelp for help on KAPPA commands. Type 'showme sun95' to browse the hypertext documentation. See the 'Release Notes' section of SUN/95 for details of the changes made for this release.  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. In [3]: 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. In [4]: collapse in=../Data/g34-3-12co_small_trim.sdf axis=VRAD estimator=max out=g34-3-12co_peakcoll.sdf WLIM=0.0   Collapsing pixel axis 3 from pixel -65 to pixel 121 inclusive...  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: In [5]: display g34-3-12co_peakcoll.sdf dev=xw mode=faint style=!  Data will be scaled from 1.20042586326599 to 22.0331649780273.  Or we can plot a png that's easy to include in the notebook session . Lets change some of the default attributes though, so it looks nicer in display. We want black text on a white background, so we will adjust the palette with palentry and lets use the warm color table with lutwarm In [6]: 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: In [7]: 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: In [8]: display g34-3-12co_peakcoll.sdf dev=peakcoll.png/PNG style=^style.dat mode=faint  Data will be scaled from 1.20042586326599 to 22.0331649780273.  Now lets mask out these points based on a 5$\sigma$detection in the peak. We could do this using the KAPPA command makesnr and thresh, but even easier is to use the errclip to do it one go. In [9]: errclip in=g34-3-12co_peakcoll.sdf limit=5 mode=SNR out=g34-3-12co_5sigpeak.sdf   Applying a lower limit on signal-to-noise ratio. 3846 pixels had signal-to-noise ratios less than 5 in "/export/data/sgraves/Tutorials/AnalysisHowTos/CubeAnalysis/g34-3-12co_peakcoll ".  In [10]: display g34-3-12co_5sigpeak.sdf dev=peakcoll_clip.png/PNG style=^style.dat mode=faint  Data will be scaled from 1.28564727306366 to 23.3510208129883.  We now want to apply this mask to our original cube. We first 'grow' this 2-D mask to be the same shape as our original cube, using manic First we get the shape of the original cube using ndftrace In [11]: ndftrace ../Data/g34-3-12co_small_trim.sdf   NDF structure ../Data/g34-3-12co_small_trim: Title: Galactic plane 30.0 Label: TA* corrected antenna temperature Units: K Shape: No. of dimensions: 3 Dimension size(s): 171 x 139 x 187 Pixel bounds : -2625:-2455, 11:149, -65:121 Total pixels : 4444803 Data Component: Type : _REAL Storage form: SIMPLE Bad pixels may be present Variance Component: Type : _REAL Storage form: SIMPLE Bad pixels may be present World Co-ordinate Systems: Number of co-ordinate Frames: 5 Current co-ordinate Frame (Frame 5): Frame title : "3-d compound coordinate system" Domain : SKY-DSBSPECTRUM First pixel centre : 34.3675, 0.0170, -30.902 Axis 1: Label : Galactic longitude Units : degrees Nominal Pixel scale: 5.96741 arc-sec Axis 2: Label : Galactic latitude Units : degrees Nominal Pixel scale: 5.98365 arc-sec Axis 3: Label : Radio velocity (LSB) Units : km/s Nominal Pixel scale: 1 km/s Extensions: FITS <_CHAR*80> History Component: Created : 2012 Oct 13 16:06:57 No. records: 28 Last update: 2019 Oct 26 00:48:52 (NDFCOPY (KAPPA 2.6-4)) Update mode: NORMAL  The part we're interested in is the pixel bounds: In [12]: ndftrace ../Data/g34-3-12co_small_trim.sdf|grep 'Pixel bounds'   Pixel bounds : -2625:-2455, 11:149, -65:121  If you were doing this in a script, you could get the values programatically using parget) as BASH arrays: In [13]: LOWER_BOUND=(parget lbound ndftrace) UPPER_BOUND=(parget ubound ndftrace)  And you can then get the third element of the array like so: In [14]: 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. In [15]: 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:

In [16]:
copybad in=../Data/g34-3-12co_small_trim.sdf ref=5sigpeakmask_3d.sdf out=g34-3-12co_5sigpeakmasked.sdf

  There were 719202 bad pixels copied to the output NDF
'/export/data/sgraves/Tutorials/AnalysisHowTos/CubeAnalysis/g34-3-12co_5sigpeak