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

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


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 )

Opening FITS files and loading the image data

I will open the FITS file and find out what it contains.


In [5]:
hdu_list = fits.open(image_file)
hdu_list.info()


Filename: /Users/adam/.astropy/cache/download/2c9202ae878ecfcb60878ceb63837f5f
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU     161   (891, 893)   int16   
1    er.mask     TableHDU        25   1600R x 4C   [F6.2, F6.2, F6.2, F6.2]   

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)


<type 'numpy.ndarray'>
(893, 891)

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]:
SIMPLE  =                    T /FITS: Compliance                                
BITPIX  =                   16 /FITS: I*2 Data                                  
NAXIS   =                    2 /FITS: 2-D Image Data                            
NAXIS1  =                  891 /FITS: X Dimension                               
NAXIS2  =                  893 /FITS: Y Dimension                               
EXTEND  =                    T /FITS: File can contain extensions               
DATE    = '2014-01-09        '  /FITS: Creation Date                            
ORIGIN  = 'STScI/MAST'         /GSSS: STScI Digitized Sky Survey                
SURVEY  = 'SERC-ER '           /GSSS: Sky Survey                                
REGION  = 'ER768   '           /GSSS: Region Name                               
PLATEID = 'A0JP    '           /GSSS: Plate ID                                  
SCANNUM = '01      '           /GSSS: Scan Number                               
DSCNDNUM= '00      '           /GSSS: Descendant Number                         
TELESCID=                    4 /GSSS: Telescope ID                              
BANDPASS=                   36 /GSSS: Bandpass Code                             
COPYRGHT= 'AAO/ROE '           /GSSS: Copyright Holder                          
SITELAT =              -31.277 /Observatory: Latitude                           
SITELONG=              210.934 /Observatory: Longitude                          
TELESCOP= 'UK Schmidt - Doubl' /Observatory: Telescope                          
INSTRUME= 'Photographic Plate' /Detector: Photographic Plate                    
EMULSION= 'IIIaF   '           /Detector: Emulsion                              
FILTER  = 'OG590   '           /Detector: Filter                                
PLTSCALE=                67.20 /Detector: Plate Scale arcsec per mm             
PLTSIZEX=              355.000 /Detector: Plate X Dimension mm                  
PLTSIZEY=              355.000 /Detector: Plate Y Dimension mm                  
PLATERA =        85.5994550000 /Observation: Field centre RA degrees            
PLATEDEC=       -4.94660910000 /Observation: Field centre Dec degrees           
PLTLABEL= 'OR14052 '           /Observation: Plate Label                        
DATE-OBS= '1990-12-22T13:49:00' /Observation: Date/Time                         
EXPOSURE=                 65.0 /Observation: Exposure Minutes                   
PLTGRADE= 'AD2     '           /Observation: Plate Grade                        
OBSHA   =             0.158333 /Observation: Hour Angle                         
OBSZD   =              26.3715 /Observation: Zenith Distance                    
AIRMASS =              1.11587 /Observation: Airmass                            
REFBETA =        66.3196420000 /Observation: Refraction Coeff                   
REFBETAP=     -0.0820000000000 /Observation: Refraction Coeff                   
REFK1   =        6423.52290000 /Observation: Refraction Coeff                   
REFK2   =       -102122.550000 /Observation: Refraction Coeff                   
CNPIX1  =                12237 /Scan: X Corner                                  
CNPIX2  =                19965 /Scan: Y Corner                                  
XPIXELS =                23040 /Scan: X Dimension                               
YPIXELS =                23040 /Scan: Y Dimension                               
XPIXELSZ=              15.0295 /Scan: Pixel Size microns                        
YPIXELSZ=              15.0000 /Scan: Pixel Size microns                        
PPO1    =       -3069417.00000 /Scan: Orientation Coeff                         
PPO2    =       0.000000000000 /Scan: Orientation Coeff                         
PPO3    =        177500.000000 /Scan: Orientation Coeff                         
PPO4    =       0.000000000000 /Scan: Orientation Coeff                         
PPO5    =        3069417.00000 /Scan: Orientation Coeff                         
PPO6    =        177500.000000 /Scan: Orientation Coeff                         
PLTRAH  =                    5 /Astrometry: Plate Centre H                      
PLTRAM  =                   42 /Astrometry: Plate Centre M                      
PLTRAS  =                23.86 /Astrometry: Plate Centre S                      
PLTDECSN= '-       '           /Astrometry: Plate Centre +/-                    
PLTDECD =                    4 /Astrometry: Plate Centre D                      
PLTDECM =                   56 /Astrometry: Plate Centre M                      
PLTDECS =                 47.9 /Astrometry: Plate Centre S                      
EQUINOX =               2000.0 /Astrometry: Equinox                             
AMDX1   =        67.1550859799 /Astrometry: GSC1 Coeff                          
AMDX2   =      0.0431478884485 /Astrometry: GSC1 Coeff                          
AMDX3   =       -292.435619180 /Astrometry: GSC1 Coeff                          
AMDX4   =  -2.68934864702E-005 /Astrometry: GSC1 Coeff                          
AMDX5   =   1.99133423290E-005 /Astrometry: GSC1 Coeff                          
AMDX6   =  -2.37011931379E-006 /Astrometry: GSC1 Coeff                          
AMDX7   =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX8   =   2.21426387429E-006 /Astrometry: GSC1 Coeff                          
AMDX9   =  -8.12841581455E-008 /Astrometry: GSC1 Coeff                          
AMDX10  =   2.48169090021E-006 /Astrometry: GSC1 Coeff                          
AMDX11  =   2.77618933926E-008 /Astrometry: GSC1 Coeff                          
AMDX12  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX13  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX14  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX15  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX16  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX17  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX18  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX19  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDX20  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY1   =        67.1593591466 /Astrometry: GSC1 Coeff                          
AMDY2   =     -0.0471363749174 /Astrometry: GSC1 Coeff                          
AMDY3   =        316.004963520 /Astrometry: GSC1 Coeff                          
AMDY4   =   2.86798151430E-005 /Astrometry: GSC1 Coeff                          
AMDY5   =  -2.00968236347E-005 /Astrometry: GSC1 Coeff                          
AMDY6   =   2.27840393227E-005 /Astrometry: GSC1 Coeff                          
AMDY7   =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY8   =   2.23885090381E-006 /Astrometry: GSC1 Coeff                          
AMDY9   =  -2.28360163464E-008 /Astrometry: GSC1 Coeff                          
AMDY10  =   2.44828851495E-006 /Astrometry: GSC1 Coeff                          
AMDY11  =  -5.76717487998E-008 /Astrometry: GSC1 Coeff                          
AMDY12  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY13  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY14  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY15  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY16  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY17  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY18  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY19  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDY20  =       0.000000000000 /Astrometry: GSC1 Coeff                          
AMDREX1 =        67.1532034737 /Astrometry: GSC2 Coeff                          
AMDREX2 =      0.0434354199559 /Astrometry: GSC2 Coeff                          
AMDREX3 =       -292.435438892 /Astrometry: GSC2 Coeff                          
AMDREX4 =   4.60919247070E-006 /Astrometry: GSC2 Coeff                          
AMDREX5 =  -3.21138058537E-006 /Astrometry: GSC2 Coeff                          
AMDREX6 =   7.23651736725E-006 /Astrometry: GSC2 Coeff                          
AMDREX7 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX8 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX9 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX10=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX11=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX12=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX13=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX14=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX15=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX16=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX17=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX18=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX19=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREX20=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY1 =        67.1522589487 /Astrometry: GSC2 Coeff                          
AMDREY2 =     -0.0481758265285 /Astrometry: GSC2 Coeff                          
AMDREY3 =        315.995683716 /Astrometry: GSC2 Coeff                          
AMDREY4 =  -7.47397531230E-006 /Astrometry: GSC2 Coeff                          
AMDREY5 =   9.55221105409E-007 /Astrometry: GSC2 Coeff                          
AMDREY6 =   7.60954485251E-006 /Astrometry: GSC2 Coeff                          
AMDREY7 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY8 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY9 =       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY10=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY11=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY12=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY13=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY14=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY15=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY16=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY17=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY18=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY19=       0.000000000000 /Astrometry: GSC2 Coeff                          
AMDREY20=       0.000000000000 /Astrometry: GSC2 Coeff                          
ASTRMASK= 'er.mask '           /Astrometry: GSC2 Mask                           
WCSAXES =                    2 /GetImage: Number WCS axes                       
WCSNAME = 'DSS               ' /GetImage: Local WCS approximation from full plat
RADESYS = 'ICRS              ' /GetImage: GSC-II calibration using ICRS system  
CTYPE1  = 'RA---TAN          ' /GetImage: RA-Gnomic projection                  
CRPIX1  =           446.000000 /GetImage: X reference pixel                     
CRVAL1  =            85.274970 /GetImage: RA of reference pixel                 
CUNIT1  = 'deg               ' /GetImage: degrees                               
CTYPE2  = 'DEC--TAN          ' /GetImage: Dec-Gnomic projection                 
CRPIX2  =           447.000000 /GetImage: Y reference pixel                     
CRVAL2  =            -2.458265 /GetImage: Dec of reference pixel                
CUNIT2  = 'deg               ' /Getimage: degrees                               
CD1_1   =        -0.0002802651 /GetImage: rotation matrix coefficient           
CD1_2   =         0.0000003159 /GetImage: rotation matrix coefficient           
CD2_1   =         0.0000002767 /GetImage: rotation matrix coefficient           
CD2_2   =         0.0002798187 /GetImage: rotation matrix coefficient           
OBJECT  = 'data              ' /GetImage: Requested Object Name                 
DATAMIN =                 3759 /GetImage: Minimum returned pixel value          
DATAMAX =                22918 /GetImage: Maximum returned pixel value          
OBJCTRA = '05 41 06.000      ' /GetImage: Requested Right Ascension (J2000)     
OBJCTDEC= '-02 27 30.00      ' /GetImage: Requested Declination (J2000)         
OBJCTX  =             12682.48 /GetImage: Requested X on plate (pixels)         
OBJCTY  =             20411.37 /GetImage: Requested Y on plate (pixels)         

SHORTCUT

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)


<type 'numpy.ndarray'>
(893, 891)

Similarly, if you only want the header, you can use fits.getheader:


In [11]:
header = fits.getheader(image_file)

Viewing the image data and getting basic statistics


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]:
<matplotlib.colorbar.Colorbar instance at 0x1092198c0>

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))


('Min:', 3759)
('Max:', 22918)
('Mean:', 9831.4816762875744)
('Stdev:', 3032.3927542049046)

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.


In [14]:
print(type(image_data.flat))


<type 'numpy.flatiter'>

In [15]:
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.


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]:
[<matplotlib.text.Text at 0x10eb72b10>,
 <matplotlib.text.Text at 0x10ec49610>,
 <matplotlib.text.Text at 0x10f10ef90>]

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.


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]:
<matplotlib.colorbar.Colorbar instance at 0x11270f878>

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.


In [22]:
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/


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)


(85.385675186882054, -2.5551867109375919)

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)


[[ 85.39407108  -2.57478146]
 [ 85.38564882  -2.58037026]
 [ 85.25946264  -2.52708506]]

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'))


<SkyCoord (FK5): equinox=J2000.000, (ra, dec) in deg
    [(85.39407107595326, -2.574781460695727),
     (85.38564881637622, -2.580370256991026),
     (85.25946263858953, -2.5270850561642217)]>
[u'05h41m34.5771s -02d34m29.2133s', u'05h41m32.5557s -02d34m49.3329s', u'05h41m02.271s -02d31m37.5062s']

You can also use WCS to get pixel coordinates for overplotting sources.


In [28]:
horsehead = coordinates.SkyCoord.from_name('Horsehead Nebula')
print(horsehead)


<SkyCoord (ICRS): ra=85.24583 deg, dec=-2.45833 deg>

In [29]:
(hxpix,hypix), = mywcs.wcs_world2pix([(horsehead.ra.deg,
                                       horsehead.dec.deg)],0)
print(hxpix,hypix)


(548.87693013963485, 445.66385271710652)

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]:
[<matplotlib.lines.Line2D at 0x11eb06e50>]

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 the world coordinate system

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)


INFO:astropy:Auto-setting vmin to  3.634e+03
INFO:astropy:Auto-setting vmax to  1.940e+04
INFO: Auto-setting vmin to  3.634e+03 [aplpy.core]
INFO: Auto-setting vmax to  1.940e+04 [aplpy.core]

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')