Calculating fluxes from SCUBA-2 maps

One common task JCMT scientists carry out is estimating the dust mass or column density from the emission in a SCUBA-2 map. To do this, it is necessary to calculate the 850um dust flux of an object/field.

Caveats

In general there are numerous scientific subtleties to which type of photometry is most appropriate for a given situation. Point sources are probably the most straight forward, followed by cases where you can integrate over an entire large area. The common (for galactic scientists) case of having extended structure throughout a large map, and hoping to break it down into separate structures lends itself to a clumpfinding approach, a filament finding approach or attempting to convert your map into hierarchical structures.

This guide will just quickly demonstrate some of the simpler cases and how to use Starlink software to calculate these. It is not a guide to which approach is appropriate for which scientific investigations, and it does not address the underlying uncertainties in these approaches.

First of all, we have to set up Starlink and KAPPA as normal. We will also create a working directory and change into it, and ensure that Starlink won't try and prompt us for any values.

(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
kappa
mkdir S2Fluxes
cd S2Fluxes/

     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.

   

Calibration units

SCUBA-2 maps are normally calibrated by our pipeline into units of either mJy/beam or mJy/arcsec$^2$. If you need more detail on the calibration process, please see Dempsey et al 2013, or our webpages at https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/calibration.

Point sources

If you want to find the flux of a point source, you would normally calibrate your map in units of Jy/beam. Each point source corresponds to one beam, and therefore the peak value in the source will tell you the total flux density of the source.

For this example, we are going to use a calibrator observation of the bright point source CRL618, in the file Data/jcmts20191022_00059_850_reduced001_nit_000.sdf. You can download this file from CADC at this link here, or you can go to the JCMT Archive at CADC and download a reduced calibrator observation of your choice. If you download a file from the archive, you'll need to convert it from FITS format to NDF (see fits2ndf in the CONVERT package), and please ensure it is a point source.

The first step is to see what units your map is currently in. We can use the KAPPA command ndftrace to do this. (Here we have also used the shell command head to only show the first 10 lines).

In [2]:
ndftrace ../Data/jcmts20191022_00059_850_reduced001_nit_000.sdf|head -n 10
   NDF structure
/export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/../Data/jcmts20191022_00
059_850_reduced001_nit_000:
      Title:  CRL618
      Label:  Flux Density
      Units:  mJy/beam

   Shape:
      No. of dimensions:  3

As you can see, this is already in mJy/beam, so we don't need to change the calibration. If we did, we would normally use the PICARD recipe CALIBRATE_SCUBA2_DATA with the recipe parameters FCF_CALTYPE=BEAM,USEFCF=True.

The nature of JCMT scan patterns mean that we normally have areas of higher noise around the map. There is also a PICARD recipe that can crop that out, CROP_SCUBA2_IMAGES. Here we will call it asking for a circle of radius 250". (Please note this requires all the input observations to a map to have been towards the same source, or it will be hard to predict where the mask will be centered.)

In [3]:
picard -log sf -nodisplay CROP_SCUBA2_IMAGES --recpars='CROP_METHOD=CIRCLE,MAP_RADIUS=250' \
    ../Data/jcmts20191022_00059_850_reduced001_nit_000.sdf
Picard Says: No display will be used
Use of uninitialized value in subroutine entry at /star/bin/oracdr/src/lib/perl5/ORAC/Basic.pm line 155.
Picard Says: Pre-starting mandatory monoliths...Done
Checking for next data file: /export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/../Data/jcmts20191022_00059_850_reduced001_nit_000.sdf
Storing: jcmts20191022_00059_850_reduced001_nit_000
Use of uninitialized value in subroutine entry at /star/bin/oracdr/src/lib/perl5/ORAC/Group.pm line 419, <DATA> line 1676.
Picard Says: Creating temporary bad observation rules file
A new group 20191022#0 has been created
Overriding PICARD instrument class to PICARD_SCUBA2_850
Sorting Groups
REDUCING: jcmts20191022_00059_850_reduced001_nit_000
Using recipe CROP_SCUBA2_IMAGES specified on command-line
Processing data for CRL618

Calling _CROP_SCUBA2_IMAGE_: trim image to specified map size
Trimming image to specified map size
Output image will be a circle of radius 250 arcsec
Masking the weights and exposure time images...
Removing temporary files...
Checking jcmts20191022_00059_850_reduced001_nit_000_tmpmask...	Removing
Checking jcmts20191022_00059_850_reduced001_nit_000_crop...	Keeping extension
Recipe took 2.657 seconds to evaluate and execute.

Picard processing complete
Processed one recipe which completed successfully
Exiting...

Picard Says: Goodbye

The output file is jcmts20191022_00059_850_reduced001_nit_000_crop. We can now take a quick look at our map to make sure it looks sensible.

In [4]:
lutwarm dev=PNG
palentry 0 white dev=PNG
palentry 1 black dev=PNG

display in=jcmts20191022_00059_850_reduced001_nit_000_crop dev=cropped.png/PNG mode=faint style=!
Data will be scaled from -176.595628266717 to 1292.29450845203.

Cropped image

Point source fluxes

You can get a quick and crude estimate of the flux of the point source by just "reading off" the peak flux in the source, using e.g. KAPPA's stats. The brightest pixel should roughly give you the flux of the source in mJy, assuming there is not a high background level.

In [5]:
stats jcmts20191022_00059_850_reduced001_nit_000_crop
   Pixel statistics for the NDF structure
/export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/jcmts20191022_00059_850_
reduced001_nit_000_crop

      Title                     : CRL618
      NDF array analysed        : DATA

         Pixel sum              : 1377520.68292101
         Pixel mean             : 7.01563882312718
         Standard deviation     : 183.610793724692
         Skewness               : 9.51239441182121
         Kurtosis               : 20.5026101428456
         Minimum pixel value    : -866.451929473621
            At pixel            : (-11, -213, 1)
            Co-ordinate         : (4:42:54.7, 36:03:22, 352.6734)
         Maximum pixel value    : 5556.1683556297
            At pixel            : (1, -2, 1)
            Co-ordinate         : (4:42:53.7, 36:06:53, 352.6734)
         Total number of pixels : 250000
         Number of pixels used  : 196350 (78.5%)
         No. of pixels excluded : 53650 (21.5%)

The peak pixel has a value of 5556.168 mJy/beam, corresponding to a source flux of 5.6 Jy. This is quite a bit higher than our website value of 4.89 ± 0.24 Jy/beam for the peak brightness (see https://www.eaobservatory.org/jcmt/instrumentation/continuum/scuba-2/calibration/calibrators/, but within the correct range.

Gaussian fit

A better estimate would be to fit a 2-D Gaussian to the source and find the peak. We can do this with the KAPPA command beamfit. We take the initial coordinates to be the position of the max pixel, as found from stats above.

In [6]:
beamfit ndf=jcmts20191022_00059_850_reduced001_nit_000_crop mode=interface pos='"4:42:53.7,36:06:53"' beams=1
    Co-ordinates    : RA(hh:mm:ss.s) Dec(ddd:mm:ss)

    Centre          : (4:42:53.7,36:06:53) +/- (0.021,0.017) arcsec
    Offset          : 0.119 +/- 0.024 arcseconds from sky reference position
    FWHM (major)    : 13.922 +/- 0.042 arcseconds
    FWHM (minor)    : 13.019 +/- 0.040 arcseconds
    Orientation     : 41.29 +/- 0.311E-01 degrees (measured from North through East)
    Amplitude       : 5544.8441824542 +/- 16.4217721767325
    Background      : 1.21604628436329 +/- 0.266401816206429
    Shape exponent  : 2.000

    RMS fit error   : 117.796970900288


From the amplitude value, this gives a flux of 5.54+/-0.02 Jy for the source.

(You can see from the output that this has a background of 1.2mJy/beam, or or 0.001 Jy, so any residual background should not be a significant effect on this value).

Integrated fluxes.

If you want to calculate the flux integrated over resolved sources you have to consider what you mean by a source. One approach for a compact object would be to use annular photometry, including a correction for the background. This is what we do to calculate our arcsecond FCF (see Dempsey et al 2013). Another approach would be to calculate the flux within an entire map or object. Here we will demonstrate calculating the integrated flux in a 2-D Clump, such as that found when running the CUPID fellwalker on a map.

Here we will use a map of G34.3, an occassional SCUBA-2 calibrator source. The file is in Data/jcmts20191018_00009_850_reduced001_nit_000.sdf

Calibration

First we need to ensure our map is in Jy/arcsec$^{2}$.

In [7]:
ndftrace ../Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
   NDF structure
/export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/../Data/jcmts20191018_00
009_850_reduced001_nit_000:
      Title:  G34.3
      Label:  Flux Density
      Units:  mJy/beam

   Shape:
      No. of dimensions:  3
      Dimension size(s):  233 x 233 x 1
      Pixel bounds     :  -121:111, -101:131, 1:1
      Total pixels     :  54289

   Data Component:
      Type        :  _DOUBLE
      Storage form:  SIMPLE
      Bad pixels may be present

   Variance Component:
      Type        :  _DOUBLE
      Storage form:  SIMPLE
      Bad pixels may be present

   Quality Component:
      Type        :  _UBYTE
      Storage form:  SIMPLE
      Bad-bits mask:  0 (binary 00000000)

   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-SPECTRUM
        First pixel centre  : 18:53:51.1, 1:08:10, 352.7287

           Axis 1:
              Label              : Right ascension
              Units              : hh:mm:ss.s
              Nominal Pixel scale: 3.99999 arc-sec

           Axis 2:
              Label              : Declination
              Units              : ddd:mm:ss
              Nominal Pixel scale: 3.99999 arc-sec

           Axis 3:
              Label              : Frequency
              Units              : GHz
              Nominal Pixel scale: 31.65514 GHz


   Extensions:
      FITS             <_CHAR*80>
      PROVENANCE       <PROVENANCE>
      SMURF            <SMURF_EXT>


   History Component:
      Created    :  2019 Oct 19 01:41:15
      No. records:  17
      Last update:  2019 Oct 19 01:41:36 (FITSMOD         (KAPPA 2.6))
      Update mode:  NORMAL

You can see that it is in mJy/beam, from the 'Unit' component shown by ndftrace. We will use the PICARD recipe CALIBRATE_SCUBA2_DATA to switch the units to mJy/arced$^2$. (This recipe is very simple: it checks what the current calibration is (if any), undoes it if required by dividing by the FCF, then multiplies by the correct standard FCF factor for the requested units.)

In [8]:
picard -log sf -nodisplay CALIBRATE_SCUBA2_DATA --recpars='FCF_CALTYPE=ARCSEC,USEFCF=True' \
../Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
Picard Says: No display will be used
Use of uninitialized value in subroutine entry at /star/bin/oracdr/src/lib/perl5/ORAC/Basic.pm line 155.
Picard Says: Pre-starting mandatory monoliths...Done
Checking for next data file: /export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/../Data/jcmts20191018_00009_850_reduced001_nit_000.sdf
Storing: jcmts20191018_00009_850_reduced001_nit_000
Use of uninitialized value in subroutine entry at /star/bin/oracdr/src/lib/perl5/ORAC/Group.pm line 419, <DATA> line 1676.
A new group 20191018#0 has been created
Overriding PICARD instrument class to PICARD_SCUBA2_850
Sorting Groups
REDUCING: jcmts20191018_00009_850_reduced001_nit_000
Using recipe CALIBRATE_SCUBA2_DATA specified on command-line
Processing data for G34.3

Calling _UNCALIBRATE_SCUBA2_DATA_: undo calibration of given data
Undoing calibration of 537000 mJy/beam/pW - output units are pW

Calling _CALIBRATE_SCUBA2_DATA_: calibrate data using standard, given or derived FCF
Picard Says: Calibrating data in mJy/arcsec**2
Multiplying jcmts20191018_00009_850_reduced001_nit_000_uncal by 2340 mJy/arcsec**2/pW
Removing temporary files...
Checking jcmts20191018_00009_850_reduced001_nit_000_uncal...	Removing
Checking jcmts20191018_00009_850_reduced001_nit_000_uncal_cal...	Keeping extension
Recipe took 1.090 seconds to evaluate and execute.

Picard processing complete
Processed one recipe which completed successfully
Exiting...

Picard Says: Goodbye

The output file name is jcmts20191022_00059_850_reduced001_nit_000_uncal_cal.sdf: this shows you that PICARD first removed the existing calibration _uncal and then applied the new calibration _cal.

We'll now run a very quick clumpfinding routine using the Fellwalker algorithm in CUPID. First we have to initialise the CUPID package:

In [9]:
cupid

     CUPID commands are now available -- (Version 2.5)

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

   

Next we will run the findclumps command with the fellwalker algorithm, but leaving most variables to their default value.

In [10]:
findclumps in=jcmts20191018_00009_850_reduced001_nit_000_uncal_cal method=FellWalker out=clumps outcat=! accept

FellWalker:
46 clumps rejected because they contain too few pixels.
4 further clumps rejected because they are smaller than the spatial beam
width.
4 further clumps rejected because they include too many bad pixels.
8 usable clumps found.


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


You can see from the output that it found 8 clumps. We can take a quick look at the data values in the first clump by examining the CUPID.CLUMPS extension within the NDF like so:

In [11]:
display in=clumps.sdf.MORE.CUPID.CLUMPS\(1\) dev=clumps.png/PNG mode=faint style=!
stats clumps.sdf.MORE.CUPID.CLUMPS\(1\)
Data will be scaled from -11.2030809988406 to 220.241529596863.

   Pixel statistics for the NDF structure
/export/data/sgraves/Tutorials/AnalysisHowTos/S2Fluxes/clumps.MORE.CUPID.CLUMPS
(1).MODEL

      Title                     : <undefined>
      NDF array analysed        : DATA

         Pixel sum              : 8828.29267215991
         Pixel mean             : 17.7274953256223
         Standard deviation     : 28.9015145519101
         Skewness               : 3.85402941028909
         Kurtosis               : 16.9699034251728
         Minimum pixel value    : 1.73469502409571
            At pixel            : (-14, 4)
            Co-ordinate         : (-14.5, 3.5)
         Maximum pixel value    : 205.394380084719
            At pixel            : (0, 0)
            Co-ordinate         : (-0.5, -0.5)
         Total number of pixels : 864
         Number of pixels used  : 498 (57.6%)
         No. of pixels excluded : 366 (42.4%)

Image of the first clump

As you can see from when we ran ndftrace above, we have 4" pixels. Therefore we need to integrate our flux over the whole region, to get a value in Jy instead of in Jy/arcsec$^{2}$. The flux in our clump is therefore the total sum of the emission, which is 8828 mJy/arcsec$^{2}$from stats, multipled by 16 square arcseconds = 141 Jy. This assumes that the background level is 0.