Links to cool pages on plotting with python

A very interesting blog on python there: http://zulko.github.io/blog/2014/09/20/vector-animations-with-python/

A package to plot interactively: http://bokeh.pydata.org/en/latest/index.html

And don't forget the impressive Mayavi package: http://mayavi.sourceforge.net/

The Seaborn library is useful to make some nice plots: http://stanford.edu/~mwaskom/software/seaborn/


A kind of summary of some of these different solutions is made there: http://pbpython.com/visualization-tools-1.html

Nice book on python, numpy, ipython etc...

A very up-to-date (august 2014) electronic book on python named "Introduction to Python for Econometrics, Statistics and Data Analysis" from Kevin Sheppard is available here:
https://www.kevinsheppard.com/images/0/09/Python_introduction.pdf
It's mainly on econometric, but most of the tools described there are also useful for astronomers.

Sending requests to MySQL and receiving the result from python, using PyMySQL

 Modern astrophysics is using every day more big databases. One of the mostly used interface to databases is MySQl (or its recent free fork MariaDB). I present in this lecture a python library to deal with MySQL databases: PyMySQL. In the lecture the examples are using access to 3MdB (https://sites.google.com/site/mexicanmillionmodels/), which requires a password. You can ask for it to me, or adapt the example to connect to other databases.
The lecture is here:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_PyMySQL.ipynb


An introduction to MySQL can be found here: https://github.com/Morisset/Python-lectures-Notebooks/blob/master/MySQL.pdf

Optimization, calling Fortran

2014 Python Lecture. Part IX

In this latest lecture of this series, I'll present some tools to optimize your code by CPU and memory profiling. It also contains some tips on using the python debugger.
The notebook is there:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Optimization.ipynb


I also give some indications on how one can call Fortran routines from within python, to accelerate the execution of some part the the code.
Here are small examples: https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Calling%20Fortran.ipynb

Object Oriented Programing. Objects, classes, etc...

2014 Python Lecture. Part VIII

In this lecture I'll introduce the basic (and some not that basic) concepts of Object Oriented Programing. I'll use an example to show how to:

• use functions to do simple jobs
• but use objects when things start to be more complex
• define classes, objects, attributes, methods, etc...
• use *args and **kwargs in functions calls
• use the class variables
• add functionalities to classes and objects
• use class inheritance
• use attributes properties

The notebook is here:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/OOP.ipynb

The astropy library

2014 Python Lecture. Part VII

The Astropy Project is a community effort to develop a single core package for Astronomy in Python and foster interoperability between Python astronomy packages. More informations here: http://www.astropy.org/

In this lecture we will see some of the facilities of the astorpy library, including:

• Constants and Units
• Data Table (a very useful one!)
• Time and Dates
• Etc...

The lecture is here:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Using_astropy.ipynb

Useful libraries

2014 Python Lecture. Part VI

This lecture will give some insights to the most useful python libraries. It is NOT exhaustive, you have to read the corresponding manual pages to find the best use you can have of them. The list of all python-included libraries is here: https://docs.python.org/2/library/

I mention in this lecture:

• time and datetime
• timeit
• os
• sys
• subprocess
• glob
• re
• urllib2

The lecture is here:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Useful_libraries.ipynb

Introduction to Scipy

2014 Python Lecture. Part V

SciPy (pronounced “Sigh Pie”) is a Python-based ecosystem of open-source software for mathematics, science, and engineering. In particular, these are some of the core packages: Numpy, Scipy library, Matplotlib, ipython, Simpy, and Panda.


In this lecture, I will show exemples covering:

• Some useful methods
• nanmean
• constants
• Integrations
• Interpolations
• 2D-interpolations
• data fitting
• multivariate estimation

The lecture is here:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Scipy.ipynb

How to make plots, images, 3D, etc, using Matplotlib

2014 Python Lecture. Part IV

This lecture is dedicated to the plotting library matplotlib. The topics are:

• Simple plot
• Controlling colors ans symbols
• Overplot
• Fixing axes limits
• Labels, titles
• Legends
• The object oriented way to use Matplotlib
• Scatter
• log plots
• Multiple plots
• Everything is object
• Error bars
• Sharing axes
• Histograms
• Boxplots
• Ticks, axes and spines
• A plot inside a plot
• Play with all the objects of a plot
• Filled regions
• 2D-histograms
• 2D data sets and images
• Contour
• 3D scatter plots
• Saving plots
• Access and clear the current figure and axe
• What's happen when not in a Notebook? plt.show() and plt.ion() commands

The lecture Notebook is there:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Matplotlib.ipynb

Interacting with files: reading writing, ascii and fits

2014 Python lecture. Part III

It's time to play with files containing data! In this lecture, we'll see how to read and write files (ascii and fits).

• Reading a simple ASCII file
• How to treat special rows (comments, header)
• classical way
• using numpy.loadtxt
• using numpy.genfromtxt

• Dealing with missing data
• Data in a fixed size format
• Writing files
• simple method
• Pickle files (python format)
• FITS files

The ipython notebook is there:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Interact%20with%20files.ipynb

Introduction to Numpy

2014 Python lecture. Part II

The introduction to Numpy can be seen here:
https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_numpy.ipynb

The topics that are presented are:

• The Array class
• create an array
• 1D, 2D 3D arrays
• creating array from scratch
• arrays share memory (views)
• random generator
• timing a command
• slicing arrays
• assignments
• using masks
• the where function
• some operations with arrays
• broadcasting
• calling scripts
• structured arrays and record arrays
• NaN other ANSI values.

Any comments are welcome.

Chris.Morisset a t Gmail.com

Python: Basics

2014 Python lecture. Part I

The introduction to Python I'm giving at IA-UNAM is accessible here:

https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/intro_Python.ipynb

I will modify this notebook during the lecture (August 2014), so reload it to have the latest version.

The topics of this first lecture are:

• Using python as a calculator
• assignments
• comments
• types
• complex numbers
• booleans
• printing strings
• strings
• Tuples, lists and dictionaries
• Blocks
• List and dictionary comprehension
• Functions, procedures
• Scripting
• Importing libraries

If you want to have an interactive session with this lecture using the ipython notebook facilities, follow the link above and download the ipynb file (download button at the right top of the web page). Save the file in a directory from where you execute the following (you must have a recent version of ipython installed):
ipython notebook
It should open a new tab in your web browser, with the list of ipynb files in the directory. Click on the one you want, will open a new tab similar to the first one, but this one is executed on YOUR computer, it means you are able to interact with the commands. You can change the commands, and execute a cell by SHIFT-ENTER. You can add comments in new cells, and save the result.

Any comments are welcome.

Brief introduction to Python

2014 Python lecture. Part 0

Back to the Python lecture, I want to share here the very quick introduction I gave before starting to play with python: https://github.com/Morisset/Python-lectures-Notebooks/blob/master/Notebooks/Intro_1.pdf
You may want to install python from Ureka from this site: http://ssb.stsci.edu/ureka/

Using ipython Notebook to teach scientific python

A very good and efficient way to teach python and python related tools, is to use ipython Notebook: http://ipython.org/notebook.html

As example of this use, the following link is a collection of lectures on python, numpy, scipy, matplotlib, use of Fortran from python, etc:
https://github.com/jrjohansson/scientific-python-lectures
Enjoy them.

Installing python and nice lectures on scientific Python.

Since more than one year without any message!... And the new one is almost nothing from me, just links to good pages.

Two easy ways to install python+ipython+numpy+matplotlib+scipy:

Anaconda from continuum (but only 64bit version for OSX, which can be a problem for MySQLdb):
Ask for the academic licence if you can.
http://continuum.io/
Once installed, you will need to add the anaconda/bin directory to your PATH and the anaconda directory to your PYTHONPATH.

Canopy from Entought:
https://www.enthought.com/products/canopy/
Once installed, you must setup the virtual environment by adding to your .tcshrc:
setenv VIRTUAL_ENV /Users/YOURNAME/Library/Enthought/Canopy_32bit/User
setenv PATH  $VIRTUAL_ENV/bin:$PATH

!!! Warning !!! When using this virtual environment, don't install using pip with the --user option. Install directly, as for example:
pip install pyfits

UPDATE: Another (better) package comes from STSCI, it's UREKA: http://ssb.stsci.edu/ureka/

Here follows good pages to learn python. Those pages are made using Notebook, which is very efficient to show/share python programs.

Here are the links:

http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-0-Scientific-Computing-with-Python.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-1-Introduction-to-Python-Programming.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-2-Numpy.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-3-Scipy.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-4-Matplotlib.ipynb
http://nbviewer.ipython.org/urls/raw.github.com/jrjohansson/scientific-python-lectures/master/Lecture-6B-HPC.ipynb

Playing with arrays: slicing, sorting, filtering, where function, etc.

You can also read the more recent post on Numpy here: http://python-astro.blogspot.mx/2014/08/introduction-to-numpy.html

Most of the data we have to manage are in form of arrays. That's why it's quite always necessary to import the numpy library to easily deal with tables, arrays, images, cubes, etc...
In this post I will discuss some of the methods and functions one may have to use in this context.

In Python, one can use lists, tuples and dictionaries to put different elements together. They even can contain elements of different types. But nothing better than numpy arrays to really "play" with the data, extract subsets, combine them using arithmetic operations.

We already discuss a little some array commands in a previous post: read-ascii-file-cont. Here we are more complete, but not exhaustive, as numpy is a whole world... Have a look at the reference guide: http://docs.scipy.org/doc/numpy/reference/ and other references at the end of this post.

Create arrays

a = np.array([1, 2, 3, 7, 5])

a

array([1, 2, 3, 7, 5])

numpy's arrays cannot contain values of different types, the values are transformed if necessary, from integer into real, and from real into string:

b = np.array([1, 2, 3, 4.])

b

array([ 1.,  2.,  3.,  4.])

c = np.array([1, 2, 3., '4'])

c

array(['1', '2', '3', '4'],
      dtype='|S1')

shapes and sizes of the arrays are obtained using method and functions:

print c.shape

(4,)

print len(c)

4

Different ways to create arrays:

a = np.arange(1, 10, 0.5)

print a

[ 1.   1.5  2.   2.5  3.   3.5  4.   4.5  5.   5.5  6.   6.5  7.   7.5  8.
  8.5  9.   9.5]

TAKE CARE, THE LATEST ELEMENT IS NOT WHAT YOU MAY EXPECT...
It's because the first element in Python is indexed by 0. The simple use of np.arange is :

a = np.arange(10)

print a

[0 1 2 3 4 5 6 7 8 9]

a starts at 0, and has 10 elements, that's ok.

Easy to create linearly spaced arrays:

c = np.linspace(0, 1, 4)

print c

[ 0.          0.33333333  0.66666667  1.        ]

or log spaced:

c = np.logspace(0, 1, 4)
c

array([  1.        ,   2.15443469,   4.64158883,  10.        ])

which is actually the same as :

c = 10**np.logspace(0, 1, 4)

To create arrays of 0. or 1.:

a = np.zeros(10)

b = np.ones((5, 4))

Here b is a 2D array.

b.shape

(5, 4)

You can use a list (or a tuple, or an array) and replicate it:

a = np.array([1, 2, 3])

b = np.tile(a, 5)

print b

[1 2 3 1 2 3 1 2 3 1 2 3 1 2 3]

c = np.tile(a, (5, 1))

print c

[[1 2 3]
 [1 2 3]
 [1 2 3]
 [1 2 3]
 [1 2 3]]

You can also create 2D arrays from 1D vectors:

x = np.linspace(-1, 1, 100)

y = x

X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel

print X.shape

(100, 100)

Outer products are obtain using... np.outer:

a = np.array([1, 2, 3])
b = np.array([1, 10, 100])

np.outer(a, b)

array([[  1,  10, 100],
       [  2,  20, 200],
       [  3,  30, 300]])
np.outer(b, a)

array([[  1,   2,   3],
       [ 10,  20,  30],
       [100, 200, 300]])

Broadcasting

Python numpy is able to add dimensions to arrays to perform operations:

a = np.array([1, 2, 3, 4])

b = np.ones((6, 4))

a * b

array([[ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.],
       [ 1.,  2.,  3.,  4.]])

Nice page on this:
http://www.scipy.org/EricsBroadcastingDoc

Indexing and Slicing

The access to elements of arrays is done using []:

a = np.arange(10)

a[5]

0

a[-1]  # last elements

9

In case of N-dims arrays, one can extract slides using :

a = np.array([1, 2, 3, 4, 5])

b = np.array([1, 10, 100, 1000])

c = np.outer(a, b)

c.shape

(5, 4)

print c 

[[   1   10  100 1000]
 [   2   20  200 2000]
 [   3   30  300 3000]
 [   4   40  400 4000]
 [   5   50  500 5000]]

print c[:,2]

[100 200 300 400 500]

print c[-1,:]

[   5   50  500 5000]

print c[-1, -1]

5000

There is a lot of methods in the array object, the best page to learn more is there: http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object

Filtering

Sometimes we need to extract elements from arrays following some criteria. There is basically two ways to do this: defining a set of indices where the condition is completed (a la WHERE in IDL), or defining a boolean mask.

Let's first define a 2D array made of 10 times 1000 random values:

a = np.random.random((10, 1000))

We want to extract the values where the 2nd and the 4th 1000-elements vectors are greater than 0.5.

Using the numpy.where function:

w1 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))

as the result is a tuple of indices, and in this case the dimension of the result is 1, better extract on the fly the array of indices from the tuple:

w2 = np.where((a[1,:] > 0.5) & (a[3,:] > 0.5))[0]

len(w2)

248  # your result may differ form this value, but must be close to 250.

We can now reduce the initial array to the desired values:

b = a[:, w2]

b.shape

(10, 248)

One can also build an array of boolean:

mask = (a[1,:] > 0.5) & (a[3,:] > 0.5)

which is actually similar to

mask2 =  np.where((a[1,:] > 0.5) & (a[3,:] > 0.5), True, False)

It can be used the same way as before:

b = a[:, mask]

b.shape

(10, 248)

These latter cases result in a 1000-element array, filled with True and False. To know the number of True values, just sum the array:

mask.sum()

248  # your result may differ form this value, but must be close to 250.

One big advantage of the mask technique is that you can combine different masks:

mask1 = a[1,:] > 0.5

mask2 = a[3,:] > 0.5

mask_total = mask1 & mask2

b = a[:, mask_total]

Sorting

a = np.array([1,2,45,2,1,46,7,-1])

b1 = np.sort(a)

ib = np.argsort(a)

b2 = a[ib]

print b1

[-1  1  1  2  2  7 45 46]

print b2 

[-1  1  1  2  2  7 45 46]

More on numpy arrays

http://scipy.org/Numpy_Example_List_With_Doc
http://scipy-lectures.github.com/intro/numpy/numpy.html#the-numpy-array-object

Plotting

Don't miss the update: "2014 Python lecture on Matplotlib" on this blog: http://python-astro.blogspot.mx/2014/09/2014-python-lecture.html

Plotting is certainly one of the most common task one would do in Python (at least for astronomers).
There is various libraries to draw plots in Python, but the mostly used and powerful is perhaps matplotlib.
The plotting part is pyplot, so before any use, one must import it, with the (widely used) alias plt:

import matplotlib.pyplot as plt

The best to see how it works is to dive into the website where there is a lot of examples of plots with the code used to generate them: http://matplotlib.sourceforge.net/, especially the gallery: http://matplotlib.sourceforge.net/gallery.html
In the following we will present a few examples to help you to start.

X-Y plot

Let's say you have 2 vectors x and y and wanna plot y vs. x. Here are some examples of how to do this.
First generate x and y (don't forget to import numpy as np):

x = np.linspace(0., 4 * np.pi, 100)

y = np.sin(x)

plt.plot(x, y)
plt.show()

This latest command is not always necessary, if you run ipython with the --pylab option.
As you can see, python automatically define the axis ranges and draw the plot using a default blue color.
If you want to overplot another function, just call plt.plot again:

plt.plot(x, y**2)

The figure can be cleaned before another plot using

plt.clf()

You may need to draw various figure in different windows at the same time. Every call to

plt.figure()

will start a new figure where the next plot(s) will be drawn. The figures (each one in a separate window) are identified with numbers, the first one being 1.
When calling plt.figure(N), a new figure is created if figure N doesn't exist, and focus will be on figure N otherwise.

The color of the line is defined using c or color keyword:

plt.plot(x, y**3, color = 'red')

plt.plot(x, y**4, c='b')

Abbreviation	Color
b	blue
g	green
r	red
c	cyan
m	magenta
y	yellow
k	black
w	white

Symbols and line style can also be easily defined:

plt.plot(x, y**3, color='r', marker='o', linestyle=':')

Symbol	Description
-	solid line
--	dashed line
-.	dash-dot line
:	dotted line
.	points
,	pixels
o	circle symbols
^	triangle up symbols
v	triangle down symbols
<	triangle left symbols
>	triangle right symbols
s	square symbols
+	plus symbols
x	cross symbols
D	diamond symbols
d	thin diamond symbols
1	tripod down symbols
2	tripod up symbols
3	tripod left symbols
4	tripod right symbols
h	hexagon symbols
H	rotated hexagon symbols
p	pentagon symbols
|	vertical line symbols
_	horizontal line symbols
steps	use gnuplot style ‘steps’ # kwarg only

Other useful line properties:

Property	Value
alpha	alpha transparency on 0-1 scale
antialiased	True or False - use antialised rendering
color	matplotlib color arg
data_clipping	whether to use numeric to clip data
label	string optionally used for legend
linestyle	one of - : -. -
linewidth	float, the line width in points
marker	one of + , o . s v x > <, etc
markeredgewidth	line width around the marker symbol
markeredgecolor	edge color if a marker is used
markerfacecolor	face color if a marker is used
markersize	size of the marker in points

If you need more control on the markers, better use scatter. For example, if one need ot change the size of the symbol according to a function:

plt.scatter(x, y**3, c=abs(y), marker='s', s=10+abs(y)*100, edgecolors='none')

The labels are set after the plot is done:

plt.xlabel('X')

plt.ylabel('Some trig function')

LaTex fans are welcome:

plt.title(r'Example of LaTex $\alpha_{\beta}$') #Notice the r before the string

One can change the axis ranges afterward:

plt.xlim((0, 10))

plt.ylim((-2, 2))

Multiple plots

To draw multiple plots of the same figure:

plt.subplot(N_y, N_x, N)

for i in np.arange(9):

    plt.subplot(3,3,i+1)

    plt.plot(x, y**i)
    plt.ylim((-2, 2))

log plots

using plt.semilogx, plt.semilogy and plt.loglog

plt.loglog(x)

plt.grid(True, which='minor')

Contours

Contour plots are done with plt.contour and plt.contourf.

x = np.linspace(-1, 1, 100)

y = x

X, Y = np.meshgrid(x, y) #this create 2D arrays containing x and y for each pixel

dist =  (X**2 + Y**2)**0.5

plt.contourf(X, Y, dist)

plt.colorbar()

CS = plt.contour(X, Y, dist, colors ='black', linewidths = 4)

plt.clabel(CS) #to print label on each contour

Saving the plot

The result of the plot can be save in PDF, EPS, JPG, BMP format, using:
plt.savefig('fig1.pdf')

More on plotting:

http://scipy-lectures.github.com/intro/matplotlib/matplotlib.html

Play with FITS files

The same lecture is now in the Notebook format, see there: http://python-astro.blogspot.mx/2014/09/interacting-with-files-reading-writing.html

What is the FITS format?

The FITS format is the most popular way to save and interchange astronomical data. The files are organized in units each of which contains a human readable header and a data. This structure is refereed as HDUs (Header/DATA Unit).

A FITS file can contain one or more HDUs, the first of which is called "primary" and the rest are called "extensions". The primary HDU usually contains 1D spectrum, 2D image or 3D data cube, although any dimension from 0 to 999 are possible. The data are 1, 2 or 4 bytes integers or 4 or 8 bytes real numbers.

The extensions can contain or arrays as in the primary HDU or ascii tables or binary tables.
If a FITS file contains only tables, it primary HDU does not contain data, but only header.

Both headers and data in a FITS file are organized in blocs of 2880 bytes. The header contain 80 bytes lines each of which consists of a keyword of 8 bytes followed in most of the cases by '= ' in the position 9 and 10 and then the value of the keyword. The rest of the line is a comment string beginning with '/'. Each header begins with the following lines

SIMPLE  =                    T / file conforms to FITS standard
BITPIX  =                   16 / number of bits per data pixel
NAXIS   =                    2 / number of data axes
NAXIS1  =                  440 / length of data axis 1
NAXIS2  =                  300 / length of data axis 2

which defines the format of the file as standard FITS, the data format and the dimensions of the stored data.

One block of 2880 bytes contains 36 lines of 80 characters per line. The header can have several blocks of 36 lines. The last block is identified by the presence of the keyword 'END' The next 2880 bytes block contains the first part of the data. The empty lines after 'END' keyword are filled with blanks and the unused bytes from the end of the data to the end of the 2880 bytes block are filled with NULLs.

Full description of the FITS format can be found at http://fits.gsfc.nasa.gov/fits_primer.html

PyFITS

As one can see, reading FITS files is not as simple as reading columns of data. Fortunately, there is a python package called PyFITS which makes the task much easier.

The package can be downloaded from

http://www.stsci.edu/institute/software_hardware/pyfits

The documentation can be found on the same web page.

Although it is possible to download and install the package it is much easier to install it with

easy_install pyfits
or
pip pyfits

Information on easy_install: HERE and on pip: HERE.
I did that on both Mac and Linux machines and it working nicely.

How to read FITS file

First, as usual, one has to import the pyfits package

import pyfits

then the file is read with

hdulist = pyfits.open('n10017o.fits')

where I used one of my FITS files from San Pedro Martir echelle spectrograph. The file can be downloaded from HERE.

The result hdulist is a list of HDU objects. In the case of a simple file, there is only one primary HDU so the list contains only one element

len(hdulist)
1

The information on what the file contains can be obtained by calling the info() method:

hdulist.info()
Filename: n10017o.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU      62   (2154, 2048)   int16

The table said that there is only a primary HDU which contains 2154 X 2048 image with data stored in 2 bytes (16 bits) integers.

As described above, the HDU (header/data unit) contains header and data. The header is a dictionary. To see what keywords were used in the header one can do:

hdulist[0].header.keys()
['SIMPLE',
 'BITPIX',
 'NAXIS',
 'NAXIS1',
 'NAXIS2',
 'EXTEND',
 'COMMENT',
 'BZERO',
 'BSCALE',
 'EXPTIME',
 'DETECTOR',
 'ORIGIN',
 'OBSERVAT',
 'TELESCOP',
 'LATITUDE',
 'LONGITUD',
 'ALTITUD',
 'SECONDAR',
 'TIMEZONE',
 'OBSERVER',
 'OBJECT',
 'INSTRUME',
 'GAINMODE',
 'FILTER',
 'IMGTYPE',
 'EQUINOX',
 'ST',
 'UT',
 'JD',
 'DATE-OBS',
 'CCDSUM',
 'RA',
 'DEC',
 'AH',
 'AIRMASS',
 'TMMIRROR',
 'TSMIRROR',
 'TAIR',
 'XTEMP',
 'HUMIDITY',
 'ATMOSBAR',
 'WIND',
 'WDATE',
 'DATE',
 'NAMPS',
 'CCDNAMPS',
 'AMPNAME',
 'CREATOR',
 'VERSION',
 'HISTORY']

and to get the value of a given keyword

hdulist[0].header['OBJECT'] 
'BD +59 363'

The header can be printed as it apears in the file by

print hdulist[0].header.ascardlist()
SIMPLE  =                    T / conforms to FITS standard                     
BITPIX  =                   16 / array data type                               
NAXIS   =                    2 / number of array dimensions                    
NAXIS1  =                 2154 / length of data axis 1                         
NAXIS2  =                 2048 / length of data axis 2                         
EXTEND  =                    T                                                 
COMMENT FITS (Flexible Image Transport System) format is defined in 'Astronomy 
COMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H  
BZERO   =                32768 / BZERO                                         
BSCALE  =                    1 / BSCALE                                        
EXPTIME =                  0.0 / Integration Time, sec.                        
DETECTOR= 'e2vm2 E2V-4240'     / CCD Type                                      
ORIGIN  = 'UNAM    '           / OAN SPM, IA-UNAM                              
OBSERVAT= 'SPM     '           / Observatory                                   
TELESCOP= '2.12m   '           / Telescope                                     
LATITUDE= '+31:02:39'          / Latitude                                      
LONGITUD= '115:27:49'          / Longitud                                      
ALTITUD =                 2800 / altitud                                       
SECONDAR=                   -1 / F/ Secondary type                             
TIMEZONE=                    8 / Time Zone                                     
OBSERVER= 'Leonid  '           / Observer's Name                               
OBJECT  = 'BD +59 363'         / Object                                        
INSTRUME= 'Echelle '           / Instrument                                    
GAINMODE=                    1 / Gain factor in the CCD                        
FILTER  = 'None    '           / Filter                                        
IMGTYPE = 'zero    '           / Image Type                                    
EQUINOX =               2011.7 / Equinox                                       
ST      = '02:19:51.2'         / Sideral Time                                  
UT      = '11:28:28'           / Universal Time                                
JD      =            2455803.5 / Julian Date                                   
DATE-OBS= '2011-08-30'         / Observation Date UTM                          
CCDSUM  = '1 1     '           / Binning [ Cols:Rows ]                         
RA      = ' 02:13:22.2'        / Right Ascension                               
DEC     = ' 52''14''44.0'      / Declination                                   
AH      = '  00:06:28.0'       / Hour Angle                                    
AIRMASS =                1.073 / Airmass                                       
TMMIRROR=                    0 / Primary Mirror Temperature (celsius degree)   
TSMIRROR=                    0 / Secundary Mirror Temperature (celsius degree) 
TAIR    =                    0 / Internal Telescope Air Temperature (celsius deg
XTEMP   =                 14.6 / Exterior Temperature (celsius degree)         
HUMIDITY=                 44.0 / % external Humidity                           
ATMOSBAR=                731.7 / Atmosferic Presure in mb                      
WIND    = 'S at 24.1 km/h'     / Wind Direction                                
WDATE   = '11:28:10, 08/30/11' / Weather Acquisition Date (Local time)         
DATE    = '2011-08-30T11:28:29' / file creation date (YYYY-MM-DDThh:mm:ss UT)  
NAMPS   =                    1 / Number of Amplifiers                          
CCDNAMPS=                    1 / Number of amplifiers used                     
AMPNAME = '1 Channel'          / Amplifier name                                
CREATOR = 'Python Oan ccds'    / Name of the software task that created the file
VERSION = '4.12D   '           / Application Software Version                  
COMMENT Visit our weather site http://www.astrossp.unam.mx/weather15           
COMMENT for complete meteorological data of your observation night             
HISTORY bin2fits V1.0                                                          
HISTORY Programmer: Enrique Colorado [ colorado@astrosen.unam.mx ]             
HISTORY Observatorio Astronomico Nacional -UNAM                                
HISTORY V1.00 By Arturo Nunez and Colorado >Ported to Python using pyfits      
HISTORY V0.50 By E. Colorado >Added interior mirrors temperatures              
HISTORY V0.49 By E. Colorado >Added BIASSEC parameter                          
HISTORY V0.48 By E. Colorado >Aditional info for autofocus calculations        
HISTORY V0.4 By E. Colorado >Now we include timezone, and remove lat. sign     
HISTORY V0.3 By E. Colorado >Now we include weather data                       
HISTORY V0.2 By E. Colorado >General OAN Working Release                       

The data in the file are accessible with

data = hdulist[0].data

and can be seen with [don't forget to import matplotlib.pyplot as plt before running this]:

plt.imshow(data)

A column from the data can be plotted with

plt.plot(data[:,1000])

where I am plotting the column number 1000. In the same way a line from the data is plotted with:

plt.plot(data[1000,:])

The data are numpy object so all manipulations are available.

Some more on displaying images

The main matplotlib function to display images is imshow. It has several parameters of which te most important are the one which controls the color scheme and the ones which control the dynamic range.

plt.imshow(data, cmap=cm.gray, vmin=1000, vmax=10000)

cmap controls the pseudocolor map (same as color= in IDL). The predefined color maps can be seen on http://matplotlib.sourceforge.net/examples/pylab_examples/show_colormaps.html

vmin and vmax controls the range of values which are mapped in [0,1] range where 0 is black (the darkest color) and 1 is white (the lightest color).

origin = 'upper' | 'lower', place the pixel with coordinates 0,0 at the upper or lower corner of the plot

extent = (xmin,xmax,ymin,ymax) - This parameter defines the numbers written on the axes. No changes on the image.

Using FITS tables

For this example I'll use a spectrum obtain with the high dispersion camera on board of IUE.
The file is opened as usual:
hdulist = pyfits.open('swp04345.mxhi')

Download the file from THERE.

but now hdulist has 2 elements (2 header/data units):

len(hdulist)
2

We can see that the primary header has dimension (), son does not contain any data. The data are in the extension.

hdulist.info()
Filename: swp04345.mxhi
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU     421   ()           uint8
1    MEHI        BinTableHDU     61   60R x 17C    [1B, 1I, 1D, 1I, 1D, 1E, 1E, 768E, 768E, 768E, 768I, 768E, 768E, 1I, 1I, 1E, 7E]

The first header contains the minimal infirmation:

print hdulist[0].header.ascardlist()[:5]
SIMPLE  =                    T / Standard FITS Format                          
BITPIX  =                    8 / Binary data                                   
NAXIS   =                    0 / Two-dimensional image                         
EXTEND  =                    T / Extensions are present                        
TELESCOP= 'IUE     '           / International Ultraviolet Explorer            

The number of axis is 0 which means there is no data block in the primary HDU. The header of the second HDU begins with the keyword XTENSION and with the specification of the data

print hdulist[1].header.ascardlist()[:5]
XTENSION= 'BINTABLE'           / Binary table extension                        
BITPIX  =                    8 / Binary data                                   
NAXIS   =                    2 / Two-dimensional table array                   
NAXIS1  =                16961 / Width of row in bytes                         
NAXIS2  =                   60 / Number of orders                              

To progress further we need to know what is in the table. As usual, the columns have names and type of the stored data. These information can be obtained using the column attribute of hdulist:

cols = hdulist[1].columns

cols.info
<bound method ColDefs.info of ColDefs(
    name = 'ORDER'; format = '1B'; unit = ' '
    name = 'NPOINTS'; format = '1I'; unit = ' '
    name = 'WAVELENGTH'; format = '1D'; unit = 'ANGSTROM'
    name = 'STARTPIX'; format = '1I'; unit = 'PIXEL'
    name = 'DELTAW'; format = '1D'; unit = 'ANGSTROM'
    name = 'SLIT HEIGHT'; format = '1E'; unit = 'PIXEL'
    name = 'LINE_FOUND'; format = '1E'; unit = 'PIXEL'
    name = 'NET'; format = '768E'; unit = 'FN'
    name = 'BACKGROUND'; format = '768E'; unit = 'FN'
    name = 'NOISE'; format = '768E'; unit = 'FN'
    name = 'QUALITY'; format = '768I'; unit = ' '
    name = 'RIPPLE'; format = '768E'; unit = 'FN'
    name = 'ABS_CAL'; format = '768E'; unit = 'ERG/CM2/S/A'
    name = 'START-BKG'; format = '1I'; unit = 'PIXEL'
    name = 'END-BKG'; format = '1I'; unit = 'PIXEL'
    name = 'SCALE_BKG'; format = '1E'; unit = ' '
    name = 'COEFF'; format = '7E'; unit = ' '
)>

the cols.info returns the names of the columns and the information of their format and units.

The data are available using (this example is NOT the right way of plotting the data, it's just an example) [and don't forget to import numpy as np to have np.arange working]:

data1 = hdulist[1].data

DTs =  data1.ABS_CAL

WLs = data1.WAVELENGTH

DWs = data1.DELTAW

for WL, DW, DT in zip(WLs, DWs, DTs):

    plot(WL + np.arange(len(DT)) * DW, DT)

Writing FITS files.

The creation of a FITS file pass through 4 steps.

1) Creation of numpy array with the data.

x = np.arange(100)
  
2) Creation of the HDU from the data.

hdu = pyfits.PrimaryHDU(x)
   
thus created, the hdu has its basic header and the data.

print hdu.header.ascardlist()
SIMPLE  =                    T / conforms to FITS standard                     
BITPIX  =                   64 / array data type                               
NAXIS   =                    1 / number of array dimensions                    
NAXIS1  =                  100                                                 
EXTEND  =                    T                                                 

3) Adding additional keywords to the header. The automatically created header contains only the required minimum of keywords. If additional keywords are needed they are added with

hdu.header.update('testkey',0.001,'some test value')

print hdu.header.ascardlist()
SIMPLE  =                    T / conforms to FITS standard                     
BITPIX  =                   64 / array data type                               
NAXIS   =                    1 / number of array dimensions                    
NAXIS1  =                  100                                                 
EXTEND  =                    T                                                 
TESTKEY =                0.001 / some test value                               

4) Once all the keywords are ready, the final HDU list have to be created and written to the file:

hdulist = pyfits.HDUList([hdu])
hdulist.writeto('new.fits')
hdulist.close()

Alternative way

Another way to deal with FITS tables is to use the ATpy library, look there:
http://atpy.github.com/index.html

References

To learn more on PyFITS: http://packages.python.org/pyfits/

Read an ascii file (cont')

In this message we'll see other ways to do the same job than in the previous post: read an ascii file.

In the previous post we read a file of the form:

0.000000 0.000000
0.095200 0.095056
0.190400 0.189251
0.285599 0.281733
0.380799 0.371662
0.475999 0.458227
0.571199 0.540641
0.666398 0.618159
0.761598 0.690079
0.856798 0.755750
...
using 2 lists, one for each column. It can be easier to manage only one 2D variable containing the whole table. Easier to apply filter on it, for example.
We first read the whole file into a single variable:

f2 = open('datas.dat', 'r')
lines = f2.readlines()
f2.close()

One can alternatively use the with formulation, which close the file automatically once used:

with open('datas.dat', 'r') as f2:

    lines = f2.readlines()

lines is not usable as it is, as it contains the lines of the file in string form, without separating the 2 values...

We now define data as a list, which will contain the lines of the file in the right format:
data = []

Now, we loop on the lines to copy the values, but after splitting the line:
for line in lines:
    p = line.split()
    data.append(p)

One can do the 4 previous commands in a single one, using list comprehension:
data = [line.split() for line in lines]

Let's have a look at the first element of the list data:

print(data[0])

['0.000000', '0.000000']

Ok, the value of x and y are not anymore in a single variable, but to really use them, we have to:

1. transform the list into a numpy array, to allow operations on the values.
2. transform this into floating point.

The 2 operations can be performed one after the other:

data2 = np.array(data)
data2 = data2.astype(float)

or in a single command (notice the f in asfarray, to transform into float):
data2 = np.asfarray(data)

The numpy arrays have a lot of properties and method. One tells us the shape of the array:
print(data2.shape)
(100, 2)

To access the distinct columns, use:

x = data2[:,0]

y = data2[:,1]

As the array is a numpy object, one can perform operations on it:

z = (data2[:,0]**2 + data2[:,1]**2)**0.5

Numpy includes functions to read ascii files and return float arrays directly. The following is then very compact and efficient:
data3 = np.loadtxt('datas.dat')


An even more complete tool to read ascii file is the following:
data4 = np.genfromtxt('datas.dat')
It allows to read formated files, deal with undefined values, etc...
One can even ask for this function to automatically put the result in 2 different variables:

x, y = np.loadtxt('datas.dat', unpack = True)

x, y = np.genfromtxt('datas.dat', unpack = True)

Once you have the data2, data3, or data4, you can plot it (we'll see more on plots latter):
plt.plot(data2[:,0], data2[:,1])

A small comment on the order of the elements in arrays in Python:
There is two ways arrays can be stored: row- or column major. It has a direct impact on the way one has to loop on the arrays. IDL is like Fortran (column major) and Python is like C (row major). It means that in Python, as you move linearly through the memory of an array, the second dimension (rightmost) changes the fastest, while in IDL the first (leftmost) dimension changes the fastest.
Consequence on the loop order in Python:

for i in range(100):
   for j in range(2):

         ... data2[i,j] ...

Another comment on boundaries in arrays:

print(data2[0:4, 0])

[ 0.        0.0952    0.1904    0.285599]
Some would expect the answer to be a 5 elements sub-array. The last index of 0:4 is not included. The rule is that a[n:m] return the elements between n and m-1.

Something quite powerful: accessing the latest elements of an array:

print(data2[-1, :])

[ 9.424778  0.      ]

print(data2[-5:-1, :])

[[ 9.043979  0.371662]
 [ 9.139179  0.281733]
 [ 9.234378  0.189251]
 [ 9.329578  0.095056]]

The np.genfromtxt method allows to read formatted column tables, as for example the following file:
#  Line      Iobs    lambda  relat_error Obs_code
H  1  4861A 1.00000    4861. 0.08000  Anabel                               
H  1  6563A 2.8667     6563. 0.19467  Anabel                               
H  1  4340A 0.4933     4340. 0.03307  Anabel                               
H  1  4102A 0.2907     4102. 0.02229  Anabel                               
H  1  3970A 0.1800     3970. 0.01253  Anabel                               
N  2  6584A 2.1681     6584. 0.08686  Anabel                               
N  2 121.7m 0.00446 1217000. 0.20000  Liu                                  
O  1  6300A 0.0147     6300. 0.00325  Anabel                               
TOTL  2326A 0.07900    2326. 0.20000  Adams                                
C  2 157.6m 0.00856 1576000. 0.20000  Liu                                  
O  1 63.17m 0.13647  631700. 0.10000  Liu                                  
O  1 145.5m 0.00446 1455000. 0.200    Liu                                  
TOTL  3727A 0.77609    3727. 0.200    Torres-Peimbert                      
S II  4070A 0.06174    4070. 0.200    Torres-Peimbert                      
S II  4078A 0.06174    4078. 0.200    Torres-Peimbert                      

read using this command:
obs  = np.genfromtxt('observations.dat', skip_header=1,
                     dtype=["a11","float","float","float","a2","a2"],
                     delimiter=[11,8,9,8,2,2],
                     names = ['label', 'i_obs', 'lambda', 'e_obs', 'na', 'observer'],
                     usecols = (0,1,2,3,5)
                     )
The skip_header tells not to read the first line. The dtype describes the type of the 6 columns. delimiter is the list of the sizes (in characters) of each column. The names are used to identify the columns. The usecols is used to specify which columns to output (here I don't need the 5th column, made of 2 spaces before the observer name)
The data are accessible using for example:
obs['i_obs']


To learn more in arrays indexing:

• http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html

To learn more about list comprehension:

• http://docs.python.org/tutorial/datastructures.html#list-comprehensions
• http://www.network-theory.co.uk/docs/pytut/ListComprehensions.html
• http://www.secnetix.de/olli/Python/list_comprehensions.hawk
  

More about np.loadtxt, np.genfromtxt:

• http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html
• http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html