import numpy as np
# Set up matplotlib and use a nicer set of plot parameters
%config InlineBackend.rc = {}
import matplotlib
matplotlib.rc_file("../../templates/matplotlibrc")
import matplotlib.pyplot as plt
%matplotlib inline
The following line is needed to download the example FITS files used here.
from astropy.utils.data import download_file
Viewing and manipulating FITS images¶
FITS files are the standard format for all astronomical data.
See details on astropy's FITS module here: https://astropy.readthedocs.org/en/latest/io/fits/index.html
from astropy.io import fits
# import pyfits
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )
Opening FITS files and loading the image data¶
I will open the FITS file and find out what it contains.
hdu_list = fits.open(image_file)
hdu_list.info()
Generally the image information is located in the PRIMARY block. The blocks are numbered and can be accessed by indexing hdu_list.
image_data = hdu_list[0].data
# image_data = hdu_list['PRIMARY'].data
You data is now stored as a 2-D numpy array. Want to know the dimensions of the image? Just look at the shape of the array.
print(type(image_data))
print(image_data.shape)
The Header is also an important component. The info above tells us it has 161 cards:
header = hdu_list['PRIMARY'].header
header
SHORTCUT¶
If you don't need to examine the FITS header, you can call fits.getdata to bypass the previous steps.
image_data = fits.getdata(image_file)
print(type(image_data))
print(image_data.shape)
Similarly, if you only want the header, you can use fits.getheader:
header = fits.getheader(image_file)
Viewing the image data and getting basic statistics¶
plt.imshow(image_data, cmap='gray')
plt.colorbar()
# To see more color maps
# http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
Let's get some basic statistics about our image. If your image contained NaNs, you'd use np.nanmin etc below
print('Min:', np.min(image_data))
print('Max:', np.max(image_data))
print('Mean:', np.mean(image_data))
print('Stdev:', np.std(image_data))
Plotting a histogram¶
To make a histogram with matplotlib.pyplot.hist(), I need to cast the data from a 2-D to array to something one dimensional.
In this case, I am using the iterable python object img_data.flat.
print(type(image_data.flat))
NBINS = 1000
histogram = plt.hist(image_data.flat, NBINS)
Displaying the image with a logarithmic scale¶
Want to use a logarithmic color scale? To do so we need to load the LogNorm object from matplotlib.
from matplotlib.colors import LogNorm
plt.imshow(image_data, cmap='gray', norm=LogNorm())
# I chose the tick marks based on the histogram above
cbar = plt.colorbar(ticks=[5.e3,1.e4,2.e4])
cbar.ax.set_yticklabels(['5,000','10,000','20,000'])
Basic image math: image stacking¶
You can perform math with the image data like any other numpy array. In this particular example, I will stack several images of M13 taken with a ~10'' telescope.
I open a series of FITS files and store the data in a list, which I've named image_concat.
url_template = 'http://data.astropy.org/tutorials/FITS-images/M13_blue_000{0}.fits'
image_list = [download_file(url_template.format(n), cache=True)
for n in ('1','2','3','4','5')]
# The long way
image_concat = []
for image in image_list:
image_concat.append(fits.getdata(image))
# The short way
#image_concat = [ fits.getdata(image) for image in IMAGE_LIST ]
Now I'll stack the images by averaging my concatenated list.
# The long way - slow, but clear
#final_image = np.zeros(shape=image_concat[0].shape)
#
#for image in image_concat:
# final_image += image
# final_image /= len(image_concat)
# The short way: numpy will turn your list of images into a 'cube'
# with axis=0 being the shortest (stacked) axis
final_image = np.mean(image_concat, axis=0)
I'm going to show the image, but I want to decide on the best stretch. To do so I'll plot a histogram of the data.
image_hist = plt.hist(final_image.flat, 1000)
I'll use the keywords vmin and vmax to set limits on the color scaling for imshow.
plt.imshow(final_image, cmap='gray', vmin=400, vmax=600)
plt.colorbar()
Writing image data to a FITS file¶
This is easy to do with the writeto() method.
You will receive an error if the file you are trying to write already exists. That's why I've set clobber=True.
Most tools in the astropy and python ecosystem use the keyword overwrite instead of clobber; clobber is inherited from IRAF.
outfile = 'stacked_M13_blue.fits'
hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, clobber=True)
Using the world coordinate system¶
"World coordinates" refer to RA/Dec or Galactic latitude/longitude.
The full documentation is available at: http://astropy.readthedocs.org/en/latest/wcs/
from astropy import wcs
Step 1. Create a WCS object
mywcs = wcs.WCS(header)
The most important tools related to a WCS are the conversion from pixels -> wcs and wcs -> pixels, accomplished with WCS.wcs_pix2world and WCS.wcs_world2pix.
The syntax is cumbersome for converting a single pixel:
xpix, ypix = [50,100]
# The '0' indicates that the pixels are indexed starting at 0, the python/c convention
# You can also use '1', which indicates the FITS/IRAF/FORTRAN convention
(ra,dec), = mywcs.wcs_pix2world([[xpix,ypix]], 0)
print(ra,dec)
But convenient for converting an array of pixels:
pix_coords = [(20,30), (50,10), (500,200)]
world_coords = mywcs.wcs_pix2world(pix_coords, 0)
print(world_coords)
The coordinates can be readily converted to astropy coordinates:
from astropy import coordinates
from astropy import units as u
my_coords = coordinates.SkyCoord(world_coords*u.deg, frame='fk5')
print(my_coords)
print(my_coords.to_string(style='hmsdms'))
You can also use WCS to get pixel coordinates for overplotting sources.
horsehead = coordinates.SkyCoord.from_name('Horsehead Nebula')
print(horsehead)
(hxpix,hypix), = mywcs.wcs_world2pix([(horsehead.ra.deg,
horsehead.dec.deg)],0)
print(hxpix,hypix)
plt.imshow(image_data, cmap='gray', norm=LogNorm())
# I chose the tick marks based on the histogram above
cbar = plt.colorbar(ticks=[5.e3,1.e4,2.e4])
cbar.ax.set_yticklabels(['5,000','10,000','20,000'])
plt.plot(hxpix, hypix, 's')
You can also create an HDU 'from scratch', but with a WCS header:
myhdu = fits.PrimaryHDU(data=image_data, header=mywcs.to_header())
Plotting the world coordinate system¶
Plotting using WCS is a more complicated task and should be left to specialized tools:
- http://wcsaxes.rtfd.org/
- http://aplpy.github.io/
import aplpy
F = aplpy.FITSFigure(hdu_list[0], figure=plt.figure(1))
F.show_grayscale()
# Can use world coordinates directly
F.show_markers(horsehead.ra.deg, horsehead.dec.deg)
from wcsaxes import WCSAxes
fig = plt.figure()
ax = WCSAxes(fig, [0.1, 0.1, 0.8, 0.8], wcs=mywcs)
fig.add_axes(ax) # note that the axes have to be added to the figure
ax.imshow(image_data, origin='lower')