Using CUPID for clumpfinding

The Starlink package CUPID contains implementations of various algorithms for Clump Finding in any number of dimensions. Here we will show you its use in 2 and 3-D, using the FellWalker algorithm.

First of all we will need to set up Starlink, initialise the Cupid and KAPPA packages, and move into our working directory. Because this is being run in a jupyter notebook, we will also turn off prompting for user input. (You should replace the /star in the first line below with the location of your Starlink installation directory: e.g. ~/star-2018A)

In [1]:
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
export ADAM_NOPROMPT=1
mkdir Clumpfinding
cd Clumpfinding
cupid
kappa


     CUPID commands are now available -- (Version 2.5)

     Type cupidhelp for help on CUPID commands.
     Type 'showme sun255' to browse the hypertext documentation.

   


     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.

Clumpfinding in 2-D

First of all we will try out some 2-D clumpfinding.

We will use the same HARP 12CO G34-3 cube from COHRS that we have been using in other tutorials.

In [2]:
cp ../Data/g34-3-12co_small_trim.sdf 12co_cube.sdf

If we integrate this over the emission, we will then have a 2-D map of integrated emission.

In [3]:
collapse in=12co_cube.sdf estimator='integ' wlim=0.0 out=12co_integ.sdf axis=VRAD low=18.0 high=120.0

   Collapsing pixel axis 3 from pixel -16 to pixel 86 inclusive...

We can use display to view this. First of all we'll set up some useful defaults for plotting, and create a quick PNG plot of the data.

In [4]:
gdclear dev=integ.png/PNG
lutwarm dev=integ.png/PNG
palentry 0 White dev=integ.png/PNG
palentry 1 Black dev=integ.png/PNG
echo 'border=1' > style.dat
echo 'color(border) = black' >>style.dat
echo 'color(numlab) = black' >>style.dat
echo 'color(ticks) = grey' >>style.dat

display in=12co_integ.sdf dev=integ.png/PNG style=^style.dat low=0.0 mode=faint

Data will be scaled from 24.4844436645508 to 177.155227661133.

12CO 3-2 integrated intensity

We will now use CUPID's findclumps with the FellWalker algorithm to identify two dimensional clumps of emission. First we will initialise cupid:

In [5]:
cupid


     CUPID commands are now available -- (Version 2.5)

     Type cupidhelp for help on CUPID commands.
     Type 'showme sun255' to browse the hypertext documentation.

We will set Fellwalker to only look at clumps with a minimum height of 5 times the RMS (Fellwalker.minheight), to only split into a different clump if there is a dip of at least 5 times the RMS between peaks(Fellwalker.MINDIP), to run 2 cleaning iterations (FellWalker.CleanIter), and to follow the clumps down to 3 times the RMS (Fellwalker.NOISE) . The RMS will be automatically calculated from the variance array.

In [6]:
findclumps in=12co_integ.sdf out=12co_integ_clumps.sdf outcat=12co_integ_clumps.FIT method=FellWalker \
config='"FellWalker.minheight=5*RMS,FellWalker.MINDIP=5*RMS,FellWalker.CleanIter=2,FellWalker.NOISE=3*RMS"' accept


FellWalker:
162 clumps rejected because they contain too few pixels.
27 usable clumps found.


Configuration parameters:
   FELLWALKER.ALLOWEDGE=1
   FELLWALKER.CLEANITER=2
   FELLWALKER.DUMPPEAK=0
   FELLWALKER.DUMPWALK=-1
   FELLWALKER.FLATSLOPE=3.5391768339269
   FELLWALKER.FWHMBEAM=2
   FELLWALKER.MAXBAD=0.05
   FELLWALKER.MAXJUMP=4
   FELLWALKER.MINDIP=17.6958841696345
   FELLWALKER.MINHEIGHT=17.6958841696345
   FELLWALKER.MINPIX=7
   FELLWALKER.NOISE=10.6175305017807
   FELLWALKER.RMS=3.5391768339269

We can use display and outline clumps to produce a map of the clumps. We are using the acucmulating postscript graphics device for KAPPA here, so that it will first display the integrated map and then contour the clumps on top of it. Alternatively, you can use Gaia and contour the output file over the input file interactively.

First we create the base image and set the current graphics device to our output file and the AVCPS device.

In [7]:
display in=12co_integ.sdf dev=integ_clumps.eps/AVCPS style=^style.dat low=0.0 mode=faint
gdset integ_clumps.eps/AVCPS

Data will be scaled from 24.4844436645508 to 177.155227661133.

All output will go to file integ_clumps.eps by default. Each subsequent
graphics command will append output to any existing file with this name.

Now we can outline the clumps using CUPID's outlineclump routine.

In [8]:
outlineclump ndf=12co_integ_clumps.sdf index=1-27

Plotting clump index 1
Alignment has occurred within the PIXEL Domain.

Plotting clump index 2
Alignment has occurred within the PIXEL Domain.

Plotting clump index 3
Alignment has occurred within the PIXEL Domain.

Plotting clump index 4
Alignment has occurred within the PIXEL Domain.

Plotting clump index 5
Alignment has occurred within the PIXEL Domain.

Plotting clump index 6
Alignment has occurred within the PIXEL Domain.

Plotting clump index 7
Alignment has occurred within the PIXEL Domain.

Plotting clump index 8
Alignment has occurred within the PIXEL Domain.

Plotting clump index 9
Alignment has occurred within the PIXEL Domain.

Plotting clump index 10
Alignment has occurred within the PIXEL Domain.

Plotting clump index 11
Alignment has occurred within the PIXEL Domain.

Plotting clump index 12
Alignment has occurred within the PIXEL Domain.

Plotting clump index 13
Alignment has occurred within the PIXEL Domain.

Plotting clump index 14
Alignment has occurred within the PIXEL Domain.

Plotting clump index 15
Alignment has occurred within the PIXEL Domain.

Plotting clump index 16
Alignment has occurred within the PIXEL Domain.

Plotting clump index 17
Alignment has occurred within the PIXEL Domain.

Plotting clump index 18
Alignment has occurred within the PIXEL Domain.

Plotting clump index 19
Alignment has occurred within the PIXEL Domain.

Plotting clump index 20
Alignment has occurred within the PIXEL Domain.

Plotting clump index 21
Alignment has occurred within the PIXEL Domain.

Plotting clump index 22
Alignment has occurred within the PIXEL Domain.

Plotting clump index 23
Alignment has occurred within the PIXEL Domain.

Plotting clump index 24
Alignment has occurred within the PIXEL Domain.

Plotting clump index 25
Alignment has occurred within the PIXEL Domain.

Plotting clump index 26
Alignment has occurred within the PIXEL Domain.

Plotting clump index 27
Alignment has occurred within the PIXEL Domain.

We'll convert that to a PNG so we can display it in the notebook (requires ImageMagick's convert to be installed -- here we use the shell command to ensure we search the general $PATH for the command convert, instead of picking up the Starlink convert package from our setup scripts.

In [9]:
command convert integ_clumps.eps -trim integ_clumps.png

Clump outlines

Each of these clumps is stored in a extension inside the output NDF from our initial findclumps command. We can copy out an individual clump using the ndfcopy command in KAPPA. If you use the hdstrace command you can see the full structure of the NDF file and any extensions within it:

In [10]:
hdstrace 12co_integ_clumps.sdf

12CO_INTEG_CLUM  <NDF>

   DATA_ARRAY     <ARRAY>         {structure}
      ORIGIN(2)      <_INTEGER>      -2625,11
      DATA(171,139)  <_INTEGER>      19,19,19,19,19,19,19,19,19,19,19,14,14,
                                     18,18,18,18,18,*,*,*,*,*,*,18,18,18,18

   TITLE          <_CHAR*19>      'Galactic plane 30.0'
   LABEL          <_CHAR*57>      'T%<s60>+%<v30>+A%^50+%<20+*%+   corrected
ant...'
   WCS            <WCS>           {structure}
      DATA(173)      <_CHAR*32>      ' Begin FrameSet',' Nframe = 5',' Curr...'
                                     ... ' End...',' End CmpMap',' End FrameSet'

   HISTORY        <HISTORY>       {structure}
      CREATED        <_CHAR*24>      '2012-OCT-13 16:06:57.666'
      CURRENT_RECORD  <_INTEGER>     30
      UPDATE_MODE    <_CHAR*8>       'NORMAL'
      RECORDS(34)    <HIST_REC>      {array of structures}

      Contents of RECORDS(1)
         COMMAND        <_CHAR*30>      'MAKECUBE        (SMURF V1.5.0)'
         DATASET        <_CHAR*60>      '/export/data/visitors/mjc/acsis/jps...'
         DATE           <_CHAR*24>      '2012-OCT-13 16:06:57.776'
         HOST           <_CHAR*6>       'ikaika'
         TEXT(20)       <_CHAR*72>      'Parameters: ALIGNSYS=FALSE AUTOGRI...'
                                        ... '   a20110607_00055_01_bl002, a2...'
         USER           <_CHAR*3>       'mjc'

   MORE           <EXT>           {structure}
      FITS(253)      <_CHAR*80>      'SIMPLE  =                    T / file...'
                                     ... 'DATASUM = '3842980826'      ...','END'
      CUPID          <CUPID_EXT>     {structure}
         CLUMPS(27)     <CLUMP>         {array of structures}

         Contents of CLUMPS(1)
            PEAK1          <_DOUBLE>       34.254821430877
            PEAK2          <_DOUBLE>       0.15326273212263
            CEN1           <_DOUBLE>       34.259453896913
            CEN2           <_DOUBLE>       0.1480484312024
            SIZE1          <_DOUBLE>       31.83049371695
            SIZE2          <_DOUBLE>       23.425359869816
            SUM            <_DOUBLE>       29404.480075836
            PEAK           <_DOUBLE>       213.23951548951
            VOLUME         <_DOUBLE>       11746.23953899
            MODEL          <NDF>           {structure}
               DATA_ARRAY     <ARRAY>         {structure}
                  DATA(26,18)    <_DOUBLE>       *,*,*,*,*,*,*,*,*,*,*,*,*,
                                                 ... *,*,*,*,*,*,*,*,*,*,*,*
                  ORIGIN(2)      <_INTEGER>      -2572,79

         QUALITY_NAMES  <QUALITY_NAMES>   {structure}
            QUAL(8)        <IRQ_QUAL>      {array of structures}

            Contents of QUAL(1)
               NAME           <_CHAR*15>      'CLUMP'
               FIXED          <_LOGICAL>      FALSE
               VALUE          <_LOGICAL>      FALSE
               BIT            <_INTEGER>      7
               COMMENT        <_CHAR*50>      'set iff a pixel is within a c...'
               READONLY       <_LOGICAL>      FALSE
               FIXED_BIT      <_LOGICAL>      FALSE

            LAST_USED      <_INTEGER>      3
            NFREE          <_INTEGER>      5
            FREE(8)        <_INTEGER>      8,7,6,5,4,0,0,0
            VALID          <_LOGICAL>      TRUE

         CONFIG(13)     <_CHAR*38>      'FELLWALKER.ALLOWEDGE=1','FELLWALKE...'
                                        'FELLWALKER.RMS=3.5391768339269'

   QUALITY        <QUALITY>       {structure}
      QUALITY        <ARRAY>         {structure}
         DATA(171,139)  <_UBYTE>        64,64,64,64,64,64,64,64,64,64,64,64,
                                        ... 128,128,128,128,128,128,96,64,64,64
         ORIGIN(2)      <_INTEGER>      -2625,11
         BAD_PIXEL      <_LOGICAL>      FALSE

End of Trace.

Here the CLUMPS are contained in an extension MORE.CUPID.CLUMPS, and each individual clump is then found as e.g. MORE.CUPID.CLUMPS\(3\) for the 3rd CLUMP. (The backslashes, \ are used here to escape the brackets so they are not interpreted by the shell we are typing this in.)

Each of these clumps is itself a valid NDF, and we can use all our normal tools to look at it. E.g., stats and display:

In [11]:
stats 12co_integ_clumps.sdf.MORE.CUPID.CLUMPS\(1\)

   Pixel statistics for the NDF structure
/export/data/sgraves/Tutorials/AnalysisHowTos/Clumpfinding/12co_integ_clumps.MO
RE.CUPID.CLUMPS(1).MODEL

      Title                     : <undefined>
      NDF array analysed        : DATA

         Pixel sum              : 29404.4800758362
         Pixel mean             : 89.3753193794412
         Standard deviation     : 30.0751660390696
         Skewness               : 1.45206166742115
         Kurtosis               : 2.43735732429543
         Minimum pixel value    : 43.6652069091797
            At pixel            : (-2554, 85)
            Co-ordinate         : (-2554.5, 84.5)
         Maximum pixel value    : 205.800048828125
            At pixel            : (-2557, 93)
            Co-ordinate         : (-2557.5, 92.5)
         Total number of pixels : 468
         Number of pixels used  : 329 (70.3%)
         No. of pixels excluded : 139 (29.7%)
In [12]:
lutwarm device=clump1.png/PNG
display 12co_integ_clumps.sdf.MORE.CUPID.CLUMPS\(1\) device=clump1.png/PNG style=^style.dat mode=faint

Data will be scaled from 59.2543419777079 to 300.222161191575.

Image of brightest clump found

The pixel axes corresponds to the pixel values in the original file it was run on. The data values are the original data values.

Clumpfinding in 3D

3D clumpfinding works exactly the same as clumpfinding in 2D. However, it can be harder to visualise the clumps. We recommend using Gaia to interactively visualise the clumps where you can scroll through the frequencies within one clump.