Python and Astronomy

An overview of Python libraries and resources for (observational) astronomy. This presentation and associated materials are available at https://github.com/AnthonyHorton/python-and-astronomy/


Python

  • General purpose, interactive, object-orientated high level programming language with a large standard library
  • Free and open source, available on many platforms
  • Large & active user community, many 3rd party libraries (Python Package Index, PyPI, hosts 79753 packages)
  • 'Default' programming language for astronomy (though C/C++, FORTRAN, etc., remain important)

In [2]:
import this

In [3]:
import antigravity

Python 2 & 3?

  • Two incompatible branches of Python currently supported, 2.x & 3.x
  • Which to use? From the Python documentation:

Short version: Python 2.x is legacy, Python 3.x is the present and future of the language

Python 3.0 was released in 2008. The final 2.x version 2.7 release came out in mid-2010, with a statement of extended support for this end-of-life release. The 2.x branch will see no new major releases after that. 3.x is under active development and has already seen over five years of stable releases, including version 3.3 in 2012, 3.4 in 2014, and 3.5 in 2015. This means that all recent standard library improvements, for example, are only available by default in Python 3.x.

  • Unless your project uses a significant amount of existing Python 2 legacy code, use Python 3

  • It is possible to write cross compatible code. Many Python 3 features backported to Python 2.7.x, some enabled by default, others can be enabled with an import:

In [4]:
from __future__ import division, print_function, absolute_import, unicode_literals
  • Python 3 code with from __future__ import at the beginning may also run in Python 2.
  • There are packages for guaranteed compatibility, e.g. six and Python-Future
  • Python-Future also has functions for automatic conversion of code from Python 2 to Python 3 or vice versa (futurize and pasteurize)

Scientific Python stack

Collection of libraries to support general scientific computing in Python

  • numpy - N-dimensional arrays, matrices and operations on them
  • Scipy - numerical algorithms including stats, optimization/fitting, image analysis, signal processing, etc.
  • Matplotlib - plotting library, includes both a MATLAB-like interactive interface (pyplot) and a Pythonic API
  • SymPy - symbolic maths library/Computer Algebra System (CAS)
  • pandas - data structures and data analysis
  • IPython - Enhanced interactive Python

An aside for MATLAB users, the numpy documentation includes a guide specifically for you that includes side by side syntax comparison tables.


In [5]:
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import rcParams

In [6]:
%matplotlib inline

In [7]:
rcParams['figure.figsize'] = 13, 8.1

IPython & Jupyter

IPython is an interactive Python environment with lots of useful features.

  • Enhanced shell with syntax highlighting, session history & persistence, autocomplete, inline help, etc.
  • Support for interactive plotting, GUI toolkits, etc.
  • Easy to use parallel computing library
  • Kernel for Jupyter

Jupyter is a web application that allows you to create 'notebooks'

  • Combine code, output, plots, rich text, $\LaTeX$ equations, embedded media in a single live document
  • Originally part of IPython but now a separate project supporting 40 languages
  • Notebooks can be embedded in webpages, exported in various document formats or as plain Python

Combination of IPython & Jupyter great for interactive calculations or data analysis, algorithm development

Astronomical libraries

There are lots! How do you decide what to use?

PyRAF

A Python interface to the venerable Image Reduction and Analysis Facility, IRAF. Discussed by Kathleen Labrie on Wednesday.

Astropy

"The Astropy Project is a community effort to develop a single core package for Astronomy in Python and foster interoperability between Python astronomy packages."

The Astropy Project encompasses the astropy core package, affiliated packages and the Astropy community.

Core package

  • Provide core functionality/common tools required across astronomy, but don’t try to do everything
  • Avoid duplication for these core/common tasks
  • Provide robust framework for building more specialised/complex tools

Affiliated packages

  • Astronomy related Python packages outside of the core package but part of Astropy Project community
  • Often (but not always) intended for eventual inclusion in core package
  • Required to be at least working towards Astropy standards (documentation, testing, etc.)

Community

Why use it?

Many other Python libraries provide astronomy functions, but...

  • Astropy is well documented, tested and actively maintained
  • Astropy is a single library providing most core astronomy functions, & has very few external dependencies (numpy, others optional) => fewer external dependencies for your code, easier to port/share/maintain
  • Astropy modules designed to work together => easier interoperability than multiple separate libraries
  • STScI software (inc. JWST analysis tools) based on Astropy & incorporated into it

astropy core package

Core data structures and transformations Module Status (as of v1.1.2)
Constants astropy.constants Reasonably stable
Units and Quantities astropy.units Reasonably stable
N-dimensional datasets astropy.nddata Actively developed
Data Tables astropy.table Reasonably stable
Time and Dates astropy.time Mature
Astronomical Coordinate Systems astropy.coordinates Reasonably stable
World Coordinate System astropy.wcs Reasonably stable
Models and Fitting astropy.modeling Actively developed
Analytic Functions astropy.analytic_functions Actively developed
Connecting up: Files and I/O Module Status (as of v1.1.2)
Unified file read/write interface
FITS File handling astropy.io.fits Mature
ASCII Tables astropy.io.ascii Mature
VOTable XML handling astropy.io.votable Mature
Miscellaneous Input/Output astropy.io.misc Mature
Astronomy computations and utilities Module Status (as of v1.1.2)
Convolution and filtering astropy.convolution Reasonably stable
Data Visualization astropy.visualization Actively developed
Cosmological Calculations astropy.cosmology Reasonably stable
Astrostatistics Tools astropy.stats Actively developed
Virtual Observatory Access astropy.vo Reasonably stable

Astropy examples

Units & constants

The Astropy units module provides the very useful Quantity class, a drop in replacement for numpy Array that includes units, checking and conversions. It's not limited to simply multiplicative conversions, through 'equivalencies' it can do more conversions including parallax, angles, spectral units, doppler shifts, spectral flux density and brightness temperature. Also included are logarithmic units, including astronomical magnitudes.

Using these features makes all sorts of code simpler to write and, more importantly, less error prone.


In [8]:
import astropy.units as u
import astropy.constants as c

In [10]:
658.28 * u.nm + 10 * u.Angstrom


Out[10]:
$659.28 \; \mathrm{nm}$

In [11]:
658.29 * u.nm + 10


---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
<ipython-input-11-bbe3a116b0ed> in <module>()
----> 1 658.29 * u.nm + 10

/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/astropy/units/quantity.py in __array_prepare__(self, obj, context)
    344                                      "argument is not a quantity (unless the "
    345                                      "latter is all zero/infinity/nan)"
--> 346                                      .format(function.__name__))
    347             except TypeError:
    348                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

A more realistic example, converting sky background surface brightness in magnitudes per square arcsecond to electrons s$^{-1}$ pixel$^{-1}$ as might be required as part of a SNR calculation.

First set up some parameters, assume imaging in $g'$ and $r'$ bands with a 3.9m telescope equipped with an efficient 1''/pixel camera and dark sky background of $g'=22.5$, $r'=21.5$ AB mag arsecond$^{-2}$


In [12]:
sky_brightness_gr = (22.5, 21.5) * u.ABmag
pivot_wavelength_gr = (4702, 6175) * u.Angstrom
band_width_gr = (138.7, 137.3) * u.nm
pixel_scale = 1.0 * u.arcsecond
sky_solid_angle = pixel_scale**2 / u.pixel
aperture_area = np.pi * (3.9**2 - 1.5**2) / 4 * u.m**2
throughput = 0.7
QE = 0.9 * u.electron / u.photon

Convert astronomical magnitudes to physical spectral flux density units


In [13]:
sky_sfd_gr = sky_brightness_gr.to(u.Watt / (u.m**2 * u.micron), \
                                equivalencies=u.equivalencies.spectral_density(pivot_wavelength_gr))
sky_sfd_gr


Out[13]:
$[4.9232904 \times 10^{-17},~7.1704561 \times 10^{-17}] \; \mathrm{\frac{W}{\mu m\,m^{2}}}$

Multiply by filter band width, telescope aperture area, throughput and solid angle to get energy flux per pixel.


In [14]:
focal_plane_flux = sky_sfd_gr * band_width_gr * aperture_area * throughput * sky_solid_angle / u.arcsecond**2
focal_plane_flux.to(u.Watt / u.pixel)


Out[14]:
$[4.8654704 \times 10^{-17},~7.0147184 \times 10^{-17}] \; \mathrm{\frac{W}{pix}}$

Characteristic photon energy for the filter bands


In [15]:
photon_energy_gr = c.h * c.c / (pivot_wavelength_gr * u.photon)
photon_energy_gr.to(u.J / u.photon)


Out[15]:
$[4.2246824 \times 10^{-19},~3.2169161 \times 10^{-19}] \; \mathrm{\frac{J}{ph}}$

Divide energy flux per pixel by photon energy to get photon flux and multiply by QE to get electron rate


In [16]:
electron_rate = focal_plane_flux / photon_energy_gr * QE
electron_rate.to(u.electron / (u.second * u.pixel))


Out[16]:
$[103.65095,~196.25151] \; \mathrm{\frac{e^{-}}{pix\,s}}$

Tables

The Astropy table module provides functions for import/export and manipulations of all sorts of tabular data. Here is a simple format conversion.


In [17]:
from astropy.table import Table

An example CSV file:


In [18]:
with open('resources/nonsense.csv') as f:
    print(f.read())


ID,RA,dec,g_mag
1,04 23 15.3,-22 36 19,26.5
2,22 17 34.7,+02 22 59,23.4
3,10 44 29.1,-48 29 51,24.0
4,09 14 23.1,-10 10 41,27.3
5,15 51 11.6,+13 43 17,25.1

Load the CSV file (format automatically detected) and display the table.


In [19]:
nonsense = Table.read('resources/nonsense.csv')
nonsense


Out[19]:
<Table length=5>
IDRAdecg_mag
int64str10str9float64
104 23 15.3-22 36 1926.5
222 17 34.7+02 22 5923.4
310 44 29.1-48 29 5124.0
409 14 23.1-10 10 4127.3
515 51 11.6+13 43 1725.1

Write the table to file in $\LaTeX$ format. There are lots of formatting options and AASTeX deluxetables are also supported.


In [20]:
nonsense.write('resources/nonsense.tex', format='ascii.latex')

Look at the output $\LaTeX$ file.


In [21]:
with open('resources/nonsense.tex') as f:
    print(f.read())


\begin{table}
\begin{tabular}{cccc}
ID & RA & dec & g_mag \\
1 & 04 23 15.3 & -22 36 19 & 26.5 \\
2 & 22 17 34.7 & +02 22 59 & 23.4 \\
3 & 10 44 29.1 & -48 29 51 & 24.0 \\
4 & 09 14 23.1 & -10 10 41 & 27.3 \\
5 & 15 51 11.6 & +13 43 17 & 25.1 \\
\end{tabular}
\end{table}

Astropy Affiliated Packages

The Astropy affiliated package list is the place to look for more specific tools that are not in the Astropy core library. To be registered packages must comply with Astropy standards:

  • Use astropy classes/functions where possible
  • Have proper documentation
  • Include a test suite
  • Connect with Astropy developer community
  • Follow Astropy coding guidelines
  • Python 3 compatibility

Registered affiliated packages

Package Stable Maintainer
Astro-SCRAPPY Yes Curtis McCully
astroplan No Brett Morris
astroquery Yes Adam Ginsburg
ccdproc Yes Steven Crawford and Matt Craig
gwcs No Nadia Dencheva
Halotools Yes Andrew Hearin
montage-wrapper Yes Thomas Robitaille
naima Yes Victor Zabalza
photutils No Larry Bradley and Brigitta Sipocz
pydl No Benjamin Alan Weaver
sncosmo Yes Kyle Barbary
spectral-cube Yes Adam Ginsburg
specutils No Wolfgang Kerzendorf
WCSAxes Yes Thomas Robitaille

Provisionally accepted packages

These don't yet meet all Astropy standards but must be working towards meeting them.

Package Stable Maintainer
APLpy Yes Thomas Robitaille and Eli Bressert
astroML Yes Jake Vanderplas
galpy Yes Jo Bovy
gammapy No Christoph Deil
ginga Yes Eric Jeschke
Glue Yes Chris Beaumont and Thomas Robitaille
imexam No Megan Sosey
MaLTPyNT Yes Matteo Bachetti
omnifit Yes Aleksi Suutarinen
pyregion Yes Jae-Joon Lee
python-cpl No Ole Streicher
PyVO No Ray Plante
reproject Yes Thomas Robitaille
spherical_geometry No Bernie Simon and Michael Droettboom

Affiliated package example

Will demonstrate some features of a few of the affliated packages by doing a simple imaging data reduction. The data are a set of images of the Orion Nebula obtained with the Huntsman Eye (prototype for the Huntsman Telephoto Array, @AstroHuntsman, Facebook) in morning twilight of 24/09/2014.


In [22]:
import os, glob, subprocess

In [23]:
import ccdproc
from astroML.plotting import hist
from matplotlib.colors import LogNorm
from astropy.wcs import WCS

Set up the data directories


In [24]:
raw_path = os.path.abspath('resources/raw/')
calib_path = os.path.abspath('resources/calibrations/')
temp_path = os.path.abspath('resources/temp/')
red_path = os.path.abspath('resources/reduced/')

Create directories for output if they don't already exist


In [25]:
os.makedirs(temp_path, mode=0o775, exist_ok=True)
os.makedirs(red_path, mode=0o775, exist_ok=True)

Change into the working directory


In [26]:
os.chdir(temp_path)

List the (pre-prepared) calibration files


In [27]:
os.listdir(calib_path)


Out[27]:
['master_unbiased_dark.fits', 'master_flat.fits', 'master_bias.fits']

Load the calibration files as CCDData objects


In [28]:
bias = ccdproc.CCDData.read(os.path.join(calib_path, 'master_bias.fits'))
dark = ccdproc.CCDData.read(os.path.join(calib_path, 'master_unbiased_dark.fits'))
flat = ccdproc.CCDData.read(os.path.join(calib_path, 'master_flat.fits'))

Simple function to create a bad pixel mask (NaNs, Infs, outside range)


In [29]:
def isbadpixel(data, min_value, max_value):
    isgoodpixel = np.logical_and(np.isfinite(data), \
                                np.logical_and(data >= min_value, data <= max_value))
    return np.logical_not(isgoodpixel)

Mask and zero bad pixels in the calibration data (there are a few)


In [30]:
bias.mask = np.where(isbadpixel(bias.data, 1,60000), True, False)
bias.data = np.where(isbadpixel(bias.data, 1,60000), 0.0, bias.data)
dark.mask = np.where(isbadpixel(dark.data, -1000, 60000), True, False)
dark.data = np.where(isbadpixel(dark.data, -1000, 60000), 0.0, dark.data)
flat.mask = np.where(isbadpixel(flat.data, 0.7, 1.4), True, False)
flat.data = np.where(isbadpixel(flat.data, 0.7, 1.4), 0.0, flat.data)


/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in greater_equal
  from ipykernel import kernelapp as app
/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ipykernel/__main__.py:2: RuntimeWarning: invalid value encountered in less_equal
  from ipykernel import kernelapp as app

Inspect the calibration data


In [31]:
plt.imshow(bias, interpolation='none', aspect='equal', cmap='cubehelix', vmin=1025, vmax=1050)
plt.colorbar()


Out[31]:
<matplotlib.colorbar.Colorbar at 0x7f64a3b8db70>

In [32]:
plt.imshow(dark, interpolation='none', aspect='equal', cmap='cubehelix', vmax=40, vmin=0)
plt.colorbar()


Out[32]:
<matplotlib.colorbar.Colorbar at 0x7f64a02ab0f0>

In [33]:
plt.imshow(flat, interpolation='none', aspect='equal', cmap='cubehelix', vmax=1.4, vmin=0.7)
plt.colorbar()


Out[33]:
<matplotlib.colorbar.Colorbar at 0x7f64a01e87f0>

List the raw data files


In [34]:
os.listdir(raw_path)


Out[34]:
['83F011167_127_light.fits',
 '83F011167_126_light.fits',
 '83F011167_129_light.fits']

Open the raw data as a list of CCDData objects


In [35]:
raw_filenames = glob.glob(os.path.join(raw_path, '*_light.fits'))

In [36]:
images = [ccdproc.CCDData.read(filename, unit='adu') for filename in raw_filenames]


WARNING: FITSFixedWarning: EPOCH = 'JNOW ' / epoch of coordinates 
a floating-point value was expected. [astropy.wcs.wcs]

Check the headers


In [37]:
images[0].header


Out[37]:
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                   16 / number of bits per data pixel                  
NAXIS   =                    2 / number of data axes                            
NAXIS1  =                 3352 / length of data axis 1                          
NAXIS2  =                 2532 / length of data axis 2                          
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 
BZERO   =                32768 / offset data range to that of unsigned short    
BSCALE  =                    1 / default scaling factor                         
EXPTIME =                 300. / exposure time (seconds)                        
TEMPERAT=                 -20. / temperature (C)                                
IMAGETYP= 'light   '           / image type                                     
FILTNUM =                    0                                                  
DATE    = '2014-09-23T18:55:31' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
SERIALNO= '83F011167'          / serial number                                  
TARGET  = 'M42     '           / target name                                    
RA      = 'Succeeded.'         / right ascension                                
DEC     = 'RA:     '           / declination                                    
EPOCH   = 'JNOW    '           / epoch of coordinates                           
OBJCTRA = 'Succeeded.'         / right ascension                                
OBJCTDEC= 'RA:     '           / declination                                    
ALTITUDE=                   5. / Altitude (deg)                                 
AZIMUTH =                  36. / Azimuth (deg)                                  

Check the exposure times


In [38]:
[print(image.header['EXPTIME']) for image in images]


300.0
10.0
60.0
Out[38]:
[None, None, None]

Inspect the images


In [39]:
plt.imshow(images[0], norm=LogNorm(vmin=1000,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


Out[39]:
<matplotlib.colorbar.Colorbar at 0x7f64a0142588>

In [40]:
plt.imshow(images[1], norm=LogNorm(vmin=250,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


Out[40]:
<matplotlib.colorbar.Colorbar at 0x7f649ff590b8>

In [41]:
plt.imshow(images[2], norm=LogNorm(vmin=10,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()


Out[41]:
<matplotlib.colorbar.Colorbar at 0x7f649fee5080>

Trim the dummy pixels from the images


In [42]:
images = [ccdproc.trim_image(image[30:,:]) for image in images]


WARNING: VerifyWarning: Keyword name 'trim_image' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

Mask bad pixels (mostly saturated ones)


In [43]:
for image in images:
    image.mask = np.where(isbadpixel(image.data, 0, 64000), True, False)

Subtract bias


In [44]:
images = [ccdproc.subtract_bias(image, bias) for image in images]


WARNING: VerifyWarning: Keyword name 'subtract_bias' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

Subtract dark


In [45]:
images = [ccdproc.subtract_dark(image, dark, exposure_time='EXPTIME', \
                                exposure_unit=u.second, scale=True) for image in images]


WARNING: VerifyWarning: Keyword name 'subtract_dark' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

Flat field


In [46]:
images = [ccdproc.flat_correct(image, flat) for image in images]


/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/astropy/nddata/mixins/ndarithmetic.py:99: RuntimeWarning: divide by zero encountered in true_divide
  data = operation(self.data * self_unit, operand.data * operand_unit)
WARNING: VerifyWarning: Keyword name 'flat_correct' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

Write results to FITS files


In [47]:
images[0].write('m42_300.fits', clobber=True)
images[1].write('m42_060.fits', clobber=True)
images[2].write('m42_010.fits', clobber=True)

Astrometric fitting using astrometry.net's solve-field.


In [46]:
subprocess.call(['solve-field', '--use-sextractor', '--no-plot', '--overwrite', \
                '--ra', '83.81', '--dec', '-5.389', '--radius', '1', \
                'm42_010.fits', 'm42_060.fits', 'm42_300.fits'])


Out[46]:
0

Extract WCS info from solve-field output


In [48]:
images[0].wcs = WCS('m42_300.wcs')
images[1].wcs = WCS('m42_060.wcs')
images[2].wcs = WCS('m42_010.wcs')


WARNING: FITSFixedWarning: The WCS transformation has more axes (2) than the image it is associated with (0) [astropy.wcs.wcs]

Reproject images onto a common WCS


In [49]:
images[1] = ccdproc.wcs_project(images[1], images[0].wcs)
images[2] = ccdproc.wcs_project(images[2], images[0].wcs)


/mnt/data/ajh/Documents/virtualenvs/python3.4.3/lib/python3.4/site-packages/ccdproc/core.py:829: RuntimeWarning: invalid value encountered in greater
  reprojected_mask = reprojected_mask > 1e-8
WARNING: VerifyWarning: Keyword name 'wcs_project' is greater than 8 characters or contains characters not allowed by the FITS standard; a HIERARCH card will be created. [astropy.io.fits.card]

Divide each image by its exposure time


In [50]:
images = [image.divide(image.header['EXPTIME']) for image in images]

Combine images with a weighted average to create a HDR image


In [51]:
combiner = ccdproc.Combiner(images)
weights = np.zeros((3, bias.data.shape[0], bias.data.shape[1]))
weights[0] = 1.0
weights[1] = 60/300.0
weights[2] = 10/300.0
combiner.weights = weights
m42 = combiner.average_combine()
m42.wcs = images[0].wcs

Inspect the result


In [52]:
plt.subplot(111, projection=m42.wcs)
plt.imshow(m42, norm=LogNorm(vmin=3,vmax=3000), cmap='cubehelix')
plt.colorbar()


Out[52]:
<matplotlib.colorbar.Colorbar at 0x7f64a017fc18>

A closer look with corrected image orientation


In [53]:
plt.subplot(111, projection=m42.wcs)
plt.imshow(m42, norm=LogNorm(vmin=3,vmax=3000), cmap='cubehelix')
plt.xlim(2400,1100)
plt.ylim(610,1640)
plt.colorbar()


Out[53]:
<matplotlib.colorbar.Colorbar at 0x7f64998140b8>

Community

All of this software is provided free for the community to use. Some of the authors will have been employed specifically to work on software but many are scientists taking time out from their personal research to work on software that benefits the community.

Cite the code you use

  • When publishing your results acknowledge or, preferably, cite the software you used.
  • This not only credits the software authors but is also good for open science/reproducibility
  • Software can be cited either through an associated paper or, increasingly, directly, e.g. software published in the Astrophysical Source Code Library (ASCL) is indexed by ADS and is citable, and code published on GitHub may be citable via a DOI from Zenodo

Contribute to community code

Publish your own code

  • If you write a piece of code that you use more than once consider publishing it, others might find it useful too

Credit

  • There is a problem with scientists not getting the credit they deserve for time spent developing research software
  • This is beginning to change, though.
  • Example is Depsy, attempt to quantify the impact of individual scientist's software contributions.