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
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
In [2]:

     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]:

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.

Image of max-collapsed 12CO data

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

Image of max-collapsed and 5sigma clipped 12CO data

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

      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

      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]:

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

You could now examine this cube in Gaia if you want to see the output cube.