Chapter 9
Advanced Analysis

 9.1 Changing co-ordinate frames
 9.2 Hybrid Data
 9.3 Position-velocity diagram
 9.4 Creating channel maps
 9.5 Clump finding

9.1 Changing co-ordinate frames

Co-ordinate systems of regridded Starlink files are described by the world co-ordinate system, (WCS). There are a number of basic co-ordinate systems (or frames) common to all NDFs; these are described by their DOMAIN (essentially their name) and are PIXEL, AXIS, GRID and FRACTION. Additional frames can be stored in the WCS component. Two common additional frames are SKY and SPECTRUM. SKY refers to two-dimensional frames (known as SkyFrames) which describe sky co-ordinates. SPECTRUM refers to one-dimensional frames (known as SpecFrames) that describe a position within a spectrum.

The choice of co-ordinates within the SkyFrame or SpecFrame is called the system. For example, for SkyFrames this may be Galactic or FK5, while for SpecFrames this may be velocity or frequency.

The current frame and the system can both be changed using the Kappa applications wcsframe and wcsattrib. The current frame for data processed with makecube (either manually or by the pipeline) typically has a DOMAIN of SKY-DSBSPECTRUM. The compound name refers to the first two axes of the frame having sky co-ordinates and the third axis having dual side-band spectral units. You can check this with ndftrace:

  % ndftrace file | grep Domain

You can change the attributes of a cube using wcsattrib. To change the SkyFrame co-ordinate system set the System attribute.

  % wcsattrib file set system galactic

To change between velocity and frequency you also use the System attribute. For SKY-DSBSPECTRUM the software knows which axis is being referred to based on the option you supply (the third in this case).

  % wcsattrib file set system freq
  % wcsattrib file set system vrad

Once a cube has been collapsed and is two-dimensional its frame becomes SKY. If you need to change frames manually use wcsframe.

  % wcsframe file sky

For offset co-ordinates you will need to set the system skyrefis to origin. The final two lines in the example below convert the offset units to arcseconds instead of radians.

  % wcsattrib file set skyrefis origin
  % wcsattrib file set ’format(1)’ ’s.*’
  % wcsattrib file set ’format(2)’ ’s.*’

9.2 Hybrid Data

To merge two hybrid observations run wcsalign on both files.

  % wcsalign in="’rawfile_01.sdf,rawfile_02.sdf’" insitu lbnd=! ubnd=! ref=!

Next trim each subsystem to remove noisy ends of the spectra in the overlap region. You can get the pixel bounds from ndftrace (here 1024). The example below trims 35 channels from the ‘left’ of rawfile_01 and 24 channels from the ‘right’ of rawfile_02.

  % ndfcopy in=’rawfile_01(35:1024,,)’ out=trimfile_01
  % ndfcopy in=’rawfile_02(0:1000,,)’  out=trimfile_02

At this point the subsystems are aligned pixel for pixel. So a subtraction with sub will generate the differences in the trimmed overlap region analysed by stats from which the median is found, which is 0.04 in the example below.

  % sub trimfile_01 trimfile_02 overlap_12
  % stats overlap_12 order
  % cadd trimfile_02 0.04 offsetfile_02

You can then apply the offset to the second subsystem with cadd. Then the first subsystem can then be merged using wcsmosaic.

  % wcsmosaic in=’"rawfile_01,offsetfile_02"’ out=rawfile12 lbnd=! method=nearest accept

You should check the final spectra to make sure there is no noise bump in the middle where they overlap.

Another technique to locate the noisy edges is apply collapse to all the spectra in an observation using the estimator=sigma option then smooth. The composite formed will show a smooth noise profile, increasing dramatically at its end. The trim limits can be estimated manually, or a linear trend fit with mfittrend will return the bounds of the flat low-noise portion of the spectrum in its aranges parameter.

  % collapse in=rawfile_01 out=profile axis=spec wlim=0.0 estimator=sigma trim
  % block in=profile out=profile_sm box=25
  % mfittrend in=profile_sm out=junk axis=1 method=region auto order=0
  % parget aranges mfittrend

In the example above, parget should return just one pair of pixel co-ordinates, which can then form part of an NDF section to extract the non-noisy part of rawfile_01 using ndfcopy. If there is more than one range you want to use the lowest and highest values.

9.3 Position-velocity diagram

You can create a position-velocity map by collapsing over the skylat or skylon axis of a reduced cube. Use the Kappa command collapse with axis=skylat:

  collapse in=reduced.sdf axis=skylat estimator=sum out=pv

You can view your pv map with Gaia or display.

To obtain a position-velocity map at an arbitrary angle there is a pvslice in the Datacube package. This can be invoked directly, or its recipe extracted from the pvslice script if you know the co-ordinates of the ends of the slice and do not need the graphics.

9.4 Creating channel maps

The Kappa application chanmap creates a two-dimensional channel map from a cube. It collapses along the nominated axis (with many of the same parameters as collapse). The number of channels to create is given by nchan and these are divided evenly between your ranges. The arrangement of the resulting panels is determined by the shape parameter. The collapsed slices are tiled with no margins to form the final image.

  % chanmap in=cube axis=3 low=-10 high=5 nchan=6 shape=3 estimator=mean

This grid of channel maps is filled from left to right, and bottom to top. It can be viewed with Gaia (see Figure 9.1) or with display.

  % display chanmap mode=faint noaxes lut=$STARLINK_DIR/bin/kappa/smooth3_lut

Channel maps can also be created using Gaia. See Section 10.2 for instructions.

Figure 9.1: Displaying a channel map with Gaia.

9.5 Clump finding

The Cupid application findclumps can be used to generate a clump catalogue. It identifies clumps of emission in one-, two- or three-dimensional NDFs. You can select from the clump-finding algorithms FellWalker, GaussClumps, ClumpFind or Reinhold. You must supply a configuration file which lists the options for whichever algorithm you chose (see the cupid manual for a full list). The result is returned as a catalogue in a text file and as a NDF pixel mask showing the clump boundaries (see Figure 9.2). See Appendix E for descriptions of the various algorithms.

  % findclumps in=map out=clumpmap outcat=clumps.FIT logfile=clumps.log \
    config=^myconfig.dat method=clumpfind rms=0.2 shape=polygon

The shape option allows findclumps to create and STC-S description (polygonal or elliptical) for each clump found. These are added as extra columns to the output catalogue.

Polygon Each polygon will have, at most, 15 vertices. For two-dimensional data the polygon is a fit to the clump’s outer boundary (the region containing all good data values). For three-dimensional data the spatial footprint of each clump is determined by rejecting the least significant 10% of spatial pixels, where “significance” is measured by the number of spectral channels that contribute to the spatial pixel. The polygon is then a fit to the outer boundary of the remaining spatial pixels.
Ellipse All data values in the clump are projected onto the spatial plane and “size” of the collapsed clump at four different position angles—all separated by 45—is found. The ellipse that generates the same sizes at the four position angles is then found and used as the clump shape.

You can plot the clump outlines over an image using Gaia. See Section 10.4 for instructions.

Figure 9.2: A velocity slice of the result of three-dimensional clump finding viewed with Gaia. This output NDF contains boundaries of the clumps with the value reflecting the number of pixels contained within that clump.