Reading images
Contents
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.93231885, -0.20554953, 1.11382749, ..., -0.48266713,
0.62491939, 0.32332385],
[-0.62965729, 0.25235816, 0.42799246, ..., 0.21027981,
-0.02531951, -0.03524356],
[-0.07380972, -1.30615523, 0.08721875, ..., -1.87320596,
-1.29376234, -0.24861751],
...,
[ 1.50038478, -2.74012182, 1.15403477, ..., 0.78683255,
-1.00991427, 0.26186583],
[ 0.43151481, 0.1382379 , 0.38706353, ..., 1.15553082,
-1.96669301, 0.11551855],
[ 1.23072078, -1.51645142, -0.26927057, ..., 1.10289169,
-0.85284388, 1.12601521]], dtype='>f8')
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.93231885, -0.20554953, 1.11382749, ..., -0.48266713,
0.62491939, 0.32332385],
[-0.62965729, 0.25235816, 0.42799246, ..., 0.21027981,
-0.02531951, -0.03524356],
[-0.07380972, -1.30615523, 0.08721875, ..., -1.87320596,
-1.29376234, -0.24861751],
...,
[ 1.50038478, -2.74012182, 1.15403477, ..., 0.78683255,
-1.00991427, 0.26186583],
[ 0.43151481, 0.1382379 , 0.38706353, ..., 1.15553082,
-1.96669301, 0.11551855],
[ 1.23072078, -1.51645142, -0.26927057, ..., 1.10289169,
-0.85284388, 1.12601521]], dtype='>f8')
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
file | simple | bitpix | naxis | naxis1 | naxis2 | imagetyp | exposure | bunit | filter | object |
---|---|---|---|---|---|---|---|---|---|---|
str13 | bool | int64 | int64 | int64 | int64 | str5 | float64 | str3 | object | object |
img-0000.fits | True | -64 | 2 | 200 | 300 | BIAS | 0.0 | adu | -- | -- |
img-0001.fits | True | -64 | 2 | 200 | 300 | BIAS | 0.0 | adu | -- | -- |
img-0002.fits | True | -64 | 2 | 200 | 300 | BIAS | 0.0 | adu | -- | -- |
img-0003.fits | True | -64 | 2 | 200 | 300 | BIAS | 0.0 | adu | -- | -- |
img-0004.fits | True | -64 | 2 | 200 | 300 | BIAS | 0.0 | adu | -- | -- |
img-0005.fits | True | -64 | 2 | 200 | 300 | DARK | 5.0 | adu | -- | -- |
img-0006.fits | True | -64 | 2 | 200 | 300 | DARK | 30.0 | adu | -- | -- |
img-0007.fits | True | -64 | 2 | 200 | 300 | DARK | 5.0 | adu | -- | -- |
img-0008.fits | True | -64 | 2 | 200 | 300 | DARK | 30.0 | adu | -- | -- |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
img-0018.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | m82 |
img-0019.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | xx cyg |
img-0020.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | m82 |
img-0021.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | xx cyg |
img-0022.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | m82 |
img-0023.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | xx cyg |
img-0024.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | m82 |
img-0025.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | xx cyg |
img-0026.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | m82 |
img-0027.fits | True | -64 | 2 | 200 | 300 | LIGHT | 30.0 | adu | V | xx 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 1.0006305082274918
In file img-0020.fits the exposure is: 30.0 with standard deviation 0.9951272812922944
In file img-0022.fits the exposure is: 30.0 with standard deviation 1.0013181677514156
In file img-0024.fits the exposure is: 30.0 with standard deviation 0.9987948739214904
In file img-0026.fits the exposure is: 30.0 with standard deviation 1.0047132884527286
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 '