Doing Astronomy with Python

One of the things that makes Python so powerful is that you can find a module for almost anything. In this article, I cover Astropy, which was originally developed by the Space Telescope Science Institute for doing astronomy calculations like image processing and observatory calculations. Because this is a Python program, you can install it with either pip or easy_install. Your Linux distribution already should have a package included. For example, in Debian-based distributions, you can install it with this:


sudo apt-get install python-astropy

There is also a separate package, python-astropy-doc, that contains extra documentation for Astropy. Because Astropy is a fairly large system, it is broken into smaller sub-packages. This should be familiar to anyone who has worked with packages like SciPy or NumPy before. So, simply using the following actually isn't very useful:


import astropy

You more likely will need to import individual sub-packages with commands like this:


from astropy.io import fits

There are sub-packages to handle file IO, cosmological calculations and coordinate systems, among many other topics. Let me provide a basic introduction to some of the available functionality so you have an idea of what you can do.

First, let's look at how to deal with data files. The common file format used in astrophysics and astronomy is the FITS file format. The PyFITS Python package was written to read and write FITS files. This code is actually the same as the code in the sub-package astropy.io.fits, so you can use it in the same way. You actually can even just drop Astropy in as a plugin with the following:


from astropy.io import fits as pyfits

This way, you can use existing file management code without having to make any changes.

The first thing you need to do is open your data file with:


from astropy.io import fits
hdulist = fits.open("My_File.fit")

This returns an object that behaves like a list. Each element of the returned object maps to a Header-Data Unit (HDU) in the file. You can get more information on the file with this command:


hdulist.info()

Each of the individual elements has a header and data portion. You can access them to see details about the data you are about to process.

Along with all of the library functions, Astropy includes a series of command-line utilities to work with FITS files. You can check the headers of a FITS file with the fitsheader utility. You can check your FITS file with fitscheck, and you even can find the differences between two files with fitsdiff.

A common computational process in astronomy is image processing. The convolution sub-package provides two categories of convolution operations: direct and FFT. You can do one-, two- and three-dimensional convolutions. The visualization sub-package handles more basic image processing like normalization and stretching. You can combine multiple transformations very easily. The + operator is overloaded to apply transformations that are "added" together in series. So, a command like this:


transform = SqrtStretch() + PercentileInterval(90.)

gives you a new function, transform, that combines the two separate transformations in a single step. This sub-package also includes a script, fits2bitmap, that can do conversions between different file formats.

A second common computational task in astronomy is doing statistics based on observations, and Astropy provides a sub-package called stats. Although the scipy.stats sub-package provides a lot of functionality, some astronomy-focused functions are missing, so the astropy.stats sub-package fills in those missing functions.

Once you have your data loaded, you can use the modeling sub-package. You can do 1D and 2D modeling with astropy.modeling. This includes curve-fitting functionality, where you can do linear and nonlinear function fitting. There are built-in functions to fit Gaussian curves and polynomials. This fitting is handled with a least-squares method. With version 1.0, you can build compound models by combining existing models with arithmetic operators.

When you are ready to start doing calculations, you will need to use constants. Instead of remembering them or using them with potential typos, Astropy includes a complete list of all the standard scientific constants that you will need when doing numerical work. You can import the entire list with this:


from astropy import constants

If you need only a few of the constants, like maybe the speed of light, you can import them individually with this:


from astropy.constants import c

The really cool thing about these constants is that they are actually "Quantity" objects. This means you can do something like change the units being used with a command like the following:


c.to('km/s')

Because it is so prevalent, you can use CGS units with c.cgs.

There are also two sub-packages to handle coordinate systems. Astronomical coordinate systems are handled by the coordinates sub-package, and world coordinate systems are handled by the wcs sub-package. In the coordinates sub-package, the core object is the SkyCoord object. There are methods of this object to handle conversions between coordinate systems or distances from one point to the origin within a given coordinate system. The wcs sub-package allows for mapping data within a FITS file onto a "real-world" coordinate system in order to analyze them correctly. This includes functionality to deal with complications, like projections onto the sphere of the sky.

You even can do cosmological calculations with Astropy. The cosmology sub-package actually includes functionality to model the evolution of the entire cosmos based on a set of initial conditions that you set. Although you can set your own cosmology, several built-in cosmologies are available. These are based on the WMAP and Planck satellite data.

Most functionality is built off a core FLRW object. This object represents a homogeneous, isotropic cosmology defined by the Friedmann-Lemaitre-Robertson-Walker metric from General Relativity. However, this class can't be used directly. You need to subclass it first. There are several included in the cosmology sub-package, such as the FlatLambdaCDM object that includes dark energy. You can look at various values, like the Hubble constant, at various points during the evolution of the cosmology. You also can include contributions to the energy density from matter, dark energy and even photons and neutrinos.

Now that you've seen a bit of what you can do with astropy, if astronomical calculations are on your radar, there is much more available. Additionally, there is the concept of affiliated packages. These are packages that basically are built on top of the core functionality provided by Astropy. Although they aren't part of Astropy, they are being built up into a community-driven environment for doing astrophysics. It definitely will be worth your while to take a look at the extended world of available packages.

Load Disqus comments