Python Tutorial

In association with Lab #2 of Phys 322, Observational Astronomy

 

Start Python:

To start Python, click on the Jupyter QTConsole icon in your Anaconda Python start menu. After a moment, a window will appear, ready for typing commands. We will be using some special libraries for data analysis and plotting, which can be loaded by typing the two commands

%pylab
from astropy.io import fits

at the command prompt.  It might take a 10 s for each of these environments to load.

Open a blank display window:

To open a blank "figure" window, simply type

figure()

at the command prompt.  You should see a blank gray window with the title Figure 1.  This is where plot or image output will go.

 

Read in an image:

We have an image of Saturn in the data repository, http://web.njit.edu/~gary/322/download. On going there, you will see near the bottom a file Saturn.fit. Right click and save it to disk. Change python's current directory setting (using the cd command) to the directory where the image was saved. The image is in a format called FITS (Flexible Image Transfer System), which is commonly used in astronomy.  Type in the following commands to change directory and read the image using the readfits procedure (change the directory name to wherever you have put the file on your system):

cd c:/users/gary/documents/
satraw = fits.open('saturn.fit')[0].data

where the argument is just the filename of the file saturn.fit.  The variable name satraw will be used to refer to the image, and we can think of it as containing the image.

 

Display the image:

To display the image, type the command:

imshow(satraw,cmap='gray')

You should see a gray-scale image of Saturn in the figure window.  The command imshow does two things—it scales the values in the image so that they cover the full 8-bit range of the display, but no more, and then it transfers the values to the figure window for display.

Use the python data cursor to find the location of Saturn:

As you move the mouse in the figure window, you will see numbers appear in the status bar at the bottom of the window showing the x and y positions of the mouse, and the intensity (in square brackets). Point as close to the center of Saturn as possible and write down the x and y numbers--you should get something like x=444, y=196. These are the pixel numbers of the center of Saturn in the image. Let's save these values in variables x and y with the following commands:

x = 444

y = 196

Extract a small portion of the image around Saturn:

You can refer to parts of an image (horizontal and vertical ranges of pixels) by using the range expression:  min:max, so that 0:10 would refer to 10 pixel range from 0 to 9 (one less than the number given).  We want to clip out a 40 by 40 pixel area around the position of Saturn and put it into a new variable called saturn.  We use the following command to do this:

saturn = satraw[y-19:y+21,x-19:x+21]

which illustrates one trick you need to know. Although you may be used to specifying coordinates as (x, y) (column, row), python expects to address array variables in the opposite order, i.e. (row, column), which means (y, x). Aside from this, you should be able to understand the above command as taking a range from columns 425:464 in the x direction and rows 177:216 in the y direction.  The geometry of the situation is shown in the figure at the right.  To verify that the small sub-image is 40 pixels by 40 pixels, use the size command, which just prints the size of arrays.  For example

saturn.shape

prints the following to the log window

(40L, 40L)

indicating that the variable saturn is of size 40 x 40 pixels.

Display the small image:

To display the small image represented by the variable saturn, just use the imagesc command as before:

imshow(saturn,cmap='gray')

which this time will show the smaller image zoomed to fill the display window.

 

Print the contents of the variable saturn:

To demonstrate that images are just arrays of numbers, print out the values in the variable saturn by typing the command

saturn

which will print some of the numbers maknig up the image.  You may notice that the numbers shown (which represent a few pixels in the corners of the image) are around 120. To find out the range of numbers over the whole image, we can use the max and min commands.  The comand

saturn.min(),saturn.max()

prints

(107, 12546)

We can display the image in numbers, and fit it all one the screen, if we scale the numbers to range from 0 to 10.  The following command does this:

for i in range(40): print (('{:2d} '*40).format(*((saturn[i]*10./saturn.max()).astype(int))))

(If you are using Python 2.7, you'll have to strip off the outer set of parens () from the above command). You should be able to see the image of Saturn in the numbers printed to the terminal window. Enlarge the window until you can see 40 numbers on each line. You can actually make out the shape of Saturn in the numbers! To understand this command, note that saturn/saturn.max() scales the numbers from 0 to 1, so we multiply these numbers by 10 to scale from 0 to 10.  The astype(int) function just converts these real numbers to integers.  The print function prints a row of 40 numbers according to format 2d, which is to say, 2 digits for each number. Since we are printing saturn[i], this is 40 numbers in each row. And finally, to print all 40 rows, we need a for loop. 

 

Display the image of Saturn in different colors

Once we have an image contained in a variable, python can display it in a nearly infinite number of ways.  Try each of the commands below, looking at the plot after each:

imshow(saturn)                   # display the image in its default color table
imshow(saturn, cmap='Spectral')  # change to "Spectral" color table
imshow(saturn, cmap='copper')    # change to "copper" color table

You can force an error to get a list of possible colormaps, by specifying an invalid colormap, e.g. imshow(saturn,' ').

 

Save (and print) a color photo of Saturn

To save an image, you simply arrange the figure window plot as you like, and then click the disk icon at the top of the plot window. After saving, you can then print the figure using a third-party program. You can add a title, axis labels, etc., e.g. title('My plot'), xlabel('RA'), ylabel('Declination') before printing.

 

A quick look at noise

Reset the display to show the satraw image in gray color table:

imshow(satraw, cmap='gray')      # Redisplay the satraw image

Although the background appears completely black, this is only because the brightest parts of the image are being scaled to the maximum of the display, which is a level of 256 (28 or 8-bit).  We can see the noise level in the image by clipping the image to a level just above the noise, say a value of 150, as follows:

imshow(satraw,cmap='gray',vmin=50,vmax=150)    # display satraw, but clip to range 50 to 150

Suddenly we can see the noisy background, and we can even see at least three moons of Saturn (zoom in to see them well). 

 

Let’s take a closer look at the noise.  Clip out the bottom section of the satraw image (the first 50 rows of the image) into a new image called noise:

noise = satraw[0:50,:]    # clip out the bottom 50 rows of satraw

where we have used a short-cut—the ‘:’ stands for the entire x range of the image, and is the same as if we had typed noise = satraw[0:50,0:765].  We can use even more of a short-cut in this case, because noise = satraw[0:50]would give the same result. Let’s plot a single row of the noise image in a new figure:

figure()       # Make a new blank figure
plot(noise[0]) # plot the bottom row of the image

The plot command just plots the array of numbers as a line plot, but we can think of it as a cut across the image.  You can see how noisy the background is.  Try plotting some other rows by changing the zero to other numbers.

 

A quantitative look at the noise

We can think of the noise as a background light level of about 120, plus random fluctuations from pixel to pixel of the CCD camera that add or subtract from the background.  How “noisy” is it?  We can quantify the noise by using the moment function as follows:

mean(noise), std(noise)    % print the mean and standard deviation

which results in two numbers being printed to the screen,

(113.086405, 5.33561)

The first number, about 113, is the mean value over all 50 lines of the noise array.  The second number is the standard deviation, s = 5.37.  This means that if the noise is distributed according to gaussian statistics, there is about a 63% chance that a particular value will be within s = 5.37 units of the mean, a better than 95% chance that it will be within 2s = 10.7 units of the mean, and a better than 99% probability that it will be within 3s = 16 units.

 

We can reduce the noise by averaging.  For example, we can average all of the rows of the noise image into a single average row (which we will name avg) by using the mean function slightly differently:

avg_noise = mean(noise,0)      % Sum over the first dimension (50 lines) of noise

which sums the first dimension (the rows) of the noise array and divides by the number of rows to obtain the average of the rows.  We could average over the columns (the second dimension) of noise by the alternative command avg2 = mean(noise,1). Let’s overplot the average of the rows onto our previous plot:

plot(avg_noise,'r')     # overplot, using a red line

You should see a red line plotted over the previous plot, which clearly has greatly reduced noise.  If we average 50 numbers, the standard deviation should drop by a factor of =7.1, to 5.37 / 7.1 = 0.76, but if we apply the mean function to our average array

std(avg_noise)

we find a standard deviation of about 2.9. That is because the background is not flat, but rather has a slope, gently increasing from left to right. This skews our measurement of the noise. To get a valid noise measurement, we should “flatten” the image to remove these slow variations across the image. The commands below demonstrate the effect of flattening on the standard deviation. The first line of code fits a straight line to the data and returns the parameters of the fit. The second line of code creates the y-values of the fitted line from the parameters. We then plot the single line of the noise array and calculate its standard deviation, and then do the same with the average noise. Now the standard deviation is 0.71, which compares well with the expected value of 0.76 above.

p = polyfit(arange(700),avg_noise[:700],1)
fit = polyval(p,arange(700))
plot(arange(700),noise[0,:700] - fit)
std(noise[0,:700] - fit)
> 4.3258231247442831
plot(arange(700),avg_noise[:700] - fit)
std(avg_noise[:700] - fit)
> 0.7102655017970132

 

Final comments

This tutorial has introduced some of the more useful commands in python, although there are many more that we could use, and the commands that we did use have other features that we have not explored.  To learn more about a command, type the command followed by '?', e.g. plot? to get on-line help.

 

In this tutorial you should have learned that images are just arrays of numbers, and mathematical manipulations are used to correct for different unwanted effects, as well as bringing out details that may not be apparent in the raw image (like the presence of the two moons of Saturn).  We will be using python to explore and analyze images for scientific measurements, as well as improving the images for aesthetic appeal.

 

Please feel free to explore python with the Saturn image or others, and indulge your curiosity.