- We are going to only skim this section, but it is very useful as a reference and to help understand some of the language constructs and what they are doing
- How to generate a random array for stellar locations or magnitude.
- Adding the PSF back into your science image using IRAF
- Making useful plots of the results.
It assumes that you already have good knowledge of aperture and PSF photometry. The PSF photometry best practices and the full process of performing psf fitting and subtraction to a cleanly subtracted image is not discussed. However, I did make a separate notebook with an outline of PSF photometry processing which you can reference. It's saved under the lecture notebook directory.
There is a good tutorial series on photometry, especially with HST and DAOPhot, given by Dave Zurek, Brad Whitmore and Harry Ferguson many moons ago. The files are still available here if you are new to psf photometry and looking for some help: http://www.stsci.edu/hst/training/events/Photometry/SSD990302/
# Import some useful modules for this notebook.
# The __future__ module allows you to use new language features
# before they become standard in a future release.
# This one we can use for opening a file and immediately closing it,
# freeing up filehandles faster. The with statement can be used for
# other things besides file handling - it's a kind of shorthand function
# that ensures whatever needs to be opened, checked or closed with the
# process is automatically executed. This has a good explanation:
# http://effbot.org/zone/python-with-statement.htm
#import with_statement - standard after python 2.6
# New print() function - print("") instead of just print ""
from __future__ import print_function
import shutil
#astropy hasn't released the formatted read fix for Daophot yet
from astropy.io import ascii
from astropy import table
#if you get an error that no login.cl is found when you import iraf
#you can either "mkiraf" in the directory the script is being imported from
#or you can make sure to have the default iraf directory setup on your home
#account: go to ~/iraf/ and "mkiraf" if you don't have it setup.
from pyraf import iraf
from iraf import noao, digiphot, daophot, ptools
from iraf import images, imcoords
from iraf import artdata
import pyfits
import numpy as np
import matplotlib.pyplot as plt
#if you didn't start your notebook with pylab inline: "ipython notebook --pylab inline"
#then call the %pylab magic command
#
#If your plots are still acting odd try plt.show()
#
#I use !open to see some of the PDF files because I'm on a MAC, you might have to use something else
If you are inside the iPython notebook and would like to start a script file, insert the magic below at the top of the cell:
%%file filename
The contents of cell following that will get saved in the specified file.
Usually Python code is saved in a filename ending in .py
, but for use only as a script, it doesn't have to be. You can call it anything you like, but unless you use it a lot, keeping .py
is a good idea.
%%file bozo.py
print("Help, I'm a clown")
Overwriting bozo.py
You can also execute that script from inside the notebook session:
!python bozo.py
Help, I'm a clown
The reason we had to specify "python" to run the code is because the file is not executable and doesn't point to the the Python interpreter. We can fix that though:
%%file sadbozo
#!/usr/bin/env python
# this points to the python installed on your PATH
print("a sad clown")
Overwriting sadbozo
# Now make it executable and run it
!chmod a+x sadbozo
!sadbozo
a sad clown
If you are working in an iPython shell (simulated in the notebook) you can use the %run
magic:
%run bozo.py
Help, I'm a clown
Be careful in this case because any variables used in the script remain in your current environment!
%%file bozo2.py
name = "Bozo the Clown"
Overwriting bozo2.py
%run bozo2.py #this runs and produces no output
name # BUT now "name" has been set in our local namespace!
'Bozo the Clown'
Important things to remember:
You can write a Python script in any text editor.
You can make the script executable by itself using a shebang (#!
) and pointing it to your environment, and make the file executable. (As shown above.)
Scripts need to import all the modules they depend on.
Running a script in a Python session will load all of its variables into your current environment. Watch for common use names! See http://docs.python.org/2/tutorial/classes.html#python-scopes-and-namespaces
For good Python coding practices, see PEP 8 (http://www.python.org/dev/peps/pep-0008/_).
def square(x):
"""
Take the square of the input value, x.
asdfja;slfsakdfj;lasjkdf;kjsad
"""
return x * x
square(4)
16
See the document string I put in the function? It's a good practice to describe what your functions do. It can be useful for yourself and other users.
Currently, Sphinx (http://sphinx-doc.org/) is the documentation tool used by SSB. It uses "reStructuredText" (http://docutils.sourceforge.net/rst.html) formatting.
square.__doc__
'\nTake the square of the input value, x.\nasdfja;slfsakdfj;lasjkdf;kjsad\n'
help(square)
Help on function square in module __main__: square(x) Take the square of the input value, x. asdfja;slfsakdfj;lasjkdf;kjsad
More on the significance of __main__
to come...
def print_me(a, b, c=2):
"""
Print the input values using new-style `print` function
available in Python 3.x.
.. note::
`c` is an optional keyword.
If not given, it defaults to 2.
"""
print(a, b, c)
print_me(2, 3, 4)
print_me(1, 2)
print_me(a=1, b=7)
2 3 4 1 2 2 1 7 2
This will crash because value for b
is not given.
It sees that we only specified a default for the c
and wants the expected value for b
.
print_me(1)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-16-5534e33b4a02> in <module>() ----> 1 print_me(1) TypeError: print_me() takes at least 2 arguments (1 given)
ERROR: TypeError: print_me() takes at least 2 arguments (1 given) [IPython.core.interactiveshell]
This will crash because there is nowhere to store 4.
print_me(1, 2, 3, 4)
help(print_me)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-17-dc27af1ec62e> in <module>() ----> 1 print_me(1, 2, 3, 4) 2 help(print_me) TypeError: print_me() takes at most 3 arguments (4 given)
ERROR: TypeError: print_me() takes at most 3 arguments (4 given) [IPython.core.interactiveshell]
def powers(x):
"""Returns multiple powers of x."""
return x, x**2, x**3, x**4
powers(3)
(3, 9, 27, 81)
Notice they are enclosed in parenthesis? This is called a tuple.
Tuples
are like lists, EXCEPT they are immutable, and always indexed the same way.
DO NOT have to contain the same data type for each element.
can contain a mutable item, like a list. That item can still be changed, but its location in the tuple cannot.
can have 1 or more elements. HOWEVER, a single-element tuple is distinquished from a number with a comma. Example: (1,)
you can create an empty tuple like this: t=tuple()
More information on tuples: http://docs.python.org/2.7/tutorial/datastructures.html#tuples-and-sequences
tup = (1, "bozo", [1, 2, 3])
print('1st item:', tup[0])
print('2nd item:', tup[1])
print('3rd item:', tup[2])
print('last element of 3rd item:', tup[2][-1]) # Tricky!
1st item: 1 2nd item: bozo 3rd item: [1, 2, 3] last element of 3rd item: 3
#want to see what type something is?
type(tup)
tuple
number = (3)
print(type(number), number)
single_tuple = (3,)
print(type(single_tuple), single_tuple)
<type 'int'> 3 <type 'tuple'> (3,)
You can also assign members of a tuple to separate variables (also known as unpacking):
a, b, c = tup
print(a)
print(b)
print(c)
1 bozo [1, 2, 3]
Now we can unpack all the returned values from our powers
function above:
first, second, third, fourth = powers(3)
print(first)
print(second)
print(third)
print(fourth)
3 9 27 81
def baddefault(inputlist, previouscases=[]):
"""Take an input list, extract negative values and
append them to previous cases. Return previous cases.
"""
for item in inputlist:
if item < 0:
previouscases.append(item)
return previouscases
# Use non-default previouscases
print(baddefault([1,2,-5], [-3,-1]))
# Use default previouscases
print(baddefault([1,2,-5]))
[-3, -1, -5] [-5]
So far, so good. Now...
# Use default previouscases, again.
print(baddefault([-1, -2]))
[-5, -1, -2]
Expected answer is
[-1, -2]
but we got
[-5, -1, -2]
where the -5 is inherited from the previous call!
Modules:
import
only imports a module once. If Python knows it is already loaded, import will not be done without an explicit reload
command.
%%file testmodule.py
def testfunction(a, b):
"""Return a sum of inputs."""
return a + b
Overwriting testmodule.py
import testmodule
testmodule.testfunction(2, 3)
5
%%file testmodule.py
def testfunction(a, b):
"""Return a multiplication of inputs."""
return a * b
Overwriting testmodule.py
testmodule.testfunction(2, 3)
5
The answer above is still the sum, not the product! We have to reload the function and try again.
reload(testmodule)
testmodule.testfunction(2, 3)
6
Layered Imports
If you have layered imports, e.g.:
a.py
imports b.py
b.py
imports c.py
When you change c.py
and do reload(a)
, it will not do reload(c)
for you as it is only looking for changes in a.py
.
So, you have reload in this order, which can become tedious:
reload(c)
reload(b)
reload(a)
Luckily, iPython has a magic dreload
(http://ipython.org/ipython-doc/dev/interactive/reference.html#dreload_) that tries to handle all those dependencies.
Avoid import *
It is generally bad form to use
from module import *
in a script because it leads to hard to understand code.
Try to be explicit:
from testmodule import testfunction
You can make your module usable as a script/executable by testing for the name "main" at the bottom:
if __name__ == "__main__":
print(testfunction(1., 2.))
This works because the code that parses the command line only runs if the module is executed as the "main" file. It can provide a nice user interface to your module or be used to execute a suite of tests to make sure it's working as expected.
%%file testmodule.py
#! /usr/bin/env python
def testfunction(a, b):
"""Return a division of inputs."""
return a / b
if __name__ == "__main__":
print(testfunction(1., 2.))
Overwriting testmodule.py
!chmod a+x testmodule.py
!testmodule.py
0.5
%%file bozo3.py
#! /usr/bin/env python
# This portion is executed when someone do this:
# import bozo3
import sys
def clown_greeting(name='Bozo'):
print("I'm a clown and my name is"), name
# This portion only runs in a script context:
if __name__ == "__main__":
# Check if shell-level command line arguments are present
if len(sys.argv) > 1:
# Ignore the first item since that is the name of this script.
# Call the function above once for each subsequent items.
for item in sys.argv[1:]:
clown_greeting(item)
else:
print("USAGE: bozo3.py name1 [name2 name3 ...]")
Overwriting bozo3.py
import bozo3
bozo3.clown_greeting()
bozo3.clown_greeting('Crusty')
I'm a clown and my name is Bozo I'm a clown and my name is Crusty
!chmod a+x bozo3.py
!bozo3.py "Ronald McDonald" It
I'm a clown and my name is Ronald McDonald I'm a clown and my name is It
Python searches the environmental variable PYTHONPATH to find the module you have specified.
Usually the current directory is part of the search path, but it might not be the FIRST thing.
You can modify the path dynamically by changing the contents of the sys.path
list (not good for standard use, but useful for experimentation sometimes) .
For more information see http://docs.python.org/2.7/tutorial/modules.html
To see what the current path is:
import sys
sys.path
['.', '', '', '/Users/sosey/ssb/Ureka/variants/common/lib/python/astropy-0.1-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/variants/common/lib/python', '/Users/sosey/ssb/Ureka/python/lib', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/distribute-0.6.28-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_core-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_ServiceManagement-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_ServerNotification-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_CoreLocation-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_AppleScriptObjC-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_ScriptingBridge-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_QTKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_PubSub-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_InterfaceBuilderKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_InstantMessage-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_InputMethodKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_DictionaryServices-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_Collaboration-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_CalendarStore-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_XgridFoundation-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_WebKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_SystemConfiguration-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_SyncServices-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_SearchKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_ScreenSaver-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_Quartz-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_PreferencePanes-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_Message-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_LaunchServices-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_LatentSemanticMapping-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_InstallerPlugins-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_FSEvents-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_ExceptionHandling-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_CoreText-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_CoreData-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_Cocoa-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_CFNetwork-2.3-py2.7-macosx-10.6-x86_64.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_Automator-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_AppleScriptKit-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pyobjc_framework_AddressBook-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/tornado-2.3-py2.7.egg', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/pip-1.1-py2.7.egg', '/Users/sosey/teaching/scientific-python-training-2012/lecture_notebooks', '/Users/sosey/ssb/Ureka/python/lib/python27.zip', '/Users/sosey/ssb/Ureka/python/lib/python2.7', '/Users/sosey/ssb/Ureka/python/lib/python2.7/plat-darwin', '/Users/sosey/ssb/Ureka/python/lib/python2.7/plat-mac', '/Users/sosey/ssb/Ureka/python/lib/python2.7/plat-mac/lib-scriptpackages', '/Users/sosey/ssb/Ureka/python/lib/python2.7/lib-tk', '/Users/sosey/ssb/Ureka/python/lib/python2.7/lib-old', '/Users/sosey/ssb/Ureka/python/lib/python2.7/lib-dynload', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/PIL', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/setuptools-0.6c11-py2.7.egg-info', '/Users/sosey/ssb/Ureka/python/lib', '/Users/sosey/ssb/Ureka/python/lib/python2.7/site-packages/IPython/extensions']
Get the value of PATH environment variable and extract the first path defined:
import os
#or maybe you want to record whose using this script?
user=os.getenv("USER")
print(("Hello %s! Thanks for using my script.")%(user))
Hello sosey! Thanks for using my script.
I wrote a script, aperphot.py
, to run aperture photometry, source finding and plotting, which will help make things run faster.
If you make the file executable you can call it from the shell command line with options:
i [my input image name]
a [what pixel aperture to perform photometry]
s [what pixel from the center to start sky measurements]
w [how wide in pixels the sky annnulus should be]
here is an example of how you might call it from the command line prompt:
aperphot.py -i image.fits -a 4 -s 8 -w 3 -p
You can also call the script from a Python session just by importing it:
import aperphot
aperphot.run(...) --> creates image.fits.phot
There is a function in the aperphot
module to plot of some of the results for quick reference.
aperphot.plotphot('image.fits.phot') --> creates image.fits.phot.pdf
I have already run the script and saved data files for this lecture.
Let's see that we got reasonable photometry back by reading the table with AstroPy reader and print out the contents:
I'll import the module and run it here
!open ../data/image.fits.phot.pdf #made by aperphot.run('../data/image.fits', aper=4, sky=8, width=3, plots=True)
Since I told it to make plots we can look at the output :)
Now that we've checked out our plots, we can play around with our detections.
Let's say we only want to keep a subset of the points we detected. We'll make a refining cut on the objects, leaving out the ones which had errors in the photometry and centering worse than 1 pixel.
reader = ascii.Daophot()
phot = reader.read('../data/image.fits.phot')
good = np.where((phot['MERR'] != 'INDEF') &
(np.abs(phot["XSHIFT"]) <= 1.) &
(np.abs(phot["YSHIFT"]) <= 1.))
phot = phot[good[0]] #notice I am overwriting the original data!
phot.pprint()
IMAGE XINIT YINIT ID ... MAG MERR PIER PERROR ------------- ------- ------- ---- ... ------ ----- ---- ------- image.fits[1] 57.106 6.779 20 ... 21.479 0.515 0 NoError image.fits[1] 122.832 4.803 22 ... 24.378 7.415 0 NoError image.fits[1] 22.615 5.853 24 ... 18.581 0.070 0 NoError image.fits[1] 68.257 6.351 25 ... 20.329 0.251 0 NoError image.fits[1] 91.009 5.895 26 ... 20.382 0.204 0 NoError image.fits[1] 214.969 5.96 28 ... 18.433 0.025 0 NoError image.fits[1] 133.468 7.586 35 ... 20.402 0.049 0 NoError image.fits[1] 16.305 10.221 40 ... 20.535 0.421 0 NoError image.fits[1] 50.934 9.889 41 ... 20.172 0.167 0 NoError image.fits[1] 228.508 10.51 45 ... 15.539 0.002 0 NoError image.fits[1] 65.875 10.536 47 ... 19.361 0.076 0 NoError ... ... ... ... ... ... ... ... ... image.fits[1] 51.908 247.496 1392 ... 20.170 0.043 0 NoError image.fits[1] 193.397 247.317 1395 ... 19.754 0.099 0 NoError image.fits[1] 220.22 246.968 1397 ... 20.175 0.042 0 NoError image.fits[1] 178.202 247.1 1401 ... 24.015 5.921 0 NoError image.fits[1] 251.794 247.987 1402 ... 19.529 0.044 0 NoError image.fits[1] 249.338 249.327 1409 ... 20.083 0.121 0 NoError image.fits[1] 19.491 249.658 1411 ... 20.652 0.176 0 NoError image.fits[1] 99.895 250.324 1413 ... 17.735 0.003 0 NoError image.fits[1] 241.12 249.704 1417 ... 18.519 0.012 0 NoError image.fits[1] 38.967 251.122 1421 ... 19.394 0.044 0 NoError image.fits[1] 130.693 252.377 1430 ... 21.260 0.705 0 NoError
Now, if we want to save these cherry picked star locations and redo the photometry, we can save the (X, Y) coordinates and rerun our phot function.
xy = phot['XCENTER', 'YCENTER']
with open('new_coords.dat', 'w') as output_file:
for line in xy.pformat(max_width=-1, max_lines=-1):
output_file.write(line + "\n")
We can rerun the photometry on our finalized list of stars.
aperphot.do_phot('../data/n8q624e8q_cal.fits', 'new_coords.dat', aper=4, sky_annulus=8, width_sky=3)
aperphot.plotphot('../data/n8q624e8q.phot')
!open ../data/n8q624e8q.phot.pdf
Not too bad... We've saved fewer stars for this run, but we have cleaner shifts on the their centering, which means they are probably a mix of bright stars and isolated stars.
For a full PSF fitting process, you will probably want to repeat the star finding and centering and photometry passes through several depths of your image as you remove the spoiler stars around your psf model stars with runs of allstar
(shown later in this notebook).
You can move on and start your PSF photometry and create a PSF using your reduced data. I'm going to skip over the steps of making a model PSF because they are fairly interactive and time consuming, but you can use the PSF photometry primer as a reference.
Let's start by just looking at a typical NICMOS PSF.
This is an image generated by TinyTim with NICMOS specifications with a normalized flux in the same filter/camera as our dataset.
tinytim = pyfits.getdata('../data/nic2_f160w_4xsub00.fits')
# Gives us the size of the image
xsize,ysize=tinytim.shape #use of tuples
print(xsize,ysize)
528 528
# Display the image in log scale so we can see the gory details.
# Everyone loves a good NICMOS PSF.
plt.imshow(np.log10(tinytim), cmap=plt.cm.gray)
<matplotlib.image.AxesImage at 0x11b6367d0>
# Let's look at a plot of the profile at Y=264
plt.plot(tinytim[164:364,264])
[<matplotlib.lines.Line2D at 0x11b7c3490>]
#you might try computing an aperture correction from the Tiny Tim image you made:
import apercor
apercor.apercor('../data/nic2_f160w_4xsub00.fits')
Setting DAOpars to NICMOS defaults... Using fwhmpsf: 3.500000 Phot aperture: 1,5,9,13,17,21,25,29,33,37,41,45,49,53,57,61,65,69,73,77,81,85,89,93,97,101,105,109,113,117,121,125,129,133,137,141,145,149,153,157,161,165,169,173,177,181,185,189,193,197,201,205,209,213,217,221,225,229,233,237,241,245,249,253,257,261 pixels Sky Annulus: 508 -> 511 pixels Removing previous photometry file Saving output files as ../data/nic2_f160w_4xsub00.fits.phot 264.941 265.019 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. saved output figure to ../data/nic2_f160w_4xsub00.fits.phot.pdf
This is a PSF generated from the drizzled NICMOS dataset using daophot.psf
, which is really just saving the model and residual information so that it can be fitted to each star in your image later using allstar
. The process to get here involves some repition on finding and removing stars to clean up your model. (see the PSF primer)
psf_model = pyfits.open('../data/dr_median.psf.1.fits')
psf_model.info() # Like catfits
Filename: ../data/dr_median.psf.1.fits No. Name Type Cards Dimensions Format 0 PRIMARY PrimaryHDU 71 (39, 39, 3) float32
There are 3 (39,39) element arrays which specify the model PSF lookup tables fit residuals and the functional fit which is stored in the primary header along with information about each star I used. Because I used a VARORDER = 1, there are 2 tables with the psf derivatives relative to the x,y location of the star in the frame.
# These are the parameters which I ended up using
psf_hdr = psf_model[0].header
psf_dat = psf_model[0].data
psf_model.close()
keys = ['VARORDER', 'FUNCTION', 'PSFRAD', 'PSFMAG','NPSFSTAR','PAR1', 'PAR2', 'PAR3']
for item in keys:
print(item, psf_hdr[item])
VARORDER 1 FUNCTION lorentz PSFRAD 9.0 PSFMAG 21.05 NPSFSTAR 7 PAR1 0.529933 PAR2 0.510016 PAR3 -0.1086451
The header says a Lorenz function for the analytic portion of the fit. We might as well take a look at the residuals from the fitting
plt.clf() #clear out the plot
plt.plot(np.abs(psf_dat[0][0]), label='Analytical fit to center') #subsampled by a factor of 2
plt.plot(np.abs(psf_dat[0][1]), label='First Deriv x')
plt.plot(np.abs(psf_dat[0][2]), label='First Deriv y')
plt.xlabel('x or y pixel')
plt.ylabel('abs value')
plt.title('Residuals for the PSF model')
plt.legend()
<matplotlib.legend.Legend at 0x11c8323d0>
Let's say we are happy with this PSF model and decide to use it to subtract some stars in our science data.
The data we used below was generated with this command:
daophot.allstar('../data/n8q624e8q_cal.fits[1]', '../data/n8q624e8q.phot',
'../data/dr_median.psf.1.fits', allstarfile="default",
rejfile="default", subimage="default", verify="no")
subimage = pyfits.getdata('../data/n8q624e8q_cal.fits1.sub.1.fits')
science_image=pyfits.getdata('../data/n8q624e8q_cal.fits','sci',1)
fig = plt.figure()
fig.set_figheight(12) #in inches
fig.set_figwidth(12) #in inches
# Shows subtracted image
ax1 = fig.add_subplot(121)
ax1.imshow(subimage,vmin=-0.1,vmax=1, cmap=plt.cm.gray)
ax1.set_title('Subtracted')
# Shows original image
ax2 = fig.add_subplot(122)
ax2.imshow(science_image,vmin=-0.1,vmax=1,cmap=plt.cm.gray)
ax2.set_title('Original')
ax2.set_xlim(0, 250)
ax2.set_ylim(0, 250)
(0, 250)
See all the junk where the stars were subtracted? This means your PSF wasn't a good fit, there are spoilers in your model star local area or the centering is a bit off. I expected some of this because I used a PSF created from a different observation set (also the same camera and filter, but drizzled at a difference scale). I just thought seeing a poor subtraction is a little more instructive in a crowded field that a good subtraction which should have featureless residuals. It's also easier to see in a talk. :)
Where do you go from here? This is where scripting comes in handy. After you've created a PSF with which you are pleased (which is an iterative process), you can script the iterative detection and centering of stars, photometry, PSF subtraction, and do a final psf subtraction to end up with a cleanly subtracted image.
Then, once you have photometry in hand with all your stars, you can go back and see how well you did by running a completeness test!
Using the same PSF model you did your photometry with, add it back into your original science image and see how well you can retrieve it. Here is a general outline of what you might do:
Create a random list of X,Y,MAG for the input artificial stars (or pull the locations or mag from a known distribution)
Add them to a copy of the science image (make a plot of the location of the stars for reference).
Detect and return photometry for all stars.
Match the input and output photometry to see how well the artificial stars were found and measured (make some plots).
Plot up a histogram of how many stars were recovered in each magnitude bin.
DONT FORGET TO TAKE CROWDING OF YOUR FIELD INTO ACCOUNT!
I won't go into much detail here, but lets look at a couple useful things.
science_image = '../data/n8q624e8q_cal.fits'
art_image = 'n8q624e8q_art.fits'
#make a copy of your science image to work with
shutil.copy(science_image, art_image)
xsize = pyfits.getval(science_image, 'naxis1', extname='sci', extver=1)
ysize = pyfits.getval(science_image, 'naxis2', extname='sci', extver=1)
Create a random list of star locations in IRAF.
starlist
in the artdata has a few distributions to choose from, and you can also let addstar
do this itself. (I like knowing ahead of time.)
artdata.starlist(starlist='n8q624e8q.starlist',
nstars=100, interactive="no", spatial="uniform",
xmin=5, xmax=xsize-5, ymin=5, ymax=ysize-1,
luminosity="uniform", minmag=18, maxmag=25)
# Take a peak at the file
!head n8q624e8q.starlist
# Fri 12:36:41 18-Jan-2013 # begin n8q624e8q.starlist # spatial uniform # xmin 5. # xmax 251. # ymin 5. # ymax 255. # luminosity uniform # minmag 18. # maxmag 25.
You could do the same thing as above in Python using the Numpy random generator.
from numpy import random
# To make sure you get the same random results every time,
# set a constant seed. Can be useful for testing.
# random.seed(1234) MENTION
#I dont want the same results every time because I will
#be calling this later to populate more stellar runs iteratively,
#so I wont seed it with the same number
# Generate random X and Y
xloc = np.asfarray(random.randint(5, xsize - 5, 100))
yloc = np.asfarray(random.randint(5, xsize - 5, 100))
# Generate random photometry between mags 18 and 25
phot = np.asfarray(random.randint(18, 25, 100)) + random.random_sample(100)
starlist = 'n8q624e8q.py_starlist'
of=open(starlist,'w')
# Save to a plain text file
for star in range(0,len(phot),1):
of.write(("%5.1f\t%5.1f\t%5.1f\t%i\n")%(xloc[star],yloc[star],phot[star],star))
of.close()
#here is another way to do the same thing
#with open(starlist, 'w') as output_file:
# for i, (xx, yy, ph) in enumerate(zip(xloc, yloc, phot), start=1):
# output_file.write('{:8.1f} {:8.1f} {:10.6f} {:4d}\n'.format(xx, yy, ph, i))
#
#And you can also use the astropy writer:
#xy = phot['XCENTER', 'YCENTER']
#xy.write('new_coords.dat', format='ascii')
!head n8q624e8q.py_starlist
71.0 87.0 19.9 0 20.0 144.0 18.2 1 83.0 141.0 18.8 2 17.0 94.0 24.8 3 135.0 81.0 20.6 4 100.0 39.0 21.1 5 63.0 147.0 22.4 6 136.0 16.0 21.7 7 80.0 41.0 23.3 8 87.0 53.0 24.6 9
Now we add the fake stars back into the science image using addstar
in IRAF.
art1='n8q624e8q_art1.fits'
#remove the old files
if os.access(art1, os.F_OK):
os.remove(art1)
os.remove(art1+".art")
iraf.addstar(image=art_image+'[1]', photfile=starlist,
psfimage='../data/dr_median.psf.1.fits',
addimage=art1, nstar=100, simple_text='yes')
Psf radius in scale units (15.): New psf radius: 15. scale units 15. pixels OUTPUT IMAGE: n8q624e8q_art1.fits Added Star: 1 - X: 71.00 Y: 87.00 Mag: 19.900 Added Star: 1 - X: 20.00 Y: 144.00 Mag: 18.200 Added Star: 2 - X: 83.00 Y: 141.00 Mag: 18.800 Added Star: 3 - X: 17.00 Y: 94.00 Mag: 24.800 Added Star: 4 - X: 135.00 Y: 81.00 Mag: 20.600 Added Star: 5 - X: 100.00 Y: 39.00 Mag: 21.100 Added Star: 6 - X: 63.00 Y: 147.00 Mag: 22.400 Added Star: 7 - X: 136.00 Y: 16.00 Mag: 21.700 Added Star: 8 - X: 80.00 Y: 41.00 Mag: 23.300 Added Star: 9 - X: 87.00 Y: 53.00 Mag: 24.600 Added Star: 10 - X: 125.00 Y: 47.00 Mag: 18.800 Added Star: 11 - X: 128.00 Y: 65.00 Mag: 19.100 Added Star: 12 - X: 58.00 Y: 67.00 Mag: 19.000 Added Star: 13 - X: 220.00 Y: 84.00 Mag: 21.400 Added Star: 14 - X: 34.00 Y: 106.00 Mag: 19.000 Added Star: 15 - X: 21.00 Y: 20.00 Mag: 20.800 Added Star: 16 - X: 156.00 Y: 22.00 Mag: 20.300 Added Star: 17 - X: 71.00 Y: 198.00 Mag: 19.200 Added Star: 18 - X: 205.00 Y: 210.00 Mag: 18.700 Added Star: 19 - X: 76.00 Y: 82.00 Mag: 23.800 Added Star: 20 - X: 194.00 Y: 166.00 Mag: 24.400 Added Star: 21 - X: 23.00 Y: 242.00 Mag: 22.300 Added Star: 22 - X: 158.00 Y: 228.00 Mag: 24.700 Added Star: 23 - X: 49.00 Y: 57.00 Mag: 21.600 Added Star: 24 - X: 239.00 Y: 245.00 Mag: 20.700 Added Star: 25 - X: 20.00 Y: 89.00 Mag: 19.300 Added Star: 26 - X: 108.00 Y: 24.00 Mag: 18.100 Added Star: 27 - X: 128.00 Y: 105.00 Mag: 24.400 Added Star: 28 - X: 231.00 Y: 106.00 Mag: 20.100 Added Star: 29 - X: 156.00 Y: 12.00 Mag: 23.300 Added Star: 30 - X: 187.00 Y: 113.00 Mag: 22.300 Added Star: 31 - X: 79.00 Y: 177.00 Mag: 24.600 Added Star: 32 - X: 34.00 Y: 45.00 Mag: 19.400 Added Star: 33 - X: 12.00 Y: 35.00 Mag: 19.800 Added Star: 34 - X: 140.00 Y: 5.00 Mag: 19.200 Added Star: 35 - X: 133.00 Y: 162.00 Mag: 24.900 Added Star: 36 - X: 22.00 Y: 201.00 Mag: 20.100 Added Star: 37 - X: 54.00 Y: 143.00 Mag: 22.300 Added Star: 38 - X: 15.00 Y: 215.00 Mag: 21.100 Added Star: 39 - X: 157.00 Y: 44.00 Mag: 21.900 Added Star: 40 - X: 18.00 Y: 93.00 Mag: 21.200 Added Star: 41 - X: 102.00 Y: 196.00 Mag: 22.900 Added Star: 42 - X: 241.00 Y: 86.00 Mag: 24.100 Added Star: 43 - X: 228.00 Y: 47.00 Mag: 21.100 Added Star: 44 - X: 198.00 Y: 225.00 Mag: 21.400 Added Star: 45 - X: 215.00 Y: 149.00 Mag: 21.400 Added Star: 46 - X: 171.00 Y: 23.00 Mag: 23.800 Added Star: 47 - X: 239.00 Y: 167.00 Mag: 22.900 Added Star: 48 - X: 214.00 Y: 180.00 Mag: 22.300 Added Star: 49 - X: 40.00 Y: 44.00 Mag: 22.000 Added Star: 50 - X: 78.00 Y: 164.00 Mag: 21.300 Added Star: 51 - X: 144.00 Y: 125.00 Mag: 20.500 Added Star: 52 - X: 180.00 Y: 170.00 Mag: 23.400 Added Star: 53 - X: 61.00 Y: 165.00 Mag: 24.700 Added Star: 54 - X: 219.00 Y: 88.00 Mag: 24.700 Added Star: 55 - X: 95.00 Y: 107.00 Mag: 19.600 Added Star: 56 - X: 144.00 Y: 208.00 Mag: 22.700 Added Star: 57 - X: 153.00 Y: 9.00 Mag: 19.400 Added Star: 58 - X: 129.00 Y: 235.00 Mag: 18.300 Added Star: 59 - X: 216.00 Y: 43.00 Mag: 22.800 Added Star: 60 - X: 166.00 Y: 71.00 Mag: 22.700 Added Star: 61 - X: 230.00 Y: 56.00 Mag: 18.300 Added Star: 62 - X: 182.00 Y: 77.00 Mag: 21.600 Added Star: 63 - X: 39.00 Y: 128.00 Mag: 19.700 Added Star: 64 - X: 12.00 Y: 103.00 Mag: 23.000 Added Star: 65 - X: 158.00 Y: 180.00 Mag: 19.600 Added Star: 66 - X: 164.00 Y: 212.00 Mag: 24.200 Added Star: 67 - X: 56.00 Y: 57.00 Mag: 19.700 Added Star: 68 - X: 181.00 Y: 206.00 Mag: 22.700 Added Star: 69 - X: 116.00 Y: 32.00 Mag: 22.800 Added Star: 70 - X: 177.00 Y: 38.00 Mag: 19.300 Added Star: 71 - X: 122.00 Y: 64.00 Mag: 21.000 Added Star: 72 - X: 167.00 Y: 192.00 Mag: 18.200 Added Star: 73 - X: 174.00 Y: 206.00 Mag: 19.700 Added Star: 74 - X: 165.00 Y: 227.00 Mag: 21.800 Added Star: 75 - X: 102.00 Y: 238.00 Mag: 21.500 Added Star: 76 - X: 72.00 Y: 20.00 Mag: 21.000 Added Star: 77 - X: 111.00 Y: 16.00 Mag: 20.100 Added Star: 78 - X: 184.00 Y: 81.00 Mag: 21.800 Added Star: 79 - X: 25.00 Y: 109.00 Mag: 20.100 Added Star: 80 - X: 150.00 Y: 166.00 Mag: 21.400 Added Star: 81 - X: 62.00 Y: 239.00 Mag: 23.900 Added Star: 82 - X: 128.00 Y: 238.00 Mag: 20.800 Added Star: 83 - X: 219.00 Y: 9.00 Mag: 24.100 Added Star: 84 - X: 11.00 Y: 199.00 Mag: 18.000 Added Star: 85 - X: 6.00 Y: 122.00 Mag: 24.200 Added Star: 86 - X: 59.00 Y: 163.00 Mag: 22.800 Added Star: 87 - X: 202.00 Y: 104.00 Mag: 24.100 Added Star: 88 - X: 84.00 Y: 204.00 Mag: 18.500 Added Star: 89 - X: 240.00 Y: 222.00 Mag: 20.000 Added Star: 90 - X: 136.00 Y: 250.00 Mag: 24.400 Added Star: 91 - X: 19.00 Y: 169.00 Mag: 21.700 Added Star: 92 - X: 139.00 Y: 133.00 Mag: 18.500 Added Star: 93 - X: 205.00 Y: 156.00 Mag: 25.000 Added Star: 94 - X: 88.00 Y: 67.00 Mag: 21.200 Added Star: 95 - X: 187.00 Y: 82.00 Mag: 21.600 Added Star: 96 - X: 211.00 Y: 178.00 Mag: 24.700 Added Star: 97 - X: 120.00 Y: 170.00 Mag: 19.600 Added Star: 98 - X: 105.00 Y: 50.00 Mag: 21.200 Added Star: 99 - X: 173.00 Y: 21.00 Mag: 24.100
Let's check out the new image to see the additions.
(Again, don't freak out. The model psf is bad to make sure you'd see them since blinking is a little hard here)
realstars = pyfits.getdata(science_image)
artstars = pyfits.getdata('n8q624e8q_art1.fits')
fig = plt.figure()
fig.set_figheight(12) #inches
fig.set_figwidth(12) #inches
# Shows original science image
ax1 = fig.add_subplot(121)
ax1.imshow(realstars,vmin=-0.1,vmax=1, cmap=plt.cm.gray)
ax1.set_title('Before')
# Shows science image with fake stars
ax2 = fig.add_subplot(122)
ax2.imshow(artstars,vmin=-0.1,vmax=1, cmap=plt.cm.gray)
ax2.set_title('After')
ax2.set_xlim(0, 250)
ax2.set_ylim(0, 250)
(0, 250)
So, now you just have to:
Perry wrote a quick matching algorithm in python which you might try as an example (it's in the scripts directory)
import match
#lets make a quick list of numbers to match
xin=np.arange(10)
yin=np.arange(10)
print(xin,yin)
[0 1 2 3 4 5 6 7 8 9] [0 1 2 3 4 5 6 7 8 9]
xfake=np.arange(5)
yfake=np.arange(5)
print(xfake,yfake)
[0 1 2 3 4] [0 1 2 3 4]
good, nomatch = match.match(xin,yin,xfake,yfake,1) #match with a tolerance of 1
print(good)
print(nomatch)
[(0, 0), (1, 1), (2, 2), (3, 3), (4, 4)] [5, 6, 7, 8, 9]
#there is also a function which matches with a secondary criteria (the mags are similar)
inmags=np.arange(10) + 20.5
print(inmags)
[ 20.5 21.5 22.5 23.5 24.5 25.5 26.5 27.5 28.5 29.5]
outmags=np.arange(5)+20.5
#now lets change the first mag so it doesn't match
outmags[0]=100.
print(outmags)
[ 100. 21.5 22.5 23.5 24.5]
match.filtermags(good,inmags,outmags,0.5)
print(good)
print(nomatch)
[(1, 1), (2, 2), (3, 3), (4, 4)] [5, 6, 7, 8, 9]
#see the first data point was removed from the list....