In [1]:
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.
In [2]:
from astropy.utils.data import download_file
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
In [3]:
from astropy.io import fits
# import pyfits
In [4]:
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )
I will open the FITS file and find out what it contains.
In [5]:
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.
In [35]:
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.
In [7]:
print(type(image_data))
print(image_data.shape)
The Header is also an important component. The info above tells us it has 161 cards:
In [8]:
header = hdu_list['PRIMARY'].header
In [9]:
header
Out[9]:
If you don't need to examine the FITS header, you can call fits.getdata to bypass the previous steps.
In [10]:
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:
In [11]:
header = fits.getheader(image_file)
In [12]:
plt.imshow(image_data, cmap='gray')
plt.colorbar()
# To see more color maps
# http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
Out[12]:
Let's get some basic statistics about our image. If your image contained NaNs, you'd use np.nanmin etc below
In [40]:
print('Min:', np.min(image_data))
print('Max:', np.max(image_data))
print('Mean:', np.mean(image_data))
print('Stdev:', np.std(image_data))
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.
In [14]:
print(type(image_data.flat))
In [15]:
NBINS = 1000
histogram = plt.hist(image_data.flat, NBINS)
Want to use a logarithmic color scale? To do so we need to load the LogNorm object from matplotlib.
In [16]:
from matplotlib.colors import LogNorm
In [17]:
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'])
Out[17]:
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.
In [18]:
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.
In [19]:
# 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.
In [20]:
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.
In [47]:
plt.imshow(final_image, cmap='gray', vmin=400, vmax=600)
plt.colorbar()
Out[47]:
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.
In [22]:
outfile = 'stacked_M13_blue.fits'
hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, clobber=True)
"World coordinates" refer to RA/Dec or Galactic latitude/longitude.
The full documentation is available at: http://astropy.readthedocs.org/en/latest/wcs/
In [23]:
from astropy import wcs
Step 1. Create a WCS object
In [24]:
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:
In [25]:
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:
In [26]:
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:
In [54]:
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.
In [28]:
horsehead = coordinates.SkyCoord.from_name('Horsehead Nebula')
print(horsehead)
In [29]:
(hxpix,hypix), = mywcs.wcs_world2pix([(horsehead.ra.deg,
horsehead.dec.deg)],0)
print(hxpix,hypix)
In [30]:
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')
Out[30]:
You can also create an HDU 'from scratch', but with a WCS header:
In [58]:
myhdu = fits.PrimaryHDU(data=image_data, header=mywcs.to_header())
Plotting using WCS is a more complicated task and should be left to specialized tools:
In [55]:
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)
In [ ]:
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')