Visualizing Astronomical Images

Visualizing Astronomical Images#

Imad Pasha, Yasmeen Asali, William Cerny, Sebastian Monzon (Yale University)

Description: Practice writing Python functions + visualization of astronomical imaging

Intended Audience: Beginner to Intermediate Undergraduate

tags: library:matplotlib, control-flow:functions, visualization, astronomical-images

Requirements: requirements.txt

Last Updated: July 23, 2024

Learning Objectives

  • Write basic functions in Python

  • Be able to appropriately scale and view FITS images

Introduction#

Images taken by telescopes (both in space and on the ground) are a crucial element of astronomical research, providing a spatially resolved view of a region (generally at a particular wavelength). These images are predominantly stored in an astronomy-specific file type known as fits — flexible image transport system (though tables can also be stored in fits files)

The astropy.io module has a useful submodule for working with fits files, which can be imported as

from astropy.io import fits

Within a fits file, information is stored in what are called extensions. These are like separate buckets/folders, and thus allow one to store multiple images, or an image + table, or multiple tables, within 1 file. But often, we have a simple file with 1 image in the 0th extension. Every extension present in a fits file has two attributes: a .header attribute, and a .data attribute. The header is like a dictionary, storing ancillary information in key-value pairs, while the data attribute stores the image, table, etc.

When we load a fits file, we usually also want to know something about the coordinates (on the sky, so RA and DEC) of the objects in our image. To do this, we need to know the WCS (world coordinate system) information for the image, which provides a translation between position in pixel (image) coordinates and sky coordinates. For many astronomical images you will work with, this information has been added to the .header attribute for the extension of the image. The astropy.wcs.WCS class can be given such a header, and compute a wcs object that will be useful when we go to plot or determine locations in our image. We can import it via

from astropy.wcs import WCS

So how do we actually open a fits file and retrieve the image and header? We’ll use a pythonic structure known as a context manager. The details of context managers aren’t super important right now, but the key point is this: context managers let us open contexts and close them automatically. For files, this means we can open a fits file, extract what we need, and close the file automatically.

Opening a Fits File#

The fits file we will be using here is an image of M33 taken by the Dragonfly Spectral Line Mapper.

You can download it here (right click and save as a .fits file).

Exercise 1

To open a simple image fits file using a context manager, we do something like the following:

from astropy.io import fits 
from astropy.wcs import WCS 

with fits.open('./M33_halpha.fits') as hdu:
    im = hdu[0].data 
    wcs = WCS(hdu[0].header)

  1. Using the M33_halpha.fits file provided, use the above snippet in your own code to open the image.

What happened in the snippet?

  • First, we said “with fits.open()” to establish it as a context manager, provided the file path, and gave the overall object the name hdu

  • Next, within the block, we set im to be the hdu indexed at 0 (0th extension) and used dot-notation to extract the .data attribute

  • Finally, we set a wcs variable to be the astropy WCS class, given the same extension but the .header attribute.

Writing an Image Opening Function#

We can, if we wish, encapsulate this operation into a function. It’s short enough that it’s not overly onerous to write out each time. But we can make these three lines into one, which might make our code easier to read when loading many different fits files throughout a script.

Recall that a function in Python is, like a mathematical function, designed to take input arguments (when needed) and return values back. So \(y = sin(x)\) is a function which takes \(x\) as an argument and returns something, which has been named \(y\) here. The default structure of a simple python function is

def function_name(arg1_name,arg2_name):
    computation = arg1_name * arg2_name 
    return computation

Here we use def to define a function (lines within the function are indented), and within the function, we assign inputs names to use when computing things, before returning some final value(s) if required.

Exercise 2

  1. Start by defining a function load_fits which has arguments fname (str), and extension (int). These are the pieces of information we needed above.

  2. Then add these lines above into the function trading out the hard coded elements for the function arguments.

  3. Finally, return im, wcs. As the majority of fits image store the image in the 0th extension, you can set the default for that argument to 0.

Warning

Remember! Don’t forget to document your function!

Visualizing the Image#

Let’s visualize the image we have loaded.

Exercise 3

  1. Use the matplotlib imshow command to plot the M33 image. To see anything useful, you will need to set the vmin and vmax arguments — a convenient choice is to compute the np.percentile of the data at [1,99] and to set those as the min and max scalings.

  2. Additionally set origin='lower' and cmap='viridis' in your imshow call.

I suggest a figure size of 10 by 10 inches. Remember we can call fig, ax = plt.subplots() to make our figure and axes objects (so here, fig,ax=plt.subplots(figsize=(10,10))). We then use ax.imshow(...) to add images to the axes.

With the right scaling, it might look something like this:

A function for Image Viewing#

The code to visualize the image above is already a few lines of code, and we often need to plot astronomical images. Let’s now make a function that can streamline some of these lines of code above for us.

Exercise 4

  1. First define a function implot. The first argument should be a variable path_or_pixels. We’re going to allow users to either provide a file name to a fits file directly, or provide an image array (as we did above). To handle this, write some logic using the isinstance() python function to check whether the input is an array or a string, and if it is a string, use the load_fits function we wrote in Exercise 1 to load the image. (Otherwise, we’ll just set our image to be the path_or_pixels input, assuming it is an array. You could use isinstance() again to check if it is a numpy array and raise an exception if not, optionally.)

  2. Include as arguments the following:

  • percentiles (list) (which will let users scale vmin and vmax, set to a default of [1,99.5]).

  • figsize (tuple): an optional argument with a default set to something like (8,8).

  • cmap (str): argument with a default of 'gray_r', passed along the imshow command, giving us control over the plotting color mapping.

Your function should return the fig and ax objects it created.

Hint: isinstance()

If you haven’t used isinstance() yet, here’s a primer.

When we want to check the type of a variable, we can use type(var_name) to see its typename. When we want to check or ensure that a variable’s type matches some desired type, we can use isinstance(variable,type) to return True or False. For example, isinstance('test',str) would return true, because the input 'test' is a string. The names for the built-in types are str,bool,list,dict,tuple,float,int, etc. For a numpy array, its “type” is np.ndarray.

  1. Try out your function to make sure it is working as expected — try inputting both our image array from above, as well as the file name. Try changing the percentiles and figsize and make sure the function responds accordingly.

Next we want to incorporate the wcs into the mix, allowing us to easily view an image with axes either in pixels or in astronomical coordinates.

Exercise 5

  1. Currently, a wcs will be accessible if the user provides the path to a file that has one. Now, add an optional argument to your function wcs with a default of None — this way, if a user has read in an image themselves (as we did above) they can provide both the array an the wcs. Add logic so that if the wcs is not None (either from loading an image, or after the function input), we add the following to our subplots command:

fig, ax = plt.subplots(figsize=figsize,subplot_kw={'projection':wcs})

What this does is turn our regular axis into a WCSAxis, which can handle plotting astronomical units along the axes of our image.

  1. Try your function again, but this time add the wcs we created when loading the image. How do your plot axes change?

  2. At the moment, the user chooses whether to provide a wcs for the pixels case, but we force a wcs plotting if a path to a fits file that contains a wcs is provided. This should be the default case, but it might be nice to give the option to turn that off. The easiest way to accomplish this is to add an ignore_wcs optional argument to your function which has a default of False. Thus, given an image file path, we can plot using coordinates by default, or the raw pixels by ignoring any wcs present. This works too if we provide a wcs but want to easily toggle it on and off.

  3. Try out your function again, using both the fits file string and the array, but now try all your toggles — i.e., fits file string, fits file string plus ‘ignore_wcs’ if you added that feature, array input, array input + wcs object input. Do you get the proper behavior? Below is an example of what the output should look like for default loading of the image.

There are now a decent number of lines of logic and checking in our function, and hopefully it is becoming clear why making such a function is so useful — we can now very quickly and easily load and plot any fits file, or any array we have on hand, with wcs on or off, with limits setting, all in one line.

Note

It is a good idea to return the fig and ax objects at the end of our function. This is useful because if we want to add to our plot (say, add a circle somewhere, or text), or change our plot (adjust the limits, etc.) we need to be able to access those objects. As we learned when we discussed scope, anything from inside a function that we want access to after the function has been run must be returned back to the larger scope of the python script we are working in.

Contents