Load FITS file and list its HDUs.
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
f = fits.open('NGC2916.V500.rscube.fits')
f.info()
print 'array shape:', f[0].data.shape
Compute wavelengths from primary header. We also have to take it to the rest frame.
def get_wave(h):
from astropy import wcs
Nl = h['NAXIS3']
w = wcs.WCS(h).sub([3])
wave = w.wcs_pix2world(np.arange(Nl), 0)
# The value returned is a list of arrays. Since we have only one, we get rid of the list.
return wave[0]
h = f[0].header
l_obs = get_wave(h)
# Redshift
v = h['MED_VEL']
c = 299792.458 # km / s
z = v / c
print 'Redshift: %.3f' % z
l_obs /= 1.0 + z
print 'Wavelengths:', l_obs
Load flux and error spectra, taking the bad pixels into account. Plotting the flux and error spectra for a given pixel.
flux_unit = 1e-16 # erg/s/cm2/AA
badpix_mask = f['badpix'].data > 0
# The 1 + z factor is the k-correction.
flux = np.ma.array(f['primary'].data, mask=badpix_mask) * (flux_unit * (1.0 + z))
error = np.ma.array(f['error'].data, mask=badpix_mask) * (flux_unit * (1.0 + z))
plt.plot(l_obs, flux[:, 32, 32], 'b-', label='flux')
plt.plot(l_obs, 10.0 * error[:, 32, 32], 'r-', label='error x 10')
plt.xlabel(r'$\lambda\ [\mathrm{\AA}]$')
plt.ylabel(r'$F_\lambda\ [\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{\AA}^{-1}]$')
plt.legend()
Taking a slice of the flux.
plt.imshow(flux[:, 31, :].T, extent=[l_obs[0], l_obs[-1], 0, flux.shape[2] - 1], aspect='auto')
plt.xlabel(r'$\lambda\ [\mathrm{\AA}]$')
plt.ylabel(r'dec. $[\mathrm{arcsec}]$')
plt.colorbar()
Taking an image of a spectral band $5635 \pm 45\,\text{\AA}$.
norm_lambda = 5635.0
spec_window = (l_obs > norm_lambda - 45.0) & (l_obs < norm_lambda + 45.0)
image = np.median(flux[spec_window], axis=0)
plt.imshow(image)
Mask that foreground object at 11 o'clock.
flux[:, 42:48, 30:36] = np.ma.masked
image = np.median(flux[spec_window], axis=0)
plt.imshow(image)
Let's begin getting fancy. Calculating the distance of each pixel to the center (brightest pixel).
y0, x0 = np.where(image == image.max())
Ny, Nx = image.shape
y, x = np.ogrid[:Ny, :Nx]
r = np.sqrt((x - x0)**2 + (y - y0)**2)
plt.imshow(r)
Now we can calculate the mean value of flux in bins of r. This is the radial profile of the flux. It may be accomplished in numpy using a little trick with the histogram function.
Usually one creates a histogram by counting the number of elements inside the bins. It is possible also to set weights to the items, so that each does not count as 1 anymore. By choosing the flux as weight, we are in fact computing the sum of flux inside the bins. By dividing the number of items (the regular histogram), we get the mean value of flux inside the bin.
# Bin edges
r_bins = np.arange(30)
# Bin centers
r_coord = np.arange(29) + 0.5
sum_image_r, _ = np.histogram(r, bins=r_bins, weights=image)
N_r, _ = np.histogram(r, bins=r_bins)
image_r = sum_image_r / N_r
plt.plot(r_coord, np.log10(image_r), 'k-')
plt.xlabel(r'$r\, [\mathrm{arcsec}]$')
plt.ylabel(r'$\log\,F_\lambda\ [\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{\AA}^{-1}]$')
A common spatial scale used in galaxies is the Half Light Radius (HLR), that is defined as the radius of a circle (or ellipse) that contains half of the light (flux, luminosity, etc). It may be calculated in a similar way to the radial profile. In fact using the sum of flux in radial bins, we only have to create a cumulative sum and find where it reaches half its maximum value. Repeating the plot with HLR as a scale.
from scipy.interpolate import interp1d
cumsum_image_r = sum_image_r.cumsum()
# We calculate r as a function of the cumsum. We can only do that because the cumsum is monotonic.
r_cs = interp1d(cumsum_image_r, r_coord)
cs_max = cumsum_image_r.max()
HLR = r_cs(cs_max / 2.0)
print 'HLR = %.1f' % HLR
plt.plot(r_coord / HLR, np.log10(image_r), 'k-')
plt.xlabel(r'$r\, [\mathrm{arcsec}]$')
plt.ylabel(r'$\log\,F_\lambda\ [\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{\AA}^{-1}]$')
Now let's do some science. Calculate D(4000) of all pixels. It is defined as the ratio between the flux in the spectral windows of [3850.0, 3950.0] and [4000.0, 4100.0] Angstroms.
# Assume evenly spaced wavelengths.
dl = l_obs[1] - l_obs[0]
blue_window = (l_obs < 3950.0) & (l_obs > 3850)
red_window = (l_obs < 4100.0) & (l_obs > 4000)
flux_b = np.trapz(flux[blue_window], l_obs[blue_window], axis=0)
flux_r = np.trapz(flux[red_window], l_obs[red_window], axis=0)
plt.imshow(flux_r / flux_b, vmin=1, vmax=2)
plt.colorbar()
The white spots in the image are pixels where the data is all masked. This may be mitigated using COMBO cubes.
Exercise: Calculate the equivalent width of $\mathrm{H}\delta$.
$$ \mathrm{EW} = \int \frac{C_\lambda - F_\lambda} {C_\lambda} \mathrm{d}\lambda $$where $C_\lambda$ is the estimated continuum inside the line. Use these limits for the measurements:
- Absorption line: 4083.500, 4122.250
- Blue continuum: 4041.600, 4079.750
- Red continuum: 4128.500, 4161.000