1.8. Reading images

Astropy provides a few ways to read in FITS images, some in the core package and others in affiliated packages.

Before exploring those, we’ll create a set of (fake) images to work with.

from pathlib import Path
from astropy.nddata import CCDData
from astropy.io import fits

1.8.1. Working with directories

The cell below contains the path to the images. In this notebook we’ll use it both to store the fake images we generate and to read images. In normal use, you wouldn’t start by writing images there, however.

If the images are in the same directory as the notebook, you can omit this or set it to an empty string ''. Having images in the same directory as the notebook is less complicated, but it’s not at all uncommon to need to work with images in a different directory.

Later, we’ll look at how to generate the full path to an image (directory plus file name) in a way that will work on any platform. One of the approaches to loading images (using ccdproc.ImageFileCollection) lets you mostly forget about this.

data_directory = 'path/to/my/images'

1.8.2. Generate some fake images

The cells below generate some fake images to use later in the notebook.

from pathlib import Path
from itertools import cycle

import numpy as np

image_path = Path(data_directory)

image_path.mkdir(parents=True, exist_ok=True)

images_to_generate = {
    'BIAS': 5,
    'DARK': 10,
    'FLAT': 3,
    'LIGHT': 10
}

exposure_times = {
    'BIAS': [0.0],
    'DARK': [5.0, 30.0],
    'FLAT': [5.0, 6.1, 7.3],
    'LIGHT': [30.0],
}

filters = {
    'FLAT': 'V',
    'LIGHT': 'V'
}

objects = {
    'LIGHT': ['m82', 'xx cyg']
}

image_size = [300, 200]

image_number = 0
for image_type, num in images_to_generate.items():
    exposures = cycle(exposure_times[image_type])
    try:
        filts = cycle(filters[image_type])
    except KeyError:
        filts = []
    
    try:
        objs = cycle(objects[image_type])
    except KeyError:
        objs = []
    for _ in range(num):
        img = CCDData(data=np.random.randn(*image_size), unit='adu')
        img.meta['IMAGETYP'] = image_type
        img.meta['EXPOSURE'] = next(exposures)
        if filts:
            img.meta['FILTER'] = next(filts)
        if objs:
            img.meta['OBJECT'] = next(objs)
        image_name = str(image_path / f'img-{image_number:04d}.fits')
        img.write(image_name)
        print(image_name)
        image_number += 1
path/to/my/images/img-0000.fits
path/to/my/images/img-0001.fits
path/to/my/images/img-0002.fits
path/to/my/images/img-0003.fits
path/to/my/images/img-0004.fits
path/to/my/images/img-0005.fits
path/to/my/images/img-0006.fits
path/to/my/images/img-0007.fits
path/to/my/images/img-0008.fits
path/to/my/images/img-0009.fits
path/to/my/images/img-0010.fits
path/to/my/images/img-0011.fits
path/to/my/images/img-0012.fits
path/to/my/images/img-0013.fits
path/to/my/images/img-0014.fits
path/to/my/images/img-0015.fits
path/to/my/images/img-0016.fits
path/to/my/images/img-0017.fits
path/to/my/images/img-0018.fits
path/to/my/images/img-0019.fits
path/to/my/images/img-0020.fits
path/to/my/images/img-0021.fits
path/to/my/images/img-0022.fits
path/to/my/images/img-0023.fits
path/to/my/images/img-0024.fits
path/to/my/images/img-0025.fits
path/to/my/images/img-0026.fits
path/to/my/images/img-0027.fits

1.8.3. Option 1: Reading a single image with astropy.io.fits

This option gives you the most flexibility but is the least adapted to CCD images specifically. What you read in is a list of FITS extensions; you must first select the one you want then access the data or header as desired.

We’ll open up the first of the fake images, img-0001.fits. To combine that with the directory name we’ll use Python 3’s pathlib, which ensures that the path combination will work on Windows too.

image_name = 'img-0001.fits'

image_path = Path(data_directory) / image_name

hdu_list = fits.open(image_path)
hdu_list.info()
Filename: path/to/my/images/img-0001.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       8   (200, 300)   float64   

The hdu_list is a list of FITS Header-Data Units. In this case there is just one, containing both the image header and data, which can be accessed as shown below.

hdu = hdu_list[0]
hdu.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  200                                                  
NAXIS2  =                  300                                                  
IMAGETYP= 'BIAS    '                                                            
EXPOSURE=                  0.0                                                  
BUNIT   = 'adu     '                                                            
hdu.data
array([[-1.32065570e+00, -8.21412051e-01,  1.56937913e-01, ...,
         8.64193546e-01, -8.20851522e-01, -4.31169944e-01],
       [-1.20840212e+00,  7.04716792e-01,  6.17997114e-01, ...,
        -1.08919523e+00, -1.09153665e+00,  1.28581634e+00],
       [-9.42112144e-01, -1.49876744e+00,  5.18880167e-02, ...,
        -2.24826273e+00,  1.77236691e+00,  1.71852854e+00],
       ...,
       [-1.38887788e-01, -6.37979471e-01,  3.63674915e-01, ...,
         1.49142010e+00, -8.68582682e-01,  1.77105366e-01],
       [ 3.26199916e-01,  9.28939953e-01, -1.12626864e+00, ...,
        -4.33055439e-01,  1.41004370e+00, -1.58166673e-01],
       [ 1.71722646e+00,  1.48466768e-03,  8.14401080e-03, ...,
         6.82543876e-01, -1.70023591e-01, -1.00105711e+00]])

The documentation for io.fits describes more of its capabilities.

1.8.4. Option 2: Use CCDData to read in a single image

Astropy contains a CCDData object for representing a single image. It’s not as flexible as using astrop.io.fits directly (for example, it assumes there is only one FITS extension and that it contains image data) but it sets up several properties that make the data easier to work with.

We’ll read in the same single image we did in the example above, img-0001.fits.

ccd = CCDData.read(image_path)

The data and header are accessed similarly to how you access it in an HDU returned by astropy.io.fits:

ccd.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  200                                                  
NAXIS2  =                  300                                                  
IMAGETYP= 'BIAS    '                                                            
EXPOSURE=                  0.0                                                  
BUNIT   = 'adu     '                                                            
ccd.data
array([[-1.32065570e+00, -8.21412051e-01,  1.56937913e-01, ...,
         8.64193546e-01, -8.20851522e-01, -4.31169944e-01],
       [-1.20840212e+00,  7.04716792e-01,  6.17997114e-01, ...,
        -1.08919523e+00, -1.09153665e+00,  1.28581634e+00],
       [-9.42112144e-01, -1.49876744e+00,  5.18880167e-02, ...,
        -2.24826273e+00,  1.77236691e+00,  1.71852854e+00],
       ...,
       [-1.38887788e-01, -6.37979471e-01,  3.63674915e-01, ...,
         1.49142010e+00, -8.68582682e-01,  1.77105366e-01],
       [ 3.26199916e-01,  9.28939953e-01, -1.12626864e+00, ...,
        -4.33055439e-01,  1.41004370e+00, -1.58166673e-01],
       [ 1.71722646e+00,  1.48466768e-03,  8.14401080e-03, ...,
         6.82543876e-01, -1.70023591e-01, -1.00105711e+00]])

There are a number of features of CCDData that make it convenient for working with WCS, slicing, and more. Some of those features will be discussed in more detail in the notebooks that follow.

1.8.5. Option 3: Working with a directory of images using ImageFileCollection

The affiliated package ccdproc provides an easier way to work with collections of images in a directory: an ImageFileCollection. The ImageFileCollection is initialized with the name of the directory containing the images.

from ccdproc import ImageFileCollection
im_collection = ImageFileCollection(data_directory)

Note that we didn’t need to worry about using pathlib to combine the directory and file name, instead we give the collection the name of the directory.

1.8.5.1. Summary of directory contents

The summary property provides an overview of the files in the directory: it’s an astropy Table, so you can access columns in the usual way.

im_collection.summary
Table masked=True length=28
filesimplebitpixnaxisnaxis1naxis2imagetypexposurebunitfilterobject
str13boolint64int64int64int64str5float64str3objectobject
img-0000.fitsTrue-642200300BIAS0.0adu----
img-0001.fitsTrue-642200300BIAS0.0adu----
img-0002.fitsTrue-642200300BIAS0.0adu----
img-0003.fitsTrue-642200300BIAS0.0adu----
img-0004.fitsTrue-642200300BIAS0.0adu----
img-0005.fitsTrue-642200300DARK5.0adu----
img-0006.fitsTrue-642200300DARK30.0adu----
img-0007.fitsTrue-642200300DARK5.0adu----
img-0008.fitsTrue-642200300DARK30.0adu----
img-0009.fitsTrue-642200300DARK5.0adu----
.................................
img-0018.fitsTrue-642200300LIGHT30.0aduVm82
img-0019.fitsTrue-642200300LIGHT30.0aduVxx cyg
img-0020.fitsTrue-642200300LIGHT30.0aduVm82
img-0021.fitsTrue-642200300LIGHT30.0aduVxx cyg
img-0022.fitsTrue-642200300LIGHT30.0aduVm82
img-0023.fitsTrue-642200300LIGHT30.0aduVxx cyg
img-0024.fitsTrue-642200300LIGHT30.0aduVm82
img-0025.fitsTrue-642200300LIGHT30.0aduVxx cyg
img-0026.fitsTrue-642200300LIGHT30.0aduVm82
img-0027.fitsTrue-642200300LIGHT30.0aduVxx cyg

1.8.5.2. Filtering and iterating over images

The great thing about ImageFileCollection is that it provides convenient ways to filter or loop over files via FITS header keyword values.

For example, looping over just the flat files is one line of code:

for a_flat in im_collection.hdus(imagetyp='FLAT'):
    print(a_flat.header['EXPOSURE'])
5.0
6.1
7.3

Instead of iterating over HDUs, as in the example above, you can iterate over just the headers (with .headers) or just the data (with .data). You can use any FITS keyword from the header as a keyword for selecting the images you want. In addition, you can return the file name while also iterating.

for a_flat, fname in im_collection.hdus(imagetyp='LIGHT', object='m82', return_fname=True):
    print(f'In file {fname} the exposure is:', a_flat.header['EXPOSURE'], 'with standard deviation ', a_flat.data.std())
In file img-0018.fits the exposure is: 30.0 with standard deviation  0.9973197037462976
In file img-0020.fits the exposure is: 30.0 with standard deviation  1.0023331489148302
In file img-0022.fits the exposure is: 30.0 with standard deviation  1.004231580886168
In file img-0024.fits the exposure is: 30.0 with standard deviation  0.9990879465700109
In file img-0026.fits the exposure is: 30.0 with standard deviation  1.0014396264899372

The documentation for ImageFileCollection describes more of its capabilities. ImageFileCollection can automatically save a copy of each image as you iterate over them, for example.

for a_flat, fname in im_collection.ccds(bunit='ADU', return_fname=True):
    print(a_flat.unit)
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
adu
a_flat.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  200                                                  
NAXIS2  =                  300                                                  
IMAGETYP= 'LIGHT   '                                                            
EXPOSURE=                 30.0                                                  
FILTER  = 'V       '                                                            
OBJECT  = 'xx cyg  '                                                            
BUNIT   = 'adu     '