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/
In [2]:
import this
In [3]:
import antigravity
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.
In [4]:
from __future__ import division, print_function, absolute_import, unicode_literals
from __future__ import
at the beginning may also run in Python 2.futurize
and pasteurize
) Collection of libraries to support general scientific computing in Python
pyplot
) and a Pythonic APIAn 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 is an interactive Python environment with lots of useful features.
Jupyter is a web application that allows you to create 'notebooks'
Combination of IPython & Jupyter great for interactive calculations or data analysis, algorithm development
There are lots! How do you decide what to use?
A Python interface to the venerable Image Reduction and Analysis Facility, IRAF. Discussed by Kathleen Labrie on Wednesday.
"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.
Many other Python libraries provide astronomy functions, but...
astropy
core packageCore 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 |
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]:
In [11]:
658.29 * u.nm + 10
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]:
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]:
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]:
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]:
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())
Load the CSV file (format automatically detected) and display the table.
In [19]:
nonsense = Table.read('resources/nonsense.csv')
nonsense
Out[19]:
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())
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:
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 |
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 |
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]:
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)
Inspect the calibration data
In [31]:
plt.imshow(bias, interpolation='none', aspect='equal', cmap='cubehelix', vmin=1025, vmax=1050)
plt.colorbar()
Out[31]:
In [32]:
plt.imshow(dark, interpolation='none', aspect='equal', cmap='cubehelix', vmax=40, vmin=0)
plt.colorbar()
Out[32]:
In [33]:
plt.imshow(flat, interpolation='none', aspect='equal', cmap='cubehelix', vmax=1.4, vmin=0.7)
plt.colorbar()
Out[33]:
List the raw data files
In [34]:
os.listdir(raw_path)
Out[34]:
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]
Check the headers
In [37]:
images[0].header
Out[37]:
Check the exposure times
In [38]:
[print(image.header['EXPTIME']) for image in images]
Out[38]:
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]:
In [40]:
plt.imshow(images[1], norm=LogNorm(vmin=250,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()
Out[40]:
In [41]:
plt.imshow(images[2], norm=LogNorm(vmin=10,vmax=64000), interpolation='none', aspect='equal', cmap='cubehelix')
plt.colorbar()
Out[41]:
Trim the dummy pixels from the images
In [42]:
images = [ccdproc.trim_image(image[30:,:]) for image in images]
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]
Subtract dark
In [45]:
images = [ccdproc.subtract_dark(image, dark, exposure_time='EXPTIME', \
exposure_unit=u.second, scale=True) for image in images]
Flat field
In [46]:
images = [ccdproc.flat_correct(image, flat) for image in images]
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]:
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')
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)
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]:
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]:
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.