In my humble opinion, astronomers could be making more use of the fantastic tools that are in the astropy. One of the obstacles to this, as identified at the Python in Astronomy workshop, was documentation and worked examples. There are some tutorials for Astropy already, as well as other introductory material for astronomers, but more never hurts.

I was using astropy to do a small task the other day and was really impressed at how the tools in astropy made this very easy. What I needed to do was take some existing FITS images of M31 and construct masks that could be used to select only pixels that were within specific ranges of M31-o-centric distance.

I wrote this as a function, because that’s how I’m used to doing things. I wouldn’t claim it is totally beautiful code, but it’s not too terrible. Here it is, hopefully commented well enough to make it understandable.

``````import numpy as np
import astropy.io.fits as fits
import astropy.units as u
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord, Angle, Distance

'''check that masks which are supposed to select only certain galactocentric
distances do this correctly
input:
output:
no return value, but masks may be written to file
'''
print 'first 2 arguments must have same length'
return

# get the limits of galactocentric distance

# set its NaN values to 0

# make an array with the x and y coordinates of each pixel
# now compute the corresponding sky coordinates for each pixel
# compute galactocentric distance using a specific set of parameters for M31
galacto_dist = correct_rgc(coords, glx_ctr=SkyCoord(10.724382145*u.deg,
41.332205842*u.deg),
glx_PA=Angle('38d'),
glx_dist=Distance(780, unit=u.kpc),
glx_incl=Angle('75.883d'))
# drop the unit on the distance
galacto_dist_kpc = galacto_dist.to(u.kpc).value

# construct a new mask with specified distances
# change only certain pixels to 1
galacto_dist_kpc < upper)]=1.0

# see if the new and old masks are equal
print '%s %s' % (mask, check)

# if the masks aren't perfectly equal, see how close they are
if not check:
print '%d matches out of %d pix (%f)' %
The only non-astropy stuff is `from gal_radii_pb import correct_rgc` — this refers to a function I wrote myself to compute deprojected galactocentric distances; you can find it here. Thanks to Jonathan Sick who fixed up the first version of it that I posted.