2.2.4 Image surface brightness limit ¶

When your science is related to extended emission (like the example here) and you are presenting your results in a scientific conference, usually the first thing that someone will ask (if you do not explicitly say it!), is the dataset’s surface brightness limit (a standard measure of the noise level), and your target’s surface brightness (a measure of the signal, either in the center or outskirts, depending on context). For more on the basics of these important concepts please see Metameasurements on full input. So in this section of the tutorial, we will measure these values for the single-exposure SDSS image of the M51 group that we downloaded in Downloading and validating input data.

Before measuring the surface brightness limit, let’s see how reliable our detection was. In other words, let’s see how “clean” our noise is (after masking all detections, as described previously in Skewness caused by signal and its measurement):

$ aststatistics det-masked.fits --quantofmean
0.5111848629

Showing that the mean is indeed very close to the median, although just about 1 quantile larger. As we saw in NoiseChisel optimization, a very small residual signal still remains in the undetected regions and this very small difference is a quantitative measure of that undetected signal. It was up to you as an exercise to improve it, so we will continue with this dataset.

The surface brightness limit of the image can be measured from the masked image and the equation in Metameasurements on full input. Let’s do it for a \(3\sigma\) surface brightness limit over an area of \(25 \rm{arcsec}^2\):

$ nsigma=3
$ zeropoint=22.5
$ areaarcsec2=25
$ std=$(aststatistics det-masked.fits --sigclip-std)
$ pixarcsec2=$(astfits det-masked.fits --pixelscale --quiet \
                       | awk '{print $3*3600*3600}')
$ astarithmetic --quiet $nsigma $std x \
                $areaarcsec2 $pixarcsec2 x \
                sqrt / $zeropoint counts-to-mag
26.0241

The customizable steps above are good for any type of mask. For example, your field of view may contain a very deep part so you need to mask all the shallow parts as well as the detections before these steps. But when your image is flat (like this), there is a much simpler method to obtain the same value through MakeCatalog (when the standard deviation image is made by NoiseChisel). NoiseChisel has already calculated the minimum (MINSTD), maximum (MAXSTD) and median (MEDSTD) standard deviation within the tiles during its processing and has stored them as FITS keywords within the SKY_STD HDU. You can see them by piping all the keywords in this HDU into grep. In Grep, each ‘.’ represents one character that can be anything so M..STD will match all three keywords mentioned above.

$ astfits r_detected.fits --hdu=SKY_STD | grep 'M..STD'

The MEDSTD value is very similar to the standard deviation derived above, so we can safely use it instead of having to mask and run Statistics.

Fortunately, MakeCatalog will use this keyword and will report the dataset’s \(n\sigma\) surface brightness limit over a certain area as keywords in the output. However, as we saw above, these types of measurements are generic to the whole image, not to any specific object. In other words, they are measurements that are used to evaluate other measurements. Therefore, in MakeCatalog, we call these “meta-measures” which are more completely described in Metameasurements on full input. In summary, to activate these types of measurements, you need to use the --meta-measures option:

$ astmkcatalog r_detected.fits -hDETECTIONS --output=sbl.fits \
               --ids --meta-measures

Before looking into the measured surface brightness limits, let’s review an important point about this call to MakeCatalog: We are only concerned with the noise (not the signal), so we do not ask for any further measurements, because they can un-necessarily slow it down. However, MakeCatalog requires at least one column, so we will only ask for the --ids column (which does not need any measurement!). The output catalog will therefore have a single row and a single column, with 1 as its value⁵⁶.

Since all the meta-measurements get written in the first (0th) HDU, with the command below you can see all the keywords that were measured with the table. Notice the group of keywords that are under the “Noise-based metameasures” title.

$ astfits sbl.fits -h0

Above the noise-based meta-measure group, you see another group of keywords under the title “Input pixel grid and value properties”. Notice how the MMSTD has the same value as NoiseChisel’s MEDSTD above: this is used directly for both the measured surface brightness limit and the noise-based magnitude limit. Using MMSTD, MakeCatalog has determined the \(n\sigma\) surface brightness limiting magnitude in these header keywords.

The multiple of \(\sigma\), or \(n\), and area used are parameters that you should set at run-time as command-line options. Therefore you should look for those under the “Option values” group of keywords. is the value of the SBL-SIGMA keyword which you can change with --sbl-sigma. The surface brightness limiting magnitude within a pixel (SBL-SIGMA) and within a shape-agnostic area of SBL-AREA arcsec\(^2\) are stored in SBL.

You will notice that the two surface brightness limiting magnitudes above have values around 3 and 4 (which is not correct!). This is because we have not given a zero point magnitude to MakeCatalog, so it uses the default value of 0. SDSS image pixel values are calibrated in units of “nanomaggy” which are defined to have a zero point magnitude of 22.5⁵⁷. So with the first command below we give the zero point value and with the second we can see the surface brightness limiting magnitudes with the correct values (around 25 and 26)

$ astmkcatalog r_detected.fits -hDETECTIONS --output=sbl.fits \
               --ids --meta-measures --zeropoint=22.5
$ astfits sbl.fits -h0 | grep SBL

As you see from SBL-SIGMA and SBL-AREA, the default multiple of sigma is 3 and the default area is 100 arcsec\(^2\) (which is commonly used in the literature). But as mentioned above, you can easily change those in every run. In the manual example above, we assumed the multiple of sigma to be 3 and the area to be 25 arcsec\(^2\), so we just need to modify the area:

$ astmkcatalog r_detected.fits -hDETECTIONS --output=sbl.fits \
               --ids --meta-measures --zeropoint=22.5 --sbl-area=25

$ astfits sbl.fits -h0 | awk '$1=="SBL" {print $3}'
26.02296

You see that the value is identical to the custom surface brightness limiting magnitude we measured above (a difference of \(0.00114\) magnitudes is negligible and hundreds of times larger than the typical errors in the zero point magnitude or magnitude measurements). But it is much more easier to have MakeCatalog do this measurement, because these values will be appended (as keywords) into your final catalog of objects within that image.

$ std=$(aststatistics masked.fits --sigclip-std)
$ astfits noisechisel.fits -hSKY_STD --update=MEDSTD,$std

We have successfully measured the image’s \(3\sigma\) surface brightness limiting magnitude over 25 arcsec\(^2\). However, as discussed in Metameasurements on full input this value is just an extrapolation of the per-pixel standard deviation. Issues like correlated noise will cause the real noise over a large area to be different. So for a more robust measurement, let’s use the upper-limit magnitude of similarly sized region. For more on the upper-limit magnitude, see the respective item in Metameasurements on full input.

In summary, the upper-limit measurements involve randomly placing the footprint of an object in undetected parts of the image many times. This results in a random distribution of brightness measurements, the standard deviation of that distribution is then converted into magnitudes. To be comparable with the results above, let’s make a circular aperture that has an area of 25 arcsec\(^2\) (thus with a radius of \(2.82095\) arcsec).

zeropoint=22.5
r_arcsec=2.82095

## Convert the radius (in arcseconds) to pixels.
r_pixel=$(astfits r_detected.fits --pixelscale -q \
                  | awk '{print '$r_arcsec'/($1*3600)}')

## Make circular aperture at pixel (100,100) position is irrelevant.
echo "1 100 100 5 $r_pixel 0 0 1 1 1" \
     | astmkprof --background=r_detected.fits \
                 --clearcanvas --mforflatpix --type=uint8 \
                 --output=lab.fits

## Do the upper-limit measurement, ignoring all NoiseChisel's
## detections as a mask for the upper-limit measurements.
astmkcatalog lab.fits -h1 --zeropoint=$zeropoint -osbl.fits \
             --sbl-area=25 --sbl-sigma=3 --meta-measures \
             --valuesfile=r_detected.fits --valueshdu=INPUT-NO-SKY \
             --upmaskfile=r_detected.fits --upmaskhdu=DETECTIONS \
             --upnsigma=3 --checkuplim=1 --upnum=1000 \
             --ids --upperlimit-sb

The sbl.fits catalog now contains the upper-limit surface brightness for a circle with an area of 25 arcsec\(^2\). You can check the value with the command below, but the great thing is that now you have both of the measured surface brightness limit in the headers (discussed above), and the upper-limit surface brightness within the table. You can also add more profiles with different shapes and sizes if necessary (just make sure they do not overlap: their position is irrelevant!). Of course, you can also use --upperlimit-sb in your actual science objects and clumps to get an object-specific or clump-specific value.

$ asttable sbl.fits -cUPPERLIMIT_SB
25.9119

You will get a slightly different value from the command above. In fact, if you run the MakeCatalog command again and look at the measured upper-limit surface brightness, it will be slightly different with your first trial! Please try exactly the same MakeCatalog command above a few times to see how it changes.

This is because of the random factor in the upper-limit measurements: every time you run it, different random points will be checked, resulting in a slightly different distribution. You can increase the random scatter by increasing the number of random checks (for example, setting --upnum=100000, compared to 1000 in the command above). But this will be slower and the results will not be exactly reproducible. The only way to ensure you get an identical result later is to fix the random number generator function and seed like the command below⁵⁸. This is a very important point for the reproducibility of any statistical process involving random numbers, please see Generating random numbers.

export GSL_RNG_TYPE=ranlxs1
export GSL_RNG_SEED=1616493518
astmkcatalog lab.fits -h1 --zeropoint=$zeropoint -osbl.fits \
             --sbl-area=25 --sbl-sigma=3 --meta-measures \
             --valuesfile=r_detected.fits --valueshdu=INPUT-NO-SKY \
             --upmaskfile=r_detected.fits --upmaskhdu=DETECTIONS \
             --upnsigma=3 --checkuplim=1 --upnum=1000 \
             --ids --upperlimit-sb --envseed

But where do all the random apertures of the upper-limit measurement fall on the image? It is good to actually inspect their location to get a better understanding for the process and also detect possible bugs/biases. When MakeCatalog is run with the --checkuplim option, it will save all the random locations and their measured brightness as an extra table HDU in the output (called CHECK-UPPERLIMIT). With the first command below you can use Gnuastro’s asttable and astscript-ds9-region to convert the successful aperture locations into a DS9 region file, and with the second can load the region file into the detections and sky-subtracted image to visually see where they are.

## Create a DS9 region file from the check table (activated
## with '--checkuplim')
$ asttable sbl.fits -hCHECK-UPPERLIMIT --noblank=RANDOM_SUM \
           | astscript-ds9-region -c1,2 --mode=img \
                                  --radius=$r_pixel \
                                  --output=apers.reg

## Have a look at the regions in relation with NoiseChisel's
## detections.
$ astscript-fits-view r_detected.fits --ds9region=apers.reg

In this example, we were looking at a single-exposure image that has no correlated noise. Because of this, the surface brightness limit and the upper-limit surface brightness are very close. They will have a bigger difference on deep datasets with stronger correlated noise (that are the result of coadding many individual exposures). As an exercise, please try measuring the upper-limit surface brightness level and surface brightness limit for the deep HST data that we used in the previous tutorial (General program usage tutorial).