Deriving a Point-Spread Function in a Crowded Field¶
following Appendix III of Peter Stetson’s User’s Manual for DAOPHOT II¶
Using pydaophot
form astwro
python package¶
All italic text here have been taken from Stetson’s manual.
The only input file for this procedure is a FITS file containing
reference frame image. Here we use sample FITS form astwro package
(NGC6871 I filter 20s frame). Below we get filepath for this image, as
well as create instances of Daophot
and Allstar
classes -
wrappers around daophot
and allstar
respectively.
One should also provide daophot.opt
, photo.opt
and
allstar.opt
in apropiriete constructors. Here default, build in,
sample, opt
files are used.
from astwro.sampledata import fits_image
frame = fits_image()
Daophot
object creates temporary working directory (runner
directory), which is passed to Allstar
constructor to share.
from astwro.pydaophot import Daophot, Allstar
dp = Daophot(image=frame)
al = Allstar(dir=dp.dir)
Daophot got FITS file in construction, which will be automatically ATTACHed.
(1) Run FIND on your frame¶
Daophot FIND
parameters Number of frames averaged, summed
are
defaulted to 1,1
, below are provided for clarity.
res = dp.FInd(frames_av=1, frames_sum=1)
Check some results returned by FIND
, every method for daophot
command returns results object.
print ("{} pixels analysed, sky estimate {}, {} stars found.".format(res.pixels, res.sky, res.stars))
9640 pixels analysed, sky estimate 12.665, 4166 stars found.
Also, take a look into runner directory
!ls $dp.dir
63d38b_NGC6871.fits allstar.opt daophot.opt i.coo
We see symlinks to input image and opt
files, and i.coo
- result
of FIND
(2) Run PHOTOMETRY on your frame¶
Below we run photometry, providing explicitly radius of aperture A1
and IS
, OS
sky radiuses.
res = dp.PHotometry(apertures=[8], IS=35, OS=50)
List of stars generated by daophot commands, can be easily get as
astwro.starlist.Starlist
being essentially pandas.DataFrame
:
stars = res.photometry_starlist
Let’s check 10 stars with least A1 error (mag_err
column).
(pandas style)
stars.sort_values('mag_err').iloc[:10]
id | x | y | mag | sky | sky_err | sky_skew | mag_err | |
---|---|---|---|---|---|---|---|---|
id | ||||||||
2631 | 2631 | 982.57 | 733.50 | 12.430 | 12.626 | 2.27 | 0.08 | 0.0012 |
2387 | 2387 | 577.37 | 666.48 | 12.118 | 15.649 | 6.55 | 0.52 | 0.0012 |
391 | 391 | 702.67 | 102.05 | 12.533 | 12.755 | 2.45 | 0.08 | 0.0012 |
697 | 697 | 502.64 | 177.66 | 12.741 | 12.794 | 2.41 | 0.09 | 0.0014 |
879 | 879 | 1091.86 | 226.61 | 12.841 | 12.902 | 2.48 | 0.10 | 0.0014 |
926 | 926 | 1107.02 | 241.15 | 12.763 | 12.866 | 2.43 | 0.11 | 0.0014 |
2277 | 2277 | 1165.50 | 636.91 | 12.742 | 12.567 | 2.36 | 0.08 | 0.0014 |
3681 | 3681 | 935.70 | 1025.92 | 13.129 | 12.528 | 2.28 | 0.07 | 0.0017 |
1753 | 1753 | 223.25 | 481.61 | 13.170 | 12.513 | 2.18 | 0.03 | 0.0017 |
2364 | 2364 | 603.22 | 662.33 | 12.908 | 16.590 | 5.20 | 0.43 | 0.0018 |
(3) SORT the output from PHOTOMETRY¶
in order of increasing apparent magnitude decreasing stellar brightness with the renumbering feature. This step is optional but it can be more convenient than not.
SORT
command of daophor
is not implemented (yet) in
pydaohot
. But we do sorting by ourself.
sorted_stars = stars.sort_values('mag')
sorted_stars.renumber()
Here we write sorted list back info photometry file at default name (overwriting existing one), because it’s convenient to use default files in next commands.
dp.write_starlist(sorted_stars, 'i.ap')
'i.ap'
(4) PICK to generate a set of likely PSF stars¶
How many stars you want to use is a function of the degree of variation you expect and the frequency with which stars are contaminated by cosmic rays or neighbor stars. […]
pick_res = dp.PIck(faintest_mag=20, number_of_stars_to_pick=40)
If no error reported, symlink to image file (renamed to i.fits
), and
all daophot output files (i.*
) are in the working directory of
runner:
ls $dp.dir
63d38b_NGC6871.fits daophot.opt i.coo
allstar.opt i.ap i.lst
One may examine and improve i.lst
list of PSF stars. Or use
astwro.tools.gapick.py
to obtain list of PSF stars optimised by
genetic algorithm.
(5) Run PSF¶
tell it the name of your complete (sorted renumbered) aperture photometry file, the name of the file with the list of PSF stars, and the name of the disk file you want the point spread function stored in (the default should be fine) […]
If the frame is crowded it is probably worth your while to generate the first PSF with the “VARIABLE PSF” option set to -1 — pure analytic PSF. That way, the companions will not generate ghosts in the model PSF that will come back to haunt you later. You should also have specified a reasonably generous fitting radius — these stars have been preselected to be as isolated as possible and you want the best fits you can get. But remember to avoid letting neighbor stars intrude within one fitting radius of the center of any PSF star.
For illustration we will set VARIABLE PSF
option, before PSf()
dp.set_options('VARIABLE PSF', 2)
psf_res = dp.PSf()
(6) Run GROUP and NSTAR or ALLSTAR on your NEI file¶
If your PSF stars have many neighbors this may take some minutes of real time. Please be patient or submit it as a batch job and perform steps on your next frame while you wait.
We use allstar
. (GROUP
and NSTAR
command are not implemented
in current version of pydaophot
). We use prepared above Allstar
object: al
operating on the same runner dir that dp
.
As parameter we set input image (we haven’t do that on constructor), and
nei
file produced by PSf()
. We do not remember name i.psf
so
use psf_res.nei_file
property.
Finally we order allstar
to produce subtracted FITS .
alls_res = al.ALlstar(image_file=frame, stars=psf_res.nei_file, subtracted_image_file='is.fits')
All result
objects, has get_buffer()
method, useful to lookup
unparsed daophot
or allstar
output:
print (alls_res.get_buffer())
63d38b_NGC6871...
Picture size: 1250 1150
File with the PSF (default 63d38b_NGC6871.psf): Input file (default 63d38b_NGC6871.ap): File for results (default i.als): Name for subtracted image (default is):
915 stars. <<
I = iteration number
R = number of stars that remain
D = number of stars that disappeared
C = number of stars that converged
I R D C
1 915 0 0 <<
2 915 0 0 <<
3 915 0 0 <<
4 724 0 191 <<
5 385 0 530 <<
6 211 0 704 <<
7 110 0 805 <<
8 67 0 848 <<
9 40 0 875 <<
10 0 0 915
Finished i
Good bye.
(8) EXIT from DAOPHOT and send this new picture to the image display¶
Examine each of the PSF stars and its environs. Have all of the PSF stars subtracted out more or less cleanly, or should some of them be rejected from further use as PSF stars? (If so use a text editor to delete these stars from the LST file.) Have the neighbors mostly disappeared, or have they left behind big zits? Have you uncovered any faint companions that FIND missed?[…]
The absolute path to subtracted file (like for most output files) is available as result’s property:
sub_img = alls_res.subtracted_image_file
We can also generate region file for psf stars:
from astwro.starlist.ds9 import write_ds9_regions
reg_file_path = dp.file_from_runner_dir('lst.reg')
write_ds9_regions(pick_res.picked_starlist, reg_file_path)
# One can run ds9 directly from notebook:
!ds9 $sub_img -regions $reg_file_path
(9) Back in DAOPHOT II ATTACH the original picture and run SUBSTAR¶
specifying the file created in step (6) or in step (8f) as the stars to subtract, and the stars in the LST file as the stars to keep.
Lookup into runner dir:
ls $al.dir
63d38b_NGC6871.fits i.ap i.nei
allstar.opt i.coo i.psf
daophot.opt i.err is.fits
i.als i.lst lst.reg
sub_res = dp.SUbstar(subtract=alls_res.profile_photometry_file, leave_in=pick_res.picked_stars_file)
You have now created a new picture which has the PSF stars still in it but from which the known neighbors of these PSF stars have been mostly removed
(10) ATTACH the new star subtracted frame and repeat step (5) to derive a new point spread function¶
(11+…) Run GROUP NSTAR or ALLSTAR¶
for i in range(3):
print ("Iteration {}: Allstar chi: {}".format(i, alls_res.als_stars.chi.mean()))
dp.image = 'is.fits'
dp.PSf()
alls_res = al.ALlstar(image_file=frame, stars='i.nei')
dp.image = frame
dp.SUbstar(subtract='i.als', leave_in='i.lst')
print ("Final: Allstar chi: {}".format(alls_res.als_stars.chi.mean()))
Iteration 0: Allstar chi: 1.14670601093
Iteration 1: Allstar chi: 1.13409726776
Iteration 2: Allstar chi: 1.1332852459
Final: Allstar chi: 1.13326229508
Check last image with subtracted PSF stars neighbours.
!ds9 $dp.SUbstar_result.subtracted_image_file -regions $reg_file_path
Once you have produced a frame in which the PSF stars and their neighbors all subtract out cleanly, one more time through PSF should produce a point-spread function you can be proud of.
dp.image = 'is.fits'
psf_res = dp.PSf()
print ("PSF file: {}".format(psf_res.psf_file))
PSF file: /var/folders/kt/1jqvm3s51jd4qbxns7dc43rw0000gq/T/pydaophot_tmpBVHrtR/i.psf