Exercise 13.3

Exercise 13.3#

Your task is to use Astropy to read in a FITS image of a galaxy emission, isolate the region of interest, remove the background noise, calibrate it and save this processed image into a new FITS file. Download the FITS image ex11.2.fits from the resource files on Vula, and then create a Python script or notebook to do the following, filling in the missing code where necessary:

Read In the FITS Image#

Read in the FITS image using the astropy.io package. What are the pixel dimensions of the image?

import numpy as np
import matplotlib.pyplot as plt

from astropy.wcs import WCS
from astropy.io import fits

##Read in FITS file
hdu = fits.open('ex11.2.fits') #header data units
hdr = hdu[0].header #Header
img = hdu[0].data #Image

print('image pixel dimensions: ', np.shape(img))

Plot the Image#

Now would be a good time to plot the image you’ve read in. As this image is two dimensional and contains light intensity as the data for each pixel, a useful matplotlib.pyplot function to use is imshow(). This produces a heat map from 2D arrays.

Firstly, we should determine the interval in which our minimum and maximum colour values should be. What are the min/max pixel values (in raw image units)?

#the min/max pixel values in img will guide our choice for the colour map max and min values
print('min/max pixel values [dn]: ', img.min(), img.max())

Now for the image itself:

#Plotting the image as a colour map
plt.imshow(img, vmin = 1, vmax=50, cmap='gist_heat')
plt.show()

Note that the values of vmin and vmax were chosen to make the plot look presentable.

Extracting the Region of Interest#

The dark regions of the image are of no interest to us. Extract a 601 by 601 pixel sub-image from the central region of img by slicing it and plot this. Follow the following steps to achieve this:

  1. Find central the pixel’s coordinates: midx and midy

  2. Extract a sub-image centred on these pixels (call it img_new). Take the region starting 300 pixels before (midx, midy) and ending 300 pixels after (midx, midy). In other words:

    • The \(x\)-range of the image should be from midx - 300 up to and including midx + 300

    • And the \(y\)-range from midy - 300 up to and including midy + 300

  3. Plot this extracted sub-image using imshow to verify that the primary region of interest is contained within the new image.

#Extracting the region of interest

# Your code here

Create a New FITS Image#

Now, to create the new FITS image from this extracted image. The purpose of the code that follows is to:

  1. Manually create the header by using astropy to create a very basic header.

  2. Get RA, Dec of midx, midy in original image, and assign those coordinates to centre pixels of new sub-image (WCS will be needed for this).

  3. Add the astronometry keywords from the ex11.2.fits header.

  4. Update BUNIT as you will be calibrating the image later.

hdu_new = fits.PrimaryHDU(img_new) 
hdulist_new = fits.HDUList([hdu_new])
hdr_new = hdulist_new[0].header 

#initialise WCSAxes:
wcs = WCS(hdr)

#get RA, Dec of midx, midy in original image, and assign those 
#coordinates to centre pixels of new sub-image:
ref_coords = wcs.wcs_pix2world(np.array([[midx,midy]]),1)

#now add the astronometry keywords from the ex11.2.fits header:
hdr_new.update(CRVAL1=ref_coords[0][0], CRVAL2=ref_coords[0][1]) 
hdr_new.update(CRPIX1=300, CRPIX2=300)
hdr_new.update(CDELT1=hdr['CDELT1']) 
hdr_new.update(CDELT2=hdr['CDELT2']) 
hdr_new.update(CTYPE1=hdr['CTYPE1']) 
hdr_new.update(CTYPE2=hdr['CTYPE2']) 

#we're going to calibrate the image, so we may as well update BUNIT:
hdr_new.update(BUNIT='app mag')
wcs_new = WCS(hdr_new)

Investigating the Background Noise#

The measurement contains background noise which should be subtracted in order to clarify the signal data. In this section you will investigate the noise distribution.

From the original full-size image, extract the 400 x 400 pixels making up the top left corner (this starts at the index (0,0) up to 400 on each axis). These pixels are mostly noise pixels. Create a histogram of the values of these pixels, and calculate the mean and standard deviation of this distribution.

#Investigating a sample of the background noise

#Your code here

The basic noise statistics are given in the FITS header:

  • The mean of the noise is given by hdr['SKY']

  • The standard deviation of the noise is given by hdr['SKYSIG']

Compare your noise sample to this header information by overlaying a normal distribution with these moments over your histogram.

#Comparing the background noise sample to the header data

#Your code here

Removing the Background Noise#

Subtract the sky background (as given in the header) from your sub-image.

sky = hdr['SKY']
img_new = img_new - sky

Once you have your sky-subtracted sub-image, apply a 2.5-sigma flux cut by setting any pixels \(< 2.5~\sigma~\) to non-values. This will get rid of most of the noise.

skysig = hdr['SKYSIG']
img_new[img_new < 2.5*skysig] = np.nan

Calibrating the New Image#

Photometrically calibrate your sky-subtracted, flux-cut sub-image using: apparent mag = MAGZP \(-2.5 \log_{10}\)(img). MAGZP is the magnitude zero point, given in the FITS header.

magzp = hdr['MAGZP'] 
img_new = magzp - 2.5*np.log10(img_new)

Saving the New FITS File#

Save the new image data to a new FITS file using fits.writeto():

fits.writeto('ex11.2_sub.fits', img_new, hdr_new, overwrite=True)

Plotting the Final Image#

Finally, using a linear colour stretch (having an image in units of mag/arcsec\(^2\) automatically provides us with a logarithmic colour stretch):

ax = plt.subplot(projection=wcs_new)
plt.imshow(img_new, cmap='gist_heat', origin='lower')
ra = ax.coords['ra'] 
dec = ax.coords['dec']
ra.set_axislabel('Right Ascension', minpad=1)
dec.set_axislabel('Declination', minpad=1)
ra.grid(color='black', alpha=0.5, linestyle='solid')
dec.grid(color='black', alpha=0.5, linestyle='solid')
plt.colorbar()

plt.show()

Make sure to include this image alongside your code in you tutorial submission.