Transcript Lecture 2

Lecture 2
• More python.
– Good books:
• “Learning Python”, Lutz & Ascher, O’Reilly (the ‘rat
book’);
• “Python in a Nutshell”, Martelli, O’Reilly (the ‘snake
book’).
• FITS files.
NASSP Masters 5003S - Computational Astronomy - 2009
OO
•
•
OO programming is the business of
constructing objects. This is done via the
class statement. You probably won’t
need to use this. However, objects
themselves are inescapable in python –
if you don’t make them, someone else
will.
In fact, almost everything in python is an
object.
OO
•
Objects have two extra features:
1. Attributes, which have names separated
from the object name by a dot, eg
fred.shape is an attribute named shape
which belongs to object type of which fred is
an example (eg a numpy array object).
2. Methods: these are attributes which are
functions. Like any function they can accept
arguments, and must still be written with
empty brackets even if there are no
arguments. Eg the append() method of lists
and the copy() method of dictionaries.
Python variables
• Ways to arrange several objects:
– A single one is called a scalar.
– A list: elements numbered*, data types can be
different. [ , 3.14, ‘fred’,…]
– A tuple: similar to a list but written with round
brackets rather than square. ( , 3.14, ‘fred’,…)
– A dictionary: elements accessed by a key,
data types can be different; curly brackets.
{‘foo’: , ‘pi’:3.14, ‘him’:’fred’…}
• These groupings are themselves objects.
*Python numbering starts from zero.
Python scalars
• Numbers
– Integers
– Floating-point
– Complex
•
•
•
•
-43, 0, 1, 219
-2.9, 4.76e6
0.2 + 5.3j
True, False
Logical
‘Fred’, ‘impossible’, ‘23’,
Strings
None
Special
Objects (includes everything)
‘True’
Python operations
• They are mostly pretty standard.
– Arithmetic binary operators:
+, -, *, /
• Note that they can have different effect on strings.
• Note that
a += 1 # etc
is short for
a = a + 1 # etc
– Logical operators:
and, or, not
– Comparison operators:
>, >=, <, <=, ==
Python traps
• The most common is really a unix trap.
Suppose you write your nice script
dothings.py. You’ve been careful with your
grammar, and have as your first line
#!/usr/bin/python
• You try to run dothings.py on the
command line, but get:
bash: ./dothings.py: Permission denied
• To fix this you have to do
chmod +x dothings.py
• This makes dothings.py ‘executable’.
NASSP Masters 5003S - Computational Astronomy - 2010
Python traps
• For some types of python object, several
names can point to a single data area.
Changing one of these then changes all
the others. This is usually not what you
want to happen.
• Since it is a common thing to want to
create a new variable name with the value
of an old one, we have to be careful of this
behaviour in python, and find workarounds for variable types where it is a
problem.
Python traps
• Compare the behavior of scalars with lists.
We can spawn a new variable name from
a scalar, which then has its own memory
area:
>>>a=4
>>>b=a
>>>a=3
>>>print b
• This should print ‘4’. But now try:
>>>aa=[1,2,3]
>>>bb=aa
>>>aa[0]=999
>>>print bb
• You’ll get ‘[999,2,3]’. Cf mutability…
Trouble with ‘change in place’
Type:
Mutable? Avoid the problem by:
Simple scalar
No
Scalar object
It
Use copy module
depends..
List
Yes
Tuple
No
Dictionary
Yes
bb=aa[:]
bb=aa.copy()
More python traps:
• The range statement:
>>>mylist = range(4)
>>>print mylist
[0, 1, 2, 3]
– Note last element is not 4, but 3.
• List slices:
>>>print mylist[:2]
[0, 1]
– Last element is not 2, but 1.
• CHECK THESE THINGS before you
accept the results of your code!
More python traps:
• Truncation during division. Try:
>>>a=4
>>>b=3
>>>print a/b
– This returns, not 1.333, but 1. This is usually
not what you want.
• Better to get into the habit of ALWAYS
dividing by a floating-point. Eg:
>>>a=4
>>>b=3
>>>print a/float(b)
– yields 1.3333. If you WANT the truncation,
use // not /.
Python control structures
if <test1>:
# do some stuff
elif <test2>: # you can have 0 or more of these
# do some other stuff
else: # you can have 0 or 1 of these
# third lot of stuff
•There is NO ‘goto’ statement in python. This is a feature.
•Anything after a # is a comment.
while <test>:
# do loop stuff
break # optional – dumps out and avoids the ‘else’.
continue # like ‘goto while’.
else: # optional - processing after normal loop exit.
# stuff to do after normal loop exit
for <item> in <list>: # etc
Some neato modules
• math
– Gives you access to constants like pi, + math
operations such as sqrt(), cos(), atan2() etc.
• os - particularly os.path
– Lets you check if files exist, lets you delete
files etc
import os
if os.path.exists(‘myfile’):
# do something. Maybe do
os.remove(‘myfile’)
NASSP Masters 5003S - Computational Astronomy - 2009
One more neato module
• sys
– sys.argv is a list of command-line arguments.
• Beware though, they are all passed as strings.
– Suppose you put the following statements in
your python script bigdogs.py:
#!/usr/bin/python
import sys
bb = sys.argv[1]
cc = int(bb)
print bb>0, cc>0
– If you then run it on the command line, with an
additional argument 0:
$ bigdogs.py 0
True False
– Because bb is a string ‘0’ and cc is an integer 0.
NASSP Masters 5003S - Computational Astronomy - 2009
Python ins and outs
• We’re going to mostly read our data from
FITS files, using a module called pyfits.
– http://www.stsci.edu/resources/software_hardwa
re/pyfits/
• We’ll crunch up the data using modules
called numpy and scipy.
– http://numpy.scipy.org/
– http://www.scipy.org/
• For graphical output we’ll use module pylab.
– http://matplotlib.sourceforge.net/index.html
Manuals and example code
• You won’t need to download manuals for
pyfits, numpy, scipy or pylab. You can get
them from my home page:
– http://www.ast.uct.ac.za/~ims/teaching/teaching.html
• Some of these manuals are huge – I
recommend you neither print or read them
in entirety, but rather:
– Read the bits you need using acroread.
– Look at the example code I’ll provide you from
time to time. There are already a few
available from my page.
NASSP Masters 5003F - Computational Astronomy - 2009
Another few cm of the python.
• Function returns. Note:
>>> def some_func(arg1,arg2):
...
a = arg1 + arg2
...
b = arg1 – arg2
...
return a, b
...
>>> c = some_func(42, 2)
>>> print c
(44, 40)
• This is a tuple. It is better to be explicit...
>>> def some_func(arg1,arg2):
...
a = arg1 + arg2
...
b = arg1 – arg2
...
return (a, b)
NASSP Masters 5003F - Computational Astronomy - 2009
Another few cm of the python.
• Also in the function call:
>>> (summ, diff) = some_func(42, 2)
– is better than
>>> c = some_func(42, 2)
• You can also force it to return a list:
>>> def anuther_func(arg1,arg2):
...
a = arg1 + arg2
...
b = arg1 – arg2
...
return [a, b]
...
>>> c = anuther_func(20, 19)
>>> print c
[39, 1]
NASSP Masters 5003F - Computational Astronomy - 2009
Another few cm of the python.
• The for statement - most of you know this
mylist = [‘peas’,’beans’,’dead dog’]
way: >>>
>>> for i in range(len(mylist)):
...
item = mylist[i]
...
print item
...
peas
beans
dead dog
• But you can also do this:
>>> mylist = [‘peas’,’beans’,’dead dog’]
>>> for item in mylist:
...
print item
...
peas
beans
dead dog
• Saves a bit of typing.
NASSP Masters 5003F - Computational Astronomy - 2009
FITS = Flexible Image Transport System
This will be the default data format for this course. See:
http://fits.gsfc.nasa.gov/fits_standard.html
Header-Data Units
“Cards”
2880-byte blocks
1
80
1
Primary
HDU
2
Header
3
4
Header
…
36
Extension 1
KEYWORD = VALUE / [UNIT] COMMENT
Unused cards of last header block filled with spaces (ASCII 32)
etc…
Binary table:
Data
Extension 2
Data
Extension 3
or
‘Image’
Data
Data
etc…
etc…
Unused bytes of last data block filled with spaces (ASCII 32)
•
•
•
•
•
WCS keywords of FITS files
WCS stands for World Coordinate System.
http://fits.gsfc.nasa.gov/fits_wcs.html
What they’re for: pixellated data – ie
samples of some quantity on a regular
grid.
WCS keywords define the mapping
between the pixel index and a world
coordinate system.
Eg: a 2d image of the sky. We want to
know which sky direction the (j,k)th pixel
corresponds to.
NASSP Masters 5003F - Computational Astronomy - 2009
Eg, projection onto a tangent plane.
WCS must encode the
relation between θ and the
pixel number.
θ
NASSP Masters 5003F - Computational Astronomy - 2009
WCS example continued
The general formula in this case is
w- wref
(p - pref) * scale = tan(w – wref).
p is the pixel coordinate and w the
world coordinate. w might eg be
right ascension or declination.
WCS must describe 4 things:
1. pref
2. wref
3. scale
4. the nature of the functional relation.
Perhaps also world units.
p- pref
Note:
(1) pref can be real-valued;
(2) By convention, p at the centre
of the 1st pixel = 1.0.
NASSP Masters 5003F - Computational Astronomy - 2009
WCS keywords for array extensions
• In what follows, n is an integer, corresponding to
one of the dimensions of the array.
–
–
–
–
CRVALn – wref.
CRPIXn – pref.
CDELTn – scale.
CTYPEn – an 8-character string encoding the function
type (eg ‘TAN---RA’). There is an agreed list of these.
– CUNITn – string encoding the unit of w (eg ‘deg’). Also an
agreed list.
• In addition, rotated coordinate systems can be
defined via either adding PCi_j keywords to the
above scheme, or replacing CDELTn by CDi_j
keywords. But I don’t want to get too deeply into
this.
• Analogous (starting with T) WCS keywords are also defined for table
columns.
NASSP Masters 5003F - Computational Astronomy - 2009
FITS interfaces and software:
• CFITSIO
– Aging, user-unfriendly interface to FITS. Most
more modern interfaces are built on top of it
though.
• FTOOLS (from HEASARC web site)
– Very useful set of tools for doing everything
you can think of to or with FITS-format data.
– But (Feb 2010: were) too many tasks bundled
into the one large, not-very-portable package.
• DS9
– A very useful, compact, portable FITS image
viewer.
NASSP Masters 5003F - Computational Astronomy - 2009
A pyfits trap:
• On the nassp computers, pyfits does not
return numpy arrays, it returns numarray
arrays. So if you do eg
>>>import numpy as nu
>>>import pyfits as pf
>>>fred = nu.zeros([3,3],nu.float)
>>>hdu = pf.PrimaryHDU(fred)
• you will get an error. A workaround is:
>>>import numarray as old_nu
>>> hdu = pf.PrimaryHDU(old_nu.array(fred))
• And similar for reading files – convert them
using numpy.array().
NASSP Masters 5003S - Computational Astronomy - 2010