Python Tutorials and Demos

Plan: Weekly 15-30 minute python tutorials, code demonstrations, and discussions Mondays at 11 AM in the library.

Organizers: Adam Ginsburg and Bernd Husemann

If you would like to give a code demonstration,
A google calendar is available
Resources from this group will be made available via a github repository
Have a request? Post it as an issue

Resources

If you want to get this notebook, you can find it on a public repository: https://github.com/ESO-python/ESOPythonTutorials

It can be viewed directly at http://nbviewer.ipython.org/github/ESO-python/ESOPythonTutorials/master/notebooks/ESOPythonDemoDay1.ipynb (or http://goo.gl/UNnpLF for short)

Astropy has very good and detailed documentation: http://www.astropy.org/ and a tutorial series (work in progress): http://www.astropy.org/astropy-tutorials/

Astropy

Astropy is a community-led effort to provide basic astronomy-related tools and libraries in python. It is meant to provide all of the functionality astrolib once provided for IDL, but much more.

Community-developed means you can (and should) contribute. All you need is a free github account.

The core functionality of astropy includes: * Units and Quantities (astropy.units) and Constants (astropy.constants) * Data Tables (astropy.table) * Astronomical Coordinate Systems (astropy.coordinates) * World Coordinate System (astropy.wcs) * FITS File handling (astropy.io.fits)

IPython and the IPython Notebook

This demonstration is written in an IPython notebook.

ipython is the best way to work with python for interactive data analysis. Notebooks are particularly useful for lectures, demonstrations, and day-to-day exploratory work.

Units

Astropy has a complete set of unit manipulation tools.

In [1]:

from astropy import units as u
from astropy import constants

Any quantity can be assigned a unit by multiplying by the appropriate unit

In [2]:

x = 6563*u.AA
x

Out[2]:

$6563 \; \mathrm{\AA}$

Units can be converted to any equivalent unit (i.e., length->length, speed->speed, etc.)

In [3]:

x.to(u.m)

Out[3]:

$6.563\times 10^{-07} \; \mathrm{m}$

In [4]:

x.to(u.pc)

Out[4]:

$2.12692\times 10^{-23} \; \mathrm{pc}$

They cannot be converted to non-equivalent units...

...unless the appropriate equivalencies are specified

In [5]:

x.to(u.THz, u.spectral())

Out[5]:

$456.792 \; \mathrm{THz}$

More complicated equivalencies, such as doppler shift (redshift) are possible. They are tied to particular reference values, though:

In [6]:

x = 6582*u.AA
ha_rest = 6562.8*u.AA
x.to(u.km/u.s, u.doppler_optical(ha_rest))

Out[6]:

$877.067 \; \mathrm{\frac{km}{s}}$

Arrays can also be given units

In [7]:

frequencies = np.linspace(1,2,50)*u.GHz
fluxes = np.random.randn(50)*u.Jy

In [8]:

import pylab as pl
pl.plot(frequencies, fluxes, drawstyle='steps-mid')
pl.xlabel("{0} ({1})".format(frequencies.unit.physical_type, frequencies.unit.to_string()))
pl.ylabel("{0} ({1})".format(fluxes.unit.physical_type, fluxes.unit.to_string()))

Out[8]:

<matplotlib.text.Text at 0x10879f990>

Coordinates

Conversion between fk5 and Galactic coordinates (and any other system) is straightforward with astropy.

In [9]:

from astropy import coordinates

In [10]:

# Query the SESAME service to get the coordinates
m31_fk5 = coordinates.FK5.from_name('M31')

In [11]:

m31_fk5

Out[11]:

<FK5 RA=10.68472 deg, Dec=41.26875 deg>

In [12]:

m31_fk5.galactic

Out[12]:

<Galactic l=121.17439 deg, b=-21.57322 deg>

In [13]:

m31_fk5.fk4

Out[13]:

<FK4 RA=10.00037 deg, Dec=40.99498 deg>

Coordinates also work with arrays

In [14]:

ra = [17.5,17.6,17.7,17.8] * u.hour
dec = [-29.2,-29.3,-29.4,-29.5] * u.deg
c = coordinates.FK5(ra,dec)

In [15]:

c, c.galactic

Out[15]:

(<FK5 RA=[ 262.5  264.   265.5  267. ] deg, Dec=[-29.2 -29.3 -29.4 -29.5] deg>,
 <Galactic l=[ 357.94945306  358.57821343  359.19077086  359.78751897] deg, b=[ 2.74153754  1.58913794  0.42968231 -0.73650041] deg>)

In [16]:

lmc = coordinates.FK5.from_name('LMC')
smc = coordinates.FK5.from_name('SMC')
m33 = coordinates.FK5.from_name('M33')

In [17]:

fig = pl.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection="mollweide")
ax.scatter(c.galactic.l.radian-2*np.pi, c.galactic.b.radian, color='k')
ax.scatter(m31_fk5.galactic.l.radian, m31_fk5.galactic.b.radian, color='orange')
ax.scatter(m33.galactic.l.radian, m33.galactic.b.radian, color='r')
ax.scatter(lmc.galactic.l.radian-2*np.pi, lmc.galactic.b.radian, color='b')
ax.scatter(smc.galactic.l.radian-2*np.pi, smc.galactic.b.radian, color='m')

Out[17]:

<matplotlib.collections.PathCollection at 0x108bbd310>

/Users/adam/virtual-python/lib/python2.7/site-packages/matplotlib-1.4.x-py2.7-macosx-10.6-intel.egg/matplotlib/projections/geo.py:489: RuntimeWarning: invalid value encountered in arcsin
  theta = np.arcsin(y / np.sqrt(2))