A First Look at an X-ray Image Dataset

Images are data. They can be 2D, from cameras, or 1D, from spectrographs, or 3D, from IFUs (integral field units). In each case, the data come packaged as an array of numbers, which we can visualize, and do calculations with.

Let's suppose we are interested in clusters of galaxies. We choose one, Abell 1835, and propose to observe it with the XMM-Newton space telescope. We are successful, we design the observations, and they are taken for us. Next: we download the data, and take a look at it.

Getting the Data

We will download our images from HEASARC, the online archive where XMM data are stored.



In [1]:

    
import astropy.io.fits as pyfits
import numpy as np
import os.path
import astropy.visualization as viz
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 10.0)

Download the example data files if we don't already have them.



In [2]:

    
targdir = 'a1835_xmm/'
!mkdir -p $targdir

imagefile = 'P0098010101M2U009IMAGE_3000.FTZ'
expmapfile = 'P0098010101M2U009EXPMAP3000.FTZ'
bkgmapfile = 'P0098010101M2X000BKGMAP3000.FTZ'

remotedir = 'http://heasarc.gsfc.nasa.gov/FTP/xmm/data/rev0/0098010101/PPS/'

for filename in [imagefile,expmapfile,bkgmapfile]:
    path = targdir + filename
    url = remotedir + filename
    if not os.path.isfile(path): # i.e. if the file does not exist already:
        !wget -nd -O $path $url 

imagefile = targdir+imagefile
expmapfile = targdir+expmapfile
bkgmapfile = targdir+bkgmapfile
        
!du -sch $targdir/*









    



4.0K	a1835_xmm//M2ptsrc.txt
0	a1835_xmm//P0098010101M2U009EXPMAP3000.FTZ
0	a1835_xmm//P0098010101M2U009IMAGE_3000.FTZ
0	a1835_xmm//P0098010101M2X000BKGMAP3000.FTZ
4.0K	total

The XMM MOS2 image

Let's find the "science" image taken with the MOS2 camera, and display it.



In [3]:

    
imfits = pyfits.open(imagefile)
imfits.info()









    



---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-3-8cfcefb6411b> in <module>()
----> 1 imfits = pyfits.open(imagefile)
      2 imfits.info()

/usr/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.py in fitsopen(name, mode, memmap, save_backup, cache, **kwargs)
    127         raise ValueError('Empty filename: %s' % repr(name))
    128 
--> 129     return HDUList.fromfile(name, mode, memmap, save_backup, cache, **kwargs)
    130 
    131 

/usr/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.py in fromfile(cls, fileobj, mode, memmap, save_backup, cache, **kwargs)
    269 
    270         return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,
--> 271                              save_backup=save_backup, cache=cache, **kwargs)
    272 
    273     @classmethod

/usr/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.py in _readfrom(cls, fileobj, data, mode, memmap, save_backup, cache, **kwargs)
    859             # raise and exception
    860             if mode in ('readonly', 'denywrite') and len(hdulist) == 0:
--> 861                 raise IOError('Empty or corrupt FITS file')
    862 
    863             # initialize/reset attributes to be used in "update/append" mode

IOError: Empty or corrupt FITS file

imfits is a FITS object, containing multiple data structures. The image itself is an array of integer type, and size 648x648 pixels, stored in the primary "header data unit" or HDU.

If we need it to be floating point for some reason, we need to cast it: im = imfits[0].data.astype('np.float32') Note that this (probably?) prevents us from using the pyfits "writeto" method to save any changes. Assuming the integer type is ok, just get a pointer to the image data.

Accessing the .data member of the FITS object returns the image data as a numpy ndarray.



In [7]:

    
im = imfits[0].data









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-14754cc8b46b> in <module>()
----> 1 im = imfits[0].data

NameError: name 'imfits' is not defined

Let's look at this with ds9.



In [10]:

    
!ds9 -log "$imagefile"

If you don't have the image viewing tool ds9, you should install it - it's very useful astronomical software. You can download it (later!) from this webpage.

We can also display the image in the notebook:



In [8]:

    
plt.imshow(viz.scale_image(im, scale='log', max_cut=40), cmap='gray', origin='lower');









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-3bdb46ca9459> in <module>()
----> 1 plt.imshow(viz.scale_image(im, scale='log', max_cut=40), cmap='gray', origin='lower');

NameError: name 'im' is not defined

Exercise

What is going on in this image?

Make a list of everything that is interesting about this image with your neighbor, and we'll discuss the features you identify in about 5 minutes time.



In [15]:

    
index = np.where(im == np.max(im))
print "image dimensions:",im.shape
print "location of maximum pixel value:",index
print "maximum pixel value: ",im[index]









    



image dimensions: (648, 648)
location of maximum pixel value: (array([348]), array([328]))
maximum pixel value:  [223]

NB. Images read in with pyfits are indexed with eg im[y,x]: ds9 shows that the maximum pixel value is at "image coordinates" x=328, y=348. pyplot knows what to do, but sometimes we may need to take the transpose of the im array. What pyplot does need to be told is that in astronomy, the origin of the image is conventionally taken to be at the bottom left hand corner, not the top left hand corner. That's what the origin=lower in the plt.imshow command was about.

We will work in image coordinates throughout this course, for simplicity. Aligning images on the sky via a "World Coordinate System" is something to be learned elsewhere.



In [ ]: