Astropy Introduction

A brief history of the Astropy project

  • Started about 4 years ago
    • Goal to standardize on core set of astronomical libraries
      • And to foster use of standard libraries in Python application and library packages
  • STScI has migrated its standard libraries to astropy and is transitioning to use it for all its released software
  • Very active community
  • Much effort going into consistency and testing
    • Testing against existing coordinate conversion libraries, for example

Table of Contents

What Astropy is

  • Open source project intended from the beginning to be maintained by a community of astronomers, astrophysicists, software developers, and any other interested parties
    • Not tied to any single institution
  • An effort to lower the barrier of entry to new contributors with good ideas; lots of documentation for new (and existing) developers for how to contribute
  • Code, documentation, and tests that keep pace with the latest best practices of the Python community at large--strong emphasis on stylistic consistency and quality
  • Most importantly, provides easy (hopefully) to use interfaces for common procedures needed to write astronomy software
    • Units and constants
    • Coordinate conversions
    • Time conversions
    • Tables
    • Convolution and filtering
    • Modeling and fitting
    • Statistics
    • VO services
    • File formats (FITS, ASCII tables, VOTable, HDF5)
    • And other sundry utilities
  • Some functionality may be too special-purpose to go into the "core" package routines, but instead can be released as an Astropy Affiliated Package. Some affiliated packages include photutils (photometry), specutils (spectroscopy), gammapy (gamma-ray data analysis), ccdproc (CCD data reduction), etc...
    • Affiliated packages may in some case also be test-beds for features to be integrated into the core package.

Some things Astropy is not

  • An "IRAF replacement"
    • Astropy does not provide a complete data reduction and analysis environment--it is a library on which data analysis tools and other more specialized libraries can be built. By using the interfaces provided by Astropy, however, all tools built with it attain a high degree of interoperability (especially important in areas like coordinates and units!)
  • An "IDL replacement"
    • Not exactly a grab bag of single-purposes functions, and more of a framework--highly object-oriented. Not all IDL routines have an exact equivalent in Astropy, though most of them could be reproduced fairly easily
    • That said, efforts are being made to make it easier to find out how to do specific tasks, and people have also been compiling lists of documentation for how to recreate popular IDL routines and IRAF tasks using Astropy (and affiliated packages)
  • A small library
    • We barely have time to scratch the surface today of what capabilities are available through Astropy
  • Unwieldy
    • By using it in one's day to day coding tasks one can get quite comfortable and efficient at using the parts of Astropy that they need


One of the most powerful, central features of Astropy is the unit conversion library, and the Quantity object that associates numerical values (and arrays) with units. Effort is made to integrate units throughout Astropy so it is important to understand their usage.

Since we may want to use a number of units in expressions, it is easiest and most concise to import the units module with:

In [190]:
from astropy import units as u

though note that this will conflict with any variable called u.

Units can then be accessed with:

In [194]:

In the notebook this just displays a LaTeX representation of the unit:

In [195]:

They all have a docstring defining them:

In [196]:
meter: base unit of length in SI

An associated physical type:

In [197]:
In [198]:

And other various associated information that can be explored through the full documentation.

To create a quantity we just multiply a value by a unit. The value can be either a scalar or a Numpy array (or an object that can be converted to an array, like a Python list):

In [199]:
1.0 * u.m
$1 \; \mathrm{m}$
In [200]:
[1.0, 2.0, 3.0] * u.m
<Quantity [ 1., 2., 3.] m>

Since u.m and u.Angstrom have the same physical type we can convert a quantity in meters to a quantity in Angstroms:

In [201]:
q = [1.0, 2.0, 3.0] * u.m
<Quantity [  1.00000000e+10,  2.00000000e+10,  3.00000000e+10] Angstrom>

Quantities can then be combined:

In [202]:
q1 = 3. * u.m
In [204]:
$3 \; \mathrm{m}$
In [203]:
q2 = 5. * / u.s / u.g**2
In [205]:
$5 \; \mathrm{\frac{cm}{s\,g^{2}}}$
In [206]:
q1 * q2
$15 \; \mathrm{\frac{cm\,m}{s\,g^{2}}}$

and converted to different units:

In [207]:
(q1 * q2).to(u.m**2 /**2 / u.s)
$150000 \; \mathrm{\frac{m^{2}}{s\,kg^{2}}}$

We can also incorporate physical constants from the astropy.constants library. For example, ff we want to compute the Gravitational force felt by a 100. kg space probe by the Sun, at a distance of 3.2au, we can do:

In [208]:
from astropy.constants import G
In [210]:
print repr(G)
<Constant name=u'Gravitational constant' value=6.67384e-11 error=8e-15 units='m3 / (kg s2)' reference=u'CODATA 2010'>
In [211]:
F = (G * 1. * u.M_sun * 100. * / (3.2 ***2
In [215]:
  Primary name | Unit definition | Aliases       
  N            | kg m / s2       | Newton, newton ,
  dyn          | 1e-05 kg m / s2 | dyne           ,
In [213]:
$0.0579271 \; \mathrm{N}$

Units don't necessarily have to have the same physical type to be converted between each other. Astropy includes a number of unit equivalencies that can be enabled (with the possibility of defining your own). For example, without equivalencies the following fails:

In [216]:
(450. * u.nm).to(u.GHz)
UnitsError                                Traceback (most recent call last)
<ipython-input-216-e4fbcb033257> in <module>()
----> 1 (450. * u.nm).to(u.GHz)

/usr/local/lib/python2.7/dist-packages/astropy/units/quantity.pyc in to(self, unit, equivalencies)
    575         unit = Unit(unit)
    576         new_val = np.asarray(
--> 577   , self.value, equivalencies=equivalencies))
    578         return self._new_view(new_val, unit)

/usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in to(self, other, value, equivalencies)
    910             If units are inconsistent
    911         """
--> 912         return self.get_converter(other, equivalencies=equivalencies)(value)
    914     def in_units(self, other, value=1.0, equivalencies=[]):

/usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in get_converter(self, other, equivalencies)
    846         except UnitsError:
    847             return self._apply_equivalences(
--> 848                 self, other, self._normalize_equivalencies(equivalencies))
    849         return lambda val: scale * _condition_arg(val)

/usr/local/lib/python2.7/dist-packages/astropy/units/core.pyc in _apply_equivalences(self, unit, other, equivalencies)
    809         raise UnitsError(
    810             "{0} and {1} are not convertible".format(
--> 811                 unit_str, other_str))
    813     def get_converter(self, other, equivalencies=[]):

UnitsError: 'nm' (length) and 'GHz' (frequency) are not convertible

But with the spectral equivalency:

In [217]:
(450. * u.nm).to(u.GHz, equivalencies=u.spectral())
$666205 \; \mathrm{GHz}$
In [218]:
(450. * u.eV).to(u.nm, equivalencies=u.spectral())
$2.7552 \; \mathrm{nm}$

Most of the Numpy functions also understand Astropy quantities (with some notable exceptions you can find in the documentation). For example, by default np.sin takes input in radians:

In [219]:
import numpy as np

But if we pass in a quantity in degrees:

In [220]:
np.sin(90 *
$1 \; \mathrm{}$

More examples and exercises with units in this notebook.


The astropy Table class provides an extension of NumPy structured arrays for storing and manipulating heterogeneous tables of data. A few notable features of this package are:

  • Initialize a table from a wide variety of input data structures and types.
  • Modify a table by adding or removing columns, changing column names, or adding new rows of data.
  • Handle tables containing missing values.
  • Include table and column metadata as flexible data structures.
  • Specify a description, units and output formatting for columns.
  • Perform operations like database joins, concatenation, and grouping.
  • Manipulate multidimensional columns.
  • Methods for Reading and writing Table objects to files
  • Integration with Astropy Units and Quantities
  • And still more capabilities there won't be time to demonstrate today...
In [221]:
from astropy.table import Table
import numpy as np

Creating tables

There is great deal of flexibility in the way that a table can be initially constructed:

  • Read an existing table from a file or web URL
  • Add columns of data one by one
  • Add rows of data one by one
  • From an existing data structure in memory:

    • List of data columns
    • Dict of data columns
    • List of row dicts
    • Numpy homgeneous array or structured array
    • List of row records (new in 0.4)

See the documentation section on Constructing a table for the gory details and plenty of examples.

The easiest (though maybe not always most useful in practice) way is to just create an empty table object and add columns to it by assigning arrays to column names, in a syntax reminiscent on dictionaries:

In [222]:
t = Table()
t['name'] = ['source 1', 'source 2', 'source 3', 'source 4']
t['flux'] = [1.2, 2.2, 3.1, 4.3]

Tables have numerous useful attributes for introspection. Just a couple examples are:

In [223]:
['name', 'flux']
In [224]:
dtype([('name', '|S8'), ('flux', '<f8')])

There are several options for viewing tables, but for a small table it will just print in full:

In [225]:
source 11.2
source 22.2
source 33.1
source 44.3
In [226]:
<open file '<fdopen>', mode 'w+b' at 0x529b300>

We can access single rows of columns, single columns, or selection sof multiple rows/columns:

In [227]:
<Column name='flux' unit=None format=None description=None>
array([ 1.2,  2.2,  3.1,  4.3])
In [228]:
<Row 1 of table
 values=('source 2', 2.2)
 dtype=[('name', '|S8'), ('flux', '<f8')]>
In [229]:
source 22.2
source 33.1

Table columns can also have units associated with them for storing quantities to tables. This is useful, for example, when we read FITS tables--if a column has an associated TUNIT keyword we can return dimensionful quantities from the table.

Masked tables are another powerful feature of the Astropy Table class:

In [230]:
t2 = Table([['x', 'y', 'z'], 
            [1.1, 2.2, 3.3]],
           names=['name', 'value'],
In [231]:
t2['value'].mask = [False, True, False]
In [232]:

Numerical operations on masked columns ignore the masked values:

In [233]:
In [234]:

We can set our own fill value which can be used when outputting the table to a format that does not natively support masked values:

In [235]:
t2['value'].fill_value = -99
In [236]:

Tables support a number of high-level operations for combining tables, including merging, concatenating, and grouping. For example, recall our first table of flux values:

In [237]:
source 11.2
source 22.2
source 33.1
source 44.3

Now say that we now got some additional flux values from a different reference for a different, but overlapping sample of sources:

In [238]:
t2 = Table()
t2['name'] = ['source 1', 'source 3', 'source 8']
t2['flux2'] = [1.4, 3.5, 8.6]

Now we can get a master table of flux measurements which are joined matching the values the name column. This includes every row from each of the two tables, which is known as an outer join:

In [239]:
from astropy.table import join
t3 = join(t, t2, keys=['name'], join_type='outer')
source 11.21.4
source 22.2--
source 33.13.5
source 44.3--
source 8--8.6

There are several other join options. For example we could choose to keep only rows where both tables had a valid measurement using an inner join:

In [240]:
t3 = join(t, t2, keys=['name'], join_type='inner')
source 11.21.4
source 33.13.5

Tables can be written to and read in from a variety of data formats:

In [241]:
t3.write('example.fits', overwrite=True)

Tables can be read back in using the interface, and then written to another format (this is an easy way to do format conversion):

In [242]:
t4 ='example.fits')
In [243]:
source 11.21.4
source 33.13.5
In [244]:
t4.write('example.vot', format='votable', overwrite=True)
In [245]:
!cat example.vot
<?xml version="1.0" encoding="utf-8"?>
<!-- Produced with version 0.4.2 -->
<VOTABLE version="1.2" xmlns="" xmlns:xsi="" xsi:noNamespaceSchemaLocation="">
 <RESOURCE type="results">
   <FIELD ID="name" arraysize="8" datatype="char" name="name"/>
   <FIELD ID="flux" datatype="double" name="flux"/>
   <FIELD ID="flux2" datatype="double" name="flux2"/>
      <TD>source 1</TD>
      <TD>source 3</TD>

Tables can also be written to a variety of ASCII formats, including custom formats. Before trying to hack together a hand-written table reader/writer look into what is possible with the

There is also an API for developing new table readers/writers for other formats, though the built-in formats cover much of the existing data out there.

In [274]:
t = Table()
t['name'] = ['source 1', 'source 2', 'source 3', 'source 4']
t['flux'] = [1.2, 2.2, 3.1, 4.3]
In [275]:
source 11.2
source 22.2
source 33.1
source 44.3
In [276]:
t.write('mytable.tex', format='latex')
In [278]:
t.write(sys.stdout, format='latex')
name & flux \\
source 1 & 1.2 \\
source 2 & 2.2 \\
source 3 & 3.1 \\
source 4 & 4.3 \\
In [279]:
text = """
name & flux \\
source 1 & 1.2 \\
source 2 & 2.2 \\
source 3 & 3.1 \\
source 4 & 4.3 \\
In [281]:'mytable.tex', format='latex')
source 11.2
source 22.2
source 33.1
source 44.3
In [277]:
import sys
<IPython.kernel.zmq.iostream.OutStream at 0x2cc3e10>

Doing more with FITS files

We briefly saw how to save and load binary tables from FITS files, but that was through a high-level interface that does not give us a great deal of access to the lower-level details of FITS files. We would also in many cases like to be able to manipulate the contents of FITS headers and the arrangement of HDUs.

And although there is a high-level generic interface for tabular data in the Table class, the equivalent for gridded data like images is still under development (see: NDData).

The full interface to FITS images that allows direct access to their image arrays and headers is via the package. If you have ever used the venerable PyFITS library, then you have already used, which is just an import of PyFITS into Astropy. By moving PyFITS into Astropy it is able to evolve new capabilities, such as the improved table interface, and improved support for units.

To open a FITS file import the fits package and use the function:

In [246]:
from import fits
In [247]:
# The following FITS file lives in the lesson repository, relative
# to this file.  However, if you get a "file not found" error you can download and open the file
# by uncommenting these lines (and commenting the last line)

# from astropy.utils import download_file
# filename = download_file('')
# hdulist =
hdulist ='data/')

The returned object, hdulist, (an instance of the HDUList class) behaves like a Python list, and each element maps to a Header-Data Unit (HDU) in the FITS file. You can view more information about the FITS file with:

In [248]:
Filename: data/
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      34   (720, 360, 30)   float32   
1    ENERGIES    BinTableHDU     19   30R x 1C     [D]   

As we can see, this file contains two HDUs. The first contains an image cube, the second a table with one column of double-precision floating point values. To access the primary HDU, which contains the main data, you can then do:

In [249]:
primary_hdu = hdulist[0]
In [250]:
primary_hdu = hdulist['PRIMARY']

To access the header use the .header attribute:

In [251]:
SIMPLE  =                    T / Written by IDL:  Thu Jan 20 07:19:05 2011      
BITPIX  =                  -32 /                                                
NAXIS   =                    3 / number of data axes                            
NAXIS1  =                  720 / length of data axis 1                          
NAXIS2  =                  360 / length of data axis 2                          
NAXIS3  =                   30 / length of data axis 3                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
FLUX    =        8.42259635886 /                                                
CRVAL1  =                   0. / Value of longitude in pixel CRPIX1             
CDELT1  =                  0.5 / Step size in longitude                         
CRPIX1  =                360.5 / Pixel that has value CRVAL1                    
CTYPE1  = 'GLON-CAR'           / The type of parameter 1 (Galactic longitude in 
CUNIT1  = 'deg     '           / The unit of parameter 1                        
CRVAL2  =                   0. / Value of latitude in pixel CRPIX2              
CDELT2  =                  0.5 / Step size in latitude                          
CRPIX2  =                180.5 / Pixel that has value CRVAL2                    
CTYPE2  = 'GLAT-CAR'           / The type of parameter 2 (Galactic latitude in C
CUNIT2  = 'deg     '           / The unit of parameter 2                        
CRVAL3  =                  50. / Energy of pixel CRPIX3                         
CDELT3  =    0.113828620540137 / log10 of step size in energy (if it is logarith
CRPIX3  =                   1. / Pixel that has value CRVAL3                    
CTYPE3  = 'photon energy'      / Axis 3 is the spectra                          
CUNIT3  = 'MeV     '           / The unit of axis 3                             
CHECKSUM= '3fdO3caL3caL3caL'   / HDU checksum updated 2009-07-07T22:31:18       
DATASUM = '2184619035'         / data unit checksum updated 2009-07-07T22:31:18 
DATE    = '2009-07-07'         /                                                
FILENAME= '$TEMPDIR/diffuse/' /File name with version number     
TELESCOP= 'GLAST   '           /                                                
INSTRUME= 'LAT     '           /                                                
ORIGIN  = 'LISOC   '           /LAT team product delivered from the LISOC       
OBSERVER= 'MICHELSON'          /Instrument PI                                   
HISTORY Scaled version of for use with P6_V11_DIFFUSE           

The Header object works very similarly to the Python dict type--the values for individual keywords can be looked up like so:

In [252]:
In [253]:

The Header interface has a wealth of other capabilities that can be read about in the full documentation.

The array data associated with an HDU can be accessed via the .data attribute, which returns a view of the data as a Numpy array. In this example the data is an image cube for a series of wavelengths. To display the first one:

In [254]:
%matplotlib inline
from matplotlib.pyplot import imshow

imshow([0], origin='lower')
<matplotlib.image.AxesImage at 0x66887d0>
In [255]:
from aplpy import FITSFigure

fig = FITSFigure('data/', slices=[0])
INFO:astropy:Auto-setting vmin to -1.114e-06
INFO:astropy:Auto-setting vmax to  1.337e-05
INFO: Auto-setting vmin to -1.114e-06 [aplpy.core]
INFO: Auto-setting vmax to  1.337e-05 [aplpy.core]

We will use this for a brief demonstration of how to write a new FITS file. Say we wanted to save just this slice to a file. We can create a new primary HDU with this data like so:

In [256]:
new_hdu = fits.PrimaryHDU([0], header=primary_hdu.header)

This copies over the header from the original file, but some of the keywords are no longer applicable, such as the WCS for the spectral axis, or the filename keyword:

In [257]:
del new_hdu.header['FILENAME']

We can use a wildcard expression to take care of the WCS:

In [258]:
del new_hdu.header['C*3']

Note, however, that we did not have to fuss around with keywords related to the basic structure of the file, such as BITPIX or the NAXIS keywords. It is not necessary to manipulate such basic structural keywords by hand--their consistency with the data is managed automatically (in the future there will be further enhancements to make WCS manipulation more transparent too):

In [259]:
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                  -32 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                  720                                                  
NAXIS2  =                  360                                                  
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
FLUX    =        8.42259635886 /                                                
CRVAL1  =                   0. / Value of longitude in pixel CRPIX1             
CDELT1  =                  0.5 / Step size in longitude                         
CRPIX1  =                360.5 / Pixel that has value CRVAL1                    
CTYPE1  = 'GLON-CAR'           / The type of parameter 1 (Galactic longitude in 
CUNIT1  = 'deg     '           / The unit of parameter 1                        
CRVAL2  =                   0. / Value of latitude in pixel CRPIX2              
CDELT2  =                  0.5 / Step size in latitude                          
CRPIX2  =                180.5 / Pixel that has value CRVAL2                    
CTYPE2  = 'GLAT-CAR'           / The type of parameter 2 (Galactic latitude in C
CUNIT2  = 'deg     '           / The unit of parameter 2                        
CHECKSUM= '3fdO3caL3caL3caL'   / HDU checksum updated 2009-07-07T22:31:18       
DATASUM = '2184619035'         / data unit checksum updated 2009-07-07T22:31:18 
DATE    = '2009-07-07'         /                                                
TELESCOP= 'GLAST   '           /                                                
INSTRUME= 'LAT     '           /                                                
ORIGIN  = 'LISOC   '           /LAT team product delivered from the LISOC       
OBSERVER= 'MICHELSON'          /Instrument PI                                   
HISTORY Scaled version of for use with P6_V11_DIFFUSE           

Writing the file out with the checksum=True option will also update the CHECKSUM and DATASUM keywords:

In [261]:
new_hdu.writeto('lat_background_model_slice.fits', checksum=True,
WARNING: Overwriting existing file 'lat_background_model_slice.fits'. []
WARNING:astropy:Overwriting existing file 'lat_background_model_slice.fits'.


This is a good place for a brief mention of the coordinates interface; specifically the SkyCoord object, a high level interface for representation and manipulation of celestial coordinates. Astropy 0.4 only includes a few common coordinate systems (ICRS, FK4, FK5, and Galactic), but future versions will include more built-in coordinate systems, and users can already define their own systems.

Creating coordinate objects is straightforward. We use units--often u.deg (a shorthand for to specify the units of the angles:

In [262]:
from astropy.coordinates import SkyCoord
from astropy import units as u
In [263]:
SkyCoord(ra=10.68458 * u.deg, dec=41.26917 * u.deg, frame='icrs')
<SkyCoord (ICRS): ra=10.68458 deg, dec=41.26917 deg>

In addition to passing angular Quantity objects, the SkyCoord class can parse coordinate strings (this can be used to pass initialize a SkyCoord object from, say, a column of coordinates read from an ASCII table with the Table interface):

In [264]:
c = SkyCoord('00h42m44.3s +41d16m9s', frame='icrs')
<SkyCoord (ICRS): ra=10.6845833333 deg, dec=41.2691666667 deg>

We can access the different components as Angle objects and convert them to other units:

In [265]:
In [266]:

...and so on.

It is also possible to convert our SkyCoord object to other coordinate frames:

In [267]:
<SkyCoord (Galactic): l=121.174305193 deg, b=-21.5728034839 deg>
In [268]:
<SkyCoord (Galactic): l=121.174305193 deg, b=-21.5728034839 deg>

Or we can create custom coordinate frames in order to adjust other parameters, such as the equinox:

In [269]:
from astropy.coordinates import FK5
In [271]:
In [272]:
frame = FK5(equinox='J2010')
In [273]:
<SkyCoord (FK5): equinox=J2010.000, ra=10.8218390257 deg, dec=41.3238612021 deg>