2.3.4 Building outer part of PSF ¶

In Saturated pixels and Segment’s clumps, and One object for the whole detection, we described how to create a Segment clump and object map, while accounting for saturated stars and not having over-fragmentation of objects in the outskirts of stars. We are now ready to start building the extended PSF.

First we will build the outer parts of the PSF, so we want the brightest stars. You will see we have several bright stars in this very large field of view, but we do not yet have a feeling how many they are, and at what magnitudes. So let’s use Gnuastro’s Query program to find the magnitudes of the brightest stars (those brighter than g-magnitude 10 in Gaia data release 3, or DR3). For more on Query, see Query.

$ astquery gaia --dataset=dr3 --overlapwith=flat/67510.fits \
           --range=phot_g_mean_mag,-inf,10 \
           --output=flat/67510-bright.fits

Now, we can easily visualize the magnitude and positions of these stars using astscript-ds9-region and the command below (for more on this script, see SAO DS9 region files from table)

$ astscript-ds9-region flat/67510-bright.fits -cra,dec \
           --namecol=phot_g_mean_mag \
           --command="ds9 flat/67510.fits -zoom to fit -zscale"

You can see that we have several stars between magnitudes 6 to 10. Let’s use astscript-psf-select-stars in the command below to select the relevant stars in the image (the brightest; with a magnitude between 6 to 10). The advantage of using this script (instead of a simple --range in Table), is that it will also check distances to nearby stars and reject those that are too close (and not good for constructing the PSF). Since we have very bright stars in this very wide-field image, we will also increase the distance to nearby neighbors with brighter or similar magnitudes (the default value is 1 arcmin). To do this, we will set --mindistdeg=0.02, which corresponds to 1.2 arcmin. The details of the options for this script are discussed in Invoking astscript-psf-select-stars.

$ mkdir outer
$ astscript-psf-select-stars flat/67510.fits \
           --magnituderange=6,10 --mindistdeg=0.02 \
           --output=outer/67510-6-10.fits

Let’s have a look at the selected stars in the image (it is very important to visually check every step when you are first discovering a new dataset).

$ astscript-ds9-region outer/67510-6-10.fits -cra,dec \
           --namecol=phot_g_mean_mag \
           --command="ds9 flat/67510.fits -zoom to fit -zscale"

Now that the catalog of good stars is ready, it is time to construct the individual stamps from the catalog above. To create stamps, first, we need to crop a fixed-size box around each isolated star in the catalog. The contaminant objects in the crop should be masked and finally, the fluxes in these cropped images should be normalized. To do these, we will use astscript-psf-stamp (for more on this script see Invoking astscript-psf-stamp).

One of the most important parameters for this script is the normalization radii --normradii. This parameter defines a ring for the flux normalization of each star stamp. The normalization of the flux is necessary because each star has a different brightness, and consequently, it is crucial for having all the stamps with the same flux level in the same region. Otherwise the final stack of the different stamps would have no sense. Depending on the PSF shape, internal reflections, ghosts, saturated pixels, and other systematics, it would be necessary to choose the --normradii appropriately.

The selection of the normalization radii is something that requires a good understanding of the data. To do that, let’s use two useful parameters that will help us in the checking of the data: --tmpdir and --keeptmp;

• With --tmpdir=checking-normradii all temporary files, including the radial profiles, will be save in that directory (instead of an internally-created name).
• With --keeptmp we will not remove the temporal files, so it is possible to have a look at them (by default the temporary directory gets deleted at the end). It is necessary to specify the --normradii even if we do not know yet the final values. Otherwise the script will not generate the radial profile.

As a consequence, in this step we put the normalization radii equal to the size of the stamps. By doing this, the script will generate the radial profile of the entire stamp. In this particular step we set it to --normradii=500,510. We also use the --nocentering option to disable sub-pixel warping in this phase (it is only relevant for the central part of the PSF). Furthermore, since there are several stars, we iterate over each row of the catalog using a while loop.

$ counter=1
$ mkdir finding-normradii
$ asttable outer/67510-6-10.fits \
           | while read -r ra dec mag; do
               astscript-psf-stamp label/67510-seg.fits \
                    --mode=wcs \
                    --nocentering \
                    --center=$ra,$dec \
                    --normradii=500,510 \
                    --widthinpix=1000,1000 \
                    --segment=label/67510-seg.fits \
                    --output=finding-normradii/$counter.fits \
                    --tmpdir=finding-normradii --keeptmp; \
               counter=$((counter+1)); \
             done

First let’s have a look at all the masked postage stamps of the cropped stars. Once they all open, feel free to zoom-in, they are all matched and locked. It is always good to check the different stamps to ensure the quality and possible two dimensional features that are difficult to detect from the radial profiles (such as ghosts and internal reflections).

$ astscript-fits-view finding-normradii/cropped-masked*.fits

If everything looks good in the image, let’s open all the radial profiles and visually check those with the command below. Note that astscript-fits-view calls the topcat graphic user interface (GUI) program to visually inspect (plot) tables. If you do not already have it, see TOPCAT.

$ astscript-fits-view finding-normradii/rprofile*.fits

After some study of this data, we could say that a good normalization ring is those pixels between R=20 and R=30 pixels. Such a ring ensures having a high number of pixels so the estimation of the flux normalization will be robust. Also, at such distance from the center the signal to noise is high and there are not obvious features that can affect the normalization. Note that the profiles are different because we are considering a wide range of magnitudes, so the fainter stars are much more noisy. However, in this tutorial we will keep these stars in order to have a higher number of stars for the outer part. In a real case scenario, we should look for stars with a much more similar brightness (smaller range of magnitudes) to not lose signal to noise as a consequence of the inclusion of fainter stars.

$ rm -r finding-normradii
$ counter=1
$ mkdir outer/stamps
$ asttable outer/67510-6-10.fits \
           | while read -r ra dec mag; do
               astscript-psf-stamp label/67510-seg.fits \
                    --mode=wcs \
                    --nocentering \
                    --center=$ra,$dec \
                    --normradii=20,30 \
                    --widthinpix=1000,1000 \
                    --segment=label/67510-seg.fits \
                    --output=outer/stamps/67510-$counter.fits; \
               counter=$((counter+1)); \
             done

After the stamps are created, we need to stack them together with a simple Arithmetic command (see Stacking operators). The stack is done using the sigma-clipped mean operator that will preserve more of the signal, while rejecting outliers (more than \(3\sigma\) with a tolerance of \(0.2\), for more on sigma-clipping see Sigma clipping). Just recall that we need to specify the number of inputs into the stacking operators, so we are reading the list of images and counting them as separate variables before calling Arithmetic.

$ numimgs=$(echo outer/stamps/*.fits | wc -w)
$ astarithmetic outer/stamps/*.fits $numimgs 3 0.2 sigclip-mean \
                -g1 --output=outer/stack.fits --wcsfile=none

Did you notice the --wcsfile=none option above? With it, the stacked image no longer has any WCS information. This is natural, because the stacked image does not correspond to any specific region of the sky any more.

Let’s compare this stacked PSF with the images of the individual stars that were used to create it. You can clearly see that the number of masked pixels is significantly decreased and the PSF is much more “cleaner”.

$ astscript-fits-view outer/stack.fits outer/stamps/*.fits

However, the saturation in the center still remains. Also, because we did not have too many images, some regions still are very noisy. If we had more bright stars in our selected magnitude range, we could have filled those outer remaining patches. In a large survey like J-PLUS (that we are using here), you can simply look into other fields that were observed soon before/after the image ID 67510 that we used here (to have a similar PSF) and get more stars in those images to add to these. In fact, the J-PLUS DR2 image ID of the field above was intentionally preserved during the steps above to show how easy it is to use images from other fields and blend them all into the output PSF.