Lecturer: Leo Singer
Jupyter Notebook Author: Leo Singer & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Summer School 2019. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-school-2019.html
Introduce the user to basic usage of Python. Includes some basic analysis of photometric data using astropy.
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with pip install <module>
. The external astromatic packages are easiest installed using package managers (e.g., rpm
, apt-get
).
None
This workshop is about doing astronomical data analysis with the Python programming language. No previous experience with Python is necessary!
Python is a powerful tool, but it comes into its own as a numerical and data analysis environment with the following packages, which you will definitely want to have:
We'll cover the basics of Python itself and then dive in to some applications to explore each of these packages.
NOTE: The purest way of interacting with Python is via its command line interpreter, which looks like this:
A relatively new but very powerful way of using Python is through the Jupyter Notebook interface, which like Mathematica allows you to intermingle computer code with generated plots. You're using one now.
from matplotlib import pyplot as plt
import numpy as np
%matplotlib inline
x = np.linspace(0, 2 * np.pi)
plt.plot(x, np.sin(x))
plt.xlabel('ppm caffeine in bloodstream')
plt.ylabel('cheeriness')
Text(0,0.5,'cheeriness')
and tables...
import astropy.table
tbl = astropy.table.Table()
tbl.add_column(astropy.table.Column(data=np.arange(5),
name='integers'))
tbl.add_column(astropy.table.Column(data=np.arange(5)**2,
name='their squares'))
tbl
integers | their squares |
---|---|
int64 | int64 |
0 | 0 |
1 | 1 |
2 | 4 |
3 | 9 |
4 | 16 |
and even notes and typeset mathematics...
And God said:
$$\nabla \cdot \mathbf{D} = \rho$$$$\nabla \cdot \mathbf{B} = 0$$$$\nabla \times \mathbf{E} = -\frac{\partial\mathbf{B}}{\partial t}$$$$\nabla \times \mathbf{H} = J + \frac{\partial\mathbf{D}}{\partial t}$$and there was light.
This is all very useful for doing interactive data analysis, so we will use the IPython Notebook interface for this tutorial. WARNING: I'm spoiling you rotten by doing this.
Python and all of the packages that we discuss in this tutorial are open source software, so there multiple options for installing them. But if you followed the instructions on the GROWTH website for downloading/installing these modules, you have already installed these dependencies. Skip to the next step.
If you have one of the common Linux/UNIX distros (for example, Ubuntu, Debian, or Fedora), then you probably already have Python and you can get Matplotlib and friends from your package manager.
For example, on Debian or Ubuntu, use:
$ sudo apt-get install jupyter-notebook python3-matplotlib python3-astropy python3-scipy
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 py37-matplotlib py37-scipy py37-jupyterlab py37-astropy
Windows does not come with Python, but popular and free builds of Python for Windows include Anaconda and Canopy. Another alternative for Windows is to set up a virtual machine with VirtualBox and then install a Linux distribution on that.
print()
function and string literals¶If this is your first time looking at Python code, the first thing that you might notice is that it is very easy to understand. For example, to print something to the screen, it's just:
print('Hello world!')
Hello world!
This is a Python statement, consisting of the built-in command print
and a string surrounded by single quotes. Double quotes are fine inside a string:
print('She said, "Hello, world!"')
She said, "Hello, world!"
But if you want single quotes inside your string, you had better delimit it with double quotes:
print("She said, 'Hello, world!'")
She said, 'Hello, world!'
If you need both single quotes and double quotes, you can use backslashes to escape characters.
print('She said, "O brave new world, that has such people in\'t!"')
She said, "O brave new world, that has such people in't!"
If you need a string that contains newlines, use triple quotes ('''
) or triple double quotes ("""
):
print("""MIRANDA
O, wonder!
How many goodly creatures are there here!
How beauteous mankind is! O brave new world
That has such people in't!""")
MIRANDA O, wonder! How many goodly creatures are there here! How beauteous mankind is! O brave new world That has such people in't!
Let's say that you need to print a few different things on the same line. Just separate them with commas, as in:
person = 'Miranda'
print("'Tis new to", person)
'Tis new to Miranda
Oops. I'm getting ahead of myself—you've now seen your first variable assignment in Python. Strings can be concatened by adding them:
'abc' + 'def'
'abcdef'
Or repeated by multiplying them:
'abcdef' * 2
'abcdefabcdef'
Python's numeric types include integers and both real and complex floating point numbers:
a = 30 # an integer
b = 0xDEADBEEF # an integer in hexadecimal
c = 3.14159 # a floating point number
d = 5.1e10 # scientific notation
e = 2.5 + 5.3j # a complex number
hungry = True # boolean literal
need_coffee = False # another boolean literal
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):
a + c
33.14159
Multiplication:
a * e
(75+159j)
Division:
a / c
9.549304651466295
*Important note: unlike C, C++, Java, etc., division of integers gives you floats*:
7 / 3
2.3333333333333335
If you want integer division, then use the double-slash //
operator:
a = 7
b = 3
7 // 3
2
The %
sign is the remainder operator:
32 % 26
6
Exponentiation is accomplished with the **
operator:
print(5 ** 3, 9**-0.5)
125 0.3333333333333333
A tuple is a sequence of values. It's just about the handiest thing since integers. A tuple is immutable: once you have created it, you cannot add items to it, remove items from it, or change items. Tuples are very handy for storing short sequences of related values or returning multiple values from a function. This is what tuples look like:
some_tuple = ('a', 'b', 'c')
another_tuple = ('caffeine', 6.674e-11, 3.14, 2.718)
nested_tuple = (5, 4, 3, 2, ('a', 'b'), 'c')
Once you have made a tuple, you might want to retrieve a value from it. You index a tuple with square brackets, *starting from zero*:
some_tuple[0]
'a'
some_tuple[1]
'b'
You can access whole ranges of values using *slice notation*:
nested_tuple[1:4]
(4, 3, 2)
Or, to count backward from the end of the tuple, use a *negative index*:
another_tuple[-1]
2.718
another_tuple[-2]
3.14
Strings can be treated just like tuples of individual charaters:
person = 'Miranda'
print(person[3:6])
and
What if you want a container like a tuple but to which you can add or remove items or alter existing items? That's a list. The syntax is almost the same, except that you create a list using square brackets []
instead of round ones ()
:
your_list = ['foo', 'bar', 'bat', 'baz']
my_list = ['xyzzy', 1, 3, 5, 7]
But you can change elements:
my_list[1] = 2
print(my_list)
['xyzzy', 2, 3, 5, 7]
Or append elements to an existing list:
my_list.append(11)
print(my_list)
['xyzzy', 2, 3, 5, 7, 11]
Or delete elements:
del my_list[0]
print(my_list)
[2, 3, 5, 7, 11]
Sometimes you need a collection of items where order doesn't necessarily matter, but each item is guaranteed to be unique. That's a set, created just like a list or tuple but with curly braces {}
:
a = {5, 6, 'foo', 7, 7, 8}
print(a)
{5, 6, 7, 8, 'foo'}
You can add items to a set:
a.add(3)
print(a)
{3, 5, 6, 7, 8, 'foo'}
Or take them away:
a.remove(3)
print(a)
{5, 6, 7, 8, 'foo'}
You also have set-theoretic intersections with the &
operator:
{1, 2, 3, 4, 5, 6} & {3, 4}
{3, 4}
And union with the |
operator:
{1, 2, 3, 4, 5, 6} | {6, 7}
{1, 2, 3, 4, 5, 6, 7}
And set difference with the -
operator:
{1, 2, 3, 4, 5, 6} - {3, 4}
{1, 2, 5, 6}
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:
my_grb = {'name': 'GRB 130702A', 'redshift': 0.145, 'ra': (14, 29, 14.78), 'dec': (15, 46, 26.4)}
my_grb
{'name': 'GRB 130702A', 'redshift': 0.145, 'ra': (14, 29, 14.78), 'dec': (15, 46, 26.4)}
You can index items in dictionaries with square braces []
, similar to tuples or lists:
my_grb['dec']
(15, 46, 26.4)
or add items to them:
my_grb['url'] = 'http://gcn.gsfc.nasa.gov/other/130702A.gcn3'
my_grb
{'name': 'GRB 130702A', 'redshift': 0.145, 'ra': (14, 29, 14.78), 'dec': (15, 46, 26.4), 'url': 'http://gcn.gsfc.nasa.gov/other/130702A.gcn3'}
or delete items from them:
del my_grb['url']
my_grb
{'name': 'GRB 130702A', 'redshift': 0.145, 'ra': (14, 29, 14.78), 'dec': (15, 46, 26.4)}
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:
{
'foods': ['chicken', 'veggie burger', 'banana'],
'cheeses': {'muenster', 'gouda', 'camembert', 'mozarella'},
(5.5, 2): 42,
'plugh': 'bat'
}
{'foods': ['chicken', 'veggie burger', 'banana'], 'cheeses': {'camembert', 'gouda', 'mozarella', 'muenster'}, (5.5, 2): 42, 'plugh': 'bat'}
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
.
a = None
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. Common choices include four spaces, two spaces, or a tab.
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:
a = 5
if a < 3:
print("i'm in the 'if' block")
messsage = 'a is less than 3'
elif a == 3:
print("i'm in the 'elif' block")
messsage = 'a is 3'
else:
print("i'm in the 'else' block")
message = 'a is greater than 3'
print(message)
i'm in the 'else' block a is greater than 3
You can chain together inequalities just like in mathematical notation:
if 0 < a <= 5:
print('a is greater than 0 but less than or equal to 5')
a is greater than 0 but less than or equal to 5
You can also combine comparison operators with the boolean and
, or
, and not
operators:
if a < 6 or a > 8:
print('yahoo!')
yahoo!
if a < 6 and a % 2 == 1:
print('a is an odd number less than 6!')
a is an odd number less than 6!
if not a == 5: # same as a != 5
print('a is not 5')
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
.
food = None
if food is None:
print('No, thanks')
else:
print('Here is your', food)
No, thanks
Likewise, there is an is not
operator:
if food is not None:
print('Yum!')
The in
and not in
operators are handy for testing for membership in a string, set, or dictionary:
if 3 in {1, 2, 3, 4, 5}:
print('indeed it is')
indeed it is
if 'i' not in 'team':
print('there is no "i" in "team"')
there is no "i" in "team"
When referring to a dictionary, the in
operator tests if the item is among the keys of the dictionary.
d = {'foo': 3, 'bar': 5, 'bat': 9}
if 'foo' in d:
print('the key "foo" is in the dictionary')
the key "foo" is in the dictionary
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.
for i in [0, 1, 2, 3]:
print(i**2)
0 1 4 9
The built-in function range
, which returns a list of numbers, is often handy here:
for i in range(4):
print(i**2)
0 1 4 9
Or you can have the range start from a nonzero value:
for i in range(-2, 4):
print(i**2)
4 1 0 1 4 9
You can iterate over the keys and values in a dictionary with .items()
:
for key, val in d.items():
print(key, '...', val**3)
foo ... 27 bar ... 125 bat ... 729
The syntax of the while
loop is similar to the if
statement:
a = 1
while a < 5:
a = a * 2
print(a)
2 4 8
Sometimes you need a loop to create one list from another. List comprehensions make this very terse. For example, the following for
loop:
a = []
for i in range(5):
a.append(i * 10)
is equivalent to this list comprehension:
a = [i * 10 for i in range(5)]
You can even incorporate conditionals into a list comprehension. The following:
a = []
for i in range(5):
if i % 2 == 0:
# i is even
a.append(i * 10)
can be written as:
a = [i * 10 for i in range(5) if i % 2 == 0]
Conditional expressions are a closely related shorthand. The following:
if 6/2 == 3:
a = 'foo'
else:
a = 'bar'
is equivalent to:
a = 'foo' if 6/2 == 3 else 'bar'
Functions are created with the def
statement. A function may either have or not have a return
statement to send back a return value.
def square(n):
return n * n
a = square(3)
print(a)
9
If you want to return multiple values from a function, return a tuple. Parentheses around the tuple are optional.
def powers(n):
return n**2, n**3
print(powers(3))
(9, 27)
If a function returns multiple values, you can automatically unpack them into multiple variables:
square, cube = powers(3)
print(square)
9
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:
def fibonacci(seed, n):
while len(seed) < n:
seed.append(seed[-1] + seed[-2])
# Note: no return statement
seed = [1, 1]
fibonacci(seed, 10)
print(seed)
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]
You can also give a function's arguments default values, such as:
def fibonacci(seed, n=6):
while len(seed) < n:
seed.append(seed[-1] + seed[-2])
# Note: no return statement
seed = [1, 1]
fibonacci(seed)
print(seed)
[1, 1, 2, 3, 5, 8]
If a function has a large number of arguments, it may be easier to read if you pass the arguments by keyword, as in:
seq = [1, 1]
fibonacci(seed=seq, n=4)
Python comes with an extensive standard library consisting of individual modules that you can opt to use with the import
statement. For example:
import math
math.sqrt(3)
1.7320508075688772
from math import pi
pi
3.141592653589793
Some particularly useful parts of the Python standard library are:
random
: random number generatorspickle
: read/write Python objects into filessqlite3
: SQLite database accesos
: operating system servicesos.path
: file path manipulationsubprocess
: launch external processesemail
: compose, parse, receive, or send e-mailpdb
: built-in debuggerre
: regular expressionshttp
: built-in lightweight web client and serveroptparse
: build pretty command-line interfacesitertools
: exotic looping constructsmultiprocessing
: parallel processingIt 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:
def sinc(x):
return math.sin(x) / x
print(sinc(0))
--------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) <ipython-input-71-e7a6fa28a489> in <module> 2 return math.sin(x) / x 3 ----> 4 print(sinc(0)) <ipython-input-71-e7a6fa28a489> in sinc(x) 1 def sinc(x): ----> 2 return math.sin(x) / x 3 4 print(sinc(0)) ZeroDivisionError: float division by zero
Oops! We know that by definition $\mathrm{sinc}(0) = 1$ , so we should catch this error:
def sinc(x):
try:
result = math.sin(x) / x
except ZeroDivisionError:
result = 1
return result
print(sinc(0))
1
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:
myfile = open('myfile.txt', 'w') # open file for writing
myfile.write("red 1\n")
myfile.write("green 2\n")
myfile.write("blue 3\n")
myfile.close()
And here is reading it:
d = {} # create empty dictionary
for line in open('myfile.txt', 'r'): # open file for reading
color, num = line.split() # break apart line by whitespace
num = int(num) # convert num to integer
d[color] = num
print(d)
{'red': 1, 'green': 2, 'blue': 3}
Numpy provides array operations and linear algebra to Python. A Numpy array is a bit like a Python list, but supports elementwise arithmetic. For example:
import numpy as np
x = np.asarray([1, 2, 3, 4, 5])
y = 2 * x
print(y)
[ 2 4 6 8 10]
Numpy arrays may have any number of dimensions:
x = np.asarray([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
x
array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
y = np.asarray([[9, 8, 7], [6, 5, 4], [3, 2, 1]])
y
array([[9, 8, 7], [6, 5, 4], [3, 2, 1]])
An array has a certain number of dimensions denoted .ndim
:
x.ndim
2
and the dimensions' individual lengths are given by .shape
:
x.shape
(3, 3)
and the total number of elements by .size
:
x.size
9
By default, multiplication is elementwise:
x * y
array([[ 9, 16, 21], [24, 25, 24], [21, 16, 9]])
To perform matrix multiplication, either convert arrays to np.matrix
or use np.dot
:
np.asmatrix(x) * np.asmatrix(y)
matrix([[ 30, 24, 18], [ 84, 69, 54], [138, 114, 90]])
np.dot(x, y)
array([[ 30, 24, 18], [ 84, 69, 54], [138, 114, 90]])
You can also perform comparison operations on arrays...
x > 5
array([[False, False, False], [False, False, True], [ True, True, True]])
Although a boolean array doesn't directly make sense in an if
statement:
if x > 5:
print('oops')
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-85-5a5904c20c71> in <module> ----> 1 if x > 5: 2 print('oops') ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
if np.any(x > 5):
print('at least some elements are greater than 5')
at least some elements are greater than 5
You can use conditional expressions like indices:
x[x > 5] = 5
x
array([[1, 2, 3], [4, 5, 5], [5, 5, 5]])
Or manipulate individual rows:
x[1, :] = -1
x
array([[ 1, 2, 3], [-1, -1, -1], [ 5, 5, 5]])
Or individual columns:
x[:, 1] += 100
x
array([[ 1, 102, 3], [ -1, 99, -1], [ 5, 105, 5]])
Other useful features include various random number generators:
from matplotlib import pyplot as plt
%matplotlib inline
# Plot histogram of 10k normal random variates
plt.hist(np.random.randn(10000))
(array([ 19., 107., 586., 1472., 2522., 2665., 1711., 710., 189., 19.]), array([-3.58212794, -2.87990784, -2.17768775, -1.47546765, -0.77324756, -0.07102746, 0.63119264, 1.33341273, 2.03563283, 2.73785292, 3.44007302]), <a list of 10 Patch objects>)
np.random.uniform(low=0, high=2*np.pi)
1.364157604885726
You've already seen a few examples of Matplotlib. If you have used MATLAB, then Matplotlib code may look familiar.
x = np.linspace(-10, 10)
y = 1 / (1 + np.exp(x))
plt.plot(x, y)
plt.annotate(
'foo bar', (x[20], y[20]), (50, 5),
textcoords='offset points',
arrowprops={'arrowstyle': '->'})
plt.grid()
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 conversionastropy.time
: Manipulation of dates and timesastropy.coordinates
: Representation of and conversion between astronomical coordinate systemsastropy.table
: Tables and gridded dataastropy.io.fits
: Manipulating FITS filesastropy.io.ascii
: Manipulating ASCII tables of many different formatsastropy.io.votable
: Virtual Observatory tablesastropy.wcs
: World Coordinate System transformationsastropy.cosmology
: Cosmological calculationsastropy.stats
: Astrostatisticsastropy.modeling
: multi-D model fitting Swiss army knifeThe Astropy project also has sevearl "affiliated packages" that have similar design but are maintained separately, including:
Let's experiment by opening up a P48 image. We'll need several modules from the Astropy package for this exercise.
import astropy.coordinates
import astropy.units as u
import astropy.io.fits
import astropy.stats
import astropy.table
import astropy.wcs
import astropy.cosmology
import scipy.optimize
import scipy.odr
I've downloaded a P48 image and put it in the data/ directory.
fits = astropy.io.fits.open('data/PTF_201307021787_i_p_scie_t041723_u016616794_f02_p003486_c11.fits')
fits
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7fc7b2e7b518>]
Let's grab the first (and only) HDU of this FITS file:
hdu = fits[0]
Then let's take a look at the contents of the header:
hdu.header
SIMPLE = T / Fits standard BITPIX = -32 / FOUR-BYTE SINGLE PRECISION FLOATING POINT NAXIS = 2 / STANDARD FITS FORMAT NAXIS1 = 2048 / STANDARD FITS FORMAT NAXIS2 = 4096 / STANDARD FITS FORMAT ORIGIN = 'Palomar Transient Factory' / Origin of these image data CREATOR = 'Infrared Processing and Analysis Center' / Creator of this FITS file TELESCOP= 'P48 ' / Name of telescope INSTRUME= 'PTF/MOSAIC' / Instrument name OBSERVER= 'KulkarniPTF' / Observer name and project CCDID = '11 ' / CCD number (0..11) DATE-OBS= '2013-07-02T04:17:23.555' / UTC shutter time YYYY-MM-DDTHH:MM:SS.SSS DATE = '2013-07-01T21:35:19' / File creation date (YYYY-MM-DDThh:mm:ss UT) REFERENC= 'http://www.astro.caltech.edu/ptf' / URL of PTF website / PROPOSAL INFORMATION PTFPRPI = 'Kulkarni' / PTF Project PI PTFPID = '30011 ' / Project type: 00000-49999 OBJECT = 'PTF_survey' / Fields object PTFFIELD= '3486 ' / PTF unique field ID PTFFLAG = '1 ' / 1 = PTF; 0 = non-PTF category / TIME AND EXPOSURE INFORMATION FILTER = 'R ' / Filter name FILTERID= '2 ' / Filter ID FILTERSL= '1 ' / Filter changer slot position EXPTIME = 60. / [s] Requested exposure time AEXPTIME= 60. / actual exposure time (sec) UTC-OBS = '2013-07-02T04:17:23.555' / UTC time shutter open YYYY-MM-DDTHH:MM:SS. OBSJD = 2456475.67874 / [day] Julian day corresponds to UTC-OBS OBSMJD = 56475.17874 / MJD corresponds to UTC-OBS (day) OBSLST = '15:11:26.46' / Mean LST corresponds to UTC-OBS 'HH:MM:SS.S' HOURANG = '0:48:14.61' / Mean HA (sHH:MM:SS.S) based on LMST at UTC-OBS HJD = 2456475.68019 / [day] Heliocentric Julian Day OBSTYPE = 'object ' / Image type (dark,science,bias,focus) IMGTYP = 'object ' / Image type (dark,science,bias,focus) / MOON AND SUN MOONRA = 30.98457 / [deg] Moon J2000.0 R.A. MOONDEC = 12.984756 / [deg] Moon J2000.0 Dec. MOONILLF= -0.306369 / [frac] Moon illuminated fraction MOONPHAS= 247.2158 / [deg] Moon phase angle MOONESB = -0. / Moon excess in sky brightness V-band MOONALT = -40.94289 / [deg] Moon altitude SUNAZ = 310.3019 / [deg] Sun azimuth SUNALT = -14.00398 / [deg] Sun altitude / PHOTOMETRY BUNIT = 'DN ' / Data number (analog-to-digital units or ADU) PHTCALEX= 1 / Was phot.-cal. module executed? PHTCALFL= 1 / Flag for image is photometric (0=N, 1=Y) PCALRMSE= 0.215766 / RMSE from (zeropoint, extinction) data fit IMAGEZPT= 23.72382 / Image magnitude zeropoint COLORTRM= -0.012228 / Image color term (g-r) ZPTSIGMA= 1.886777 / Robust dispersion of SEx-SDSS magnitudes IZPORIG = 'SDSS ' / Photometric-calibration origin ZPRULE = 'DIRECT ' / Photometric-calibration method MAGZPT = 23.77773 / Magnitude zeropoint at airmass=1 EXTINCT = -0.195707 / Extinction APSFILT = 'r ' / SDSS filter used in abs phot cal APSCOL = 'r-i ' / SDSS color used in abs phot cal APRMS = 0.05618604 / RMS in mag of final abs phot cal APBSRMS = 0.0344094 / RMS in mag of final abs phot cal for bright sta APNSTDI1= 177080 / Number of standard stars in first iteration APNSTDIF= 160107 / Number of standard stars in final iteration APCHI2 = 771280.24477198 / Chi2 of final abs phot cal APDOF = 160096. / Dof of chi2 of final abs phot cal APMEDJD = 2456475.80246722 / Median JD used in abs phot cal APPN01 = 'ZeroPoint' / Name of parameter abs phot cal 01 APPAR01 = 23.68493517 / Value of parameter abs phot cal 01 APPARE01= 0.0026937 / Error of parameter abs phot cal 01 APPN02 = 'ColorTerm' / Name of parameter abs phot cal 02 APPAR02 = -0.03174945 / Value of parameter abs phot cal 02 APPARE02= 0.00383537 / Error of parameter abs phot cal 02 APPN03 = 'AirMassTerm' / Name of parameter abs phot cal 03 APPAR03 = -0.28459543 / Value of parameter abs phot cal 03 APPARE03= 0.00226667 / Error of parameter abs phot cal 03 APPN04 = 'AirMassColorTerm' / Name of parameter abs phot cal 04 APPAR04 = 0.24261075 / Value of parameter abs phot cal 04 APPARE04= 0.00313212 / Error of parameter abs phot cal 04 APPN05 = 'TimeTerm' / Name of parameter abs phot cal 05 APPAR05 = 0.19911481 / Value of parameter abs phot cal 05 APPARE05= 0.00219937 / Error of parameter abs phot cal 05 APPN06 = 'Time2Term' / Name of parameter abs phot cal 06 APPAR06 = 1.65698453 / Value of parameter abs phot cal 06 APPARE06= 0.02684713 / Error of parameter abs phot cal 06 APPN07 = 'XTerm ' / Name of parameter abs phot cal 07 APPAR07 = 0.02703671 / Value of parameter abs phot cal 07 APPARE07= 0.00053429 / Error of parameter abs phot cal 07 APPN08 = 'YTerm ' / Name of parameter abs phot cal 08 APPAR08 = -0.01528394 / Value of parameter abs phot cal 08 APPARE08= 0.00131972 / Error of parameter abs phot cal 08 APPN09 = 'Y2Term ' / Name of parameter abs phot cal 09 APPAR09 = 0.00938432 / Value of parameter abs phot cal 09 APPARE09= 0.00208283 / Error of parameter abs phot cal 09 APPN10 = 'Y3Term ' / Name of parameter abs phot cal 10 APPAR10 = 0.03443254 / Value of parameter abs phot cal 10 APPARE10= 0.00825142 / Error of parameter abs phot cal 10 APPN11 = 'XYTerm ' / Name of parameter abs phot cal 11 APPAR11 = 0.01025838 / Value of parameter abs phot cal 11 APPARE11= 0.00188474 / Error of parameter abs phot cal 11 / ASTROMETRY WCSAXES = 2 / Number of axes in world coordinate system CRVAL1 = 217.309490953823 / [deg] RA of reference point CRVAL2 = 16.6581466832508 / [deg] DEC of reference point CRPIX1 = 1497.193 / [pix] Image reference point CRPIX2 = 609.3031 / [pix] Image reference point CTYPE1 = 'RA---TAN-SIP' / TAN (gnomic) projection + SIP distortions CTYPE2 = 'DEC--TAN-SIP' / TAN (gnomic) projection + SIP distortions CUNIT1 = 'deg ' / Image axis-1 celestial-coordinate units CUNIT2 = 'deg ' / Image axis-2 celestial-coordinate units CRTYPE1 = 'deg ' / Data units of CRVAL1 CRTYPE2 = 'deg ' / Data units of CRVAL2 CD1_1 = 0.000280377484830748 / Transformation matrix CD1_2 = -1.71606271442079E-06 CD2_1 = -1.61073877029747E-06 CD2_2 = -0.000280943199913645 OBJRA = '14:22:34.454' / Requested field J2000.0 Ra. OBJDEC = '+16:52:30.00' / Requested field J2000.0 Dec. OBJRAD = 215.64356 / [deg] Requested field RA (J2000.0) OBJDECD = 16.875 / [deg] Requested field Dec (J2000.0) PIXSCALE= 1.01 / [arcsec/pix] Pixel scale EQUINOX = 2000. / [yr] Equatorial coordinates definition LONPOLE = 180. LATPOLE = 0. / IMAGE QUALITY SEEING = 1.71 / [pix] Seeing FWHM PEAKDIST= 0.389898704794061 / [pix] Mean dist brightest pixel-centroid pixel ELLIP = 0.107 / Mean image ellipticity A/B ELLIPPA = -35.01 / [deg] Mean image ellipticity PA FBIAS = 689.3355 / [DN] Floating bias of the image SATURVAL= 53000. / [DN] Saturation value of the CCD array FWHMSEX = 1.85 / [arcsec] SExtractor SEEING estimate MDSKYMAG= 20.39761 / [mag/s-arcsec^2] Median sky obsolete MSMAPCZP= 20.05711 / [mag/s-arcsec^2] Median sky abs. phot. cal. LIMITMAG= 21.74218 / [mag/s-arcsec^2] Limiting magnitude obsolete LMGAPCZP= 21.40168 / [mag/s-arcsec^2] Limiting mag. abs. phot. cal. MEDFWHM = 2.608781 / [arcsecond] Median FWHM MEDELONG= 1.168383 / [dimensionless] Median elongation STDELONG= 0.6377625 / [dimensionless] Std. dev. of elongation MEDTHETA= -7.846824 / [deg] Atan(median sin(theta)/median cos(theta)) STDTHETA= 64.7667 / [deg] Atan(stddev sin(theta)/stddev cos(theta)) MEDDLMAG= 30.16257 / [mag/s-arcsec^2] Median (MU_MAX-MAG_AUTO) STDDLMAG= 1.823812 / [mag/s-arcsec^2] Stddev of (MU_MAX-MAG_AUTO) / OBSERVATORY AND TCS OCS_TIME= '2013-07-02T04:17:23.461' / UTC Date for OCS calc time-dep params OPERMODE= 'OCS ' / Mode of operation: OCS | Manual | N/A SOFTVER = '1.1.1.1 ' / Softwere version (TCS.Camera.OCS.Sched) OCS_VER = '1 ' / OCS software version and date TCS_VER = '1 ' / TCS software version and date SCH_VER = '1 ' / OCS-Scheduler software version and date MAT_VER = '7.7.0.471' / Matlab version HDR_VER = '1 ' / Header version TRIGGER = 'N/A ' / trigger ID for TOO, e.g. VOEVENT-Nr TCSMODE = 'Star ' / TCS fundamental mode TCSSMODE= 'Active ' / TCS fundamental submode TCSFMODE= 'Pos ' / TCS focus mode TCSFSMOD= 'On-Target' / TCS focus submode TCSDMODE= 'Stop ' / TCS dome mode TCSDSMOD= 'N/A ' / TCS dome submode TCSWMODE= 'Slave ' / TCS windscreen mode TCSWSMOD= 'N/A ' / TCS windscreen submode OBSLAT = 33.3574 / [deg] Telescope geodetic latitude in WGS84 OBSLON = -116.8599 / [deg] Telescope geodetic longitude in WGS84 OBSALT = 1703.2 / [m] Telescope geodetic altitude in WGS84 DEFOCUS = 0. / [mm] Focus position - nominal focus FOCUSPOS= 1.3785 / [mm] Exposures focusPos DOMESTAT= 'open ' / Dome status at begining of exposure TRACKRA = 10.6 / [arcsec/hr] Track speed RA rel to sidereal TRACKDEC= 1.3 / [arcsec/hr] Track speed Dec rel to sidereal AZIMUTH = 216.6835 / [deg] Telescope Azimuth ALTITUDE= 70.18806 / [deg] Telescope altitude AIRMASS = 1.062788 / Telescope airmass TELRA = 215.7994 / [deg] Telescope ap equinox of date RA TELDEC = 16.8139 / [deg] Telescope ap equinox of date Dec TELHA = 12.0595 / [deg] Telescope ap equinox of date HA DOMEAZ = 215.3933 / [deg] Dome azimuth WINDSCAL= 10.1089 / [deg] Wind screen altitude WINDDIR = 1.6 / [deg] Azimuth of wind direction WINDSPED= 9.1656 / Wind speed (km/hour) OUTTEMP = 23.33333 / [C] Outside temperature OUTRELHU= 0.409 / [frac] Outside relative humidity OUTDEWPT= 9.277778 / [C] Outside dew point / INSTRUMENT TELEMETRY PANID = '_p48s ' / PAN identification DHSID = '_p48s ' / DHS identification CCDSEC = '[1:2048,1:4096]' / CCD section CCDSIZE = '[1:2048,1:4096]' / CCD size DATASEC = '[1:2048,1:4096]' / Data section DETSEC = '[1:2048,1:4096]' / Detector section ROISEC = '[1:2048,1:4096]' / ROI section FPA = 'P48MOSAIC' / Focal plan array CCDNAME = 'W7C1 ' / Detector mfg serial number CHECKSUM= ' ' / Image header unit checksum DATASUM = ' ' / Image data unit checksum DHEINF = 'SDSU, Gen-III' / Controller info DHEFIRM = '/usr/src/dsp/20090618/tim_m.lod' / DSP software CAM_VER = '20090615.1.3.100000' / Camera server date.rev.cfitsio LV_VER = '8.5 ' / LabVIEW software version PCI_VER = '2.0c ' / Astropci software version DETID = 'PTF/MOSAIC' / Detector ID AUTHOR = 'PTF/OCS/TCS/Camera' / Source for header information DATAMIN = 0. / Minimum value for array ROISTATE= 'ROI ' / ROI State (FULL | ROI) LEDBLUE = 'OFF ' / 470nm LED state (ON | OFF) LEDRED = 'OFF ' / 660nm LED state (ON | OFF) LEDNIR = 'OFF ' / 880nm LED state (ON | OFF) CCD9TEMP= 175.002 / [K] 0x0 servo temp sensor on CCD09 HSTEMP = 150.25 / [K] 0x1 heat spreader temp DHE0TEMP= 303.015 / [K] 0x2 detector head electronics temp, master DHE1TEMP= 304.935 / [K] 0x3 detector head electronics temp, slave DEWWTEMP= 292.787 / [K] 0x4 dewar wall temp HEADTEMP= 141.038 / [K] 0x5 cryo cooler cold head temp CCD5TEMP= 175.312 / [K] 0x6 temp sensor on CCD05 CCD11TEM= 176.263 / [K] 0x7 temp sensor on CCD11 CCD0TEMP= 169.376 / [K] 0x8 temp sensor on CCD00 RSTEMP = 241.756 / [K] 0x9 temp sensor on radiation shield DEWPRESS= 2.3 / [milli-torr] Dewar pressure DETHEAT = 28.5 / [%] Detector focal plane heater power NAMPSXY = '6 2 ' / Number of amplifiers in x y CCDSUM = '1 1 ' / [pix] Binning in x and y MODELFOC= 'N/A ' / MODELFOC EXPCKSUM= '4VHn7S9l4SEl4S9l' / Primary header unit checksum EXPDTSUM= ' 0' / Primary data unit checksum GAIN = 1.5 / [e-/D.N.] Gain of detector. READNOI = 5.2 / [e-] Read noise of detector. DARKCUR = 0.1 / [e-/s] Dark current of detector / SCAMP DISTORTION KEYWORDS RADECSYS= 'ICRS ' / Astrometric system PV1_0 = 0. / Projection distortion parameter PV1_1 = 1. / Projection distortion parameter PV1_2 = 0. / Projection distortion parameter PV1_4 = -0.00173789497161513 / Projection distortion parameter PV1_5 = 6.00083003595371E-05 / Projection distortion parameter PV1_6 = -0.00054585707358405 / Projection distortion parameter PV1_7 = -0.00133490210858131 / Projection distortion parameter PV1_8 = 0.000173741003278833 / Projection distortion parameter PV1_9 = -0.000547160675264968 / Projection distortion parameter PV1_10 = 0.000251433351740168 / Projection distortion parameter PV1_12 = -0.00328774077188989 / Projection distortion parameter PV1_13 = -0.000280710215945655 / Projection distortion parameter PV1_14 = 0.000714951620183291 / Projection distortion parameter PV1_15 = -0.000275634808657675 / Projection distortion parameter PV1_16 = 0.000192684172156633 / Projection distortion parameter PV2_0 = 0. / Projection distortion parameter PV2_1 = 1. / Projection distortion parameter PV2_2 = 0. / Projection distortion parameter PV2_4 = 0.00025061777603867 / Projection distortion parameter PV2_5 = -0.00125376481783182 / Projection distortion parameter PV2_6 = -0.000118547045320473 / Projection distortion parameter PV2_7 = -0.000709046387896288 / Projection distortion parameter PV2_8 = 0.000780095509401425 / Projection distortion parameter PV2_9 = -0.000924905030225836 / Projection distortion parameter PV2_10 = -0.000525520061118934 / Projection distortion parameter PV2_12 = -0.000283472622900895 / Projection distortion parameter PV2_13 = 0.000738353372557653 / Projection distortion parameter PV2_14 = 0.000164958172111789 / Projection distortion parameter PV2_15 = -0.00144743155797026 / Projection distortion parameter PV2_16 = -0.000422174086212961 / Projection distortion parameter FGROUPNO= 1 / SCAMP field group label ASTIRMS1= 0. / Astrom. dispersion RMS (intern., high S/N) ASTIRMS2= 0. / Astrom. dispersion RMS (intern., high S/N) ASTRRMS1= 2.60633E-05 / Astrom. dispersion RMS (ref., high S/N) ASTRRMS2= 2.679676E-05 / Astrom. dispersion RMS (ref., high S/N) ASTINST = 1 / SCAMP astrometric instrument label FLXSCALE= 0. / SCAMP relative flux scale MAGZEROP= 0. / SCAMP zero-point PHOTIRMS= 0. / mag dispersion RMS (internal, high S/N) RA_RMS = 0.4284936 / [arcsec] RMS of SCAMP fit from 2MASS matching DEC_RMS = 0.3289862 / [arcsec] RMS of SCAMP fit from 2MASS matching ASTROMN = 1003 / Number of stars in SCAMP astrometric solution SCAMPPTH= '/ptf/pos/archive/fallbackcal/scamp/11/' / SCAMP catalog path SCAMPFIL= 'PTF_201102044118_c_e_sdss_t095258_u003860515_f01_p003486_c11.fits' / SIP DISTORTION KEYWORDS A_ORDER = 4 / Distortion order for A A_0_2 = -1.53991686981733E-07 / Projection distortion parameter A_0_3 = -1.9961153583439E-11 / Projection distortion parameter A_0_4 = 4.28220501859733E-15 / Projection distortion parameter A_1_1 = -1.48140435828959E-08 / Projection distortion parameter A_1_2 = -4.37451084272745E-11 / Projection distortion parameter A_1_3 = 6.11763923987906E-15 / Projection distortion parameter A_2_0 = -4.87160559968282E-07 / Projection distortion parameter A_2_1 = -1.27084757090404E-11 / Projection distortion parameter A_2_2 = 1.57779409644097E-14 / Projection distortion parameter A_3_0 = -1.04765134994984E-10 / Projection distortion parameter A_3_1 = 7.95789862286398E-15 / Projection distortion parameter A_4_0 = -7.23704568938464E-14 / Projection distortion parameter A_DMAX = 2.44291521352919 / Projection distortion parameter B_ORDER = 4 / Distortion order for B B_0_2 = -6.7373693835656E-08 / Projection distortion parameter B_0_3 = -5.54765729401967E-11 / Projection distortion parameter B_0_4 = 6.16120047283346E-15 / Projection distortion parameter B_1_1 = -3.5264356557202E-07 / Projection distortion parameter B_1_2 = -6.12610617364653E-11 / Projection distortion parameter B_1_3 = 1.64880764414739E-14 / Projection distortion parameter B_2_0 = 3.39463080480267E-08 / Projection distortion parameter B_2_1 = -7.40973776817625E-11 / Projection distortion parameter B_2_2 = -2.8629870353495E-15 / Projection distortion parameter B_3_0 = 4.14105136847929E-11 / Projection distortion parameter B_3_1 = -3.22124747225646E-14 / Projection distortion parameter B_4_0 = 9.51825972208521E-15 / Projection distortion parameter B_DMAX = 2.95057287643024 / Projection distortion parameter AP_ORDER= 4 / Distortion order for AP AP_0_1 = 2.05105623663409E-08 / Projection distortion parameter AP_0_2 = 1.54047990787523E-07 / Projection distortion parameter AP_0_3 = 1.99841723199859E-11 / Projection distortion parameter AP_0_4 = -4.27260819246964E-15 / Projection distortion parameter AP_1_0 = -7.88989114580002E-08 / Projection distortion parameter AP_1_1 = 1.49024754213322E-08 / Projection distortion parameter AP_1_2 = 4.39773302986539E-11 / Projection distortion parameter AP_1_3 = -6.1078457561122E-15 / Projection distortion parameter AP_2_0 = 4.87351326761802E-07 / Projection distortion parameter AP_2_1 = 1.27034030399091E-11 / Projection distortion parameter AP_2_2 = -1.56781266527055E-14 / Projection distortion parameter AP_3_0 = 1.05609755183277E-10 / Projection distortion parameter AP_3_1 = -7.88269146660008E-15 / Projection distortion parameter AP_4_0 = 7.2479268759029E-14 / Projection distortion parameter BP_ORDER= 4 / Distortion order for BP BP_0_1 = -3.61828763934933E-08 / Projection distortion parameter BP_0_2 = 6.74317721340757E-08 / Projection distortion parameter BP_0_3 = 5.55651238033361E-11 / Projection distortion parameter BP_0_4 = -6.14956183405539E-15 / Projection distortion parameter BP_1_0 = -5.01841161647254E-08 / Projection distortion parameter BP_1_1 = 3.52809293826258E-07 / Projection distortion parameter BP_1_2 = 6.13665672637537E-11 / Projection distortion parameter BP_1_3 = -1.64300062518365E-14 / Projection distortion parameter BP_2_0 = -3.4054112344702E-08 / Projection distortion parameter BP_2_1 = 7.45298070253847E-11 / Projection distortion parameter BP_2_2 = 2.92717865788252E-15 / Projection distortion parameter BP_3_0 = -4.14183265768535E-11 / Projection distortion parameter BP_3_1 = 3.23323470749901E-14 / Projection distortion parameter BP_4_0 = -9.48831790286268E-15 / Projection distortion parameter / DATA FLOW ORIGNAME= '/data/PTF_default_56753.fits' / Filename as written by the camera FILENAME= 'PTF201307021787_2_o_56753.fits' / Filename of delivered camera image PROCORIG= 'IPAC-PTF pipelines' / Processing origin PROCDATE= 'Wed Jul 3 19:07:12 2013' / Processing date/time (Pacific time) PTFVERSN= 5. / Version of PTFSCIENCEPIPELINE program PMASKPTH= '/ptf/pos/archive/fallbackcal/pmasks/' / Pathname of pixel mask PMASKFIL= '70sOn35s_pixmask_chip11.trimmed.v4.fits' / Filename of pixel mask SFLATPTH= '/ptf/pos/sbx2/2013/07/02/f2/c11/cal/p4/cId98684/' / Pathname of super SFLATFIL= 'PTF_201307020000_i_s_flat_t120000_u000098684_f02_p000000_c11.fits' SBIASPTH= '/ptf/pos/sbx2/2013/07/02/f2/c11/cal/p1/cId98678/' / Pathname of super SBIASFIL= 'PTF_201307020000_i_s_bias_t120000_u000098678_f00_p000000_c11.fits' DBNID = 1550 / Database night ID DBEXPID = 342625 / Database exposure ID DBRID = 5743856 / Database raw-image ID DBPID = 16616794 / Database processed-image ID DBFID = 2 / Database filter ID DBPIID = 1 / Database P.I. ID DBPRID = 17 / Database project ID DBFIELD = 342625 / Database field ID DBSVID = 51 / Database software-version ID DBCVID = 57 / Database config-data-file ID INFOBITS= 0 / Database infobits (2^2 and 2^3 excluded)
Now let's plot the image data. But let's use sigma-clipping to pick a nice scale for the image.
mean, median, std = astropy.stats.sigma_clipped_stats(hdu.data)
plt.figure(figsize=(20, 10))
plt.imshow(hdu.data, vmin=mean-std, vmax=mean+3*std, cmap='binary')
plt.xlabel('pixel $x$')
plt.ylabel('pixel $y$')
Text(0, 0.5, 'pixel $y$')