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.


In [1]:
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.


In [2]:
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.


In [3]:
res = dp.FInd(frames_av=1, frames_sum=1)

Check some results returned by FIND, every method for daophot command returns results object.


In [4]:
print ("{} pixels analysed, sky estimate {}, {} stars found.".format(res.pixels, res.sky, res.stars))


9640 pixels analysed, sky estimate 12.665, 4166 stars found.

In [ ]:

Also, take a look into runner directory


In [5]:
!ls -lt $dp.dir


total 536
lrwxr-xr-x  1 michal  staff      60 Jun 26 18:25 63d38b_NGC6871.fits -> /Users/michal/projects/astwro/astwro/sampledata/NGC6871.fits
lrwxr-xr-x  1 michal  staff      65 Jun 26 18:25 allstar.opt -> /Users/michal/projects/astwro/astwro/pydaophot/config/allstar.opt
lrwxr-xr-x  1 michal  staff      65 Jun 26 18:25 daophot.opt -> /Users/michal/projects/astwro/astwro/pydaophot/config/daophot.opt
-rw-r--r--  1 michal  staff  258438 Jun 26 18:25 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.


In [6]:
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:


In [7]:
stars = res.photometry_starlist

Let's check 10 stars with least A1 error (mag_err column). (pandas style)


In [8]:
stars.sort_values('mag_err').iloc[:10]


Out[8]:
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.


In [9]:
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.


In [10]:
dp.write_starlist(sorted_stars, 'i.ap')


Out[10]:
'i.ap'

In [11]:
!head -n20 $dp.PHotometry_result.photometry_file


 NL    NX    NY  LOWBAD HIGHBAD  THRESH     AP1  PH/ADU  RNOISE    FRAD 
  2  1250  1150    -3.9 31000.0    5.81    8.00    9.00    1.70    6.00


      1  577.370  666.480   12.118
        15.649  6.55  0.52  0.0012

      2  982.570  733.500   12.430
        12.626  2.27  0.08  0.0012

      3  702.670  102.050   12.533
        12.755  2.45  0.08  0.0012

      4  603.270  675.390   12.727
        16.515  7.82  0.58  0.0020

      5  502.640  177.660   12.741
        12.794  2.41  0.09  0.0014

      6 1165.500  636.910   12.742

In [12]:
dp.PHotometry_result.photometry_file


Out[12]:
'/var/folders/kt/1jqvm3s51jd4qbxns7dc43rw0000gq/T/pydaophot_tmpDu5p8c/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. [...]


In [13]:
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:


In [14]:
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()


In [15]:
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 .


In [16]:
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:


In [17]:
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:


In [18]:
sub_img = alls_res.subtracted_image_file

We can also generate region file for psf stars:


In [19]:
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)

In [20]:
# 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:


In [21]:
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

In [22]:
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


In [23]:
for i in range(3):
    print ("Iteration {}: Allstar chi: {}".format(i, alls_res.als_stars.chi.mean()))
    dp.image = 'is.fits'
    respsf = dp.PSf()
    print ("Iteration {}: PSF chi: {}".format(i, respsf.chi))
    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 0: PSF chi: 0.0249
Iteration 1: Allstar chi: 1.13409726776
Iteration 1: PSF chi: 0.0249
Iteration 2: Allstar chi: 1.1332852459
Iteration 2: PSF chi: 0.0249
Final:       Allstar chi: 1.13326229508

In [37]:
alls_res.als_stars


Out[37]:
id x y mag mag_err sky iter chi sharp
id
540 540 1151.854 75.195 17.2302 0.0179 12.584 4.0 0.989 0.012
3499 3499 145.037 57.111 20.2348 0.1469 12.660 4.0 1.100 0.057
928 928 208.005 103.241 18.0819 0.0308 12.590 4.0 1.105 -0.048
3 3 702.673 102.077 12.4965 0.0069 12.565 4.0 1.502 0.031
992 992 688.449 120.159 18.3611 0.0300 12.682 4.0 0.891 0.098
2814 2814 744.407 108.595 20.1617 0.1286 12.609 4.0 1.027 -0.108
657 657 1152.035 115.646 18.1858 0.0258 12.469 4.0 0.862 0.054
2913 2913 1098.630 133.994 20.6351 0.2112 12.630 4.0 1.138 -0.085
930 930 1093.865 150.934 18.4099 0.0426 12.491 4.0 1.252 0.067
3512 3512 718.902 125.976 20.0723 0.1243 12.619 4.0 1.076 0.112
566 566 484.442 134.669 17.3278 0.0223 12.617 4.0 1.161 -0.006
2863 2863 1153.645 137.794 20.4459 0.1453 12.598 4.0 0.917 -0.156
346 346 707.580 145.109 16.8294 0.0159 12.625 4.0 1.031 -0.046
2269 2269 466.273 146.019 19.8062 0.1285 12.611 4.0 1.380 -0.582
257 257 455.804 147.555 16.4950 0.0244 12.654 4.0 1.805 -0.038
1857 1857 547.820 173.936 19.2825 0.0681 12.670 4.0 1.101 0.132
5 5 502.664 177.672 12.7116 0.0049 12.532 4.0 1.005 0.027
2323 2323 1203.448 150.666 19.5926 0.0862 12.567 4.0 1.093 -0.021
323 323 1120.644 173.486 16.7282 0.0186 12.588 4.0 1.274 0.066
2528 2528 481.494 172.071 20.6795 0.2050 12.547 4.0 1.012 -0.781
878 878 453.996 189.749 18.3533 0.0354 12.577 4.0 1.068 -0.094
1791 1791 1173.956 198.837 18.7949 0.0448 12.656 4.0 1.014 -0.081
1191 1191 1132.810 198.778 19.1139 0.0559 12.580 4.0 1.024 0.109
95 95 1121.696 200.301 15.3751 0.0240 12.584 4.0 2.573 0.004
2704 2704 1132.567 209.790 19.8789 0.1062 12.549 4.0 1.089 0.092
1789 1789 1120.242 216.888 19.6838 0.1006 12.586 4.0 1.201 0.429
2174 2174 1100.309 200.454 19.1560 0.0643 12.570 4.0 1.135 -0.168
345 345 549.124 208.800 16.8098 0.0179 12.799 4.0 1.166 0.028
7 7 1107.026 241.171 12.7259 0.0051 12.561 4.0 1.049 0.032
2502 2502 771.897 265.143 20.2399 0.1426 12.701 4.0 1.076 0.143
... ... ... ... ... ... ... ... ... ...
533 533 1130.590 171.811 18.1224 0.0774 12.663 10.0 2.684 0.939
2643 2643 465.721 190.288 22.0037 0.7000 12.611 10.0 1.102 -0.308
134 134 1121.832 207.707 19.8053 0.1211 12.669 10.0 1.236 -0.163
529 529 767.248 281.278 20.1183 0.1376 12.712 10.0 1.035 0.022
59 59 796.890 307.708 18.0256 0.0318 12.608 10.0 1.027 0.201
70 70 58.475 323.981 19.0317 0.0668 12.486 10.0 1.037 0.166
448 448 727.441 328.938 20.5069 0.1767 12.540 10.0 1.048 -0.107
3559 3559 485.145 395.435 20.6073 0.1848 12.420 10.0 0.951 -0.353
94 94 440.254 444.983 19.8242 0.1087 12.626 10.0 0.997 -0.158
215 215 691.102 476.562 19.8805 0.1053 13.857 10.0 1.010 -0.342
2254 2254 706.034 485.949 21.2466 0.2909 13.770 10.0 0.879 -0.285
65 65 603.545 500.642 18.6272 0.0415 13.576 10.0 1.014 0.334
444 444 628.324 507.129 19.3640 0.2163 14.205 10.0 3.192 1.826
465 465 695.796 531.668 19.0811 0.1558 13.996 10.0 2.946 1.061
202 202 686.225 532.391 18.9456 0.3005 14.103 10.0 6.208 1.546
310 310 625.648 597.763 19.7984 0.1080 14.116 10.0 1.034 -0.468
3613 3613 1133.147 621.466 21.9036 0.5948 12.357 10.0 1.064 -2.376
188 188 421.221 643.606 18.8318 0.0584 12.907 10.0 1.172 0.373
252 252 619.800 649.275 19.4646 0.2096 15.520 10.0 2.810 1.647
491 491 636.751 657.274 19.9176 0.1930 14.739 10.0 1.767 1.809
9 9 602.951 662.130 13.5473 0.1348 17.268 10.0 20.256 0.141
4 4 603.273 675.388 12.6545 0.0107 18.295 10.0 2.246 -0.012
1766 1766 754.723 724.891 20.5705 0.1741 12.958 10.0 0.988 0.620
987 987 788.989 745.930 20.4957 0.1565 12.984 10.0 0.875 0.082
50 50 254.283 794.702 19.3301 0.0730 12.371 10.0 0.945 -0.103
385 385 1137.059 809.156 20.1285 0.1247 12.600 10.0 1.010 -0.243
702 702 1161.335 821.090 20.0755 0.1252 12.567 10.0 1.078 -0.022
216 216 797.846 1042.250 20.8959 0.2613 12.607 10.0 1.060 1.278
372 372 844.859 1073.881 20.8228 0.2350 12.406 10.0 1.102 -0.004
824 824 792.384 1079.963 18.3105 0.0383 12.673 10.0 1.123 0.153

915 rows × 9 columns

Check last image with subtracted PSF stars neighbours.


In [24]:
!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.


In [25]:
dp.image = 'is.fits'
psf_res = dp.PSf()
print ("PSF file: {}".format(psf_res.psf_file))


PSF file: /var/folders/kt/1jqvm3s51jd4qbxns7dc43rw0000gq/T/pydaophot_tmpDu5p8c/i.psf

In [ ]: