2.3.7 Subtracting the PSF ¶

Previously (in Uniting the different PSF components) we constructed a full PSF, from the central pixel to a radius of 500 pixels. Now, let’s use the PSF to subtract the scattered light from each individual star in the image.

By construction, the pixel values of the PSF came from the normalization of the individual stamps (that were created for stars of different magnitudes). As a consequence, it is necessary to compute a scale factor to fit that PSF image to each star. This is done with the same astscript-psf-scale-factor command that we used previously in Uniting the different PSF components. The difference is that now we are not aiming to join two different PSF parts but looking for the necessary scale factor to match the star with the PSF. Afterwards, we will use astscript-psf-subtract for placing the PSF image at the desired coordinates within the same pixel grid as the image. Finally, once the stars have been modeled by the PSF, we will subtract it.

First, let’s start with a single star. Later, when the basic idea has been explained, we will generalize the method for any number of stars. With the following command we obtain the coordinates (RA and DEC) and magnitude of the brightest star in the image (which is on the top edge of the image):

$ mkdir single-star
$ center=$(asttable flat/67510-bright.fits --sort phot_g_mean_mag \
                    --column=ra,dec --head 1 \
                    | awk '{printf "%s,%s", $1, $2}')
$ echo $center

With the center position of that star, let’s obtain the flux factor using the same normalization ring we used for the creation of the outer part of the PSF:

$ scale=$(astscript-psf-scale-factor label/67510-seg.fits \
                   --mode=wcs --quiet \
                   --psf=psf.fits \
                   --center=$center \
                   --normradii=10,15 \
                   --segment=label/67510-seg.fits)

Now we have all the information necessary to model the star using the PSF: the position on the sky and the flux factor. Let’s use this data with the script astscript-psf-subtract for modeling this star and have a look with DS9.

$ astscript-psf-subtract label/67510-seg.fits \
           --mode=wcs \
           --psf=psf.fits \
           --scale=$scale \
           --center=$center \
           --output=single-star/subtracted.fits

$ astscript-fits-view label/67510-seg.fits single-star/subtracted.fits \
           --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 4"

You will notice that there is something wrong with this “subtraction”! The box of the extended PSF is clearly visible! The sky noise under the box is clearly larger than the rest of the noise in the image. Before reading on, please try to think about the cause of this yourself.

To understand the cause, let’s look at the scale factor, the number of stamps used to build the outer part (and its square root):

$ echo $scale
$ ls outer/stamps/*.fits | wc -l
$ ls outer/stamps/*.fits | wc -l | awk '{print sqrt($1)}'

You see that the scale is almost 19! As a result, the PSF has been multiplied by 19 before being subtracted. However, the outer part of the PSF was created with only a handful of star stamps. When you stack \(N\) images, the stack’s signal-to-noise ratio (S/N) improves by \(\sqrt{N}\). We had 8 images for the outer part, so the S/N has only improved by a factor of just under 3! When we multiply the final stacked PSF with 19, we are also scaling up the noise by that same factor (most importantly: in the outer most regions where there is almost no signal). So the stacked image’s noise-level is \(19/3=6.3\) times larger than the noise of the input image. This terrible noise-level is what you clearly see as the footprint of the PSF.

To confirm this, let’s use the commands below to subtract the faintest of the bright-stars catalog (note the use of --tail when finding the central position). You will notice that the scale factor (\(\sim1.3\)) is now smaller than 3. So when we multiply the PSF with this factor, the PSF’s noise level is lower than our input image and we should not see any footprint like before. Note also that we are using a larger zoom factor, because this star is smaller in the image.

$ center=$(asttable flat/67510-bright.fits --sort phot_g_mean_mag \
                    --column=ra,dec --tail 1 \
                    | awk '{printf "%s,%s", $1, $2}')

$ scale=$(astscript-psf-scale-factor label/67510-seg.fits \
                   --mode=wcs --quiet \
                   --psf=psf.fits \
                   --center=$center \
                   --normradii=10,15 \
                   --segment=label/67510-seg.fits)
$ echo $scale

$ astscript-psf-subtract label/67510-seg.fits \
           --mode=wcs \
           --psf=psf.fits \
           --scale=$scale \
           --center=$center \
           --output=single-star/subtracted.fits

$ astscript-fits-view label/67510-seg.fits single-star/subtracted.fits \
           --ds9center=$center --ds9mode=wcs --ds9extra="-zoom 10"

In a large survey like J-PLUS, it is easy to use more and more bright stars from different pointings (ideally with similar FWHM and similar telescope properties⁵⁹) to improve the S/N of the PSF. As explained before, we designed the output files of this tutorial with the 67510 (which is this image’s pointing label in J-PLUS) where necessary so you see how easy it is to add more pointings to use in the creation of the PSF.

Let’s consider now more than one single star. We should have two things in mind:

• The brightest (subtract-able, see the point below) star should be the first star to be subtracted. This is because of its extended wings which may affect the scale factor of nearby stars. So we should sort the catalog by magnitude and come down from the brightest.
• We should only subtract stars where the scale factor is less than the S/N of the PSF (in relation to the data).

Since it can get a little complex, it is easier to implement this step as a script (that is heavily commented for you to easily understand every step; especially if you put it in a good text editor with color-coding!). You will notice that script also creates a .log file, which shows which star was subtracted and which one was not (this is important, and will be used below!).

#!/bin/bash

# Abort the script on first error.
set -e

# ID of image to subtract stars from.
imageid=67510

# Get S/N level of the final PSF in relation to the actual data:
snlevel=$(ls outer/stamps/*.fits | wc -l | awk '{print sqrt($1)}')

# Put a copy of the image we want to subtract the PSF from in the
# final file (this will be over-written after each subtraction).
subtracted=subtracted/$imageid.fits
cp label/$imageid-seg.fits $subtracted

# Name of log-file to keep status of the subtraction of each star.
logname=subtracted/$imageid.log
echo "# Column 1: RA   [deg, f64] Right ascension of star." >  $logname
echo "# Column 2: Dec  [deg, f64] Declination of star."     >> $logname
echo "# Column 3: Stat [deg, f64] Status (1: subtracted)"   >> $logname

# Go over each item in the bright star catalog:
asttable flat/67510-bright.fits -cra,dec --sort phot_g_mean_mag  \
    | while read -r ra dec; do

    # Put a comma between the RA/Dec to pass to options.
    center=$(echo $ra $dec | awk '{printf "%s,%s", $1, $2}')

    # Calculate the scale value
    scale=$(astscript-psf-scale-factor label/67510-seg.fits \
                   --mode=wcs --quiet\
                   --psf=psf.fits \
                   --center=$center \
                   --normradii=10,15 \
                   --segment=label/67510-seg.fits)

    # Subtract this star if the scale factor is less than the S/N
    # level calculated above.
    check=$(echo $snlevel $scale \
                | awk '{if($1>$2) c="good"; else c="bad"; print c}')
    if [ $check = good ]; then

        # A temporary file to subtract this star.
        subtmp=subtracted/$imageid-tmp.fits

        # Subtract this star from the image where all previous stars
        # were subtracted.
        astscript-psf-subtract $subtracted \
                 --mode=wcs \
                 --psf=psf.fits \
                 --scale=$scale \
                 --center=$center \
                 --output=$subtmp

        # Rename the temporary subtracted file to the final one:
        mv $subtmp $subtracted

        # Keep the status for this star.
        status=1
    else
        # Let the user know this star did not work, and keep the status
        # for this star.
        echo "$center: $scale is larger than $snlevel"
        status=0
    fi

    # Keep the status in a log file.
    echo "$ra $dec $status" >> $logname
done

Copy the contents above into a file called subtract-psf-from-cat.sh and run the following commands. Just note that in the script above, we assumed the output is written in the subtracted/, directory, so we will first make that.

$ mkdir subtracted
$ chmod +x subtract-psf-from-cat.sh
$ ./subtract-psf-from-cat.sh

$ astscript-fits-view label/67510-seg.fits subtracted/67510.fits

Can you visually find the stars that have been subtracted? Its a little hard, is not it? This shows that you done a good job this time (the sky-noise is not significantly affected)! So let’s subtract the actual image from the PSF-subtracted image to see the scattered light field of the subtracted stars. With the second command below we will zoom into the brightest subtracted star, but of course feel free to zoom-out and inspect the others also.

$ astarithmetic label/67510-seg.fits subtracted/67510.fits - \
                --output=scattered-light.fits -g1

$ center=$(asttable subtracted/67510.log --equal=Stat,1 --head=1 \
                    -cra,dec | awk '{printf "%s,%s", $1, $2}')

$ astscript-fits-view label/67510-seg.fits subtracted/67510.fits \
                      scattered-light.fits \
                      --ds9center=$center --ds9mode=wcs \
                      --ds9extra="-scale limits -0.5 1.5 -match scale" \
                      --ds9extra="-lock scale yes -zoom 10" \
                      --ds9extra="-tile mode column"

## We can always make it easily, so let's remove this.
$ rm scattered-light.fits

You will probably have noticed that in the scattered light field there are some patches that correspond to the saturation of the stars. Since we obtained the scattered light field by subtracting PSF-subtracted image from the original image, it is natural that we have such saturated regions. To solve such inconvenience, this script also has an option to not make the subtraction of the PSF but to give as output the modeled star. For doing that, it is necessary to run the script with the option --modelonly. We encourage the reader to obtain such scattered light field model. In some scenarios it could be interesting having such way of correcting the PSF. For example, if there are many faint stars that can be modeled at the same time because their flux do not affect each other. In such situation, the task could be easily parallelized without having to wait to model the brighter stars before the fainter ones. At the end, once all stars have been modeled, a simple Arithmetic command could be used to sum the different modeled-PSF stamps to obtain the entire scattered light field.

In general you see that the subtraction has been done nicely and almost all the extended wings of the PSF have been subtracted. The central regions of the stars are not perfectly subtracted:

• Some may get too dark at the center. This may be due to the non-linearity of the CCD counting (as discussed previously in Uniting the different PSF components).
• Others may have a strong gradient: one side is too positive and one side is too negative (only in the very central few pixels). This is due to the non-accurate positioning: most probably this happens because of imperfect astrometry.

Note also that during this process we assumed that the PSF does not vary with the CCD position or any other parameter. In other words, we are obtaining an averaged PSF model from a few star stamps that are naturally different, and this also explains the residuals on each subtracted star.

We let as an interesting exercise the modeling and subtraction of other stars, for example, the non saturated stars of the image. By doing this, you will notice that in the core region the residuals are different compared to the residuals of brighter stars that we have obtained.

In general, in this tutorial we have showed how to deal with the most important challenges for constructing an extended PSF. Each image or dataset will have its own particularities that you will have to take into account when constructing the PSF.