Images and Image Plotting

Learning Objectives

• Read a FITS file
• Load and display an image
• Use array slicing operations to display part of a cube
• Use SunPy to explore a multi-dimensional array

Arrays as images

All photographic images represent a measurement of how much light hits the receiver. For instance, the Hubble image below is obtained by measuring the brightnesses of distant stars:

With traditional optical cameras, this measurement results in an image which is continuous, as it is projected directly onto paper or some other object. In order to store images digitally, they need to be divided into discrete chunks, pixels, each of which contains the value of the measurement in that small portion of the image. In this representation, an image is simply a grid of numbers, which allows it to be easily stored as an array with a shape equal to the resolution of the image.

SunPy has some example data that we’ll be using to demontrate this principle.

import matplotlib.pyplot as plt
from astropy.io import fits

# Get the example data for this lesson
from astropy.utils.data import download_file
cube_file = download_file("http://data.sunpy.org/sunkit-sst/CRISP_TXY_Cube.fits.gz", cache=True)

The above command will download the example data if it isn’t present, and store it in a FITS file. We will read the data from this file, but for this lesson we will not be considering any coordinate information in the header. We will be covering that tomorrow.

# Check the file exists
cube_file

'/home/stuart/.astropy/cache/download/py3/a2460f944e47f8d2bca6f8e366f2004c'

# Load the data
cube = fits.getdata(cube_file)

This file is a time, solar-x, solar-y cube. To start off with we will be using the first image.

# Check what dimensions our data have so we know how to slice the array
cube.shape

(50, 966, 980)

sunspot = cube[0]
sunspot.shape

(966, 980)

Once read in as an array, the image can be processed in the same ways as any other array. For instance, we can easily find the highest, lowest and mean values of the image, the type of the variables stored in the array, and the resolution of the image:

# Output the image minimum, mean and maximum.
print('Image min:', sunspot.min(), '; Image mean:', sunspot.mean(), '; Image max: ', sunspot.max())

# Output the array dtype.
print('Data type:', sunspot.dtype)

# Output image size.
print('Image size:', sunspot.shape)

Image min: 0.0 ; Image mean: 1875.7292 ; Image max:  8677.0
Data type: >f4
Image size: (966, 980)

Plotting images

While storing an image as a grid of numbers is very useful for analysis, we still need to be able to visually inspect the image. This can be achieved with plt.imshow(), which allocates a colour to every element in the array according to its value.

# Display image array with imshow()
plt.imshow(sunspot)

<matplotlib.image.AxesImage at 0x7fb671537898>

When plotting an image in this way, you will often need to know what actual values correspond to the colours. To find this out, we can draw a colour bar alongside the image which indicates the mapping of values to colours:

plt.imshow(sunspot)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fb671152da0>

Fortunately, matplotlib provides a large variety of colour maps which are suitable for various different purposes (more on this later). plt.imshow() has a cmap keyword argument which can be passed a string defining the desired colour map.

# Display the image with a more monochrome colour map.
plt.imshow(sunspot, cmap='gray')
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fb6710837f0>

The full list of available colour maps (for matplotlib 2.2.3) and some discussion on how to use them can be found here.

# 1

plt.imshow(cube[:,720,:].T, aspect='auto')

<matplotlib.image.AxesImage at 0x7fb670ff2198>

# 2

# Display my image
plt.imshow(cube[:,720,:].T, aspect='auto', cmap='magma')
plt.colorbar()
plt.xlabel("Time")
plt.ylabel("y pixels")

Text(0,0.5,'y pixels')

Value limits

The default behaviour of imshow() in terms of colour mapping is that the colours cover the full range of the data so that the lower end (blue, in the plots above) represents the smallest value in the array, and the higher end (red) represents the greatest value.

This is fine if rest of the values are fairly evenly spaced between these extremes. However, if we have a very low minimum or very high maximum compared to the rest of the image, this default scaling is unhelpful. To deal with this problem, imshow() allows you to set the minimum and maximum values used for the scaling with the vmin and vmax keywords.

plt.imshow(sunspot, cmap='gray', vmin=75, vmax=6000)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fb67162a978>

As you can see, this allows us to increase the contrast of the image at the cost of discounting extreme values, or we can include a broader range of values but see less detail. Similar effects can also be achieved with the norm keyword, which allows you to set how imshow() scales values in order to map them to colours (linear or logarithmic scaling, for example).

Axes

You will notice in the above plots that the axes are labelled with the pixel coordinates of the image. You will also notice that the origin of the axes is in the top left corner rather than the bottom left. This is a convention in image-drawing, but can be changed if necessary by setting the origin keyword to 'lower' when calling imshow():

plt.imshow(sunspot, cmap='gray', origin='lower')

<matplotlib.image.AxesImage at 0x7fb670e5b8d0>

imshow() also allows you to change the upper and lower values of each axis, and the appropriate tick labels will be drawn. This feature can be used to apply physical spatial scales to the image (if you know them) rather than going purely on pixel positions, which may be less useful. This is done with the extent keyword, which takes a list of values corresponding to lower and upper x values and the lower and upper y values (in that order).

plt.imshow(sunspot, cmap='gray', origin='lower', extent=[413, 470, 223, 277])

<matplotlib.image.AxesImage at 0x7fb670e34240>

# 1 

# Display the coins image with adjusted value range
plt.imshow(sunspot, cmap='magma', vmin=60, vmax=5000)
plt.colorbar()

<matplotlib.colorbar.Colorbar at 0x7fb66e7466d8>

Interpolation

In order to display a smooth image, imshow() automatically interpolates to find what values should be displayed between the given data points. The default interpolation scheme is 'linear', which interpolates linearly between points, as you might expect. The interpolation can be changed with yet another keyword in imshow(). Here are a few examples:

fig, ax = plt.subplots(2, 2, figsize=(10, 10))
smallim = cube[10, 250:750, :]
ax[0, 0].imshow(smallim) # Default (linear) interpolation
ax[0, 1].imshow(smallim, interpolation='bicubic') # Bicubic interpolation
ax[1, 0].imshow(smallim, interpolation='nearest') # Nearest-neighbour interpolation
ax[1, 1].imshow(smallim, interpolation='none') # No interpolation

<matplotlib.image.AxesImage at 0x7fb66e630c88>

This can be a useful way to change how the image appears. For instance, if the exact values of the data are extremely important, little or no interpolation may be appropriate so the original values are easier to discern, whereas a high level of interpolation can be used if the smoothness of the image is more important than the actual numbers.

Note that that 'none' in the imshow() call above is NOT the same as None. None tells imshow() you are not passing it a variable for the interpolation keyword, so it uses the default, whereas 'none' explicitly tells it not to interpolate.

SunPy Cube Explorer

SunPy has an ImageAnimator class that will create sliders for each axis not used for plotting the image. By default it uses the last two axes of the array to plot the image.

from sunpy.visualization.imageanimator import ImageAnimator

ImageAnimator(cube, colorbar=True, cmap="plasma")

<sunpy.visualization.imageanimator.ImageAnimator at 0x7fb66dbd1cf8>