Observational Astronomy Primer Exercises: 3.) Reducing SPIRE Data, Aperture Photometry

Tasks/exercises for SPIRE data in HIPE

As an example, work through one supernova remnant from the Hi-Gal survey, G332.4+0.1.  You can find more information about this object at, e.g., SIMBAD.

  • Try to download the Hi-Gal data that has already been calibrated and mapped.  Check the data quality to determine if you need to re-reduce or make the map again (for example, if there is some residual striping from the cross-scans).
  • If you decide you need to re-reduce it, you can follow the links I’ve listed below.
  • You will need the OBSID – you can look for this with the HSA Tool at  http://www.cosmos.esa.int/web/herschel/science-archive
    –> Input the object name or coordinates and it will return a list of observations and IDs
    There will be many observations with PACS/SPIRE/HIFI. Choose the ones you want to see.  The SPIRE large map obs for G332.4+0.1 were OBSID 1342192055
  • MAKE SURE you are using a recent version of HIPE if you plan to reprocess the data.  You will need LOTS of RAM to reduce Herschel data – aim for a machine with at least 32GB.



Data reduction pipeline in HIPE


==> See tutorials on NASA/IPAC’s site:
And the DRG for reprocessing http://herschel.esac.esa.int/hcss-doc-15.0/load/spire_drg/html/photometer_launchpad.html#d0e571
The tutorial for the SPIRE large map pipeline reprocessing:
That one is quite good – much better than was available when I started learning how to reduce SPIRE data, with many picture examples of things to look out for in your final images
They have also made available a large number of web video tutorials over the past few years:



SPIRE photometry ‘recipe’

http://herschel.esac.esa.int/hcss-doc-15.0/load/spire_drg/html/ch06s09.html  –> section
Also see the HIPE Data Analyis Guide for a description of basically any analysis task you would care to do in HIPE:

Overall summary of commands from the SPIRE ‘recipe’, downloading level-2 data from the HSA:

obsid = 1342192055  # Specify the observation ID from the HSA  (Kes 32)
alpha = 2   #For a source with spectrum S(ν) proportional to Sα (SPIRE default pipeline assumes α = -1)
array = "PSW"  # Run for an individual SPIRE band: "PSW", "PMW", "PLW"
obs = getObservation(obsid, useHsa=True, instrument='SPIRE') # Loading an observation of Gamma Dra from the HSA
# obs = getObservation(obsid, poolName='mypool', instrument='SPIRE') # Alternative, for observation from your own local pool
mapExtd = obs.level2.refs["extd"+array].product  #Extract  Extended (MJy/sr) calibrated maps from the Observation Context

cal = spireCal() #Load the calibration tree
#  --> if that doesn't work: cal = spireCal(calTree="spire_cal", saveTree=1)
beamCorrTable  = cal.phot.refs["ColorCorrBeam"].product
kCorrExtdTable = cal.phot.colorCorrKList.refs[0].product
beamArea  = beamCorrTable.meta["beamPipeline"+array.title()+"Arc"].double
kCorrExtd = kCorrExtdTable.getAlphaCorrection(alpha, array)

mapExtended = convertImageUnit(image=mapExtd, newUnit='Jy/pixel') # convert maps from MJy/sr to Jy/pix
mapExtendedCorrected = imageMultiply(image1=mapExtended, scalar=kCorrExtd) #Color Correction

ra  = '244.252833'   # Target RA and string values:  16h 17m 00.68s
dec = '-50.799300'   # Target Dec and string values: -50d 47' 57.48''
photrad = 200.0      #photometry source radius, in arcsec  (I'm just making these up - you will have to choose an appropriate region)
phot_ann_in = 300.0  #photometry annulus inner radius, in arcsec
phot_ann_out = 300.0 #photometry annulus outer radius, in arcsec
# Carry out circular/annulus Aperture Photometry
annularPSW = annularSkyAperturePhotometry(image= mapExtendedCorrected, \
    centerRA=ra, centerDec=dec, fractional=1, \
    radiusArcsec=photrad, innerArcsec=phot_ann_in, outerArcsec=phot_ann_out)
# annularSkyAperturePhotometry() #Alternatively, define the photometry regions manually

flux = annularPSW.getTargetTotal()  # Final target brightness in Jy/pixel
print 'PSW flux = %5.3f Jy'%(flux)


Observational Astronomy Primer Exercises: 2.) Python plotting / fitting / etc.

You may find some of these tasks useful as we move forward doing photometry and running HIPE.  This is just some useful material that you can work on over the next few weeks or months.

You can write & run python scripts in a few different ways, and the best way for you is really just up to your preference.  I like to write my code in gedit (Pluma in Ubuntu Mate), and interactively test/run in ipython in the terminal.  But there are also several IDEs such as Spyder, Eclipse, Geany, PyCharm (similar to those used by Matlab or IDL).

There are many great tutorials on the web for learning python.  I started to make a few very simple tutorials elsewhere on this site ( http://research.endlessfernweh.com/tutorials/ ) but there are a variety of better tutorials out there:
Practical Python For Astronomers  https://python4astronomers.github.io/  is a great comprehensive one tailored to astronomy, but I find that it’s a little bit hard to jump right in
Here is a page with powerpoint-style tutorials on many things in python (again, targeted at astronomers): https://astro.uni-bonn.de/~rschaaf/Python2008/
Astrobetter is a blog that often has very useful posts about coding and other important things in astronomy.  Here is one post on python: http://www.astrobetter.com/blog/2014/07/07/astropy-tutorials-learn-how-to-do-common-astro-tasks-with-astropy-and-python/

For plotting (see below as well), I recommend using the pywcsgrid2 module, but you can also use APLpy and kapteyn to make great  images.  Here are a couple pages for using pywcsgrid2:
Kapteyn has many other useful tools for reprojection, image smoothing, etc.  https://www.astro.rug.nl/software/kapteyn/

Here is a list of useful goals & concepts to work on.

  • importing modules, syntax
  • Create some functions (try some useful ones for astronomy, such as a Gaussian, blackbody, Modified BB)
  • Load tables from text files, or DS9 regs.
  • Save data to text files.
  • Fit function parameters – see my web tutorial for a good overview
  • Plot single function and 2D (heatmap, maybe try a 2D gaussian?) – see my web tutorial
    • include labels – title, axis labels (modify size, color, etc.)
    • include legend (1D), colorbar (2D)
  • Plot single images – maybe start with Kepler’s SNR because we already downloaded those fits files
    • Can use APLpy, pywcsgrid2, kapteyn… whatever package you like
    • Include labels
    • include colorbar & units
    • play with colormaps, intensity limits, scaling
    • maybe include size scale bar, compass, beam size…
    • draw a region (e.g., ellipse)
  • Plot two images side by side (say, Kepler infrared and radio; follow the pywcsgrid2 tutorial)
    • match coordinates
    • make contours


Observational Astronomy Primer Exercises: 1.) Viewing FITS Images

To help you become familiar with astronomical data files – .FITS format – I would like to give you a small task: please experiment with some supernova images in DS9.   Here is what I’m thinking for some small goals:

  • Download data for one astronomical target from two different telescopes and display them together, side-by-side.
  • Stretch the image scaling to appropriate levels and experiment with color maps.
  • Match the coordinates.
  • Make labels.
  • Create regions.
  • Output/print an image (png, jpg, etc.) of your results & save ‘backup’ file.


I discuss photometry in another post & how to do simple flux extraction with DS9 and HIPE later.

Try Kepler’s SNR (G4.5+6.8) in IR  http://chandra.harvard.edu/photo/openFITS/multiwavelength_data.html   and try to find one more image from another band that is not listed on that page.  I will discuss some good online databases tomorrow, but one good option would be a radio image from the VLA – you can download VLA data from the NRAO image archive retrieval tool  https://archive.nrao.edu/archive/archiveimage.html




Notes on image headers

A good .fits file will include necessary WCS (pointing) and brightness units.  And beam size information if needed.


The main header ‘card’ to look for is ‘BUNIT’, which should be something like Jy/beam, MJy/sr, W/m2/pixel, etc.  Sometimes images you download (such as HST, GALEX…) will still be in units of counts/sec or other system that will need to be converted to physical units before you can integrate fluxes.  If you need to perform these conversions, try looking at the instrument webpages or papers in the literature to get the conversion factors, otherwise you can try to contact the PI who owns the data.

To convert the beam area – that is, to calculate the beam size as pixels/beam or sr/beam etc – our job is easy if the PSF can be approximated by a 2D Gaussian.  Properly-formed radio, sub-mm and IR headers should have the beam parameters listed as BMAJ, BMIN, and BPA – the beam major axis FWHM, minor axis FWHM, and position angle, respectively.  All 3 are normally in units of degrees.  The PA in the FITS standard is defined as counter-clockwise from ‘image north’ or the upper y-axis, as opposed to the mathematical standard right x-axis assumed by python and other programming languages.

WCS (World Coordinate System)

The information about which pixel corresponds to what RA/DEC on the sky is given in the header as well.  Here is the NASA dictionary of standard header cards:  http://heasarc.gsfc.nasa.gov/docs/fcg/standard_dict.html
Basically, all the header cards beginning with C (CTYPE, CDELT, …) have to do with the the axis and pointing information.  Here is a brief description of some:
CTYPE gives the projection type – SIN, TAN, …
CRPIX gives the reference pixel number (starting from 1)
CRVAL gives the pointing coordinate of the reference pixel
CUNIT gives the units of the reference coordinate values (default is degrees)
CDELT gives the width of each pixel in sky coordinates
CROTA gives the image rotation angle, if applicable.

The number after the main keyword denotes which image axis – 1 is for the x-axis (RA), 2 is for the y-axis (DEC), 3 is for the z-axis (velocity/frequency/wavelength) in a 3D data cube, 4 is for the Stokes axis.  So, CRPIX1 is the reference pixel for the x-axis, CDELT2 is the pixel width in the y-axis, CUNIT3=’Hz’ means that the spectral axis channels are stepped in frequency, etc.

Instead of CRDELT/CROTA etc. you may find your image using an alternate grid reference system based on the cards CD1_1, CD2_1, CD1_2 and CD2_2.  These four parameters define a matrix grid of the pixel scale and rotation.