Table Of Contents

Previous topic


Next topic

Python Installation and Understanding Packages

This Page

More Examples

This page contains further examples of applications of Python scripts to Astronomy beyond those in the Quick tour of Python:

Reading a table and plotting with asciitable

The Fermi Gamma-ray satellite has a nice catalog of AGN available through HEASARC. The script below will read in the catalog data using the asciitable module, do some basic filtering with NumPy, and make a couple of plots with matplotlib

import asciitable   # Make external package available

# Read table.
# ==> dat[column_name] and dat[row_number] both valid <==
dat ='fermi_agn.dat')

redshift = dat['redshift']    # array of values from 'redshift' column
flux = dat['photon_flux']
gamma = dat['spectral_index']

# Select rows that have a measured redshift
with_z = (redshift != -999)

semilogx(flux, gamma, '.b', label='All')  # First plot!
semilogx(flux[with_z], gamma[with_z], 'or', label='With Z')
xlabel('Flux (photon/cm$^2$/s)')   # latex works
ylabel('Spectral index $\Gamma$')

# Select low- and high-z samples
lowz = with_z & (redshift < 0.8)
highz = with_z & (redshift >= 0.8)

bins = arange(1.2, 3.0, 0.1)    # values from 1.2 to 3.0 by 0.1
hist(gamma[lowz], bins, color='b', alpha=0.5, label='z < 0.8')
hist(gamma[highz], bins, color='r', alpha=0.5, label='z > 0.8')
xlabel('Spectral index $\Gamma$')
title('$\Gamma$ for low-z and high-z samples')
legend(loc='upper left')

asciitable.write(dat[with_z], 'fermi_agn_with_z.dat')
../_images/scatter.png ../_images/hist.png

Curve fitting with SciPy

SciPy provides curve_fit, a simple and useful implementation of the Levenburg-Marquardt non-linear minimization algorithm. This example shows a code to generate a fake dataset and then fit with a gaussian, returning the covariance matrix for parameter uncertainties.

from scipy.optimize import curve_fit

# Create a function
# ==> First encounter with *whitespace* in Python <==
def gaussian(x, a, b, c):
    val = a * exp(-(x - b)**2 / c**2)
    return val

# Generate fake data.
# Note: functions in random package, array arithmetic (exp)
n = 100
x = random.uniform(-10., 10., n)
y = exp(-(x - 3.)**2 / 4) * 10. + random.normal(0., 2., n)
e = random.uniform(0.1, 1., n)

# Fit
popt, pcov = curve_fit(gaussian, x, y, sigma=e)

# Print results
print "Scale =  %.3f +/- %.3f" % (popt[0], sqrt(pcov[0, 0]))
print "Offset = %.3f +/- %.3f" % (popt[1], sqrt(pcov[1, 1]))
print "Sigma =  %.3f +/- %.3f" % (popt[2], sqrt(pcov[2, 2]))

# Plot data
errorbar(x, y, yerr=e, linewidth=1, color='black', fmt=None)

# Plot model
xm = linspace(-10., 10., 100)  # 100 evenly spaced points
plot(xm, gaussian(xm, popt[0], popt[1], popt[2]))

# Save figure

The plotted fit result is as shown below:


Intermission: NumPy, Matplotlib, and SciPy

These three packages are the workhorses of scientific Python.

  • NumPy is the fundamental package for scientific computing in Python [NumPy Reference]
  • Matplotlib is one of many plotting packages. Started as a Matlab clone.
  • SciPy is a collection of mathematical algorithms and convenience functions [SciPy Reference]

Synthetic images

This example demonstrates how to create a synthetic image of a cluster, including convolution with a Gaussian filter and the addition of noise.

import pyfits
from scipy.ndimage import gaussian_filter

# Create empty image
nx, ny = 512, 512
image = zeros((ny, nx))

# Set number of stars
n = 10000

# Generate random positions
r = random.random(n) * nx
theta = random.uniform(0., 2. * pi, n)

# Generate random fluxes
f = random.random(n) ** 2

# Compute position
x = nx / 2 + r * cos(theta)
y = ny / 2 + r * sin(theta)

# Add stars to image
# ==> First for loop and if statement <==
for i in range(n):
    if x[i] >= 0 and x[i] < nx and y[i] >= 0 and y[i] < ny:
        image[y[i], x[i]] += f[i]

# Convolve with a gaussian
image = gaussian_filter(image, 1)

# Add noise
image += random.normal(3., 0.01, image.shape)

# Write out to FITS image
pyfits.writeto('cluster.fits', image, clobber=True)

The simulated cluster image is below:


Running existing compiled codes

In addition to just doing computations and plotting, Python is great for gluing together other codes and doing system type tasks.

import os
import asciitable

smoothing = 30  # Smoothing window length
freqs = [2, 4]  # Frequency values for making data
noises = [1, 5] # Noise amplitude inputs


# Loop over freq and noise values, running standalone code to create noisy data
# and smooth it.  Get the data back into Python and plot.
plot_num = 1
for freq in freqs:
    for noise in noises:
        # Run the compiled code "make_data" to make data as a list of x, y, y_smooth
        cmd = 'make_data %s %s %s' % (freq, noise, smoothing)
        print 'Running', cmd
        out = os.popen(cmd).read()
        # out now contains the output from <cmd> as a single string

        # Write the output to a file
        filename = 'data_%s_%s' % (freq, noise)
        open(filename, 'w').write(out)

        # Parse the output string as a table
        dat =

        # Make a plot
        subplot(2, 2, plot_num)
        plot(dat['x'], dat['y'])
        plot(dat['x'], dat['y_smooth'], linewidth=3, color='r')

        plot_num += 1

And much much more...