Hands-on Exercise 0: Introduction to Python & Astropy

by Leo Singer (2014) and Eric Bellm (2015-2016)

Introduction

Our hands-on exercises will use the Python programming language. No previous experience with Python is necessary!

Python is a powerful tool out of the box, but its strengh comes from the huge number of external packages. The following are vital for data analysis:

• Matplotlib: plotting interactive or publication-quality figures
• Numpy: vectorized arithmetic and linear algebra
• Scipy: curated collection of algorithms for root finding, interpolation, integration, signal processing, statistics, linear algebra, and much more
• IPython Notebook, the Mathematica-like interface that you are using now, and last but not least
• Astropy, a community library for astronomy. See also Robitaille et al. (2013).

We'll cover the basics of Python itself and then dive in to some applications to explore each of these packages.

We are using the Jupyter Notebook interface, which like Mathematica allows you to intermingle computer code with generated plots and output.

The notebook contains two main types of cells: Markdown cells for text, like this one...

Double-click in this text cell to edit it:

WRITE YOUR NAME HERE

and then press Shift-Enter to execute it.

How to get Python/Matplotlib/Numpy/Scipy

Python and all of the packages that we discuss in this tutorial are open source software, so there are multiple options for installing them.

You should have already installed Python before our tutorial, but this section is here for your future reference.

The most current version of Python is 3.5, but there is a large existing codebase of Python 2.7 code, so you may see both in the wild. Unfortunately there are small differences between the two versions. Here is one guide to the differences.

Recommended on all platforms: Anaconda, Canopy

Anaconda and Canopy are free, self-contained Python distributions for Linux, Mac, Windows. Both come with an excellent selection of up-to-date scientific Python libraries, including Astropy.

TIP These are both great Python environments for everyday use. If you need the latest version of Python in a headless environment like a cluster or web server, consider the lightweight Miniconda distribution.

You can install multiple versions of python alongside each other using conda environments: for example, to run python 2.7 as well as 3.5. Documentation is here; here is the key idea:

conda create --name summerschool python=3.5 anaconda

source activate summerschool

Advanced: Debian/Ubuntu

Most operating systems already come with Python. Release schedules vary, but generally the 'testing' versions of Debian and Ubuntu come with very up-to-date versions of Python. Debian Jessie and Ubuntu 13.10 (Saucy Salamander) both come with the latest version of Astropy. On these systems, you can set up a pretty complete Python environment with:

$ sudo apt-get install ipython-notebook python-matplotlib python-scipy python-pip python-astropy

Advanced: MacPorts

Every version of Mac OS comes with a Python interpreter, but it's slightly easier to obtain Matplotlib and Numpy if you use a package manager such as MacPorts, HomeBrew, or Fink. I use MacPorts (and contribute to it, too), so that's what I suggest. Install MacPorts and then do:

$ sudo port install py35-pip py35-matplotlib py35-scipy py35-astropy py35-ipython +notebook

For Python 2.7, replace 35 above with 27.

How to install Python packages

The Python community has developed a huge collection of packages for various tasks. The best place to find Python libraries is the Python Package Index (PyPI, https://pypi.python.org/).

No matter how you installed Python, you can install any Python package from PyPI using pip. For example:

$ pip install astroquery

or

$ pip install --user astroquery

In Anaconda, you can also install some popular packages using conda. For example:

$ conda install astropy

On Debian or Ubuntu, you can get many Python packages with apt-get and on MacPorts with port.

Python basics

We're going to give you a whirlwind tour of Python--just enough to get you up to speed for the rest of this week's hands-on exercises. Because Python is widely used outside of astronomy, there are lots of great resources to learn more when you are ready!

AstroBetter resources list

The Official Python Tutorial

Python as a calculator

Python does what you expect when you type in math expressions:

Data Types

Python has several useful data types. We've already met a couple--Python integers (int) and floating points (float) are just numbers!

Numeric and boolean literals

Python's numeric types include integers and both real and complex floating point numbers:

The expressions above assign values to the variables, so we can use (or modify) the variables instead of the values.

By the way, all of the text on a given line after the trailing hash sign (#) is a comment, ignored by Python.

The arithmetic operators in Python are similar to C, C++, Java, and so on. There is addition (and subtraction):

Multiplication:

Division:

But beware that, like in C, C++, Java, etc., in python 2, division of integers gives you integers. In Python 3, division of integers is true division.

If you want true division in python 2, convert one or both of the operands to floating point. Or even better, use the future module to make your python 2 code compatible with python 3!

If you are using Python 2.2 or higher, you can enable true division by putting the following statement at the top of your script:

from __future__ import division

However, in both Python 2 and Python 3, the double-slash operator // represents integer division:

The % sign is the remainder (modulus) operator:

Exponentiation is accomplished with the ** operator:

Exercise

Calculate the reciprocal of the cube root of 18:

Strings

Strings are text--they are indicated by pairs of quotes:

In the iPython notebook, executing a cell with an expression in it returns that expression, evaluated:

To display strings in functions, etc., we will print them:

This is a Python statement, consisting of the built-in command print and a string surrounded by single quotes.

Note: print changed from a built-in statement to a function in Python 3. If you are using Python 3, just put the arguments in parentheses. To get this behavior in previous versions of Python, put this at the top of your script:

from __future__ import print_function

If you need single or double quotes inside a string, use the opposite kind outside the string:

If you need a string that contains newlines, use triple quotes (''') or triple double quotes ("""):

Strings are concatenated (joined) by "adding" them:

Part of the magic of Python is that every datatype in python is an object. This means that, along with data there are methods that operate on that data. Strings have many useful methods built-in. Let's use a filename as an example:

If we want to check if a file ends with a certain suffix, we can use the .endswith() method:

We can split the filename into pieces (a list--see the next section!) based on certain characters:

Note that we can get help on Python functions with help.

Or if you are typing them in a notebook, try pressing Shift-Tab after you type the function name!

The dir function lists the methods of a given python object. Methods with names surrounded by __double underscores__ are internal; you can ignore them for now.

Exercise

Replace "PTF" in the filename variable with your name, and make the whole string uppercase.

Lists

Lists are flexible containers which can contain python objects of any type. You can add or remove items or alter existing items. Lists are surrounded by square brackets and have items separated by commas.

Once you have made a list, you might want to retrieve a value from it. You index a list with square brackets, starting from zero:

You can access whole ranges of values using slice notation:

Or, to count backward from the end of the list, use a negative index:

Lists can be nested:

len gives the number of items in a list, but note that nested lists only count as one object:

Strings can be treated just like lists of individual charaters:

You can change elements:

Or append elements to an existing list:

Or delete elements:

Exercise

Replace the second element of another_list with 7, remove the nested list, and append 12.

Dictionaries

Sometimes, you want a collection that is like a list, but whose indices are strings or other Python values. That's a dictionary. Dictionaries are handy for any type of database-like operation, or for storing mappings from one set of values to another. You create a dictionary by enclosing a list of key-value pairs in curly braces:

You can index items in dictionaries with square braces [], similar to lists, but you provide the name of the mapping element rather than a numerical position. Dictionaries may store items in any order!

or add items to them:

or delete items from them:

Dictionary keys can be any immutable kind of Python object: tuples, strings, integers, and floats are all fine. Values in a dictionary can be any Python value at all, including lists or other dictionaries.

Exercise

Add two vegetable choices to the menu_options dictionary and replace your least-favorite entree with something better. Add another guest. Then print a string listing your choices for cheese, vegetable, entree, and dessert separated by commas.

The None object

Sometimes you need to represent the absence of a value, for instance, if you have a gap in a dataset. You might be tempted to use some special value like -1 or 99 for this purpose, but don't! Use the built-in object None.

Conditionals

In Python, control flow statements such as conditionals and loops have blocks indicated with indentation. Any number of spaces or tabs is fine, as long as you are consistent within a block. Four spaces is the recommended standard.

You can use the if...elif...else statement to have different bits of code run depending on the truth or falsehood of boolean expressions. For example:

You can chain together inequalities just like in mathematical notation:

You can also combine comparison operators with the boolean and, or, and not operators:

Exercise

Write a conditional expression that tests if a is an odd number less than 6:

The comparison operator is tests whether two Python values are not only equal, but represent the same object. Since there is only one None object, the is operator is particularly useful for detecting None.

Likewise, there is an is not operator:

The in and not in operators are handy for testing for membership in a string, list, or dictionary:

When referring to a dictionary, the in operator tests if the item is among the keys of the dictionary.

The for and while loops

In Python, there are just two types of loops: for and while. for loops are useful for repeating a set of statements for each item in a collection (tuple, set, list, dictionary, or string). while loops are not as common, but can be used to repeat a set of statements until a boolean expression becomes false.

The built-in function range, which returns a list of numbers, is often handy here:

Or you can have the range start from a nonzero value:

You can iterate over the keys and values in a dictionary with .items():

The syntax of the while loop is similar to the if statement:

Exercise

Print the keys in our menu_options dictionary in uppercase, followed by the number of items in the value. You will get an error: how can you work around it?

List comprehensions

Sometimes you need a loop to create one list from another. List comprehensions make this very terse. For example, the following for loop:

is equivalent to this list comprehension:

You can even incorporate conditionals into a list comprehension. The following:

can be written as:

Exercise

Return a list of the lengths (in characters--use len()) of all of the keys in menu_options except number_of_guests. Hint: menu_options.keys() returns a list of keys.

Functions

Functions are created with the def statement. A function may either have or not have a return statement to send back a return value.

If you want to return multiple values from a function, separate them by commas:

If a function returns multiple values, you can automatically unpack them into multiple variables:

If you pass a mutable value such as a list to a function, then the function may modify that value. For example, you might implement the Fibonacci sequence like this:

You can also give a function's arguments default values, such as:

If a function has a large number of arguments, it may be easier to read if you pass the arguments by keyword, as in:

Exercise

Write a function that plays "FizzBuzz":

given an input number n, it prints the numbers from 1-n, unless the number is a multiple of 3, in which case print "Fizz". If the number is a multiple of 5, print "Buzz". If both, print "FizzBuzz".

run your function for n=20.

Python standard library

Python comes with an extensive standard library consisting of individual modules that you can opt to use with the import statement.

You can import specific functions or variables with the from x import y syntax:

or import the whole module in a "namespace" with import x or import long_module as y. Then you can use all of the functions and variables by calling x.func() or y.variable.

Some particularly useful parts of the Python standard library are:

• random: random number generators
• pickle: read/write Python objects into files
• os: operating system services
• os.path: file path manipulation
• sqlite3: SQLite database acces
• subprocess: launch external processes
• email: compose, parse, receive, or send e-mail
• pdb: built-in debugger
• re: regular expressions
• SimpleHTTPServer: built-in lightweight web server
• optparse: build pretty command-line interfaces
• itertools: exotic looping constructs
• multiprocessing: parallel processing

help() and dir() are useful for finding what commands are available in the standard library.

Exercise:

Use the os module to print the current working directory. Then list its contents. (Look here for functions.)

Error handling

It can be important for your code to be able to handle error conditions. For example, let's say that you are implementing a sinc function:

Oops! We know that by definition $\mathrm{sinc}(0) = 1$ , so we should catch this error:

Reading and writing files

The built-in open function opens a file and returns a file object that you can use to read or write data. Here's an example of writing data to a file:

And here is reading it:

Numpy & Matplotlib

Numpy provides array operations and linear algebra to Python. A Numpy array is a bit like a Python list, but arrays usually contain only numbers and support elementwise arithmetic. For example:

By default, multiplication is elementwise:

To perform matrix multiplication, either convert arrays to np.matrix or use np.dot:

Numpy arrays may have any number of dimensions:

An array has a certain number of dimensions denoted .ndim:

and the dimensions' individual lengths are given by .shape:

and the total number of elements by .size:

You can also perform comparison operations on arrays...

Although a boolean array doesn't directly make sense in an if statement:

You can use conditional expressions like indices:

Or manipulate individual rows:

Or individual columns:

Other useful features include various random number generators:

And statistical functions like mean and standard devation:

Matplotlib is the standard python plotting package. Line plots, histograms, scatter plots, and error bar plots are the most useful:

Exercise

Generate x = 100 random numbers between 0 and 2 pi, and make a scatter plot of x vs y = sin(x).

Extra challenge After you successfully plot that, color all the points with absolute value greater than 0.8 red and the others black.

Astropy

Astropy is a core Python package for astronomy. It is formed from the merger of a number of other Python astronomy packages, but also contains a lot of original code. Core features include:

• astropy.constants, astropy.units: Physical constants, units, and unit conversion
• astropy.time: Manipulation of dates and times
• astropy.coordinates: Representation of and conversion between astronomical coordinate systems
• astropy.table: Tables and gridded data
• astropy.io.fits: Manipulating FITS files
• astropy.io.ascii: Manipulating ASCII tables of many different formats
• astropy.io.votable: Virtual Observatory tables
• astropy.wcs: World Coordinate System transformations
• astropy.cosmology: Cosmological calculations
• astropy.stats: Astrostatistics
• astropy.modeling: multi-D model fitting Swiss army knife

The Astropy project also has sevearl "affiliated packages" that have similar design but are maintained separately, including:

• APLPy: High-level astronomical map making with Matplotlib
• WCSAxes: Nuts-and-bolts astronomical mapmaking for Matplotlib experts
• Photutils: Aperture photometry
• Astroquery: Query astronomical databases

Let's experiment by opening up a P48 image. We'll need several modules from the Astropy package for this exercise.

We just have time to scratch the surface of what astropy is capable of. astropy.constants and astropy.units are both quite useful:

Unit conversion is as simple as adding .to(unit):

Exercise

Use astropy.constants to calculate the mean density of Earth in g/cm^3.

Let's work with one of the PTF coadded images we'll use later in the hands-on sessions.

Let's grab the first HDU ("header data unit") of this FITS file:

Then let's take a look at the contents of the header:

Now let's plot the image data. But let's use sigma-clipping to pick a nice scale for the image.

Let's look at a star that will quickly become our favorite: it has coordinates $\alpha_\mathrm{J2000}, \delta_\mathrm{J2000} = (312.503802, -0.706603)$.

We'll open the catalog file from IPAC and find this object.

The astropy SkyCoord code enables efficient crossmatching between positions.

Let's create a coordinates object to represent the target that we are searching for:

Now, let's match this position against the source catalog:

In uncrowded PTF fields like this one, we typically use a match radius of $1.5^{\prime\prime}$. separation is an astropy Angle, which has units, so we can easily test if the matched object is the correct one:

Here is the closest-matching row:

Let's look at where these catalog sources appear on the image. We'll need to convert the RA, Dec in the table to pixel coordinates. That's where the World Coordinate System (WCS) transformation comes in.

For more sophisticated plotting of astronomical images, WCSAxes or APLPy are useful modules.

Exercise

The catalog file above is for R-band. Load the g-band catalog also and crossmatch both catalogs against each other using a $1.5^{\prime\prime}$ match radius. Scatter plot g-r color against r-band magnitude. (You'll want to adjust the x & y limits to zoom in on the data, i.e. ignore extreme outliers - sources with MAG_AUTO = 99.)