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
)
export STARLINK_DIR=/star
source $STARLINK_DIR/etc/profile
export ADAM_NOPROMPT=1
mkdir Clumpfinding
cd Clumpfinding
cupid
kappa
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.
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.
collapse in=12co_cube.sdf estimator='integ' wlim=0.0 out=12co_integ.sdf axis=VRAD low=18.0 high=120.0
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.
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
We will now use CUPID's findclumps
with the FellWalker algorithm to identify two dimensional clumps of emission. First we will initialise cupid:
cupid
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.
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
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.
display in=12co_integ.sdf dev=integ_clumps.eps/AVCPS style=^style.dat low=0.0 mode=faint
gdset integ_clumps.eps/AVCPS
Now we can outline the clumps using CUPID's outlineclump
routine.
outlineclump ndf=12co_integ_clumps.sdf index=1-27
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.
command convert integ_clumps.eps -trim integ_clumps.png
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:
hdstrace 12co_integ_clumps.sdf
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
:
stats 12co_integ_clumps.sdf.MORE.CUPID.CLUMPS\(1\)
lutwarm device=clump1.png/PNG
display 12co_integ_clumps.sdf.MORE.CUPID.CLUMPS\(1\) device=clump1.png/PNG style=^style.dat mode=faint
The pixel axes corresponds to the pixel values in the original file it was run on. The data values are the original data values.
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.