gnuastro (0.22)

Browse

Build Log

Usage

This is gnuastro.info, produced by makeinfo version 7.1 from
gnuastro.texi.

This book documents version 0.22 of the GNU Astronomy Utilities
(Gnuastro). Gnuastro provides various programs and libraries for
astronomical data manipulation and analysis.

Permission is granted to copy, distribute and/or modify this
document under the terms of the GNU Free Documentation License,
Version 1.3 or any later version published by the Free Software
Foundation; with no Invariant Sections, no Front-Cover Texts, and
no Back-Cover Texts. A copy of the license is included in the
section entitled "GNU Free Documentation License".
INFO-DIR-SECTION Astronomy
START-INFO-DIR-ENTRY
* Gnuastro: (gnuastro). GNU Astronomy Utilities.
* libgnuastro: (gnuastro)Gnuastro library. Full Gnuastro library doc.

* help-gnuastro: (gnuastro)help-gnuastro mailing list. Getting help.

* bug-gnuastro: (gnuastro)Report a bug. How to report bugs

* Arithmetic: (gnuastro)Arithmetic. Arithmetic operations on pixels.
* astarithmetic: (gnuastro)Invoking astarithmetic. Options to Arithmetic.

* BuildProgram: (gnuastro)BuildProgram. Compile and run programs using Gnuastro's library.
* astbuildprog: (gnuastro)Invoking astbuildprog. Options to BuildProgram.

* ConvertType: (gnuastro)ConvertType. Convert different file types.
* astconvertt: (gnuastro)Invoking astconvertt. Options to ConvertType.

* Convolve: (gnuastro)Convolve. Convolve an input file with kernel.
* astconvolve: (gnuastro)Invoking astconvolve. Options to Convolve.

* CosmicCalculator: (gnuastro)CosmicCalculator. For cosmological params.
* astcosmiccal: (gnuastro)Invoking astcosmiccal. Options to CosmicCalculator.

* Crop: (gnuastro)Crop. Crop region(s) from image(s).
* astcrop: (gnuastro)Invoking astcrop. Options to Crop.

* Fits: (gnuastro)Fits. View and manipulate FITS extensions and keywords.
* astfits: (gnuastro)Invoking astfits. Options to Fits.

* MakeCatalog: (gnuastro)MakeCatalog. Make a catalog from labeled image.
* astmkcatalog: (gnuastro)Invoking astmkcatalog. Options to MakeCatalog.

* MakeProfiles: (gnuastro)MakeProfiles. Make mock profiles.
* astmkprof: (gnuastro)Invoking astmkprof. Options to MakeProfiles.

* Match: (gnuastro)Match. Match two separate catalogs.
* astmatch: (gnuastro)Invoking astmatch. Options to Match.

* NoiseChisel: (gnuastro)NoiseChisel. Detect signal in noise.
* astnoisechisel: (gnuastro)Invoking astnoisechisel. Options to NoiseChisel.

* Segment: (gnuastro)Segment. Segment detections based on signal structure.
* astsegment: (gnuastro)Invoking astsegment. Options to Segment.

* Query: (gnuastro)Query. Access remote databases for downloading data.
* astquery: (gnuastro)Invoking astquery. Options to Query.

* Statistics: (gnuastro)Statistics. Get image Statistics.
* aststatistics: (gnuastro)Invoking aststatistics. Options to Statistics.

* Table: (gnuastro)Table. Read and write FITS binary or ASCII tables.
* asttable: (gnuastro)Invoking asttable. Options to Table.

* Warp: (gnuastro)Warp. Warp a dataset to a new grid.
* astwarp: (gnuastro)Invoking astwarp. Options to Warp.

* astscript: (gnuastro)Installed scripts. Gnuastro's installed scripts.
* astscript-ds9-region: (gnuastro)Invoking astscript-ds9-region. Options to this script
* astscript-fits-view: (gnuastro)Invoking astscript-fits-view. Options to this script
* astscript-pointing-simulate: (gnuastro)Invoking astscript-pointing-simulate. Options to this script
* astscript-psf-scale-factor: (gnuastro)Invoking astscript-psf-scale-factor. Options to this script
* astscript-psf-select-stars: (gnuastro)Invoking astscript-psf-select-stars. Options to this script
* astscript-psf-stamp: (gnuastro)Invoking astscript-psf-stamp. Options to this script
* astscript-psf-subtract: (gnuastro)Invoking astscript-psf-subtract. Options to this script
* astscript-psf-unite: (gnuastro)Invoking astscript-psf-unite. Options to this script
* astscript-radial-profile: (gnuastro)Invoking astscript-radial-profile. Options to this script
* astscript-sort-by-night: (gnuastro)Invoking astscript-sort-by-night. Options to this script
* astscript-zeropoint: (gnuastro)Invoking astscript-zeropoint. Options to this script
END-INFO-DIR-ENTRY

File: gnuastro.info, Node: Citing and acknowledging Gnuastro, Prev: Writing scripts to automate the steps, Up: General program usage tutorial

2.1.23 Citing and acknowledging Gnuastro
----------------------------------------

In conclusion, we hope this extended tutorial has been a good starting
point to help in your exciting research. If this book or any of the
programs in Gnuastro have been useful for your research, please cite the
respective papers, and acknowledge the funding agencies that made all of
this possible. Without citations, we will not be able to secure future
funding to continue working on Gnuastro or improving it, so please take
software citation seriously (for all the scientific software you use,
not just Gnuastro).

To help you in this, all Gnuastro programs have a ‘--cite’ option to
facilitate the citation and acknowledgment. Just note that it may be
necessary to cite additional papers for different programs, so please
try it out on all the programs that you used, for example:

$ astmkcatalog --cite
$ astnoisechisel --cite

File: gnuastro.info, Node: Detecting large extended targets, Next: Building the extended PSF, Prev: General program usage tutorial, Up: Tutorials

2.2 Detecting large extended targets
====================================

The outer wings of large and extended objects can sink into the noise
very gradually and can have a large variety of shapes (for example, due
to tidal interactions). Therefore separating the outer boundaries of
the galaxies from the noise can be particularly tricky. Besides causing
an under-estimation in the total estimated brightness of the target,
failure to detect such faint wings will also cause a bias in the noise
measurements, thereby hampering the accuracy of any measurement on the
dataset. Therefore even if they do not constitute a significant
fraction of the target's light, or are not your primary target, these
regions must not be ignored. In this tutorial, we will walk you through
the strategy of detecting such targets using *note NoiseChisel::.

*Do not start with this tutorial:* If you have not already completed
*note General program usage tutorial::, we strongly recommend going
through that tutorial before starting this one. Basic features like
access to this book on the command-line, the configuration files of
Gnuastro's programs, benefiting from the modular nature of the programs,
viewing multi-extension FITS files, or using NoiseChisel's outputs are
discussed in more detail there.

We will try to detect the faint tidal wings of the beautiful M51
group(1) in this tutorial. We will use a dataset/image from the public
Sloan Digital Sky Survey (http://www.sdss.org/), or SDSS. Due to its
more peculiar low surface brightness structure/features, we will focus
on the dwarf companion galaxy of the group (or NGC 5195).

* Menu:

* Downloading and validating input data:: How to get and check the input data.
* NoiseChisel optimization:: Detect the extended and diffuse wings.
* Skewness caused by signal and its measurement:: How signal changes the distribution.
* Image surface brightness limit:: Standards to quantify the noise level.
* Achieved surface brightness level:: Calculate the outer surface brightness.
* Extract clumps and objects:: Find sub-structure over the detections.

---------- Footnotes ----------

(1) <https://en.wikipedia.org/wiki/M51_Group>

File: gnuastro.info, Node: Downloading and validating input data, Next: NoiseChisel optimization, Prev: Detecting large extended targets, Up: Detecting large extended targets

2.2.1 Downloading and validating input data
-------------------------------------------

To get the image, you can use the simple field search
(https://dr12.sdss.org/fields) tool of SDSS. As long as it is covered by
the SDSS, you can find an image containing your desired target either by
providing a standard name (if it has one), or its coordinates. To
access the dataset we will use here, write ‘NGC5195’ in the "Object
Name" field and press "Submit" button.

*Type the example commands:* Try to type the example commands on your
terminal and use the history feature of your command-line (by pressing
the "up" button to retrieve previous commands). Do not simply copy and
paste the commands shown here. This will help simulate future
situations when you are processing your own datasets.

You can see the list of available filters under the color image. For
this demonstration, we will use the r-band filter image. By clicking on
the "r-band FITS" link, you can download the image. Alternatively, you
can just run the following command to download it with GNU Wget(1). To
keep things clean, let's also put it in a directory called ‘ngc5195’.
With the ‘-O’ option, we are asking Wget to save the downloaded file
with a more manageable name: ‘r.fits.bz2’ (this is an r-band image of
NGC 5195, which was the directory name).

$ mkdir ngc5195
$ cd ngc5195
$ topurl=https://dr12.sdss.org/sas/dr12/boss/photoObj/frames
$ wget $topurl/301/3716/6/frame-r-003716-6-0117.fits.bz2 -Or.fits.bz2

When you want to reproduce a previous result (a known analysis, on a
known dataset, to get a known result: like the case here!) it is
important to verify that the file is correct: that the input file has
not changed (on the remote server, or in your own archive), or there was
no downloading problem. Otherwise, if the data have changed in your
server/archive, and you use the same script, you will get a different
result, causing a lot of confusion!

One good way to verify the contents of a file is to store its
_Checksum_ in your analysis script and check it before any other
operation. The _Checksum_ algorithms look into the contents of a file
and calculate a fixed-length string from them. If any change (even in a
bit or byte) is made within the file, the resulting string will change,
for more see Wikipedia (https://en.wikipedia.org/wiki/Checksum). There
are many common algorithms, but a simple one is the SHA-1 algorithm
(https://en.wikipedia.org/wiki/SHA-1) (Secure Hash Algorithm 1) that you
can calculate easily with the command below (the second line is the
output, and the checksum is the first/long string: it is independent of
the file name)

$ sha1sum r.fits.bz2
5fb06a572c6107c72cbc5eb8a9329f536c7e7f65 r.fits.bz2

If the checksum on your computer is different from this, either the
file has been incorrectly downloaded (most probable), or it has changed
on SDSS servers (very unlikely(2)). To get a better feeling of
checksums open your favorite text editor and make a test file by writing
something in it. Save it and calculate the text file's SHA-1 checksum
with ‘sha1sum’. Try renaming that file, and you'll see the checksum has
not changed (checksums only look into the contents, not the
name/location of the file). Then open the file with your text editor
again, make a change and re-calculate its checksum, you'll see the
checksum string has changed.

Its always good to keep this short checksum string with your
project's scripts and validate your input data before using them. You
can do this with a shell conditional like this:

filename=r.fits.bz2
expected=5fb06a572c6107c72cbc5eb8a9329f536c7e7f65
sum=$(sha1sum $filename | awk '{print $1}')
if [ $sum = $expected ]; then
echo "$filename: validated"
else
echo "$filename: wrong checksum!"
exit 1
fi

Now that we know you have the same data that we wrote this tutorial
with, let's continue. The SDSS server keeps the files in a Bzip2
compressed file format (that have a ‘.bz2’ suffix). So we will first
decompress it with the following command to use it as a normal FITS
file. By convention, compression programs delete the original file
(compressed when uncompressing, or uncompressed when compressing). To
keep the original file, you can use the ‘--keep’ or ‘-k’ option which is
available in most compression programs for this job. Here, we do not
need the compressed file any more, so we will just let ‘bunzip’ delete
it for us and keep the directory clean.

$ bunzip2 r.fits.bz2

---------- Footnotes ----------

(1) To make the command easier to view on screen or in a page, we
have defined the top URL of the image as the ‘topurl’ shell variable.
You can just replace the value of this variable with ‘$topurl’ in the
‘wget’ command.

(2) If your checksum is different, try uncompressing the file with
the ‘bunzip2’ command after this, and open the resulting FITS file. If
it opens and you see the image of M51 and NGC5195, then there was no
download problem, and the file has indeed changed on the SDSS servers!
In this case, please contact us at ‘bug-gnuastro@gnu.org’.

File: gnuastro.info, Node: NoiseChisel optimization, Next: Skewness caused by signal and its measurement, Prev: Downloading and validating input data, Up: Detecting large extended targets

2.2.2 NoiseChisel optimization
------------------------------

In *note Detecting large extended targets:: we downloaded the single
exposure SDSS image. Let's see how NoiseChisel operates on it with its
default parameters:

$ astnoisechisel r.fits -h0

As described in *note NoiseChisel and Multi-Extension FITS files::,
NoiseChisel's default output is a multi-extension FITS file. Open the
output ‘r_detected.fits’ file and have a look at the extensions, the
0-th extension is only meta-data and contains NoiseChisel's
configuration parameters. The rest are the Sky-subtracted input, the
detection map, Sky values and Sky standard deviation.

$ ds9 -mecube r_detected.fits -zscale -zoom to fit

Flipping through the extensions in a FITS viewer, you will see that
the first image (Sky-subtracted image) looks reasonable: there are no
major artifacts due to bad Sky subtraction compared to the input. The
second extension also seems reasonable with a large detection map that
covers the whole of NGC5195, but also extends towards the bottom of the
image where we actually see faint and diffuse signal in the input image.

Now try flipping between the ‘DETECTIONS’ and ‘SKY’ extensions. In
the ‘SKY’ extension, you'll notice that there is still significant
signal beyond the detected pixels. You can tell that this signal
belongs to the galaxy because the far-right side of the image (away from
M51) is dark (has lower values) and the brighter parts in the Sky image
(with larger values) are just under the detections and follow a similar
pattern.

The fact that signal from the galaxy remains in the ‘SKY’ HDU shows
that NoiseChisel can be optimized for a much better result. The ‘SKY’
extension must not contain any light around the galaxy. Generally, any
time your target is much larger than the tile size and the signal is
very diffuse and extended at low signal-to-noise values (like this
case), this _will_ happen. Therefore, when there are large objects in
the dataset, *the best place* to check the accuracy of your detection is
the estimated Sky image.

When dominated by the background, noise has a symmetric distribution.
However, signal is not symmetric (we do not have negative signal).
Therefore when non-constant(1) signal is present in a noisy dataset, the
distribution will be positively skewed. For a demonstration, see Figure
1 of Akhlaghi and Ichikawa 2015 (https://arxiv.org/abs/1505.01664).
This skewness is a good measure of how much faint signal we have in the
distribution. The skewness can be accurately measured by the difference
in the mean and median (assuming no strong outliers): the more distant
they are, the more skewed the dataset is. This important concept will
be discussed more extensively in the next section (*note Skewness caused
by signal and its measurement::).

However, skewness is only a proxy for signal when the signal has
structure (varies per pixel). Therefore, when it is approximately
constant over a whole tile, or sub-set of the image, the constant
signal's effect is just to shift the symmetric center of the noise
distribution to the positive and there will not be any skewness (major
difference between the mean and median). This positive(2) shift that
preserves the symmetric distribution is the Sky value. When there is a
gradient over the dataset, different tiles will have different constant
shifts/Sky-values, for example, see Figure 11 of Akhlaghi and Ichikawa
2015 (https://arxiv.org/abs/1505.01664).

To make this very large diffuse/flat signal detectable, you will
therefore need a larger tile to contain a larger change in the values
within it (and improve number statistics, for less scatter when
measuring the mean and median). So let's play with the tessellation a
little to see how it affects the result. In Gnuastro, you can see the
option values (‘--tilesize’ in this case) by adding the ‘-P’ option to
your last command. Try running NoiseChisel with ‘-P’ to see its default
tile size.

You can clearly see that the default tile size is indeed much smaller
than this (huge) galaxy and its tidal features. As a result,
NoiseChisel was unable to identify the skewness within the tiles under
the outer parts of M51 and NGC 5159 and the threshold has been
over-estimated on those tiles. To see which tiles were used for
estimating the quantile threshold (no skewness was measured), you can
use NoiseChisel's ‘--checkqthresh’ option:

$ astnoisechisel r.fits -h0 --checkqthresh

Did you see how NoiseChisel aborted after finding and applying the
quantile thresholds? When you call any of NoiseChisel's ‘--check*’
options, by default, it will abort as soon as all the check steps have
been written in the check file (a multi-extension FITS file). This
allows you to focus on the problem you wanted to check as soon as
possible (you can disable this feature with the ‘--continueaftercheck’
option).

To optimize the threshold-related settings for this image, let's play
with this quantile threshold check image a little. Do not forget that
"_Good statistical analysis is not a purely routine matter, and
generally calls for more than one pass through the computer_" (Anscombe
1973, see *note Science and its tools::). A good scientist must have a
good understanding of her tools to make a meaningful analysis. So do
not hesitate in playing with the default configuration and reviewing the
manual when you have a new dataset (from a new instrument) in front of
you. Robust data analysis is an art, therefore a good scientist must
first be a good artist. So let's open the check image as a
multi-extension cube:

$ ds9 -mecube r_qthresh.fits -zscale -cmap sls -zoom to fit

The first extension (called ‘CONVOLVED’) of ‘r_qthresh.fits’ is the
convolved input image where the threshold(s) is(are) defined (and later
applied to). For more on the effect of convolution and thresholding,
see Sections 3.1.1 and 3.1.2 of Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664). The second extension
(‘QTHRESH_ERODE’) has a blank/white value for all the pixels of any tile
that was identified as having significant signal. The other tiles have
the measured threshold over them. The next two extensions
(‘QTHRESH_NOERODE’ and ‘QTHRESH_EXPAND’) are the other two quantile
thresholds that are necessary in NoiseChisel's later steps. Every step
in this file is repeated on the three thresholds.

Play a little with the color bar of the ‘QTHRESH_ERODE’ extension,
you clearly see how the non-blank tiles around NGC 5195 have a gradient.
As one line of attack against discarding too much signal below the
threshold, NoiseChisel rejects outlier tiles. Go forward by three
extensions to ‘VALUE1_NO_OUTLIER’ and you will see that many of the
tiles over the galaxy have been removed in this step. For more on the
outlier rejection algorithm, see the latter half of *note Quantifying
signal in a tile::.

Even though much of the galaxy's footprint has been rejected as
outliers, there are still tiles with signal remaining: play with the DS9
color-bar and you still see a gradient near the outer tidal feature of
the galaxy. Before trying to correct this, let's look at the other
extensions of this check image. We will use a ‘*’ as a wild-card that
can be 1, 2 or 3. In the ‘THRESH*_INTERP’ extensions, you see that all
the blank tiles have been interpolated using their nearest neighbors
(the relevant option here is ‘--interpnumngb’). In the following
‘THRESH*_SMOOTH’ extensions, you can see the tile values after smoothing
(configured with ‘--smoothwidth’ option). Finally, in
‘QTHRESH-APPLIED’, you see the thresholded image: pixels with a value of
1 will be eroded later, but pixels with a value of 2 will pass the
erosion step un-touched.

Let's get back to the problem of optimizing the result. You have two
strategies for detecting the outskirts of the merging galaxies: 1)
Increase the tile size to get more accurate measurements of skewness.
2) Strengthen the outlier rejection parameters to discard more of the
tiles with signal (primarily by increasing ‘--outliernumngb’).
Fortunately in this image we have a sufficiently large region on the
right side of the image that the galaxy does not extend to. So we can
use the more robust first solution. In situations where this does not
happen (for example, if the field of view in this image was shifted to
the left to have more of M51 and less sky) you are limited to a
combination of the two solutions or just to the second solution.

*Skipping convolution for faster tests:* The slowest step of NoiseChisel
is the convolution of the input dataset. Therefore when your dataset is
large (unlike the one in this test), and you are not changing the input
dataset or kernel in multiple runs (as in the tests of this tutorial),
it is faster to do the convolution separately once (using *note
Convolve::) and use NoiseChisel's ‘--convolved’ option to directly feed
the convolved image and avoid convolution. For more on ‘--convolved’,
see *note NoiseChisel input::.

To better identify the skewness caused by the flat NGC 5195 and M51
tidal features on the tiles under it, we have to choose a larger tile
size. Let's try a tile size of 100 by 100 pixels and inspect the check
image.

$ astnoisechisel r.fits -h0 --tilesize=100,100 --checkqthresh
$ ds9 -mecube r_qthresh.fits -zscale -cmap sls -zoom to fit

You can clearly see the effect of this increased tile size: the tiles
are much larger and when you look into ‘VALUE1_NO_OUTLIER’, you see that
all the tiles are nicely grouped on the right side of the image (the
farthest from M51, where we do not see a gradient in ‘QTHRESH_ERODE’).
Things look good now, so let's remove ‘--checkqthresh’ and let
NoiseChisel proceed with its detection.

$ astnoisechisel r.fits -h0 --tilesize=100,100
$ ds9 -mecube r_detected.fits -zscale -cmap sls -zoom to fit

The detected pixels of the ‘DETECTIONS’ extension have expanded a
little, but not as much. Also, the gradient in the ‘SKY’ image is
almost fully removed (and does not fall over M51 anymore). However, on
the bottom-right of the m51 detection, we see many holes gradually
increasing in size. This hints that there is still signal out there.
Let's check the next series of detection steps by adding the
‘--checkdetection’ option this time:

$ astnoisechisel r.fits -h0 --tilesize=100,100 --checkdetection
$ ds9 -mecube r_detcheck.fits -zscale -cmap sls -zoom to fit

The output now has 16 extensions, showing every step that is taken by
NoiseChisel. The first and second (‘INPUT’ and ‘CONVOLVED’) are clear
from their names. The third (‘THRESHOLDED’) is the thresholded image
after finding the quantile threshold (last extension of the output of
‘--checkqthresh’). The fourth HDU (‘ERODED’) is new: it is the
name-stake of NoiseChisel, or eroding pixels that are above the
threshold. By erosion, we mean that all pixels with a value of ‘1’
(above the threshold) that are touching a pixel with a value of ‘0’
(below the threshold) will be flipped to zero (or "carved" out)(3). You
can see its effect directly by going back and forth between the
‘THRESHOLDED’ and ‘ERODED’ extensions.

In the fifth extension (‘OPENED-AND-LABELED’) the image is "opened",
which is a name for eroding once, then dilating (dilation is the inverse
of erosion). This is good to remove thin connections that are only due
to noise. Each separate connected group of pixels is also given its
unique label here. Do you see how just beyond the large M51 detection,
there are many smaller detections that get smaller as you go more
distant? This hints at the solution: the default number of erosions is
too much. Let's see how many erosions take place by default (by adding
‘-P | grep erode’ to the previous command)

$ astnoisechisel r.fits -h0 --tilesize=100,100 -P | grep erode

We see that the value of ‘erode’ is ‘2’. The default NoiseChisel
parameters are primarily targeted to processed images (where there is
correlated noise due to all the processing that has gone into the
warping and stacking of raw images, see *note NoiseChisel optimization
for detection::). In those scenarios 2 erosions are commonly necessary.
But here, we have a single-exposure image where there is no correlated
noise (the pixels are not mixed). So let's see how things change with
only one erosion:

$ astnoisechisel r.fits -h0 --tilesize=100,100 --erode=1 \
--checkdetection
$ ds9 -mecube r_detcheck.fits -zscale -cmap sls -zoom to fit

Looking at the ‘OPENED-AND-LABELED’ extension again, we see that the
main/large detection is now much larger than before. While the
immediately-outer connected regions are still present, they have
decreased dramatically, so we can pass this step.

After the ‘OPENED-AND-LABELED’ extension, NoiseChisel goes onto
finding false detections using the undetected pixels. The process is
fully described in Section 3.1.5. (Defining and Removing False
Detections) of Akhlaghi and Ichikawa 2015
(https://arxiv.org/pdf/1505.01664.pdf). Please compare the extensions
to what you read there and things will be very clear. In the last HDU
(‘DETECTION-FINAL’), we have the final detected pixels that will be used
to estimate the Sky and its Standard deviation. We see that the main
detection has indeed been detected very far out, so let's see how the
full NoiseChisel will estimate the Sky and its standard deviation (by
removing ‘--checkdetection’):

$ astnoisechisel r.fits -h0 --tilesize=100,100 --erode=1
$ ds9 -mecube r_detected.fits -zscale -cmap sls -zoom to fit

The ‘DETECTIONS’ extension of ‘r_detected.fits’ closely follows what
the ‘DETECTION-FINAL’ of the check image (looks good!). If you go ahead
to the ‘SKY’ extension, things still look good. But it can still be
improved.

Look at the ‘DETECTIONS’ again, you will see the right-ward edges of
M51's detected pixels have many "holes" that are fully surrounded by
signal (value of ‘1’) and the signal stretches out in the noise very
thinly (the size of the holes increases as we go out). This suggests
that there is still undetected signal and that we can still dig deeper
into the noise.

With the ‘--detgrowquant’ option, NoiseChisel will "grow" the
detections in to the noise. Its value is the ultimate limit of the
growth in units of quantile (between 0 and 1). Therefore
‘--detgrowquant=1’ means no growth and ‘--detgrowquant=0.5’ means an
ultimate limit of the Sky level (which is usually too much and will
cover the whole image!). See Figure 2 of Akhlaghi 2019
(https://arxiv.org/pdf/1909.11230.pdf) for more on this option. Try
running the previous command with various values (from 0.6 to higher
values) to see this option's effect on this dataset. For this
particularly huge galaxy (with signal that extends very gradually into
the noise), we will set it to ‘0.75’:

$ astnoisechisel r.fits -h0 --tilesize=100,100 --erode=1 \
--detgrowquant=0.75
$ ds9 -mecube r_detected.fits -zscale -cmap sls -zoom to fit

Beyond this level (smaller ‘--detgrowquant’ values), you see many of
the smaller background galaxies (towards the right side of the image)
starting to create thin spider-leg-like features, showing that we are
following correlated noise for too much. Please try it for yourself by
changing it to ‘0.6’ for example.

When you look at the ‘DETECTIONS’ extension of the command shown
above, you see the wings of the galaxy being detected much farther out,
But you also see many holes which are clearly just caused by noise.
After growing the objects, NoiseChisel also allows you to fill such
holes when they are smaller than a certain size through the
‘--detgrowmaxholesize’ option. In this case, a maximum area/size of
10,000 pixels seems to be good:

$ astnoisechisel r.fits -h0 --tilesize=100,100 --erode=1 \
--detgrowquant=0.75 --detgrowmaxholesize=10000
$ ds9 -mecube r_detected.fits -zscale -cmap sls -zoom to fit

When looking at the raw input image (which is very "shallow": less
than a minute exposure!), you do not see anything so far out of the
galaxy. You might just think to yourself that "this is all noise, I
have just dug too deep and I'm following systematics"! If you feel like
this, have a look at the deep images of this system in Watkins 2015
(https://arxiv.org/abs/1501.04599), or a 12 hour deep image of this
system (with a 12-inch telescope):
<https://i.redd.it/jfqgpqg0hfk11.jpg>(4). In these deeper images you
clearly see how the outer edges of the M51 group follow this exact
structure, below in *note Achieved surface brightness level::, we will
measure the exact level.

As the gradient in the ‘SKY’ extension shows, and the deep images
cited above confirm, the galaxy's signal extends even beyond this. But
this is already far deeper than what most (if not all) other tools can
detect. Therefore, we will stop configuring NoiseChisel at this point
in the tutorial and let you play with the other options a little more,
while reading more about it in the papers: Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664) and 2019
(https://arxiv.org/abs/1909.11230)) and *note NoiseChisel::. When you
do find a better configuration feel free to contact us for feedback. Do
not forget that good data analysis is an art, so like a sculptor, master
your chisel for a good result.

To avoid typing all these options every time you run NoiseChisel on
this image, you can use Gnuastro's configuration files, see *note
Configuration files::. For an applied example of setting/using them,
see *note Option management and configuration files::.

*This NoiseChisel configuration is NOT GENERIC:* Do not use the
configuration derived above, on another instrument's image _blindly_.
If you are unsure, just use the default values. As you saw above, the
reason we chose this particular configuration for NoiseChisel to detect
the wings of the M51 group was strongly influenced by the noise
properties of this particular image. Remember *note NoiseChisel
optimization for detection::, where we looked into the very deep XDF
image which had strong correlated noise?

As long as your other images have similar noise properties (from the
same data-reduction step of the same instrument), you can use your
configuration on any of them. But for images from other instruments,
please follow a similar logic to what was presented in these tutorials
to find the optimal configuration.

*Smart NoiseChisel:* As you saw during this section, there is a clear
logic behind the optimal parameter value for each dataset. Therefore,
we plan to add capabilities to (optionally) automate some of the choices
made here based on the actual dataset, please join us in doing this if
you are interested. However, given the many problems in existing
"smart" solutions, such automatic changing of the configuration may
cause more problems than they solve. So even when they are implemented,
we would strongly recommend quality checks for a robust analysis.

---------- Footnotes ----------

(1) by constant, we mean that it has a single value in the region we
are measuring.

(2) In processed images, where the Sky value can be over-estimated,
this constant shift can be negative.

(3) Pixels with a value of ‘2’ are very high signal-to-noise pixels,
they are not eroded, to preserve sharp and bright sources.

(4) The image is taken from this Reddit discussion:
<https://www.reddit.com/r/Astronomy/comments/9d6x0q/12_hours_of_exposure_on_the_whirlpool_galaxy/>

File: gnuastro.info, Node: Skewness caused by signal and its measurement, Next: Image surface brightness limit, Prev: NoiseChisel optimization, Up: Detecting large extended targets

2.2.3 Skewness caused by signal and its measurement
---------------------------------------------------

In the previous section (*note NoiseChisel optimization::) we showed how
to customize NoiseChisel for a single-exposure SDSS image of the M51
group. During the customization, we also discussed the skewness caused
by signal. In the next section (*note Image surface brightness
limit::), we will use this to measure the surface brightness limit of
the image. However, to better understand NoiseChisel and also, the
image surface brightness limit, understanding the skewness caused by
signal, and how to measure it properly are very important. Therefore
now that we have separated signal from noise, let's pause for a moment
and look into skewness, how signal creates it, and find the best way to
measure it.

Let's start masking all the detected pixels found at the end of the
previous section (*note NoiseChisel optimization::) and having a look at
the noise distribution with Gnuastro's Arithmetic and Statistics
programs as shown below (while visually inspecting the masked image with
DS9 in the middle).

$ astarithmetic r_detected.fits -hINPUT-NO-SKY set-in \
r_detected.fits -hDETECTIONS set-det \
in det nan where -odet-masked.fits
$ ds9 det-masked.fits
$ aststatistics det-masked.fits

You will see that Gnuastro's Statistics program prints an ASCII
histogram when no option is given (it is shown below). This is done to
give you a fast and easy view of the distribution of values in the
dataset (pixels in an image, or rows in a table's column).

-------
Input: det-masked.fits (hdu: 1)
-------
Number of elements: 903920
Minimum: -0.113543
Maximum: 0.130339
Median: -0.00216306
Mean: -0.0001893073877
Standard deviation: 0.02569057188
-------
Histogram:
| ** *
| * ** * *
| ** ** * *
| * ** ** ** *
| ** ** ** ** * **
| ** ** ** ** * ** *
| * ** ** ** ** * ** **
| ** ** ** ** **** ** ** *
| ** ** ** ** ** **** ** ** ** *
| ** ** ** ** ** ** ******* ** ** ** *
|*********** ** ** ** ******************* ** ** ** ** ***** ** ***** **
|----------------------------------------------------------------------

This histogram shows a roughly symmetric noise distribution, so let's
have a look at its skewness. The most commonly used definition of
skewness is known as the "Pearson's first skewness coefficient". It
measures the difference between the mean and median, in units of the
standard deviation (STD):

$$\rm{Skewness}\equiv\frac{(\rm{mean}-\rm{median})}{\rm{STD}}$$

The logic behind this definition is simple: as more signal is added
to the same pixels that originally only have raw noise (skewness is
increased), the mean shifts to the positive faster than the median, so
the distance between the mean and median should increase. Let's measure
the skewness (as defined above) over the image without any signal. Its
very easy with Gnuastro's Statistics program (and piping the output to
AWK):

$ aststatistics det-masked.fits --mean --median --std \
| awk '{print ($1-$2)/$3}'
0.0768279

We see that the mean and median are only $0.08\sigma$ (rounded) away
from each other (which is very close)! All pixels with significant
signal are masked, so this is expected, and everything is fine. Now,
let's check the pixel distribution of the sky-subtracted input (where
pixels with significant signal remain, and are not masked):

$ ds9 r_detected.fits
$ aststatistics r_detected.fits -hINPUT-NO-SKY
-------
Input: r_detected.fits (hdu: INPUT-NO-SKY)
Unit: nanomaggy
-------
Number of elements: 3049472
Minimum: -0.113543
Maximum: 159.25
Median: 0.0241158
Mean: 0.1057885317
Standard deviation: 0.698167489
-------
Histogram:
|*
|*
|*
|*
|*
|*
|*
|*
|*
|*
|******************************************* *** ** **** * * * * *
|----------------------------------------------------------------------

Comparing the distributions above, you can see that the _minimum_ value
of the image has not changed because we have not masked the minimum
values. However, as expected, the _maximum_ value of the image has
changed (from $0.13$ to $159.25$). This is clearly evident from the
ASCII histogram: the distribution is very elongated because the galaxy
inside the image is extremely bright.

Now, let's limit the displayed information with the ‘--lessthan=0.13’
option of Statistics as shown below (to only use values less than 0.13;
the maximum of the image where all signal is masked).

$ aststatistics r_detected.fits -hINPUT-NO-SKY --lessthan=0.13
-------
Input: r_detected.fits (hdu: INPUT-NO-SKY)
Range: up to (exclusive) 0.13.
Unit: nanomaggy
-------
Number of elements: 2531949
Minimum: -0.113543
Maximum: 0.126233
Median: 0.0137138
Mean: 0.01735551527
Standard deviation: 0.03590550597
-------
Histogram:
| *
| * ** **
| * * ** ** **
| * * ** ** ** *
| * ** * ** ** ** *
| ** ** * ** ** ** * *
| * ** ** * ** ** ** * *
| ** ** ** * ** ** ** ** * ** *
| * ** ** **** ** ** ** **** ** ** **
| * ** ** ** **** ** ** ** ******* ** ** ** * ** ** **
|***** ** ********** ** ** ********** ** ********** ** ************* **
|----------------------------------------------------------------------

The improvement is obvious: the ASCII histogram better shows the
pixel values near the noise level. We can now compare with the
distribution of ‘det-masked.fits’ that we found earlier. The ASCII
histogram of ‘det-masked.fits’ was approximately symmetric, while this
is asymmetric in this range, especially in outer (to the right, or
positive) direction. The heavier right-side tail is a clear visual
demonstration of skewness that is caused by the signal in the un-masked
image.

Having visually confirmed the skewness, let's quantify it with
Pearson's first skewness coefficient. Like before, we can simply use
Gnuastro's Statistics and AWK for the measurement and calculation:

$ aststatistics r_detected.fits --mean --median --std \
| awk '{print ($1-$2)/$3}'
0.116982

The difference between the mean and median is now approximately
$0.12\sigma$. This is larger than the skewness of the masked image
(which was approximately $0.08\sigma$). At a glance (only looking at
the numbers), it seems that there is not much difference between the two
distributions. However, visually looking at the non-masked image, or
the ASCII histogram, you would expect the quantified skewness to be much
larger than that of the masked image, but that has not happened! Why is
that?

The reason is that the presence of signal does not only shift the
mean and median, it _also_ increases the standard deviation! To see
this for yourself, compare the standard deviation of ‘det-masked.fits’
(which was approximately $0.025$) to ‘r_detected.fits’ (without
‘--lessthan’; which was approximately $0.699$). The latter is almost 28
times larger!

This happens because the standard deviation is defined only in a
symmetric (and Gaussian) distribution. In a non-Gaussian distribution,
the standard deviation is poorly defined and is not a good measure of
"width". Since Pearson's first skewness coefficient is defined in units
of the standard deviation, this very large increase in the standard
deviation has hidden the much increased distance between the mean and
median after adding signal.

We therefore need a better unit or scale to quantify the distance
between the mean and median. A unit that is less affected by skewness
or outliers. One solution that we have found to be very useful is the
quantile units or quantile scale. The quantile scale is defined by
first sorting the dataset (which has $N$ elements). If we want the
quantile of a value $V$ in a distribution, we first find the nearest
data element to $V$ in the sorted dataset. Let's assume the nearest
element is the $i$-th element, counting from 0, after sorting. The
quantile of V in that distribution is then defined as $i/(N-1)$ (which
will have a value between 0 and 1).

The quantile of the median is obvious from its definition: 0.5. This
is because the median is defined to be the middle element of the
distribution after sorting. We can therefore define skewness as the
quantile of the mean ($q_m$). If $q_m\sim0.5$ (the median), then the
distribution (of signal blended in noise) is symmetric (possibly
Gaussian, but the functional form is irrelevant here). A larger value
for $|q_m-0.5|$ quantifies a more skewed the distribution. Furthermore,
a $q_m>0.5$ signifies a positive skewness, while $q_m<0.5$ signifies a
negative skewness.

Let's put this definition to a test on the same two images we have
already created. Fortunately Gnuastro's Statistics program has the
‘--quantofmean’ option to easily calculate $q_m$ for you. So testing is
easy:

$ aststatistics det-masked.fits --quantofmean
0.51295636

$ aststatistics r_detected.fits -hINPUT-NO-SKY --quantofmean
0.8105163158

The two quantiles of mean are now very distinctly different ($0.51$
and $0.81$): differing by about $0.3$ (on a scale of 0 to 1)! Recall
that when defining skewness with Pearson's first skewness coefficient,
their difference was negligible ($0.04\sigma$)! You can now better
appreciate why we discussed quantile so extensively in *note NoiseChisel
optimization::. In case you would like to know more about the usage of
the quantile of the mean in Gnuastro, please see *note Quantifying
signal in a tile::, or watch this video demonstration:
<https://peertube.stream/w/35b7c398-9fd7-4bcf-8911-1e01c5124585>.

File: gnuastro.info, Node: Image surface brightness limit, Next: Achieved surface brightness level, Prev: Skewness caused by signal and its measurement, Up: Detecting large extended targets

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
*note Quantifying measurement limits::). So in this section of the
tutorial, we will measure these values for this image and this target.

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 *note 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 *note 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 *note Quantifying measurement limits::.
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. In fact, MakeCatalog also uses this keyword and will report
the dataset's $n\sigma$ surface brightness limit as keywords in the
output (not as measurement columns, since it is related to the noise,
not labeled signal):

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

Before looking into the measured surface brightness limits, let's review
some important points about this call to MakeCatalog first:
• 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(1).
• If we do not ask for any noise-related column (for example, the
signal-to-noise ratio column with ‘--sn’, among other noise-related
columns), MakeCatalog is not going to read the noise standard
deviation image (again, to speed up its operation when it is
redundant). We are thus using the ‘--forcereadstd’ option (short
for "force read standard deviation image") here so it is ready for
the surface brightness limit measurements that are written as
keywords.

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 "Surface brightness limit (SBL)" title.

$ astfits sbl.fits -h1

Since all the keywords of interest here start with ‘SBL’, we can get a
more cleaner view with this command.

$ astfits sbl.fits -h1 | grep ^SBL

Notice how the ‘SBLSTD’ has the same value as NoiseChisel's ‘MEDSTD’
above. Using ‘SBLSTD’, MakeCatalog has determined the $n\sigma$ surface
brightness limiting magnitude in these header keywords. The multiple of
$\sigma$, or $n$, is the value of the ‘SBLNSIG’ keyword which you can
change with the ‘--sfmagnsigma’. The surface brightness limiting
magnitude within a pixel (‘SBLNSIG’) and within a pixel-agnostic area of
‘SBLAREA’ arcsec$^2$ are stored in ‘SBLMAG’.

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(2). 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 --zeropoint=22.5 \
--output=sbl.fits --forcereadstd --ids
$ astfits sbl.fits -h1 | grep ^SBL

As you see from ‘SBLNSIG’ and ‘SBLAREA’, the default multiple of
sigma is 1 and the default area is 1 arcsec$^2$. Usually higher values
are used for these two parameters. Following the manual example we did
above, you can ask for the multiple of sigma to be 3 and the area to be
25 arcsec$^2$:

$ astmkcatalog r_detected.fits -hDETECTIONS --zeropoint=22.5 \
--output=sbl.fits --sfmagarea=25 --sfmagnsigma=3 \
--forcereadstd --ids
$ astfits sbl.fits -h1 | awk '/^SBLMAG /{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.

*Custom STD for MakeCatalog's Surface brightness limit:* You can
manually change/set the value of the ‘MEDSTD’ keyword in your input STD
image with *note Fits:::

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

With this change, MakeCatalog will use your custom standard deviation
for the surface brightness limit. This is necessary in scenarios where
your image has multiple depths and during your masking, you also mask
the shallow regions (as well as the detections of course).

We have successfully measured the image's $3\sigma$ surface
brightness limiting magnitude over 25 arcsec$^2$. However, as discussed
in *note Quantifying measurement limits:: 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 *note Quantifying measurement
limits::.

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 \
--sfmagarea=25 --sfmagnsigma=3 --forcereadstd \
--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 surface brightness limiting magnitude 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. 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
decrease 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(3). This is a very important point regarding any statistical
process involving random numbers, please see *note Generating random
numbers::.

export GSL_RNG_TYPE=ranlxs1
export GSL_RNG_SEED=1616493518
astmkcatalog lab.fits -h1 --zeropoint=$zeropoint -osbl.fits \
--sfmagarea=25 --sfmagnsigma=3 --forcereadstd \
--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 print all the random locations and their measured brightness as a
table in a file with the suffix ‘_upcheck.fits’. 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 lab_upcheck.fits --noblank=RANDOM_SUM \
| astscript-ds9-region -c1,2 --mode=img \
--radius=$r_pixel

## Have a look at the regions in relation with NoiseChisel's
## detections.
ds9 r_detected.fits[INPUT-NO-SKY] -regions load ds9.reg
ds9 r_detected.fits[DETECTIONS] -regions load ds9.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 stacking 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 (*note General program usage tutorial::).

---------- Footnotes ----------

(1) Recall that NoiseChisel's output is a binary image: 0-valued
pixels are noise and 1-valued pixel are signal. NoiseChisel does not
identify sub-structure over the signal, this is the job of Segment, see
*note Extract clumps and objects::.

(2) From <https://www.sdss.org/dr12/algorithms/magnitudes>

(3) You can use any integer for the seed. One recommendation is to
run MakeCatalog without ‘--envseed’ once and use the randomly generated
seed that is printed on the terminal.

File: gnuastro.info, Node: Achieved surface brightness level, Next: Extract clumps and objects, Prev: Image surface brightness limit, Up: Detecting large extended targets

2.2.5 Achieved surface brightness level
---------------------------------------

In *note NoiseChisel optimization:: we customized NoiseChisel for a
single-exposure SDSS image of the M51 group and in *note Image surface
brightness limit:: we measured the surface brightness limit and the
upper-limit surface brightness level (which are both measures of the
noise level). In this section, let's do some measurements on the
outer-most edges of the M51 group to see how they relate to the noise
measurements found in the previous section.

For this measurement, we will need to estimate the average flux on
the outer edges of the detection. Fortunately all this can be done with
a few simple commands using *note Arithmetic:: and *note MakeCatalog::.
First, let's separate each detected region, or give a unique
label/counter to all the connected pixels of NoiseChisel's detection map
with the command below. Recall that with the ‘set-’ operator, the
popped operand will be given a name (‘det’ in this case) for easy usage
later.

$ astarithmetic r_detected.fits -hDETECTIONS set-det \
det 2 connected-components -olabeled.fits

You can find the label of the main galaxy visually (by opening the
image and hovering your mouse over the M51 group's label). But to have
a little more fun, let's do this automatically (which is necessary in a
general scenario). The M51 group detection is by far the largest
detection in this image, this allows us to find its ID/label easily. We
will first run MakeCatalog to find the area of all the labels, then we
will use Table to find the ID of the largest object and keep it as a
shell variable (‘id’):

# Run MakeCatalog to find the area of each label.
$ astmkcatalog labeled.fits --ids --geo-area -h1 -ocat.fits

## Sort the table by the area column.
$ asttable cat.fits --sort=AREA_FULL

## The largest object, is the last one, so we will use '--tail'.
$ asttable cat.fits --sort=AREA_FULL --tail=1

## We only want the ID, so let's only ask for that column:
$ asttable cat.fits --sort=AREA_FULL --tail=1 --column=OBJ_ID

## Now, let's put this result in a variable (instead of printing)
$ id=$(asttable cat.fits --sort=AREA_FULL --tail=1 --column=OBJ_ID)

## Just to confirm everything is fine.
$ echo $id

We can now use the ‘id’ variable to reject all other detections:

$ astarithmetic labeled.fits $id eq -oonly-m51.fits

Open the image and have a look. To separate the outer edges of the
detections, we will need to "erode" the M51 group detection. So in the
same Arithmetic command as above, we will erode three times (to have
more pixels and thus less scatter), using a maximum connectivity of 2
(8-connected neighbors). We will then save the output in ‘eroded.fits’.

$ astarithmetic labeled.fits $id eq 2 erode 2 erode 2 erode \
-oeroded.fits

In ‘labeled.fits’, we can now set all the 1-valued pixels of
‘eroded.fits’ to 0 using Arithmetic's ‘where’ operator added to the
previous command. We will need the pixels of the M51 group in
‘labeled.fits’ two times: once to do the erosion, another time to find
the outer pixel layer. To do this (and be efficient and more readable)
we will use the ‘set-i’ operator (to give this image the name '‘i’').
In this way we can use it any number of times afterwards, while only
reading it from disk and finding M51's pixels once.

$ astarithmetic labeled.fits $id eq set-i i \
i 2 erode 2 erode 2 erode 0 where -oedge.fits

Open the image and have a look. You'll see that the detected edge of
the M51 group is now clearly visible. You can use ‘edge.fits’ to mark
(set to blank) this boundary on the input image and get a visual feeling
of how far it extends:

$ astarithmetic r.fits -h0 edge.fits nan where -oedge-masked.fits

To quantify how deep we have detected the low-surface brightness
regions (in units of signal to-noise ratio), we will use the command
below. In short it just divides all the non-zero pixels of ‘edge.fits’
in the Sky subtracted input (first extension of NoiseChisel's output) by
the pixel standard deviation of the same pixel. This will give us a
signal-to-noise ratio image. The mean value of this image shows the
level of surface brightness that we have achieved. You can also break
the command below into multiple calls to Arithmetic and create temporary
files to understand it better. However, if you have a look at *note
Reverse polish notation:: and *note Arithmetic operators::, you should
be able to easily understand what your computer does when you run this
command(1).

$ astarithmetic edge.fits -h1 set-edge \
r_detected.fits -hSKY_STD set-skystd \
r_detected.fits -hINPUT-NO-SKY set-skysub \
skysub skystd / edge not nan where meanvalue --quiet

We have thus detected the wings of the M51 group down to roughly
1/3rd of the noise level in this image which is a very good achievement!
But the per-pixel S/N is a relative measurement. Let's also measure the
depth of our detection in absolute surface brightness units; or
magnitudes per square arc-seconds (see *note Brightness flux
magnitude::). We will also ask for the S/N and magnitude of the full
edge we have defined. Fortunately doing this is very easy with
Gnuastro's MakeCatalog:

$ astmkcatalog edge.fits -h1 --valuesfile=r_detected.fits \
--zeropoint=22.5 --ids --sb --sn --magnitude
$ asttable edge_cat.fits
1 25.6971 55.2406 15.8994

We have thus reached an outer surface brightness of $25.70$
magnitudes/arcsec$^2$ (second column in ‘edge_cat.fits’) on this single
exposure SDSS image! This is very similar to the surface brightness
limit measured in *note Image surface brightness limit:: (which is a big
achievement!). But another point in the result above is very
interesting: the total S/N of the edge is $55.24$ with a total edge
magnitude(2) of 15.90!!! This is very large for such a faint signal
(recall that the mean S/N per pixel was 0.32) and shows a very important
point in the study of galaxies: While the per-pixel signal in their
outer edges may be very faint (and invisible to the eye in noise), a lot
of signal hides deeply buried in the noise.

In interpreting this value, you should just have in mind that
NoiseChisel works based on the contiguity of signal in the pixels.
Therefore the larger the object, the deeper NoiseChisel can carve it out
of the noise (for the same outer surface brightness). In other words,
this reported depth, is the depth we have reached for this object in
this dataset, processed with this particular NoiseChisel configuration.
If the M51 group in this image was larger/smaller than this (the field
of view was smaller/larger), or if the image was from a different
instrument, or if we had used a different configuration, we would go
deeper/shallower.

---------- Footnotes ----------

(1) ‘edge.fits’ (extension ‘1’) is a binary (0 or 1 valued) image.
Applying the ‘not’ operator on it, just flips all its pixels (from ‘0’
to ‘1’ and vice-versa). Using the ‘where’ operator, we are then setting
all the newly 1-valued pixels (pixels that are not on the edge) to
NaN/blank in the sky-subtracted input image (‘r_detected.fits’,
extension ‘INPUT-NO-SKY’, which we call ‘skysub’). We are then dividing
all the non-blank pixels (only those on the edge) by the sky standard
deviation (‘r_detected.fits’, extension ‘SKY_STD’, which we called
‘skystd’). This gives the signal-to-noise ratio (S/N) for each of the
pixels on the boundary. Finally, with the ‘meanvalue’ operator, we are
taking the mean value of all the non-blank pixels and reporting that as
a single number.

(2) You can run MakeCatalog on ‘only-m51.fits’ instead of ‘edge.fits’
to see the full magnitude of the M51 group in this image.

File: gnuastro.info, Node: Extract clumps and objects, Prev: Achieved surface brightness level, Up: Detecting large extended targets

2.2.6 Extract clumps and objects (Segmentation)
-----------------------------------------------

In *note NoiseChisel optimization:: we found a good detection map over
the image, so pixels harboring signal have been differentiated from
those that do not. For noise-related measurements like the surface
brightness limit, this is fine. However, after finding the pixels with
signal, you are most likely interested in knowing the sub-structure
within them. For example, how many star forming regions (those bright
dots along the spiral arms) of M51 are within this image? What are the
colors of each of these star forming regions? In the outer most wings
of M51, which pixels belong to background galaxies and foreground stars?
And many more similar questions. To address these questions, you can
use *note Segment:: to identify all the "clumps" and "objects" over the
detection.

$ astsegment r_detected.fits --output=r_segmented.fits
$ ds9 -mecube r_segmented.fits -cmap sls -zoom to fit -scale limits 0 2

Open the output ‘r_segmented.fits’ as a multi-extension data cube
with the second command above and flip through the first and second
extensions, zoom-in to the spiral arms of M51 and see the detected
clumps (all pixels with a value larger than 1 in the second extension).
To optimize the parameters and make sure you have detected what you
wanted, we recommend to visually inspect the detected clumps on the
input image.

For visual inspection, you can make a simple shell script like below.
It will first call MakeCatalog to estimate the positions of the clumps,
then make an SAO DS9 region file and open ds9 with the image and region
file. Recall that in a shell script, the numeric variables (like ‘$1’,
‘$2’, and ‘$3’ in the example below) represent the arguments given to
the script. But when used in the AWK arguments, they refer to column
numbers.

To create the shell script, using your favorite text editor, put the
contents below into a file called ‘check-clumps.sh’. Recall that
everything after a ‘#’ is just comments to help you understand the
command (so read them!). Also note that if you are copying from the PDF
version of this book, fix the single quotes in the AWK command.

#! /bin/bash
set -e # Stop execution when there is an error.
set -u # Stop execution when a variable is not initialized.

# Run MakeCatalog to write the coordinates into a FITS table.
# Default output is `$1_cat.fits'.
astmkcatalog $1.fits --clumpscat --ids --ra --dec

# Use Gnuastro's Table and astscript-ds9-region to build the DS9
# region file (a circle of radius 1 arcseconds on each point).
asttable $1"_cat.fits" -hCLUMPS -cRA,DEC \
| astscript-ds9-region -c1,2 --mode=wcs --radius=1 \
--output=$1.reg

# Show the image (with the requested color scale) and the region file.
ds9 -geometry 1800x3000 -mecube $1.fits -zoom to fit \
-scale limits $2 $3 -regions load all $1.reg

# Clean up (delete intermediate files).
rm $1"_cat.fits" $1.reg

Finally, you just have to activate the script's executable flag with the
command below. This will enable you to directly/easily call the script
as a command.

$ chmod +x check-clumps.sh

This script does not expect the ‘.fits’ suffix of the input's
filename as the first argument. Because the script produces
intermediate files (a catalog and DS9 region file, which are later
deleted). However, we do not want multiple instances of the script (on
different files in the same directory) to collide (read/write to the
same intermediate files). Therefore, we have used suffixes added to the
input's name to identify the intermediate files. Note how all the ‘$1’
instances in the commands (not within the AWK command(1)) are followed
by a suffix. If you want to keep the intermediate files, put a ‘#’ at
the start of the last line.

The few, but high-valued, bright pixels in the central parts of the
galaxies can hinder easy visual inspection of the fainter parts of the
image. With the second and third arguments to this script, you can set
the numerical values of the color map (first is minimum/black, second is
maximum/white). You can call this script with any(2) output of Segment
(when ‘--rawoutput’ is _not_ used) with a command like this:

$ ./check-clumps.sh r_segmented -0.1 2

Go ahead and run this command. You will see the intermediate
processing being done and finally it opens SAO DS9 for you with the
regions superimposed on all the extensions of Segment's output. The
script will only finish (and give you control of the command-line) when
you close DS9. If you need your access to the command-line before
closing DS9, add a ‘&’ after the end of the command above.

While DS9 is open, slide the dynamic range (values for black and
white, or minimum/maximum values in different color schemes) and zoom
into various regions of the M51 group to see if you are satisfied with
the detected clumps. Do not forget that through the "Cube" window that
is opened along with DS9, you can flip through the extensions and see
the actual clumps also. The questions you should be asking yourself are
these: 1) Which real clumps (as you visually _feel_) have been missed?
In other words, is the _completeness_ good? 2) Are there any clumps
which you _feel_ are false? In other words, is the _purity_ good?

Note that completeness and purity are not independent of each other,
they are anti-correlated: the higher your purity, the lower your
completeness and vice-versa. You can see this by playing with the
purity level using the ‘--snquant’ option. Run Segment as shown above
again with ‘-P’ and see its default value. Then increase/decrease it
for higher/lower purity and check the result as before. You will see
that if you want the best purity, you have to sacrifice completeness and
vice versa.

One interesting region to inspect in this image is the many bright
peaks around the central parts of M51. Zoom into that region and
inspect how many of them have actually been detected as true clumps. Do
you have a good balance between completeness and purity? Also look out
far into the wings of the group and inspect the completeness and purity
there.

An easier way to inspect completeness (and only completeness) is to
mask all the pixels detected as clumps and visually inspecting the rest
of the pixels. You can do this using Arithmetic in a command like
below. For easy reading of the command, we will define the shell
variable ‘i’ for the image name and save the output in ‘masked.fits’.

$ in="r_segmented.fits -hINPUT-NO-SKY"
$ clumps="r_segmented.fits -hCLUMPS"
$ astarithmetic $in $clumps 0 gt nan where -oclumps-masked.fits

Inspecting ‘clumps-masked.fits’, you can see some very diffuse peaks
that have been missed, especially as you go farther away from the group
center and into the diffuse wings. This is due to the fact that with
this configuration, we have focused more on the sharper clumps. To put
the focus more on diffuse clumps, you can use a wider convolution
kernel. Using a larger kernel can also help in detecting the existing
clumps to fainter levels (thus better separating them from the
surrounding diffuse signal).

You can make any kernel easily using the ‘--kernel’ option in *note
MakeProfiles::. But note that a larger kernel is also going to wash-out
many of the sharp/small clumps close to the center of M51 and also some
smaller peaks on the wings. Please continue playing with Segment's
configuration to obtain a more complete result (while keeping reasonable
purity). We will finish the discussion on finding true clumps at this
point.

The properties of the clumps within M51, or the background objects
can then easily be measured using *note MakeCatalog::. To measure the
properties of the background objects (detected as clumps over the
diffuse region), you should not mask the diffuse region. When measuring
clump properties with *note MakeCatalog:: and using the ‘--clumpscat’,
the ambient flux (from the diffuse region) is calculated and subtracted.
If the diffuse region is masked, its effect on the clump brightness
cannot be calculated and subtracted.

To keep this tutorial short, we will stop here. See *note
Segmentation and making a catalog:: and *note Segment:: for more on
using Segment, producing catalogs with MakeCatalog and using those
catalogs.

---------- Footnotes ----------

(1) In AWK, ‘$1’ refers to the first column, while in the shell
script, it refers to the first argument.

(2) Some modifications are necessary based on the input dataset:
depending on the dynamic range, you have to adjust the second and third
arguments. But more importantly, depending on the dataset's world
coordinate system, you have to change the region ‘width’, in the AWK
command. Otherwise the circle regions can be too small/large.

File: gnuastro.info, Node: Building the extended PSF, Next: Sufi simulates a detection, Prev: Detecting large extended targets, Up: Tutorials

2.3 Building the extended PSF
=============================

Deriving the extended PSF of an image is very important in many aspects
of the analysis of the objects within it. Gnuastro has a set of
installed scripts, designed to simplify the process following the recipe
of Infante-Sainz et al. 2020 (https://arxiv.org/abs/1911.01430); for
more, see *note PSF construction and subtraction::. An overview of the
process is given in *note Overview of the PSF scripts::.

* Menu:

* Preparing input for extended PSF:: Which stars should be used?
* Saturated pixels and Segment's clumps:: Saturation is a major hurdle!
* One object for the whole detection:: Avoiding over-segmentation in objects.
* Building outer part of PSF:: Building the outermost PSF wings.
* Inner part of the PSF:: Going towards the PSF center.
* Uniting the different PSF components:: Merging all the components into one PSF.
* Subtracting the PSF:: Having the PSF, we now want to subtract it.

File: gnuastro.info, Node: Preparing input for extended PSF, Next: Saturated pixels and Segment's clumps, Prev: Building the extended PSF, Up: Building the extended PSF

2.3.1 Preparing input for extended PSF
--------------------------------------

We will use an image of the M51 galaxy group in the r (SDSS) band of the
Javalambre Photometric Local Universe Survey (J-PLUS) to extract its
extended PSF. For more information on J-PLUS, and its unique features
visit: <http://www.j-plus.es>.

First, let's download the image from the J-PLUS web page using
‘wget’. But to have a generalize-able, and easy to read command, we
will define some base variables (in all-caps) first. After the download
is complete, open the image with SAO DS9 (or any other FITS viewer you
prefer!) to have a feeling of the data (and of course, enjoy the beauty
of M51 in such a wide field of view):

$ urlend="jplus-dr2/get_fits?id=67510"
$ urlbase="http://archive.cefca.es/catalogues/vo/siap/"
$ mkdir jplus-dr2
$ wget $urlbase$urlend -O jplus-dr2/67510.fits.fz
$ astscript-fits-view jplus-dr2/67510.fits.fz

After enjoying the large field of view, have a closer look at the
edges of the image. Please zoom in to the corners. You will see that
on the edges, the pixel values are either zero or with significantly
different values than the main body of the image. This is due to the
dithering pattern that was used to make this image and happens in all
imaging surveys(1). To avoid potential issues or problems that these
regions may cause, we will first crop out the main body of the image
with the command below. To keep the top-level directory clean, let's
also put the crop in a directory called ‘flat’.

$ mkdir flat
$ astcrop jplus-dr2/67510.fits.fz --section=225:9275,150:9350 \
--mode=img -oflat/67510.fits
$ astscript-fits-view flat/67510.fits

Please zoom into the edges again, you will see that they now have the
same noise-level as the rest of the image (the problematic parts are now
gone).

---------- Footnotes ----------

(1) Recall the cropping in a previous tutorial for a similar reason
(varying "depth" across the image): *note Dataset inspection and
cropping::.

File: gnuastro.info, Node: Saturated pixels and Segment's clumps, Next: One object for the whole detection, Prev: Preparing input for extended PSF, Up: Building the extended PSF

2.3.2 Saturated pixels and Segment's clumps
-------------------------------------------

A constant-depth (flat) image was created in the previous section (*note
Preparing input for extended PSF::). As explained in *note Overview of
the PSF scripts::, an important step when building the PSF is to mask
other sources in the image. Therefore, before going onto selecting
stars, let's detect all significant signal, and identify the clumps of
background objects over the wings of the extended PSF.

There is a problem however: the saturated pixels of the bright stars
are going to cause problems in the segmentation phase. To see this
problem, let's make a $1000\times1000$ crop around a bright star to
speed up the test (and its solution). Afterwards we will apply the
solution to the whole image.

$ astcrop flat/67510.fits --mode=wcs --widthinpix --width=1000 \
--center=203.3916736,46.7968652 --output=saturated.fits
$ astnoisechisel saturated.fits --output=sat-nc.fits
$ astsegment sat-nc.fits --output=sat-seg.fits
$ astscript-fits-view sat-seg.fits

Have a look at the ‘CLUMPS’ extension. You will see that instead of
a single clump at the center of the bright star, we have many clumps!
This has happened because of the saturated pixels! When saturation
occurs, the sharp peak of the profile is lost (like cutting off the tip
of a mountain to build a telescope!) and all saturated pixels get a
noisy value close to the saturation level. To see this saturation noise
run the last command again and in SAO DS9, set the "Scale" to "min max"
and zoom into the center. You will see the noisy saturation pixels at
the center of the star in red.

This noise-at-the-peak disrupts Segment's assumption to expand clumps
from a local maxima: each noisy peak is being treated as a separate
local maxima and thus a separate clump. For more on how Segment defines
clumps, see Section 3.2.1 and Figure 8 of Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664). To have the center identified as a
single clump, we should mask these saturated pixels in a way that suites
Segment's non-parametric methodology.

First we need to find the saturation level! The saturation level is
usually fixed for any survey or input data that you receive from a
certain database, so you will usually have to do this only once (the
first time you get data from that database). Let's make a smaller crop
of $50\times50$ pixels around the star with the first command below.
With the next command, please look at the crop with DS9 to visually
understand the problem. You will see the saturated pixels as the noisy
red pixels in the center of the image. A non-saturated star will have a
single pixel as the maximum and will not have such a large area covered
by a noisy constant value (find a few stars in the image and see for
yourself). Visual and qualitative inspection of the process is very
important for understanding the solution.

$ astcrop saturated.fits --mode=wcs --widthinpix --width=50 \
--center=203.3916736,46.7968652 --output=sat-center.fits
$ astscript-fits-view sat-center.fits --ds9scale=minmax

To quantitatively identify the saturation level in this image, let's
have a look at the distribution of pixels with a value larger than 100
(above the noise level):

$ aststatistics sat-center.fits --greaterequal=100
Histogram:
|*
|*
|*
|*
|* *
|** *
|*** **
|**** **
|****** ****
|********** * * * ******
|************************* ************ * *** ******* *** ************
|----------------------------------------------------------------------

The peak you see in the right end (larger values) of the histogram
shows the saturated pixels (a constant level, with some scatter due to
the large Poisson noise). If there was no saturation, the number of
pixels should have decreased at increasing values; until reaching the
maximum value of the profile in one pixel. But that is not the case
here. Please try this experiment on a non-saturated (fainter) star to
see what we mean.

If you still have not experimented on a non-saturated star, please
stop reading this tutorial! Please open ‘flat/67510.fits’ in DS9,
select a fainter/smaller star and repeat the last three commands (with a
different center). After you have confirmed the point above (visually,
and with the histogram), please continue with the rest of this tutorial.

Finding the saturation level is easy with Statistics (by using the
‘--lessthan’ option until the histogram becomes as expected: only
decreasing). First, let's try ‘--lessthan=3000’:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=3000
-------
Histogram:
|*
|*
|*
|*
|*
|**
|*** *
|**** *
|******* **
|*********** * * * * * * * ****
|************************* * ***** ******* ***** ** ***** * ********
|----------------------------------------------------------------------

We still see an increase in the histogram around 3000. Let's try a
threshold of 2500:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=2500
-------
Histogram:
|*
|*
|**
|**
|**
|**
|****
|*****
|*********
|************* * * * *
|********************************* ** ** ** *** ** * **** ** *****
|----------------------------------------------------------------------

The peak at the large end of the histogram has gone! But let's have
a closer look at the values (the resolution of an ASCII histogram is
limited!). To do this, we will ask Statistics to save the histogram
into a table with the ‘--histogram’ option, then look at the last 20
rows:

$ aststatistics sat-center.fits --greaterequal=100 --lessthan=2500 \
--histogram --output=sat-center-hist.fits
$ asttable sat-center-hist.fits --tail=20
2021.1849112701 1
2045.0495397186 0
2068.9141681671 1
2092.7787966156 1
2116.6434250641 0
2140.5080535126 0
2164.3726819611 0
2188.2373104095 0
2212.101938858 1
2235.9665673065 1
2259.831195755 2
2283.6958242035 0
2307.560452652 0
2331.4250811005 1
2355.289709549 1
2379.1543379974 1
2403.0189664459 2
2426.8835948944 1
2450.7482233429 2
2474.6128517914 2

Since the number of points at the extreme end are increasing (from 1
to 2), We therefore see that a value 2500 is still above the saturation
level (the number of pixels has started to increase)! A more reasonable
saturation level for this image would be 2200! As an exercise, you can
try automating this selection with AWK.

Therefore, we can set the saturation level in this image(1) to be
2200. Let's mask all such pixels with the command below:

$ astarithmetic saturated.fits set-i i i 2200 gt nan where \
--output=sat-masked.fits
$ astscript-fits-view sat-masked.fits --ds9scale=minmax

Please see the peaks of several bright stars, not just the central
very bright star. Zoom into each of the peaks you see. Besides the
central very bright one that we were looking at closely until now, only
one other star is saturated (its center is NaN, or Not-a-Number). Try
to find it.

But we are not done yet! Please zoom-in to that central bright star
and have another look on the edges of the vertical "bleeding" saturated
pixels, there are strong positive/negative values touching it (almost
like "waves"). These will also cause problems and have to be masked!
So with a small addition to the previous command, let's dilate the
saturated regions (with 2-connectivity, or 8-connected neighbors) four
times and have another look:

$ astarithmetic saturated.fits set-i i i 2200 gt \
2 dilate 2 dilate 2 dilate 2 dilate \
nan where --output=sat-masked.fits
$ astscript-fits-view sat-masked.fits --ds9scale=minmax

Now that saturated pixels (and their problematic neighbors) have been
masked, we can convolve the image (recall that Segment will use the
convolved image for identifying clumps) with the command below.
However, we will use the Spatial Domain convolution which can account
for blank pixels (for more on the pros and cons of spatial and frequency
domain convolution, see *note Spatial vs. Frequency domain::). We will
also create a Gaussian kernel with $\rm{FWHM}=2$ pixels, truncated at
$5\times\rm{FWHM}$.

$ astmkprof --kernel=gaussian,2,5 --oversample=1 -okernel.fits
$ astconvolve sat-masked.fits --kernel=kernel.fits --domain=spatial \
--output=sat-masked-conv.fits
$ astscript-fits-view sat-masked-conv.fits --ds9scale=minmax

Please zoom-in to the star and look closely to see how after
spatial-domain convolution, the problematic pixels are still NaN. But
Segment requires the profile to start with a maximum value and decrease.
So before feeding into Segment, let's fill the blank values with the
maximum value of the neighboring pixels in both the input and convolved
images (see *note Interpolation operators::):

$ astarithmetic sat-masked.fits 2 interpolate-maxofregion \
--output=sat-fill.fits
$ astarithmetic sat-masked-conv.fits 2 interpolate-maxofregion \
--output=sat-fill-conv.fits
$ astscript-fits-view sat-fill* --ds9scale=minmax

Have a closer look at the opened images. Please zoom-in (you will
notice that they are already matched and locked, so they will both
zoom-in together). Go to the centers of the saturated stars and confirm
how they are filled with the largest non-blank pixel. We can now feed
this image to NoiseChisel and Segment as the convolved image:

$ astnoisechisel sat-fill.fits --convolved=sat-fill-conv.fits \
--output=sat-nc.fits
$ astsegment sat-nc.fits --convolved=sat-fill-conv.fits \
--output=sat-seg.fits --rawoutput
$ ds9 -mecube sat-seg.fits -zoom to fit -scale limits -1 1

See the ‘CLUMPS’ extension. Do you see how the whole center of the star
has indeed been identified as a single clump? We thus achieved our aim
and did not let the saturated pixels harm the identification of the
center!

If the issue was only clumps (like in a normal deep image
processing), this was the end of Segment's special considerations.
However, in the scenario here, with the very extended wings of the
bright stars, it usually happens that background objects become "clumps"
in the outskirts and will rip the bright star outskirts into separate
"objects". In the next section (*note One object for the whole
detection::), we will describe how you can modify Segment to avoid this
issue.

---------- Footnotes ----------

(1) In raw exposures, this value is usually around 65000 (close to
$2^{16}$, since most CCDs have 16-bit pixels; see *note Numeric data
types::). But that is not the case here, because this is a
processed/stacked image that has been calibrated.

File: gnuastro.info, Node: One object for the whole detection, Next: Building outer part of PSF, Prev: Saturated pixels and Segment's clumps, Up: Building the extended PSF

2.3.3 One object for the whole detection
----------------------------------------

In *note Saturated pixels and Segment's clumps::, we described how you
can run Segment such that saturated pixels do not interfere with its
clumps. However, due to the very extended wings of the PSF, the default
definition of "objects" should also be modified for the scenario here.
To better see the problem, let's inspect now the ‘OBJECTS’ extension,
focusing on those objects with a label between 50 to 150 (which include
the main star):

$ astscript-fits-view sat-seg.fits -hOBJECTS --ds9scale="limits 50 150"

We can see that the detection corresponding to the star has been
broken into different objects. This is not a good object segmentation
image for our scenario here. Since those objects in the outer wings of
the bright star's detection harbor a lot of the extended PSF. We want to
keep them with the same "object" label as the star (we only need to mask
the "clumps" of the background sources). To do this, we will make the
following changes to Segment's options (see *note Segmentation options::
for more on this options):
• Since we want the extended diffuse flux of the PSF to be taken as a
single object, we want all the grown clumps to touch. Therefore,
it is necessary to decrease ‘--gthresh’ to very low values, like
$-10$. Recall that its value is in units of the input standard
deviation, so ‘--gthresh=-10’ corresponds to $-10\sigma$. The
default value is not for such extended sources that dominate all
background sources.
• Since we want all connected grown clumps to be counted as a single
object in any case, we will set ‘--objbordersn=0’ (its smallest
possible value).

Let's make these changes and check if the star has been kept as a single
object in the ‘OBJECTS’ extension or not:

$ astsegment sat-nc.fits --convolved=sat-fill-conv.fits \
--gthresh=-10 --objbordersn=0 \
--output=sat-seg.fits --rawoutput
$ astscript-fits-view sat-seg.fits -hOBJECTS --ds9scale="limits 50 150"

Now we can extend these same steps to the whole image. To detect
signal, we can run NoiseChisel using the command below. We modified the
default value to two of the options, below you can see the reason for
these changes. See *note Detecting large extended targets:: for more on
optimizing NoiseChisel.
• Since the image is so large, we have increased ‘--outliernumngb’ to
get better outlier statistics on the tiles. The default value is
primarily for small images, so this is usually the first thing you
should do when running NoiseChisel on a real/large image.
• Since the image is not too deep (made from few exposures), it does
not have strong correlated noise, so we will decrease
‘--detgrowquant’ and increase ‘--detgrowmaxholesize’ to better
extract signal.
Furthermore, since both NoiseChisel and Segment need a convolved image,
we will do the convolution before and feed it to both (to save running
time). But in the first command below, let's delete all the temporary
files we made above.

Since the image is large (+300 MB), to avoid wasting storage, any
temporary file that is no longer necessary for later processing is
deleted after it is used. You can visually check each of them with DS9
before deleting them (or not delete them at all!). Generally, within a
pipeline it is best to remove such large temporary files, because space
runs out much faster than you think (for example, once you get good
results and want to use more fields).

$ rm *.fits
$ mkdir label
$ astmkprof --kernel=gaussian,2,5 --oversample=1 \
-olabel/kernel.fits
$ astarithmetic flat/67510.fits set-i i i 2200 gt \
2 dilate 2 dilate 2 dilate 2 dilate nan where \
--output=label/67510-masked-sat.fits
$ astconvolve label/67510-masked-sat.fits --kernel=label/kernel.fits \
--domain=spatial --output=label/67510-masked-conv.fits
$ rm label/kernel.fits
$ astarithmetic label/67510-masked-sat.fits 2 interpolate-maxofregion \
--output=label/67510-fill.fits
$ astarithmetic label/67510-masked-conv.fits 2 interpolate-maxofregion \
--output=label/67510-fill-conv.fits
$ rm label/67510-masked-conv.fits
$ astnoisechisel label/67510-fill.fits --outliernumngb=100 \
--detgrowquant=0.8 --detgrowmaxholesize=100000 \
--convolved=label/67510-fill-conv.fits \
--output=label/67510-nc.fits
$ rm label/67510-fill.fits
$ astsegment label/67510-nc.fits --output=label/67510-seg-raw.fits \
--convolved=label/67510-fill-conv.fits --rawoutput \
--gthresh=-10 --objbordersn=0
$ rm label/67510-fill-conv.fits
$ astscript-fits-view label/67510-seg-raw.fits

We see that the saturated pixels have not caused any problem and the
central clumps/objects of bright stars are now a single clump/object.
We can now proceed to estimating the outer PSF. But before that, let's
make a "standard" segment output: one that can safely be fed into
MakeCatalog for measurements and can contain all necessary outputs of
this whole process in a single file (as multiple extensions).

The main problem is again the saturated pixels: we interpolated them
to be the maximum of their nearby pixels. But this will cause problems
in any measurement that is done over those regions. To let MakeCatalog
know that those pixels should not be used, the first extension of the
file given to MakeCatalog should have blank values on those pixels. We
will do this with the commands below:

## First HDU of Segment (Sky-subtracted input)
$ astarithmetic label/67510-nc.fits -hINPUT-NO-SKY \
label/67510-masked-sat.fits isblank nan where \
--output=label/67510-seg.fits
$ astfits label/67510-seg.fits --update=EXTNAME,INPUT-NO-SKY

## Second and third HDUs: CLUMPS and OBJECTS
$ astfits label/67510-seg-raw.fits --copy=CLUMPS --copy=OBJECTS \
--output=label/67510-seg.fits

## Fourth HDU: Sky standard deviation (from NoiseChisel):
$ astfits label/67510-nc.fits --copy=SKY_STD \
--output=label/67510-seg.fits

## Clean up all the un-necessary files:
$ rm label/67510-masked-sat.fits label/67510-nc.fits \
label/67510-seg-raw.fits

You can now simply run MakeCatalog on this image and be sure that
saturated pixels will not affect the measurements. As one example, you
can use MakeCatalog to find the clumps containing saturated pixels:
recall that the ‘--area’ column only calculates the area of non-blank
pixels, while ‘--geo-area’ calculates the area of the label (independent
of their blank-ness in the values image):

$ astmkcatalog label/67510-seg.fits --ids --ra --dec --area \
--geo-area --clumpscat --output=cat.fits

The information of the clumps that have been affected by saturation
can easily be found by selecting those with a differing value in the
‘AREA’ and ‘AREA_FULL’ columns:

## With AWK (second command, counts the number of rows)
$ asttable cat.fits -hCLUMPS | awk '$5!=$6'
$ asttable cat.fits -hCLUMPS | awk '$5!=$6' | wc -l

## Using Table arithmetic (compared to AWK, you can use column
## names, save as FITS, and be faster):
$ asttable cat.fits -hCLUMPS -cRA,DEC --noblankend=3 \
-c'arith AREA AREA AREA_FULL eq nan where'

## Remove the table (which was just for a demo)
$ rm cat.fits

We are now ready to start building the outer parts of the PSF in *note
Building outer part of PSF::.

File: gnuastro.info, Node: Building outer part of PSF, Next: Inner part of the PSF, Prev: One object for the whole detection, Up: Building the extended PSF

2.3.4 Building outer part of PSF
--------------------------------

In *note Saturated pixels and Segment's clumps::, and *note 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 *note
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 *note 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 *note 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
*note 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 *note 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 *note 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 *note 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.

File: gnuastro.info, Node: Inner part of the PSF, Next: Uniting the different PSF components, Prev: Building outer part of PSF, Up: Building the extended PSF

2.3.5 Inner part of the PSF
---------------------------

In *note Building outer part of PSF::, we were able to create a stack of
the outer-most behavior of the PSF in a J-PLUS survey image. But the
central part that was affected by saturation and non-linearity is still
remaining, and we still do not have a "complete" PSF! In this section,
we will use the same steps before to make stacks of more inner regions
of the PSF to ultimately unite them all into a single PSF in *note
Uniting the different PSF components::.

For the outer PSF, we selected stars in the magnitude range of 6 to
10. So let's have a look and see how many stars we have in the
magnitude range of 12 to 13 with a more relaxed condition on the minimum
distance for neighbors, ‘--mindistdeg=0.01’ (36 arcsec, since these
stars are fainter), and use the ds9 region script to visually inspect
them:

$ mkdir inner
$ astscript-psf-select-stars flat/67510.fits \
--magnituderange=12,13 --mindistdeg=0.01 \
--output=inner/67510-12-13.fits

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

We have 41 stars, but if you zoom into their centers, you will see
that they do not have any major bleeding-vertical saturation any more.
Only the very central core of some of the stars is saturated. We can
therefore use these stars to fill the strong bleeding footprints that
were present in the outer stack of ‘outer/stack.fits’. Similar to
before, let's build ready-to-stack crops of these stars. To get a
better feeling of the normalization radii, follow the same steps of
*note Building outer part of PSF:: (setting ‘--tmpdir’ and ‘--keeptmp’).
In this case, since the stars are fainter, we can set a smaller size for
the individual stamps, ‘--widthinpix=500,500’, to speed up the
calculations:

$ counter=1
$ mkdir inner/stamps
$ asttable inner/67510-12-13.fits \
| while read -r ra dec mag; do
astscript-psf-stamp label/67510-seg.fits \
--mode=wcs \
--normradii=5,10 \
--center=$ra,$dec \
--widthinpix=500,500 \
--segment=label/67510-seg.fits \
--output=inner/stamps/67510-$counter.fits; \
counter=$((counter+1)); \
done

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

We are now ready to unite the two stacks we have constructed: the outer
and the inner parts.

File: gnuastro.info, Node: Uniting the different PSF components, Next: Subtracting the PSF, Prev: Inner part of the PSF, Up: Building the extended PSF

2.3.6 Uniting the different PSF components
------------------------------------------

In *note Building outer part of PSF:: we built the outer part of the
extended PSF and the inner part was built in *note Inner part of the
PSF::. The outer part was built with very bright stars, and the inner
part using fainter stars to not have saturation in the core of the PSF.
The next step is to join these two parts in order to have a single PSF.
First of all, let's have a look at the two stacks and also to their
radial profiles to have a good feeling of the task. Note that you will
need to have TOPCAT to run the last command and plot the radial profile
(see *note TOPCAT::).

$ astscript-fits-view outer/stack.fits inner/stack.fits
$ astscript-radial-profile outer/stack.fits -o outer/profile.fits
$ astscript-radial-profile inner/stack.fits -o inner/profile.fits
$ astscript-fits-view outer/profile.fits inner/profile.fits

From the visual inspection of the images and the radial profiles, it
is clear that we have saturation in the center for the outer part. Note
that the absolute flux values of the PSFs are meaningless since they
depend on the normalization radii we used to obtain them. The uniting
step consists in scaling up (or down) the inner part of the PSF to have
the same flux at the junction radius, and then, use that flux-scaled
inner part to fill the center of the outer PSF. To get a feeling of the
process, first, let's open the two radial profiles and find the factor
manually first:

1. Run this command to open the two tables in *note TOPCAT:::
$ astscript-fits-view outer/profile.fits inner/profile.fits
2. On the left side of the screen, under "Table List", you will see
the two imported tables. Click on the first one (profile of the
outer part) so it is shown first.
3. Under the "Graphics" menu item, click on "Plane plot". A new
window will open with the plot of the first two columns: ‘RADIUS’
on the horizontal axis and ‘MEAN’ on the vertical. The rest of the
steps are done in this window.
4. In the bottom settings, within the left panel, click on the "Axes"
item. This will allow customization of the plot axes.
5. In the bottom-right panel, click on the box in front of "Y Log" to
make the vertical axis logarithmic-scaled.
6. On the "Layers" menu, select "Add Position Control" to allow adding
the profile of the inner region. After it, you will see that a new
red-blue scatter plot icon opened on the bottom-left menu (with a
title of ‘<no table>’).
7. On the bottom-right panel, in the drop-down menu in front of
‘Table:’, select ‘2: profile.fits’. Afterwards, you will see the
radial profile of the inner stack as the newly added blue plot.
Our goal here is to find the factor that is necessary to multiply
with the inner profile so it matches the outer one.
8. On the bottom-right panel, in front of ‘Y:’, you will see ‘MEAN’.
Click in the white-space after it, and type this: ‘*100’. This
will display the ‘MEAN’ column of the inner profile, after
multiplying it by 100. Afterwards, you will see that the inner
profile (blue) matches more cleanly with the outer (red);
especially in the smaller radii. At larger radii, it does not drop
like the red plot. This is because of the extremely low
signal-to-noise ratio at those regions in the fainter stars used to
make this stack.
9. Take your mouse cursor over the profile, in particular over the
bump around a radius of 100 pixels. Scroll your mouse down-ward to
zoom-in to the profile and up-ward to zoom-out. You can
click-and-hold any part of the profile and if you move your cursor
(while still holding the mouse-button) to look at different parts
of the profile. This is particular helpful when you have zoomed-in
to the profile.
10. Zoom-in to the bump around a radius of 100 pixels until the
horizontal axis range becomes around 50 to 130 pixels.
11. You clearly see that the inner stack (blue) is much more noisy
than the outer (red) stack. By "noisy", we mean that the scatter
of the points is much larger. If you further zoom-out, you will
see that the shallow slope at the larger radii of the inner (blue)
profile has also affected the height of this bump in the inner
profile. This is a _very important_ point: this clearly shows that
the inner profile is too noisy at these radii.
12. Click-and-hold your mouse to see the inner parts of the two
profiles (in the range 0 to 80). You will see that for radii less
than 40 pixels, the inner profile (blue) points loose their scatter
(and thus have a good signal-to-noise ratio).
13. Zoom-in to the plot and follow the profiles until smaller radii
(for example, 10 pixels). You see that for each radius, the inner
(blue) points are consistently above the outer (red) points. This
shows that the $\times100$ factor we selected above was too much.
14. In the bottom-right panel, change the ‘100’ to ‘80’ and zoom-in to
the same region. At each radius, the blue points are now below the
red points, so the scale-factor 80 is not enough. So let's
increase it and try ‘90’. After zooming-in, you will notice that
in the inner-radii (less than 30 pixels), they are now very
similar. The ultimate aim of the steps below is to find this
factor automatically.
15. But before continuing, let's focus on another important point
about the central regions: non-linearity and saturation. While you
are zoomed-in (from the step above), follow (click-and-drag) the
profile towards smaller radii. You will see that smaller than a
radius of 10, they start to diverge. But this time, the outer
(red) profile is getting a shallower slope and diverges
significantly from about the radius of 8. We had masked all
saturated pixels before, so this divergence for radii smaller than
10 shows the effect of the CCD's non-linearity (where the number of
electrons will not be linearly correlated with the number of
incident photons). This is present in all CCDs and pixels beyond
this level should not be used in measurements (or properly
corrected).

The items above were only listed so you get a good mental/visual
understanding of the logic behind the operation of the next script (and
to learn how to tune its parameters where necessary):
‘astscript-psf-scale-factor’. This script is more general than this
particular problem, but can be used for this special case also. Its job
is to take a model of an object (PSF, or inner stack in this case) and
the position of an instance of that model (a star, or the outer stack in
this case) in a larger image.

Instead of dealing with radial profiles (that enforce a certain
shape), this script will put the centers of the inner and outer PSFs
over each other and divide the outer by the inner. Let's have a look
with the command below. Just note that we are running it with
‘--keeptmp’ so the temporary directory with all the intermediate files
remain for further clarification:

$ astscript-psf-scale-factor outer/stack.fits \
--psf=inner/stack.fits --center=501,501 \
--mode=img --normradii=10,15 --keeptmp
$ astscript-fits-view stack_psfmodelscalefactor/cropped-*.fits \
stack_psfmodelscalefactor/for-factor-*.fits

With the second command, you see the four steps of the process: the
first two images show the cropped outer and inner stacks (to same width
image). The third shows the radial position of each pixel (which is
used to only keep the pixels within the desired radial range). The
fourth shows the per-pixel division of the outer by the inner within the
requested radii. The sigma-clipped median of these pixels is finally
reported. Unlike the radial profile method (which averages over a
circular/elliptical annulus for each radius), this method imposes no
a-priori shape on the PSF. This makes it very useful for complex PSFs
(like the case here).

To continue, let's remove the temporary directory and re-run the
script but with ‘--quiet’ mode so we can put the output in a shell
variable.

$ rm -r stack_psfmodelscalefactor
$ scale=$(astscript-psf-scale-factor outer/stack.fits \
--psf=inner/stack.fits --center=501,501 \
--mode=img --normradii=10,15 --quiet)
$ echo $scale

Now that we know the scaling factor, we are ready to unite the outer
and the inner part of the PSF. To do that, we will use the script
‘astscript-psf-unite’ with the command below (for more on this script,
see *note Invoking astscript-psf-unite::). The basic parameters are the
inner part of the PSF (given to ‘--inner’), the inner part's scale
factor (‘--scale’), and the junction radius (‘--radius’). The inner
part is first scaled, and all the pixels of the outer image within the
given radius are replaced with the pixels of the inner image. Since the
flux factor was computed for a ring of pixels between 10 and 15 pixels,
let's set the junction radius to be 12 pixels (roughly in between 10 and
15):

$ astscript-psf-unite outer/stack.fits \
--inner=inner/stack.fits --radius=12 \
--scale=$scale --output=psf.fits

Let's have a look at the outer stack and the final PSF with the command
below. Since we want several other DS9 settings to help you directly
see the main point, we are using ‘--ds9extra’. After DS9 is opened, you
can see that the center of the PSF has now been nicely filled. You can
click on the "Edit" button and then the "Colorbar" and hold your cursor
over the image and move it. You can see that besides filling the inner
regions nicely, there is also no major discontinuity in the 2D image
around the union radius of 12 pixels around the center.

$ astscript-fits-view outer/stack.fits psf.fits --ds9scale=minmax \
--ds9extra="-scale limits 0 22000 -match scale" \
--ds9extra="-lock scale yes -zoom 4 -scale log"

Nothing demonstrates the effect of a bad analysis than actually
seeing a bad result! So let's choose a bad normalization radial range
(50 to 60 pixels) and unite the inner and outer parts based on that.
The last command will open the two PSFs together in DS9, you should be
able to immediately see the discontinuity in the union radius.

$ scale=$(astscript-psf-scale-factor outer/stack.fits \
--psf=inner/stack.fits --center=501,501 \
--mode=img --normradii=50,60 --quiet)

$ astscript-psf-unite outer/stack.fits \
--inner=inner/stack.fits --radius=55 \
--scale=$scale --output=psf-bad.fits

$ astscript-fits-view psf-bad.fits psf.fits --ds9scale=minmax \
--ds9extra="-scale limits 0 50 -match scale" \
--ds9extra="-lock scale yes -zoom 4 -scale log"

As you see, the selection of the normalization radii and the unite
radius are very important. The first time you are trying to build the
PSF of a new dataset, it has to be explored with a visual inspection of
the images and radial profiles. Once you have found a good
normalization radius for a certain part of the PSF in a survey, you can
generally use it comfortably without change. But for a new survey, or a
different part of the PSF, be sure to repeat the visual checks above to
choose the best radii. As a summary, a good junction radius is one
that:

• Is large enough to not let saturation and non-linearity (from the
outer profile) into the inner region.
• Is small enough to have a sufficiently high signal to noise ratio
(from the inner profile) to avoid adding noise in the union radius.

Now that the complete PSF has been obtained, let's remove that
bad-looking PSF, and stick with the nice and clean PSF for the next step
in *note Subtracting the PSF::.

$ rm -rf psf-bad.fits

File: gnuastro.info, Node: Subtracting the PSF, Prev: Uniting the different PSF components, Up: Building the extended PSF

2.3.7 Subtracting the PSF
-------------------------

Previously (in *note 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
*note 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(1)) 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 *note
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.

---------- Footnotes ----------

(1) For example, in J-PLUS, the baffle of the secondary mirror was
adjusted in 2017 because it produced extra spikes in the PSF. So all
images after that date have a PSF with 4 spikes (like this one), while
those before it have many more spikes.

File: gnuastro.info, Node: Sufi simulates a detection, Next: Detecting lines and extracting spectra in 3D data, Prev: Building the extended PSF, Up: Tutorials

2.4 Sufi simulates a detection
==============================

It is the year 953 A.D. and Abd al-rahman Sufi (903 - 986 A.D.)(1) is in
Shiraz as a guest astronomer. He had come there to use the advanced 123
centimeter astrolabe for his studies on the ecliptic. However,
something was bothering him for a long time. While mapping the
constellations, there were several non-stellar objects that he had
detected in the sky, one of them was in the Andromeda constellation.
During a trip he had to Yemen, Sufi had seen another such object in the
southern skies looking over the Indian ocean. He was not sure if such
cloud-like non-stellar objects (which he was the first to call 'Sahābi'
in Arabic or 'nebulous') were real astronomical objects or if they were
only the result of some bias in his observations. Could such diffuse
objects actually be detected at all with his detection technique?

He still had a few hours left until nightfall (when he would continue
his studies on the ecliptic) so he decided to find an answer to this
question. He had thoroughly studied Claudius Ptolemy's (90 - 168 A.D)
Almagest and had made lots of corrections to it, in particular in
measuring the brightness. Using his same experience, he was able to
measure a magnitude for the objects and wanted to simulate his
observation to see if a simulated object with the same brightness and
size could be detected in simulated noise with the same detection
technique. The general outline of the steps he wants to take are:

1. Make some mock profiles in an over-sampled image. The initial mock
image has to be over-sampled prior to convolution or other forms of
transformation in the image. Through his experiences, Sufi knew
that this is because the image of heavenly bodies is actually
transformed by the atmosphere or other sources outside the
atmosphere (for example, gravitational lenses) prior to being
sampled on an image. Since that transformation occurs on a
continuous grid, to best approximate it, he should do all the work
on a finer pixel grid. In the end he can resample the result to
the initially desired grid size.

2. Convolve the image with a point spread function (PSF, see *note
PSF::) that is over-sampled to the same resolution as the mock
image. Since he wants to finish in a reasonable time and the PSF
kernel will be very large due to oversampling, he has to use
frequency domain convolution which has the side effect of dimming
the edges of the image. So in the first step above he also has to
build the image to be larger by at least half the width of the PSF
convolution kernel on each edge.

3. With all the transformations complete, the image should be
resampled to the same size of the pixels in his detector.

4. He should remove those extra pixels on all edges to remove
frequency domain convolution artifacts in the final product.

5. He should add noise to the (until now, noise-less) mock image.
After all, all observations have noise associated with them.

Fortunately Sufi had heard of GNU Astronomy Utilities from a
colleague in Isfahan (where he worked) and had installed it on his
computer a year before. It had tools to do all the steps above. He had
used MakeProfiles before, but was not sure which columns he had chosen
in his user or system-wide configuration files for which parameters, see
*note Configuration files::. So to start his simulation, Sufi runs
MakeProfiles with the ‘-P’ option to make sure what columns in a catalog
MakeProfiles currently recognizes, and confirm the output image
parameters. In particular, Sufi is interested in the recognized columns
(shown below).

$ astmkprof -P

[[[ ... Truncated lines ... ]]]

# Output:
type float32 # Type of output: e.g., int16, float32, etc.
mergedsize 1000,1000 # Number of pixels along first FITS axis.
oversample 5 # Scale of oversampling (>0 and odd).

[[[ ... Truncated lines ... ]]]

# Columns, by info (see `--searchin'), or number (starting from 1):
ccol 2 # Coord. columns (one call for each dim.).
ccol 3 # Coord. columns (one call for each dim.).
fcol 4 # sersic (1), moffat (2), gaussian (3), point
# (4), flat (5), circumference (6), distance
# (7), custom-prof (8), azimuth (9),
# custom-img (10).
rcol 5 # Effective radius or FWHM in pixels.
ncol 6 # Sersic index or Moffat beta.
pcol 7 # Position angle.
qcol 8 # Axis ratio.
mcol 9 # Magnitude.
tcol 10 # Truncation in units of radius or pixels.

[[[ ... Truncated lines ... ]]]

In Gnuastro, column counting starts from 1, so the columns are ordered
such that the first column (number 1) can be an ID he specifies for each
object (and MakeProfiles ignores), each subsequent column is used for
another property of the profile. It is also possible to use column
names for the values of these options and change these defaults, but
Sufi preferred to stick to the defaults. Fortunately MakeProfiles has
the capability to also make the PSF which is to be used on the mock
image and using the ‘--prepforconv’ option, he can also make the mock
image to be larger by the correct amount and all the sources to be
shifted by the correct amount.

For his initial check he decides to simulate the nebula in the
Andromeda constellation. The night he was observing, the PSF had
roughly a FWHM of about 5 pixels, so as the first row (profile) in the
table below, he defines the PSF parameters. Sufi sets the radius column
(‘rcol’ above, fifth column) to ‘5.000’, he also chooses a Moffat
function for its functional form. Remembering how diffuse the nebula in
the Andromeda constellation was, he decides to simulate it with a mock
Sérsic index 1.0 profile. He wants the output to be 499 pixels by 499
pixels, so he can put the center of the mock profile in the central
pixel of the image which is the 250th pixel along both dimensions (note
that an even number does not have a "central" pixel).

Looking at his drawings of it, he decides a reasonable effective
radius for it would be 40 pixels on this image pixel scale (second row,
5th column below). He also sets the axis ratio (0.4) and position angle
(-25 degrees) to approximately correct values too, and finally he sets
the total magnitude of the profile to 3.44 which he had measured. Sufi
also decides to truncate both the mock profile and PSF at 5 times the
respective radius parameters. In the end he decides to put four stars
on the four corners of the image at very low magnitudes as a visual
scale. While he was preparing the catalog, one of his students
approached him and was also following the steps.

As described above, the catalog of profiles to build will be a table
(multiple columns of numbers) like below:

0 0.000 0.000 2 5 4.7 0.0 1.0 30.0 5.0
1 250.0 250.0 1 40 1.0 -25 0.4 3.44 5.0
2 50.00 50.00 4 0 0.0 0.0 0.0 6.00 0.0
3 450.0 50.00 4 0 0.0 0.0 0.0 6.50 0.0
4 50.00 450.0 4 0 0.0 0.0 0.0 7.00 0.0
5 450.0 450.0 4 0 0.0 0.0 0.0 7.50 0.0

This contains all the "data" to build the profile, and you can easily
pass it to Gnuastro's MakeProfiles: since Sufi already knows the columns
and expected values very good, he has placed the information in the
proper columns. However, when the student sees this, he just sees a
mumble-jumble of numbers! Generally, Sufi explains to the student that
even if you know the number positions and values very nicely today, in a
couple of months you will forget! It will then be very hard for you to
interpret the numbers properly. So you should never use naked data (or
data without any extra information).

Data (or information) that describes other data is called "metadata"!
One common example is column names (the name of a column is itself a
data element, but data that describes the lower-level data within that
column: how to interpret the numbers within it). Sufi explains to his
student that Gnuastro has a convention for adding metadata within a
plain-text file; and guides him to *note Gnuastro text table format::.
Because we do not want metadata to be confused with the actual data, in
a plain-text file, we start lines containing metadata with a '‘#’'. For
example, see the same data above, but this time with metadata for every
column:

# Column 1: ID [counter, u8] Identifier
# Column 2: X [pix, f32] Horizontal position
# Column 3: Y [pix, f32] Vertical position
# Column 4: PROFILE [name, u8] Radial profile function
# Column 5: R [pix, f32] Effective radius
# Column 6: N [n/a, f32] Sersic index
# Column 7: PA [deg, f32] Position angle
# Column 8: Q [n/a, f32] Axis ratio
# Column 9: MAG [log, f32] Magnitude
# Column 10: TRUNC [n/a, f32] Truncation (multiple of R)
0 0.000 0.000 2 5 4.7 0.0 1.0 30.0 5.0
1 250.0 250.0 1 40 1.0 -25 0.4 3.44 5.0
2 50.00 50.00 4 0 0.0 0.0 0.0 6.00 0.0
3 450.0 50.00 4 0 0.0 0.0 0.0 6.50 0.0
4 50.00 450.0 4 0 0.0 0.0 0.0 7.00 0.0
5 450.0 450.0 4 0 0.0 0.0 0.0 7.50 0.0

The numbers now make much more sense for the student! Before
continuing, Sufi reminded the student that even though metadata may not
be strictly/technically necessary (for the computer programs), metadata
are critical for human readers! Therefore, a good scientist should
never forget to keep metadata with any data that they create, use or
archive.

To start simulating the nebula, Sufi creates a directory named
‘simulationtest’ in his home directory. Note that the ‘pwd’ command
will print the "parent working directory" of the current directory (its
a good way to confirm/check your current location in the full file
system: it always starts from the root, or '‘/’').

$ mkdir ~/simulationtest
$ cd ~/simulationtest
$ pwd
/home/rahman/simulationtest

It is possible to use a plain-text editor to manually put the catalog
contents above into a plain-text file. But to easily automate catalog
production (in later trials), Sufi decides to fill the input catalog
with the redirection features of the command-line (or shell). Sufi's
student was not familiar with this feature of the shell! So Sufi
decided to do a fast demo; giving the following explanations while
running the commands:

Shell redirection allows you to "re-direct" the "standard output" of
a program (which is usually printed by the program on the command-line
during its execution; like the output of ‘pwd’ above) into a file. For
example, let's simply "echo" (or print to standard output) the line
"This is a test.":

$ echo "This is a test."
This is a test.

As you see, our statement was simply "echo"-ed to the standard output!
To redirect this sentence into a file (instead of simply printing it on
the standard output), we can simply use the ‘>’ character, followed by
the name of the file we want it to be dumped in.

$ echo "This is a test." > test.txt

This time, the ‘echo’ command did not print anything in the terminal.
Instead, the shell (command-line environment) took the output, and
"re-directed" it into a file called ‘test.txt’. Let's confirm this with
the ‘ls’ command (‘ls’ is short for "list" and will list all the files
in the current directory):

$ ls
test.txt

Now that you confirm the existence of ‘test.txt’, you can see its
contents with the ‘cat’ command (short for "concatenation"; because it
can also merge multiple files together):

$ cat test.txt
This is a test.

Now that we have written our first line in ‘test.txt’, let's try adding
a second line (do not forget that our final catalog of objects to
simulate will have multiple lines):

$ echo "This is my second line." > test.txt
$ cat test.txt
This is my second line.

As you see, the first line that you put in the file is no longer
present! This happens because '‘>’' always starts dumping content to a
file from the start of the file. In effect, this means that any
possibly pre-existing content is over-written by the new content! To
append new lines (or dumping new content at the end of existing
content), you can use '‘>>’'. for example, with the commands below,
first we will write the first sentence (using '‘>’'), then use '‘>>’' to
add the second and third sentences. Finally, we will print the contents
of ‘test.txt’ to confirm that all three lines are preserved.

$ echo "My first sentence." > test.txt
$ echo "My second sentence." >> test.txt
$ echo "My third sentence." >> test.txt
$ cat test.txt
My first sentence.
My second sentence.
My third sentence.

The student thanked Sufi for this explanation and now feels more
comfortable with redirection. Therefore Sufi continues with the main
project. But before that, he deletes the temporary test file:

$ rm test.txt

To put the catalog of profile data and their metadata (that was
described above) into a file, Sufi uses the commands below. While Sufi
was writing these commands, the student complained that "I could have
done in this in a text editor". Sufi reminded the student that it is
indeed possible; but it requires manual intervention. The advantage of
a solution like below is that it can be automated (for example, adding
more rows; for more profiles in the final image).

$ echo "# Column 1: ID [counter, u8] Identifier" > cat.txt
$ echo "# Column 2: X [pix, f32] Horizontal position" >> cat.txt
$ echo "# Column 3: Y [pix, f32] Vertical position" >> cat.txt
$ echo "# Column 4: PROF [name, u8] Radial profile function" \
>> cat.txt
$ echo "# Column 5: R [pix, f32] Effective radius" >> cat.txt
$ echo "# Column 6: N [n/a, f32] Sersic index" >> cat.txt
$ echo "# Column 7: PA [deg, f32] Position angle" >> cat.txt
$ echo "# Column 8: Q [n/a, f32] Axis ratio" >> cat.txt
$ echo "# Column 9: MAG [log, f32] Magnitude" >> cat.txt
$ echo "# Column 10: TRUNC [n/a, f32] Truncation (multiple of R)" \
>> cat.txt
$ echo "0 0.000 0.000 2 5 4.7 0.0 1.0 30.0 5.0" >> cat.txt
$ echo "1 250.0 250.0 1 40 1.0 -25 0.4 3.44 5.0" >> cat.txt
$ echo "2 50.00 50.00 4 0 0.0 0.0 0.0 6.00 0.0" >> cat.txt
$ echo "3 450.0 50.00 4 0 0.0 0.0 0.0 6.50 0.0" >> cat.txt
$ echo "4 50.00 450.0 4 0 0.0 0.0 0.0 7.00 0.0" >> cat.txt
$ echo "5 450.0 450.0 4 0 0.0 0.0 0.0 7.50 0.0" >> cat.txt

To make sure if the catalog's content is correct (and there was no typo
for example!), Sufi runs '‘cat cat.txt’', and confirms that it is
correct.

Now that the catalog is created, Sufi is ready to call MakeProfiles
to build the image containing these objects. He looks into his records
and finds that the zero point magnitude for that night, and that
particular detector, was 18 magnitudes. The student was a little
confused on the concept of zero point, so Sufi pointed him to *note
Brightness flux magnitude::, which the student can study in detail
later. Sufi therefore runs MakeProfiles with the command below:

$ astmkprof --prepforconv --mergedsize=499,499 --zeropoint=18.0 cat.txt
MakeProfiles 0.22 started on Sat Oct 6 16:26:56 953
- 6 profiles read from cat.txt
- Random number generator (RNG) type: ranlxs1
- Basic RNG seed: 1652884540
- Using 12 threads.
---- row 3 complete, 5 left to go
---- row 4 complete, 4 left to go
---- row 6 complete, 3 left to go
---- row 5 complete, 2 left to go
---- ./0_cat_profiles.fits created.
---- row 1 complete, 1 left to go
---- row 2 complete, 0 left to go
- ./cat_profiles.fits created. 0.092573 seconds
-- Output: ./cat_profiles.fits
MakeProfiles finished in 0.293644 seconds

Sufi encourages the student to read through the printed output. As
the statements say, two FITS files should have been created in the
running directory. So Sufi ran the command below to confirm:

$ ls
0_cat_profiles.fits cat_profiles.fits cat.txt

The file ‘0_cat_profiles.fits’ is the PSF Sufi had asked for, and
‘cat_profiles.fits’ is the image containing the main objects in the
catalog. Sufi opened the main image with the command below (using SAO
DS9):

$ astscript-fits-view cat_profiles.fits --ds9scale=95

The student could clearly see the main elliptical structure in the
center. However, the size of ‘cat_profiles.fits’ was surprising for the
student, instead of 499 by 499 (as we had requested), it was 2615 by
2615 pixels (from the command below):

$ astfits cat_profiles.fits
Fits (GNU Astronomy Utilities) 0.22
Run on Sat Oct 6 16:26:58 953
-----
HDU (extension) information: 'cat_profiles.fits'.
Column 1: Index (counting from 0, usable with '--hdu').
Column 2: Name ('EXTNAME' in FITS standard, usable with '--hdu').
Column 3: Image data type or 'table' format (ASCII or binary).
Column 4: Size of data in HDU.
-----
0 MKPROF-CONFIG no-data 0
1 Mock profiles float32 2615x2615

So Sufi explained why oversampling is important in modeling, especially
for parts of the image where the flux change is significant over a
pixel. Recall that when you oversample the model (for example, by 5
times), for every desired pixel, you get 25 pixels ($5\times5$). Sufi
then explained that after convolving (next step below) we will
down-sample the image to get our originally desired size/resolution.

After seeing the image, the student complained that only the large
elliptical model for the Andromeda nebula can be seen in the center. He
could not see the four stars that we had also requested in the catalog.
So Sufi had to explain that the stars are there in the image, but the
reason that they are not visible when looking at the whole image at
once, is that they only cover a single pixel! To prove it, he centered
the image around the coordinates 2308 and 2308, where one of the stars
is located in the over-sampled image [you can do this in ‘ds9’ by
selecting "Pan" in the "Edit" menu, then clicking around that position].
Sufi then zoomed in to that region and soon, the star's non-zero pixel
could be clearly seen.

Sufi explained that the stars will take the shape of the PSF (cover
an area of more than one pixel) after convolution. If we did not have
an atmosphere and we did not need an aperture, then stars would only
cover a single pixel with normal CCD resolutions. So Sufi convolved the
image with this command:

$ astconvolve --kernel=0_cat_profiles.fits cat_profiles.fits \
--output=cat_convolved.fits
Convolve started on Sat Oct 6 16:35:32 953
- Using 8 CPU threads.
- Input: cat_profiles.fits (hdu: 1)
- Kernel: 0_cat_profiles.fits (hdu: 1)
- Input and Kernel images padded. 0.075541 seconds
- Images converted to frequency domain. 6.728407 seconds
- Multiplied in the frequency domain. 0.040659 seconds
- Converted back to the spatial domain. 3.465344 seconds
- Padded parts removed. 0.016767 seconds
- Output: cat_convolved.fits
Convolve finished in: 10.422161 seconds

When convolution finished, Sufi opened ‘cat_convolved.fits’ and the four
stars could be easily seen now:

$ astscript-fits-view cat_convolved.fits --ds9scale=95

It was interesting for the student that all the flux in that single
pixel is now distributed over so many pixels (the sum of all the pixels
in each convolved star is actually equal to the value of the single
pixel before convolution). Sufi explained how a PSF with a larger FWHM
would make the points even wider than this (distributing their flux in a
larger area). With the convolved image ready, they were prepared to
resample it to the original pixel scale Sufi had planned [from the ‘$
astmkprof -P’ command above, recall that MakeProfiles had over-sampled
the image by 5 times]. Sufi explained the basic concepts of warping the
image to his student and ran Warp with the following command:

$ astwarp --scale=1/5 --centeroncorner cat_convolved.fits
Warp started on Sat Oct 6 16:51:59 953
Using 8 CPU threads.
Input: cat_convolved.fits (hdu: 1)
matrix:
0.2000 0.0000 0.4000
0.0000 0.2000 0.4000
0.0000 0.0000 1.0000

$ astfits cat_convolved_scaled.fits --quiet
0 WARP-CONFIG no-data 0
1 Warped float32 523x523

‘cat_convolved_scaled.fits’ now has the correct pixel scale. However,
the image is still larger than what we had wanted, it is $523\times523$
pixels (not our desired $499\times499$). The student is slightly
confused, so Sufi also resamples the PSF with the same scale by running

$ astwarp --scale=1/5 --centeroncorner 0_cat_profiles.fits
$ astfits 0_cat_profiles_scaled.fits --quiet
0 WARP-CONFIG no-data 0
1 Warped float32 25x25

Sufi notes that $25=12+12+1$ and that $523=499+12+12$. He goes on to
explain that frequency space convolution will dim the edges and that is
why he added the ‘--prepforconv’ option to MakeProfiles above. Now that
convolution is done, Sufi can remove those extra pixels using Crop with
the command below. Crop's ‘--section’ option accepts coordinates
inclusively and counting from 1 (according to the FITS standard), so the
crop region's first pixel has to be 13, not 12.

$ astcrop cat_convolved_scaled.fits --section=13:*-12,13:*-12 \
--mode=img --zeroisnotblank
Crop started on Sat Oct 6 17:03:24 953
- Read metadata of 1 image. 0.001304 seconds
---- ...nvolved_scaled_cropped.fits created: 1 input.
Crop finished in: 0.027204 seconds

To fully convince the student, Sufi checks the size of the output of the
crop command above:

$ astfits cat_convolved_scaled_cropped.fits --quiet
0 n/a no-data 0
1 n/a float32 499x499

Finally, ‘cat_convolved_scaled_cropped.fits’ is $499\times499$ pixels
and the mock Andromeda galaxy is centered on the central pixel. This is
the same dimensions as Sufi had desired in the beginning. All this
trouble was certainly worth it because now there is no dimming on the
edges of the image and the profile centers are more accurately sampled.

The final step to simulate a real observation would be to add noise
to the image. Sufi set the zero point magnitude to the same value that
he set when making the mock profiles and looking again at his
observation log, he had measured the background flux near the nebula had
a _per-pixel_ magnitude of 7 that night. For more on how the background
value determines the noise, see *note Noise basics::. So using these
values he ran Arithmetic's ‘mknoise-sigma-from-mean’ operator, and with
the second command, he visually inspected the image. The
‘mknoise-sigma-from-mean’ operator takes the noise standard deviation in
linear units, not magnitudes (which are logarithmic). Therefore within
the same Arithmetic command, he has converted the sky background
magnitude to counts using Arithmetic's ‘counts-to-mag’ operator.

$ astarithmetic cat_convolved_scaled_cropped.fits \
7 18 mag-to-counts mknoise-sigma-from-mean \
--output=out.fits

$ astscript-fits-view out.fits

The ‘out.fits’ file now contains the noised image of the mock catalog
Sufi had asked for. The student had not observed the nebula in the sky,
so when he saw the mock image in SAO DS9 (with the second command
above), he understood why Sufi was dubious: it was very diffuse!

Seeing how the ‘--output’ option allows the user to specify the name
of the output file, the student was confused and wanted to know why Sufi
had not used it more regularly before? Sufi explained that for
intermediate steps, you can rely on the automatic output of the programs
(see *note Automatic output::). Doing so will give all the intermediate
files a similar basic name structure, so in the end you can simply
remove them all with the Shell's capabilities, and it will be familiar
for other users. So Sufi decided to show this to the student by making
a shell script from the commands he had used before.

The command-line shell has the capability to read all the separate
input commands from a file. This is useful when you want to do the same
thing multiple times, with only the names of the files or minor
parameters changing between the different instances. Using the shell's
history (by pressing the up keyboard key) Sufi reviewed all the commands
and then he retrieved the last 5 commands with the ‘$ history 5’
command. He selected all those lines he had input and put them in a
text file named ‘mymock.sh’. Then he defined the ‘edge’ and ‘base’
shell variables for easier customization later, and before every
command, he added some comments (lines starting with <#>) for future
readability. Finally, Sufi pointed the student to Gnuastro's *note
General program usage tutorial::, which has a full section on *note
Writing scripts to automate the steps::.

#!/bin/bash

edge=12
base=cat

# Stop running next commands if one fails.
set -e

# Remove any (possibly) existing output (from previous runs)
# before starting.
rm -f out.fits

# Run MakeProfiles to create an oversampled FITS image.
astmkprof --prepforconv --mergedsize=499,499 --zeropoint=18.0 \
"$base".txt

# Convolve the created image with the kernel.
astconvolve "$base"_profiles.fits \
--kernel=0_"$base"_profiles.fits \
--output="$base"_convolved.fits

# Scale the image back to the intended resolution.
astwarp --scale=1/5 --centeroncorner "$base"_convolved.fits

# Crop the edges out (dimmed during convolution). '--section'
# accepts inclusive coordinates, so the start of the section
# must be one pixel larger than its end.
st_edge=$(( edge + 1 ))
astcrop "$base"_convolved_scaled.fits --zeroisnotblank \
--mode=img --section=$st_edge:*-$edge,$st_edge:*-$edge

# Add noise to the image.
$ astarithmetic "$base"_convolved_scaled_cropped.fits \
7 18 mag-to-counts mknoise-sigma-from-mean \
--output=out.fits

# Remove all the temporary files.
rm 0*.fits "$base"*.fits

He used this chance to remind the student of the importance of
comments in code or shell scripts! Just like metadata in a dataset,
when writing the code, you have a good mental picture of what you are
doing, so writing comments might seem superfluous and excessive.
However, in one month when you want to re-use the script, you have lost
that mental picture and remembering it can be time-consuming and
frustrating. The importance of comments is further amplified when you
want to share the script with a friend/colleague. So it is good to
accompany any step of a script, or code, with useful comments while you
are writing it (create a good mental picture of why you are doing
something: do not just describe the command, but its purpose).

Sufi then explained to the eager student that you define a variable
by giving it a name, followed by an ‘=’ sign and the value you want.
Then you can reference that variable from anywhere in the script by
calling its name with a ‘$’ prefix. So in the script whenever you see
‘$base’, the value we defined for it above is used. If you use advanced
editors like GNU Emacs or even simpler ones like Gedit (part of the
GNOME graphical user interface) the variables will become a different
color which can really help in understanding the script. We have put
all the ‘$base’ variables in double quotation marks (‘"’) so the
variable name and the following text do not get mixed, the shell is
going to ignore the ‘"’ after replacing the variable value. To make the
script executable, Sufi ran the following command:

$ chmod +x mymock.sh

Then finally, Sufi ran the script, simply by calling its file name:

$ ./mymock.sh

After the script finished, the only file remaining is the ‘out.fits’
file that Sufi had wanted in the beginning. Sufi then explained to the
student how he could run this script anywhere that he has a catalog if
the script is in the same directory. The only thing the student had to
modify in the script was the name of the catalog (the value of the
‘base’ variable in the start of the script) and the value to the ‘edge’
variable if he changed the PSF size. The student was also happy to hear
that he will not need to make it executable again when he makes changes
later, it will remain executable unless he explicitly changes the
executable flag with ‘chmod’.

The student was really excited, since now, through simple shell
scripting, he could really speed up his work and run any command in any
fashion he likes allowing him to be much more creative in his works.
Until now he was using the graphical user interface which does not have
such a facility and doing repetitive things on it was really frustrating
and some times he would make mistakes. So he left to go and try
scripting on his own computer. He later reminded Sufi that the second
tutorial in the Gnuastro book as more complex commands in data analysis,
and a more advanced introduction to scripting (see *note General program
usage tutorial::).

Sufi could now get back to his own work and see if the simulated
nebula which resembled the one in the Andromeda constellation could be
detected or not. Although it was extremely faint(2). Therefore, Sufi
ran Gnuastro's detection software (*note NoiseChisel::) to see if this
object is detectable or not. NoiseChisel's output (‘out_detected.fits’)
is a multi-extension FITS file, so he used Gnuastro's
‘astscript-fits-view’ program in the second command to see the output:

$ astnoisechisel out.fits

$ astscript-fits-view out_detected.fits

In the "Cube" window (that was opened with DS9), if Sufi clicked on
the "Next" button to see the pixels that were detected to contain
significant signal. Fortunately the nebula's shape was detectable and
he could finally confirm that the nebula he kept in his notebook was
actually observable. He wrote this result in the draft manuscript that
would later become "Book of fixed stars"(3).

He still had to check the other nebula he saw from Yemen and several
other such objects, but they could wait until tomorrow (thanks to the
shell script, he only has to define a new catalog). It was nearly
sunset and they had to begin preparing for the night's measurements on
the ecliptic.

---------- Footnotes ----------

(1) In Latin Sufi is known as Azophi. He was an Iranian astronomer.
His manuscript "Book of fixed stars" contains the first recorded
observations of the Andromeda galaxy, the Large Magellanic Cloud and
seven other non-stellar or 'nebulous' objects.

(2) The brightness of a diffuse object is added over all its pixels
to give its final magnitude, see *note Brightness flux magnitude::. So
although the magnitude 3.44 (of the mock nebula) is orders of magnitude
brighter than 6 (of the stars), the central galaxy is much fainter. Put
another way, the brightness is distributed over a large area in the case
of a nebula.

(3) <https://en.wikipedia.org/wiki/Book_of_Fixed_Stars>

File: gnuastro.info, Node: Detecting lines and extracting spectra in 3D data, Next: Color images with full dynamic range, Prev: Sufi simulates a detection, Up: Tutorials

2.5 Detecting lines and extracting spectra in 3D data
=====================================================

3D data cubes are an increasingly common format of data products in
observational astronomy. As opposed to 2D images (where each 2D
"picture element" or "pixel" covers an infinitesimal area on the surface
of the sky), 3D data cubes contain "volume elements" or "voxels" that
are also connected in a third dimension.

The most common case of 3D data in observational astrophysics is when
the first two dimensions are spatial (RA and Dec on the sky), and the
third dimension is wavelength. This type of data is generically (also
outside of astronomy) known as Hyperspectral imaging(1). For example
high-level data products of Integral Field Units (IFUs) like MUSE(2) in
the optical, ACIS(3) in the X-ray, or in the radio where most data are
3D cubes.

In this tutorial, we'll use a small crop of a reduced deep MUSE cube
centered on the Abell 370 (https://en.wikipedia.org/wiki/Abell_370)
galaxy cluster from the Pilot-WINGS survey; see Lagattuta et al. 2022
(https://arxiv.org/abs/2202.04663). Abell 370 has a spiral galaxy in
its background that is stretched due to the cluster's gravitational
potential to create a beautiful arch. If you haven't seen it yet, have
a look at some of its images in the Wikipedia link above before
continuing.

The Pilot-WINGS survey data are available in its webpage(4). The
cube of the _core_ region is 10.2GBs. This can be prohibitively large
to download (and later process) on many networks and smaller computers.
Therefore, in this demonstration we won't be using the full cube. We
have prepared a small crop(5) of the full cube that you can download
with the first command below. The randomly selected crop is centered on
(RA,Dec) of (39.96769,-1.58930), with a width of about 27 arcseconds.

$ mkdir tutorial-3d
$ cd tutorial-3d
$ wget http://akhlaghi.org/data/a370-crop.fits # Downloads 287 MB

In the sections below, we will first review how you can visually
inspect a 3D datacube in DS9 and interactively see the spectra of any
region. We will then subtract the continuum emission, detect the
emission-lines within this cube and extract their spectra. We will
finish with creating pseudo narrow-band images optimized for some of the
emission lines.

* Menu:

* Viewing spectra and redshifted lines:: Interactively see the spectra of an object
* Sky lines in optical IFUs:: How to see sky lines in a cube.
* Continuum subtraction:: Removing the continuum from a data cube.
* 3D detection with NoiseChisel:: Finding emission-lines and their spectra.
* 3D measurements and spectra:: Measuring 3d properties including spectra.
* Extracting a single spectrum and plotting it:: Extracting a single vector row.
* Pseudo narrow-band images:: Collapsing the third dimension into a 2D image.

---------- Footnotes ----------

(1) <https://en.wikipedia.org/wiki/Hyperspectral_imaging>

(2) <https://en.wikipedia.org/wiki/Multi-unit_spectroscopic_explorer>

(3) <https://en.wikipedia.org/wiki/Advanced_CCD_Imaging_Spectrometer>

(4) <https://astro.dur.ac.uk/~hbpn39/pilot-wings.html>

(5) You can download the full cube and create the crop your self with
the commands below. Due to the decompression of the +10GB file that is
necessary for the compressed downloaded file (note that its suffix is
‘.fits.gz’), the Crop command will take a little long.
$ wget https://astro.dur.ac.uk/~hbpn39/pilotWINGS/A370_PilotWINGS_CORE.fits.gz
$ astcrop A370_PilotWINGS_CORE.fits.gz -hDATA --mode=img \
--section=200:300,100:200 -oa370-crop.fits --metaname=DATA
$ astcrop A370_PilotWINGS_CORE.fits.gz -hSTAT --mode=img --append \
--section=200:300,100:200 -oa370-crop.fits --metaname=STAT

File: gnuastro.info, Node: Viewing spectra and redshifted lines, Next: Sky lines in optical IFUs, Prev: Detecting lines and extracting spectra in 3D data, Up: Detecting lines and extracting spectra in 3D data

2.5.1 Viewing spectra and redshifted lines
------------------------------------------

In *note Detecting lines and extracting spectra in 3D data:: we
downloaded a small crop from the Pilot-WINGS survey of Abell 370
cluster; observed with MUSE. In this section, we will review how you can
visualize/inspect a datacube using that example. With the first command
below, we'll open DS9 such that each 2D slice of the cube (at a fixed
wavelength) is seen as a single image. If you move the slider in the
"Cube" window (that also opens), you can view the same field at
different wavelengths. We are ending the first command with a '‘&’' so
you can continue viewing DS9 while using the command-line (press one
extra ‘ENTER’ to see the prompt). With the second command, you can see
that the spacing between each slice is $1.25\times10^{-10}$ meters (or
1.25 Angstroms).

$ astscript-fits-view a370-crop.fits -h1 --ds9scale="limits -5 20" &

$ astfits a370-crop.fits --pixelscale
Basic info. for --pixelscale (remove info with '--quiet' or '-q')
Input: a370-crop.fits (hdu 1) has 3 dimensions.
Pixel scale in each FITS dimension:
1: 5.55556e-05 (deg/pixel) = 0.2 (arcsec/pixel)
2: 5.55556e-05 (deg/pixel) = 0.2 (arcsec/pixel)
3: 1.25e-10 (m/slice)
Pixel area (on each 2D slice) :
3.08642e-09 (deg^2) = 0.04 (arcsec^2)
Voxel volume:
3.85802e-19 (deg^2*m) = 5e-12 (arcsec^2*m) = 0.05 (arcsec^2*A)

In the DS9 "Cube" window, you will see two numbers on the two sides
of the scroller. The left number is the wavelength in meters (WCS
coordinate in 3rd dimension) and the right number is the slice number
(slice number or array coordinates in 3rd dimension). You can manually
edit any of these numbers and press ENTER to go to that slice in any
coordinate system. If you want to go one-by-one, simply press the
"Next" button. The first few slides are very noisy, but in the rest the
noise level decreases and the galaxies are more obvious.

As you slide between the different wavelengths, you see that the
noise-level is not constant and in some slices, the sky noise is very
strong (for example, go to slice 3201 and press the "Next" button). We
will discuss these issues below (in *note Sky lines in optical IFUs::).
To view the spectra of a region in DS9 take the following steps:

1. Click somewhere on the image (to make sure DS9 receives your
keyboard inputs), then press ‘Ctrl+R’ to activate regions and click
on the brightest galaxy of this cube (center-right, at RA, Dec of
39.9659175 and -1.5893075).
2. A thin green circle will show up; this is called a "region" in DS9.
3. Double-click on the region, and you will see a "Circle" window.
4. Within the "Circle" window, click on the "Analysis" menu and select
"Plot 3D".
5. A second "Circle" window will open that shows the spectra within
your selected region. This is just the sum of values on each slice
within the region.
6. Don't close the second "circle" window (that shows the spectrum).
Click and hold the region in DS9, and move it to other objects
within the cube. You will see that the spectrum changes as you
move the region, and you can see that different objects have very
different spectra. You can even see the spectra of only one part
of a galaxy, not the whole galaxy.
7. Take the region back to the first (brightest) galaxy that we
originally started with.
8. Slide over different wavelengths in the "Cube" window, you will see
the light-blue line moving through the spectrum as you slide to
different wavelengths. This line shows the wavelength of the
displayed image in the main window over the spectra.
9. The strongest emission line in this galaxy appears to be around
8500 Angstroms or $8.5\times10^{-7}$ meters. From the position of
the Balmer break (https://en.wikipedia.org/wiki/Balmer_jump)
(blue-ward of 5000 Angstroms for this galaxy), the strong seems to
be H-alpha.
10. To confirm that this is H-alpha, you can select the "Edit" menu in
the spectrum window and select "Zoom".
11. Double-click and hold (for next step also) somewhere before the
strongest line and slightly above the continuum (for example at
‘8E-07’ in the horizontal and
$60\times10^{-20}$erg/Angstrom/cm$^2$/s on the vertical). As you
move your cursor (while holding), you will see a rectangular box
getting created.
12. Move the bottom-left corner of the box to somewhere after the
strongest line and below the continuum. For example at ‘9E-07’ and
$20\times10^{-20}$erg/Angstrom/cm$^2$/s.
13. Once you remove your finger from the mouse/touchpad, it will
zoom-in to that part of the spectrum.
14. To zoom out to the full spectrum, just press the right mouse
button over the spectra (or tap with two fingers on a touchpad).
15. Select that zoom-box again to see the brightest line much more
clearly. You can also see the two lines of the Nitrogen II doublet
that sandwich H-alpha. Beside its relative position to the Balmer
break, this is further evidence that the strongest line is H-alpha.
16. Let's have a look at the galaxy in its best glory: right over the
H-alpha line: Move the wavelength slider accurately (by pressing
the "Previous" or "Next" buttons) such that the blue line falls in
the middle of the H-alpha line. We see that the wavelength at this
slice is ‘8.56593e-07’ meters or 8565.93 Angstroms. Please compare
the image of the galaxy at this wavelength with the wavelengths
before (by pressing "Next" or "Previous"). You will also see that
it is much more extended and brighter than other wavelengths!
H-alpha shows the un-obscured star formation of the galaxy!

*Automaticly going to next slice:* When you want to get a general
feeling of the cube, pressing the "Next" button many times is annoying
and slow. To automatically shift between the slices, you can press the
"Play" button in the DS9 "Cube" window. You can adjust the time it
stays on each slice by clicking on the "Interval" menu and selecting
lower values.

Knowing that this is H-alpha at 8565.93 Angstroms, you can get the
redshift of the galaxy with the first command below and the location of
all other expected lines in Gnuastro's spectral line database with the
second command. Because there are many lines in the second command
(more than 200!), with the third command, we'll only limit it to the
Balmer series (that start with ‘H-’) using ‘grep’. The output of the
second command prints the metadata on the top (that is not shown any
more in the third command due to the ‘grep’ call). To be complete, the
first column is the observed wavelength of the given line in the given
redshift and the second column is the name of the line.

# Redshift where H-alpha falls on 8565.93.
$ astcosmiccal --obsline=H-alpha,8565.93 --usedredshift
0.305221

# Wavelength of all lines in Gnuastro's database at this redshift
$ astcosmiccal --obsline=H-alpha,8565.93 --listlinesatz

# Only the Balmer series (Lines starting with 'H-'; given to Grep).
$ astcosmiccal --obsline=H-alpha,8565.93 --listlinesatz | grep H-
4812.13 H-19
4818.29 H-18
4825.61 H-17
4834.36 H-16
4844.95 H-15
4857.96 H-14
4874.18 H-13
4894.79 H-12
4921.52 H-11
4957.1 H-10
5006.03 H-9
5076.09 H-8
5181.83 H-epsilon
5353.68 H-delta
5665.27 H-gamma
6345.11 H-beta
8565.93 H-alpha
4758.84 H-limit

Zoom-out to the full spectrum and move the displayed slice to the
location of the first emission line that is blue-ward (at shorter
wavelengths) of H-alpha (at around 6300 Angstroms) and follow the
previous steps to confirm that you are on its center. You will see that
it falls exactly on $6.34468\times10^{-7}$ m or 6344.68 Angstroms. Now,
have a look at the Balmer lines above. You have found the H-beta line!

The rest of the Balmer series
(https://en.wikipedia.org/wiki/Balmer_series) that you see in the list
above (like H-gamma, H-delta and H-epsilon) are visible only as
absorption lines. Please check their location by moving the blue line
on the wavelengths above and confirm the spectral absorption lines with
the ones above. The Balmer break is caused by the fact that these
stronger Balmer absorption lines become too close to each other.

Looking back at the full spectrum, you can also confirm that the only
other relatively strong emission line in this galaxy, that is on the
blue side of the spectrum is the weakest OII line that is approximately
located at 4864 Angstroms in the observed spectra of this galaxy. The
numbers after the various OII emission lines show their rest-frame
wavelengths ("OII" can correspond to many electron transitions, so we
should be clear about which one we are talking about).

$ astcosmiccal --obsline=H-alpha,8565.93 --listlinesatz | grep O-II-
4863.3 O-II-3726
4866.93 O-II-3728
5634.82 O-II-4317
5762.42 O-II-4414
9554.21 O-II-7319
9568.22 O-II-7330

Please stop here and spend some time on doing the exercise above on
other galaxies in the this cube to get a feeling of types of galaxy
spectral features (and later on the full/large cube). You will notice
that only star-forming galaxies have such strong emission lines! If you
enjoy it, go get the full non-cropped cube and investigate the spectra,
redshifts and emission/absorption lines of many more galaxies.

But going into those higher-level details of the physical meaning of
the spectra (as intriguing as they are!) is beyond the scope of this
tutorial. So we have to stop at this stage unfortunately. Now that you
have a relatively good feeling of this small cube, let's start doing
some analysis to extract the spectra of the objects in this cube.

File: gnuastro.info, Node: Sky lines in optical IFUs, Next: Continuum subtraction, Prev: Viewing spectra and redshifted lines, Up: Detecting lines and extracting spectra in 3D data

2.5.2 Sky lines in optical IFUs
-------------------------------

As we were visually inspecting the cube in *note Viewing spectra and
redshifted lines::, we noticed some slices with very bad noise. They
will later affect our detection within the cube, so in this section
let's have a fast look at them here. We'll start by looking at the two
cubes within the downloaded FITS file:

$ astscript-fits-view a370-crop.fits

The cube on the left is the same cube we studied before. The cube on
the right (which is called ‘STAT’) shows the variance of each voxel. Go
to slice 3195 and press "Next" to view the subsequent slices. Initially
(for the first 5 or 6 slices), the noise looks reasonable. But as you
pass slice 3206, you will see that the noise becomes very bad in both
cubes. It stays like this until about slice 3238! As you go through
the whole cube, you will notice that these slices are much more frequent
in the reddest wavelengths.

These slices are affected by the emission lines from our own
atmosphere! The atmosphere's emission in these wavelengths
significantly raises the background level in these slices. As a result,
the Poisson noise also increases significantly (see *note Photon
counting noise::). During the data reduction, the excess background
flux of each slice is removed as the Sky (or the mean of undetected
pixels, see *note Sky value::). However, the increased Poisson noise
(scatter of pixel values) remains!

To see spectrum of the sky emission lines, simply put a region
somewhere in the ‘STAT’ cube and generate its spectrum (as we did in
*note Viewing spectra and redshifted lines::). You will clearly see the
comb-like shape of atmospheric emission lines and can use this to know
where to expect them.

File: gnuastro.info, Node: Continuum subtraction, Next: 3D detection with NoiseChisel, Prev: Sky lines in optical IFUs, Up: Detecting lines and extracting spectra in 3D data

2.5.3 Continuum subtraction
---------------------------

In *note Viewing spectra and redshifted lines::, we visually inspected
some of the most prominent emission lines of the brightest galaxy of the
demo MUSE cube (see *note Detecting lines and extracting spectra in 3D
data::). Here, we will remove the "continuum" flux from under the
emission lines to see them more distinctly.

Within a spectra, the continuum is the local "background" flux in the
third/wavelength dimension. In other words, it is the flux that would
be present at that wavelength if the emission line didn't exist.
Therefore, to accurately measure the flux of the emission line, we first
need to subtract the continuum. One crude way of estimating the
continuum flux at every slice is to use the sigma-clipped median value
of that same pixel in the $\pm{N/2}$ slides around it (for more on
sigma-clipping, see *note Sigma clipping::).

In this case, $N=100$ should be a good first approximate (since it is
much larger than any of the absorption or emission lines). With the
first command below, let's use Arithmetic's filtering operators for
estimating the sigma-clipped median only along the third dimension for
every pixel in every slice (see *note Filtering operators::). With the
second command, have a look at the filtered cube and spectra. Note that
the first command is computationally expensive and may take a minute or
so.

$ astarithmetic a370-crop.fits set-i --output=filtered.fits \
3 0.2 1 1 100 i filter-sigclip-median

$ astscript-fits-view filtered.fits -h1 --ds9scale="limits -5 20"

Looking at the filtered cube above, and sliding through the different
wavelengths, you will see the noise in each slice has been significantly
reduced! This is expected because each pixel's value is now calculated
from 100 others (along the third dimension)! Using the same steps as
*note Viewing spectra and redshifted lines::, plot the spectra of the
brightest galaxy. Then, have a look at its spectra. You see that the
emission lines have been significantly smoothed out to become almost(1)
invisible.

You can now subtract this "continuum" cube from the input cube to
create the emission-line cube. In fact, as you see below, we can do it
in a single Arithmetic command (blending the filtering and subtraction
in one command). Note how the only difference with the previous
Arithmetic command is that we added an ‘i’ before the ‘3’ and a ‘-’
after ‘filter-sigclip-median’. For more on Arithmetic's powerful
notation, see *note Reverse polish notation::. With the second command
below, let's view the input _and_ continuum-subtracted cubes together:

$ astarithmetic a370-crop.fits set-i --output=no-continuum.fits \
i 3 0.2 1 1 100 i filter-sigclip-median -

$ astscript-fits-view a370-crop.fits no-continuum.fits -h1 \
--ds9scale="limits -5 20"

Once the cubes are open, slide through the different wavelengths.
Comparing the left (input) and right (continuum-subtracted) slices, you
will rarely see any galaxy in the continuum-subtracted one! As its name
suggests, the continuum flux is continuously present in all the
wavelengths (with gradual change)! But the continuum has been
subtracted now; so in the right-side image, you don't see anything on
wavelengths that don't contain a spectral emission line. Some dark
regions also appear; these are absorption lines! Please spend a few
minutes sliding through the wavelengths and seeing how the emission
lines pop-up and disappear again. It is almost like scuba diving, with
fish appearing out of nowhere and passing by you.

Let's go to slice 3046 (corresponding to 8555.93 Angstroms; just
before the H-alpha line for the brightest galaxy in *note Viewing
spectra and redshifted lines::). Now press the "Next" button to change
slices one by one until there is no more emission in the brightest
galaxy. As you go to redder slices, you will see that not only does the
brightness increase, but the position of the emission also changes.
This is the Doppler effect
(https://en.wikipedia.org/wiki/Doppler_effect) caused by the rotation of
the galaxy: the side that rotating towards us gets blue-shifted to bluer
slices and the one that is going away from us gets redshifted to redder
slices. If you go to the emission lines of the other galaxies, you will
see that they move with a different angle! We can use this to derive
the galaxy's rotational properties and kinematics (Gnuastro doesn't have
this feature yet).

To see the Doppler shift in the spectrum, plot the spectrum over the
top-side of the galaxy (which is visible in slice 3047). Then Zoom-in
to the H-alpha line (as we did in *note Viewing spectra and redshifted
lines::) and press "Next" until you reach the end of the H-alpha
emission-line. You see that by the time H-alpha disappears in the
spectrum, within the cube, the emission shifts in the vertical axis by
about 15 pixels! Then, move the region across the same path that the
emission passed. You will clearly see that the H-alpha and Nitrogen II
lines also move with you, in the zoomed-in spectra. Again, try this for
several other emission lines, and several other galaxies to get a good
feeling of this important concept when using hyper-spectral 3D data.

---------- Footnotes ----------

(1) For more on why Sigma-clipping is only a crude solution to
background removal, see Akhlaghi and Ichikawa 2015
(https://arxiv.org/abs/1505.01664).

File: gnuastro.info, Node: 3D detection with NoiseChisel, Next: 3D measurements and spectra, Prev: Continuum subtraction, Up: Detecting lines and extracting spectra in 3D data

2.5.4 3D detection with NoiseChisel
-----------------------------------

In *note Continuum subtraction:: we subtracted the continuum emission,
leaving us with only noise and the absorption and emission lines. The
absorption lines are negative and will be missed by detection methods
that look for a positive skewness(1) (like *note NoiseChisel::). So we
will focus on the detection and extraction of emission lines here.

The first step is to extract the voxels that contain emission signal.
To do that, we will be using *note NoiseChisel::. NoiseChisel and *note
Segment:: operate on 2D images or 3D cubes. But by default, they are
configured for 2D images (some parameters like tile size take a
different number of values based on the dimensionality). Therefore, to
do 3D detection, the first necessary step is to run NoiseChisel with the
default 3D configuration file.

To see where Gnuastro's programs are installed, you can run the
following command (the printed output is the default location when you
install Gnuastro from source, but if you used another installation
method or manually set a different location, you will see a different
output, just use that):

$ which astnoisechisel
/usr/local/bin/astnoisechisel

As you see, the compiled binary programs (like NoiseChisel) are
installed in the ‘bin/’ sub-directory of the install path (‘/usr/local’
in the example above, may be different on your system). The
configuration files are in the ‘etc/’ sub-directory of the install path
(here only showing NoiseChisel's configuration files):

$ ls /usr/local/etc/astnoisechisel*.conf
/usr/local/etc/astnoisechisel-3d.conf
/usr/local/etc/astnoisechisel.conf

We should therefore call NoiseChisel with the 3D configuration file like
below (please change ‘/usr/local’ to any directory that you find from
the ‘which’ command above):

$ astnoisechisel --config=/usr/local/etc/astnoisechisel-3d.conf \
no-continuum.fits --output=det.fits

But having to add this long ‘--config’ option is annoying and makes
the command hard to read! To simplify the calling of NoiseChisel in 3D,
let's first make a shell alias called ‘astnoisechisel-3d’ using the
‘alias’ command. Afterwards, we can just use the alias. Afterwards (in
the second command below), we are calling the alias, producing the same
output as above. Finally (with the last command), let's have a look at
NoiseChisel's output:

$ alias astnoisechisel-3d="astnoisechisel \
--config=/usr/local/etc/astnoisechisel-3d.conf"

$ astnoisechisel-3d no-continuum.fits --output=det.fits

$ astscript-fits-view det.fits

Similar to its 2D outputs, NoiseChisel's output contains four
extensions/HDUs (see *note NoiseChisel output::). For a multi-extension
file with 3D data, ‘astscript-fits-view’ shows each cube as a separate
DS9 "Frame". In this way, as you slide through the wavelengths, you see
the same slice in all the cubes. The third and fourth extensions are
the Sky and Sky standard deviation, which are not relevant here, so you
can close them. To do that, press on the "Frame" button (in the top row
of buttons), then press "delete" two times in the second row of buttons.

As a final preparation, manually set the scale of ‘INPUT-NO-SKY’ cube
to a fixed range so the changing flux/noise in each slice doesn't
interfere with visually comparing the data in the slices as you move
around:
1. Click on the ‘INPUT-NO-SKY’ cube, so it is selected.
2. Click on the "Scale" menu, then the "Scale Parameters".
3. For the "Low" value set -2 and for the "High" value set 5.
4. In the "Cube" window, slide between the slices to confirm that the
noise level is visually fixed.
5. Go back to the first slice for the next steps. Note that the first
and last couple of slices have much higher noise, don't worry about
those.

As you press the "Next" button in the first few slides, you will
notice that the ‘DETECTION’ cube is fully black: showing that nothing
has been detected. The first detection pops up in the 55th slice for
the galaxy on the top of this cube. As you press "Next" you will see
that the detection fades away and other detections pop up. Spend a few
minutes shifting between the different slices and comparing the detected
voxels with the emission lines in the continuum-subtracted cube (the
‘INPUT-NO-SKY’ extension).

Go ahead to slice 2815 and press "Next" a few times. You will notice
that the detections suddenly start covering the whole slice and until
slice 2859 where the detection map becomes normal (no extra
detections!). This is the effect of the sky lines we mentioned before
in *note Sky lines in optical IFUs::. The increased noise makes the
reduction very hard and as a result, a lot of artifacts appear. To
reduce the effect of sky lines, we can divide the cube by its standard
deviation (the square root of the variance or ‘STAT’ extension; see
*note Sky lines in optical IFUs::) and run NoiseChisel afterwards.

$ astarithmetic no-continuum.fits -h1 a370-crop.fits -hSTAT sqrt / \
--output=sn.fits

$ astnoisechisel-3d sn.fits --output=det.fits

$ astscript-fits-view det.fits

After the new detection map opens have another look at the specific
slices mentioned above (from slice 2851 to 2859). You will see that
there are no more detection maps that cover the whole field of view.
Scroll the slide counter across the whole cube, you will rarely see such
effects by Sky lines any more. But this is just a crude solution and
doesn't remove all sky line artifacts. For example go to slide 650 and
press "Next". You will see that the artifacts caused by this sky line
are so strong that the solution above wasn't successful. For these very
strong emission lines, we need to improve the reduction. But generally,
since the number of sky-line affected slices has significantly
decreased, we can go ahead.

---------- Footnotes ----------

(1) But if you want to detect the absorption lines, just multiply the
cube by $-1$ and repeat the same steps here (the noise is symmetric
around 0).

File: gnuastro.info, Node: 3D measurements and spectra, Next: Extracting a single spectrum and plotting it, Prev: 3D detection with NoiseChisel, Up: Detecting lines and extracting spectra in 3D data

2.5.5 3D measurements and spectra
---------------------------------

In the context of optical IFUs or radio IFUs in astronomy, a "Spectrum"
is defined as separate measurements on each 2D slice of the 3D cube.
Each 2D slice is defined by the first two FITS dimensions: the first
FITS dimension is the horizontal axis and the second is the vertical
axis. As with the tutorial on 2D image analysis (in *note Segmentation
and making a catalog::), let's run Segment to see how it works in 3D.
Like NoiseChisel above, to simplify the commands, let's make an alias
(*note 3D detection with NoiseChisel::):

$ alias astsegment-3d="astsegment \
--config=/usr/local/etc/astsegment-3d.conf"

$ astsegment-3d det.fits --output=seg.fits

$ astscript-fits-view seg.fits

You see that we now have 3D clumps and 3D objects. So we can go
ahead to do measurements. MakeCatalog can do single-valued measurements
(as in 2D) on 3D datasets also. For example, with the command below,
let's get the flux-weighted center (in the three dimensions) and sum of
pixel values. There isn't usually a standard name for the third WCS
dimension (unlike Ra/Dec). So in Gnuastro, we just call it ‘--w3’.
With the second command, we are having a look at the first 5 rows. Note
that we are not using ‘-Y’ with ‘asttable’ anymore because it the
wavelength column would only be shown as zero (since it is in meters!).

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum --output=cat.fits

$ asttable cat.fits -h1 -O --txtf64p=5 --head=5
# Column 1: OBJ_ID [counter ,i32,] Object identifier.
# Column 2: RA [deg ,f64,] Flux weighted center (WCS axis 1).
# Column 3: DEC [deg ,f64,] Flux weighted center (WCS axis 2).
# Column 4: AWAV [m ,f64,] Flux weighted center (WCS axis 3).
# Column 5: SUM [input-units,f32,] Sum of sky subtracted values.
1 3.99677e+01 -1.58660e+00 4.82994e-07 7.311189e+02
2 3.99660e+01 -1.58927e+00 4.86411e-07 7.872681e+03
3 3.99682e+01 -1.59141e+00 4.90609e-07 1.314548e+03
4 3.99677e+01 -1.58666e+00 4.90816e-07 7.798024e+02
5 3.99659e+01 -1.58930e+00 4.93657e-07 3.255210e+03

Besides the single-valued measurements above (that are shared with 2D
inputs), on 3D cubes, MakeCatalog can also do per-slice measurements.
The options for these measurements are formatted as ‘--*in-slice’. With
the command below, you can check their list:

$ astmkcatalog --help | grep in-slice
--area-in-slice [3D input] Number of labeled in each slice.
--area-other-in-slice [3D input] Area of other lab. in projected area.
--area-proj-in-slice [3D input] Num. voxels in '--sum-proj-in-slice'.
--sum-err-in-slice [3D input] Error in '--sum-in-slice'.
--sum-in-slice [3D input] Sum of values in each slice.
--sum-other-err-in-slice [3D input] Area in '--sum-other-in-slice'.
--sum-other-in-slice [3D input] Sum of other lab. in projected area.
--sum-proj-err-in-slice [3D input] Error of '--sum-proj-in-slice'.
--sum-proj-in-slice [3D input] Sum of projected area in each slice.

For every label and measurement, these options will give many values
in a vector column (see *note Vector columns::). Let's have a look by
asking for the sum of values and area of each label in each slice
associated to each label with the command below. There is just one
important point: in *note 3D detection with NoiseChisel::, we ran
NoiseChisel on the signal-to-noise image, not the continuum-subtracted
image! So the values to use for the measurement of each label should
come from the ‘no-continuum.fits’ file (not ‘seg.fits’).

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum \
--area-in-slice --sum-in-slice --output=cat.fits \
--valuesfile=no-continuum.fits

$ asttable -i cat.fits
--------
seg_cat.fits (hdu: 1)
------- ----- ---- -------
No.Name Units Type Comment
------- ----- ---- -------
1 OBJ_ID counter int32 Object identifier.
2 RA deg float64 Flux wht center (WCS 1).
3 DEC deg float64 Flux wht center (WCS 2).
4 AWAV m float64 Flux wht center (WCS 3).
5 SUM input-units float32 Sum of sky-subed values.
6 AREA-IN-SLICE counter int32(3681) Number of pix. in each slice.
7 SUM-IN-SLICE input-units float32(3681) Sum of values in each slice.
--------
Number of rows: 211
--------

You can see that the new ‘AREA-IN-SLICE’ and ‘SUM-IN-SLICE’ columns
have a ‘(3681)’ in their types. This shows that unlike the
single-valued columns before them, in these columns, each row has 3681
values (a "vector" column). If you are not already familiar with vector
columns, please take a few minutes to read *note Vector columns::.
Since a MUSE data cube has 3681 slices, this is effectively the spectrum
of each object.

Let's find the object that corresponds to the H-alpha emission of the
brightest galaxy (that we found in *note Viewing spectra and redshifted
lines::). That emission line was around 8565.93 Angstroms, so let's
look for the objects within $\pm5$ Angstroms of that value (between 8560
to 8570 Angstroms):

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -cobj_id,ra,dec -Y
198 39.965897 -1.589279

From the command above, we see that at this wavelength, there was
only one object. Let's extract its spectrum by asking for the
‘sum-in-slice’ column:

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 \
-carea-in-slice,sum-in-slice

If you look into the outputs, you will see that it is a single line!
It contains a long list of 0 values at the start and ‘nan’ values in the
end. If you scroll slowly, in the middle of each you will see some
non-zero and non-NaN numbers. To help interpret this more easily, let's
transpose these vector columns (so each value of the vector column
becomes a row in the output). We will use the ‘--transpose’ option of
Table for this (just note that since transposition changes the number of
rows, it can only be used when your table only has vector columns and
they all have the same number of elements (as in this case, for more):

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 \
-carea-in-slice,sum-in-slice --transpose

We now see the measurements on each slice printed in a separate line
(making it much more easier to visually read). However, without a
counter, it is very hard to interpret them. Let's pipe the output to a
new Table command and use column arithmetic's ‘counter’ operator for
displaying the slice number (see *note Size and position operators::).
Note that since we are piping the output, we also added ‘-O’ so the
column metadata are also passed to the new instance of Table:

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -O \
-carea-in-slice,sum-in-slice --transpose \
| asttable -c'arith $1 counter swap',2
...[[truncated]]...
3040 0 nan
3041 0 nan
3042 0 nan
3043 0 nan
3044 1 4.311140e-01
3045 18 3.936019e+00
3046 161 -5.800080e+00
3047 360 2.967184e+02
3048 625 1.912855e+03
3049 823 5.140487e+03
3050 945 7.174101e+03
3051 999 6.967604e+03
3052 1046 6.468591e+03
3053 1025 6.457354e+03
3054 996 6.599119e+03
3055 966 6.762280e+03
3056 873 5.014052e+03
3057 649 2.003334e+03
3058 335 3.167579e+02
3059 131 1.670975e+01
3060 25 -2.953789e+00
3061 0 nan
3062 0 nan
3063 0 nan
3064 0 nan
...[[truncated]]...

$ astscript-fits-view seg.fits

After DS9 opens with the last command above, go to slice 3044 (which
is the first non-NaN slice in the spectrum above). In the ‘OBJECTS’
extension of this slice, you see several non-zero pixels. The few
non-zero pixels on the bottom have a label of 197 and the single
non-zero pixel at a higher Y axis position has a label of 198 (which as
we saw above, was the label of the H-alpha emission of this galaxy).
The few 197 labeled pixels in this slice are the last voxels of the NII
emission that is just blue-ward of H-alpha.

The single pixel you see in slice 3044 is why you see a value of 1 in
the ‘AREA-IN-SLICE’ column. As you go to the next slices, if you count
the pixels, you will see they add up to the same number you see in that
column. The values in the ‘SUM-IN-SLICE’ are the sum of values in the
continuum-subtracted cube for those same voxels. You should now be able
to understand why the ‘--sum-in-slice’ column has NaN values in all
other slices: because this label doesn't exist in any other slice!
Also, within slices that contain label 198, this column only uses the
voxels that have the label. So as you see in the second column above,
the area that is used in each changes.

Therefore ‘--sum-in-slice’ or ‘area-in-slice’ are the raw 3D spectrum
of each 3D emission-line. This is a different concept from the
traditional "spectrum" where the same area is used over all the slices.
To get that you should use the ‘--sumprojinslice’ column of MakeCatalog.
All the ‘--*in-slice’ options that contain a ‘proj’ in their name are
measurements over the fixed "projection" of the 3D volume on the 2D
surface of each slice. To see the effect, let's also ask MakeCatalog to
measure this projected sum column:

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum \
--area-in-slice --sum-in-slice --sum-proj-in-slice \
--output=cat.fits --valuesfile=no-continuum.fits
$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -O \
-carea-in-slice,sum-in-slice,sum-proj-in-slice \
--transpose \
| asttable -c'arith $1 counter swap',2,3
...[[truncated]]...
3040 0 nan 8.686357e+02
3041 0 nan 4.384907e+02
3042 0 nan 4.994813e+00
3043 0 nan -1.595918e+02
3044 1 4.311140e-01 -2.793141e+02
3045 18 3.936019e+00 -3.251023e+02
3046 161 -5.800080e+00 -2.709914e+02
3047 360 2.967184e+02 1.049625e+02
3048 625 1.912855e+03 1.841315e+03
3049 823 5.140487e+03 5.108451e+03
3050 945 7.174101e+03 7.149740e+03
3051 999 6.967604e+03 6.913166e+03
3052 1046 6.468591e+03 6.442184e+03
3053 1025 6.457354e+03 6.393185e+03
3054 996 6.599119e+03 6.572642e+03
3055 966 6.762280e+03 6.716916e+03
3056 873 5.014052e+03 4.974084e+03
3057 649 2.003334e+03 1.870787e+03
3058 335 3.167579e+02 1.057906e+02
3059 131 1.670975e+01 -2.415764e+02
3060 25 -2.953789e+00 -3.534623e+02
3061 0 nan -3.745465e+02
3062 0 nan -2.532008e+02
3063 0 nan -2.372232e+02
3064 0 nan -2.153670e+02
...[[truncated]]...

As you see, in the new ‘SUM-PROJ-IN-SLICE’ column, we have a
measurement in each slice: including slices that do not have the label
of 198 at all. Also, the area used to measure this sum is the same in
all slices (similar to a classical spectrometer's output).

However, there is a big problem: have a look at the sums in slices
3040 and 3041: the values are increasing! This is because of the
emission in the NII line that also falls over the projected area of
H-alpha. This shows the power of IFUs as opposed to classical
spectrometers: we can distinguish between individual lines based on
spatial position and do measurements in 3D!

Finally, in case you want the spectrum with the continuum, you just
have to change the file given to ‘--valuesfile’:

$ astmkcatalog seg.fits --ids --ra --dec --w3 --sum \
--area-in-slice --sum-in-slice --sum-proj-in-slice \
--valuesfile=a370-crop.fits \
--output=cat-with-continuum.fits

File: gnuastro.info, Node: Extracting a single spectrum and plotting it, Next: Pseudo narrow-band images, Prev: 3D measurements and spectra, Up: Detecting lines and extracting spectra in 3D data

2.5.6 Extracting a single spectrum and plotting it
--------------------------------------------------

In *note 3D measurements and spectra:: we measured the spectra of all
the objects with the MUSE data cube of this demonstration tutorial.
Let's now write the resulting spectra for our object 198 into a file to
view our measured spectra in TOPCAT for a more visual inspection. But
we don't want slice numbers (which are specific to MUSE), we want the
horizontal axis to be in Angstroms. To do that, we can use the WCS
information:

‘CRPIX3’
The "Coordinate Reference PIXel" in the 3rd dimension (or slice
number of reference) Let's call this $s_r$.
‘CRVAL3’
The "Coordinate Reference VALue" in the 3rd dimension (the WCS
coordinate of the slice in ‘CRPIX3’. Let's call this $\lambda_r$
‘CDELT3’
The "Coordinate DELTa" in the 3rd dimension, or how much the WCS
changes with every slice. Let's call this $\delta$.

To find the $\lambda$ (wavelength) of any slice with number $s$, we can
simply use this equation:

$$\lambda=\lambda_r+\delta(s-s_r)$$

Let's extract these three values from the FITS WCS keywords as shell
variables to automatically do this within Table's column arithmetic.
Here we are using the technique that is described in *note Separate
shell variables for multiple outputs::.

$ eval $(astfits seg.fits --keyvalue=CRPIX3,CRVAL3,CDELT3 -q \
| xargs printf "sr=%s; lr=%s; d=%s;")

## Just for a check:
$ echo $sr
1.000000e+00
$ echo $lr
4.749679687500000e-07
$ echo $d
1.250000000000000e-10

Now that we have the necessary constants, we can simply convert the
equation above into *note Reverse polish notation:: and use column
arithmetic to convert the slice counter into wavelength in the command
of *note 3D measurements and spectra::.

$ asttable cat.fits --range=AWAV,8.560e-7,8.570e-7 -O \
-carea-in-slice,sum-in-slice,sum-proj-in-slice \
--transpose \
| asttable -c'arith $1 counter '$sr' - '$d' x '$lr' + f32 swap' \
-c2,3 --output=spectrum-obj-198.fits \
--colmetadata=1,WAVELENGTH,m,"Wavelength of slice." \
--colmetadata=2,"AREA-IN-SLICE",voxel,"No. of voxels."

$ astscript-fits-view spectrum-obj-198.fits

Once TOPCAT opens, take the following steps:

1. In the "Graphics" menu, select "Plane plot".
2. Change ‘AREA-IN-SLICE’ to ‘SUM-PROJ-IN-SLICE’.
3. Select the "Form" tab.
4. Click on the button with the large green "+" button and select "Add
line".
5. Un-select the "Mark" item that was originally selected.

Of course, the table in ‘spectrum-obj-198.fits’ can be plotted using any
other plotting tool you prefer to use in your scientific papers.

File: gnuastro.info, Node: Pseudo narrow-band images, Prev: Extracting a single spectrum and plotting it, Up: Detecting lines and extracting spectra in 3D data

2.5.7 Pseudo narrow-band images
-------------------------------

In *note Continuum subtraction:: we subtracted/separated the continuum
from the emission/absorption lines of our galaxy in the MUSE cube.
Let's visualize the morphology of the galaxy at some of the spectral
lines to see how it looks. To do this, we will create pseudo
narrow-band 2D images by collapsing the cube along the third dimension
within a certain wavelength range that is optimized for that flux.

Let's find the wavelength range that corresponds to H-alpha emission
we studied in *note Extracting a single spectrum and plotting it::.
Fortunately MakeCatalog can calculate the minimum and maximum position
of each label along each dimension like the command below. If you
always need these values, you can include these columns in the same
MakeCatalog with ‘--sum-proj-in-slice’. Here we are running it
separately to help you follow the discussion there.

$ astmkcatalog seg.fits --output=cat-ranges.fits \
--ids --min-x --max-x --min-y --max-y --min-z --max-z

Let's extract the minimum and maximum positions of this particular
object with the first command and with the second, we'll write them into
different shell variables. With the second command, we are writing
those six values into a single string in the format of Crop's *note Crop
section syntax::. For more on the ‘eval’-based shell trick we used
here, see *note Separate shell variables for multiple outputs::.
Finally, we are running Crop and viewing the cropped 3D cube.

$ asttable cat-ranges.fits --equal=OBJ_ID,198 \
-cMIN_X,MAX_X,MIN_Y,MAX_Y,MIN_Z,MAX_Z
56 101 11 61 3044 3060

$ eval $(asttable cat-ranges.fits --equal=OBJ_ID,198 \
-cMIN_X,MAX_X,MIN_Y,MAX_Y,MIN_Z,MAX_Z \
| xargs printf "section=%s:%s,%s:%s,%s:%s; ")

$ astcrop no-continuum.fits --mode=img --section=$section \
--output=crop-no-continuum.fits

$ astscript-fits-view crop-no-continuum.fits

Go through the slices and you will only see this particular region of
the full cube. We can now collapse the third dimension of this image
into a 2D pseudo-narrow band image with Arithmetic's *note
Dimensionality changing operators:::

$ astarithmetic crop-no-continuum.fits 3 collapse-sum \
--output=collapsed-all.fits

$ astscript-fits-view collapsed-all.fits

During the collapse, used all the pixels in each slice. This is not
good for the faint outskirts in the peak of the emission line: the noise
of the slices with less signal decreases the over-all signal-to-noise
ratio in the pseudo-narrow band image. So let's set all the pixels that
aren't labeled with this object as NaN, then collapse. To do that, we
first need to crop the ‘OBJECT’ cube in ‘seg.fits’. With the second
command, please have a look to confirm how the labels change as a
function of wavelength.

$ astcrop seg.fits -hOBJECTS --mode=img --section=$section \
--output=crop-obj.fits

$ astscript-fits-view crop-obj.fits

Let's use Arithmetic to first set all the pixels that are not equal
to 198 in ‘collapsed-obj.fits’ to be NaN in ‘crop-no-continuum.fits’.
With the second command, we are opening the two collapsed images
together:

$ astarithmetic crop-no-continuum.fits set-i \
crop-obj.fits set-o \
i o 198 ne nan where 3 collapse-sum \
--output=collapsed-obj.fits

$ astscript-fits-view collapsed-all.fits collapsed-obj.fits \
--ds9extra="-lock scalelimits yes -blink"

Let it blink a few times and focus on the outskirts: you will see
that the diffuse flux in the outskirts has indeed been preserved better
in the object-based collapsed narrow-band image. But this is a little
hard to appreciate in the 2D image. To see it better practice, let's
get the two radial profiles. We will approximately assume a position
angle of -80 and axis ratio of 0.7(1). With the final command below, we
are opening both radial profiles in TOPCAT to visualize them. We are
also undersampling the radial profile to have better signal-to-noise
ratio in the outer radii:

$ astscript-radial-profile collapsed-all.fits \
--position-angle=-80 --axis-ratio=0.7 \
--undersample=2 --output=collapsed-all-rad.fits

$ astscript-radial-profile collapsed-obj.fits \
--position-angle=-80 --axis-ratio=0.7 \
--undersample=2 --output=collapsed-obj-rad.fits

To view the difference, let's merge the two profiles (the ‘MEAN’
column) into one table and simply print the two profiles beside each
other. We will then pipe the resulting table containing both columns to
a second call to Gnuastro's Table and use column arithmetic to subtract
the two mean values and divide them by the optimized one (to get the
fractional difference):

$ asttable collapsed-all-rad.fits --catcolumns=MEAN -O \
--catcolumnfile=collapsed-obj-rad.fits \
| asttable -c1,2,3 -c'arith $3 $2 - $3 /' \
--colmetadata=2,MEAN-ALL \
--colmetadata=3,MEAN-OBJ \
--colmetadata=4,DIFF,frac,"Fractional diff." -YO
# Column 1: RADIUS [pix ,f32,] Radial distance
# Column 2: MEAN-ALL [input-units,f32,] Mean of sky subtracted values.
# Column 3: MEAN-OBJ [input-units,f32,] Mean of sky subtracted values.
# Column 4: DIFF [frac ,f32,] Fractional diff.
0.000 436.737 450.256 0.030
2.000 371.880 384.071 0.032
4.000 313.429 320.138 0.021
6.000 275.744 280.102 0.016
8.000 152.214 154.470 0.015
10.000 59.311 62.207 0.047
12.000 18.466 20.396 0.095
14.000 6.940 8.671 0.200
16.000 3.052 4.256 0.283
18.000 1.590 2.848 0.442
20.000 1.430 2.550 0.439
22.000 0.838 1.975 0.576

As you see, beyond a radius of 10, the last fractional difference column
becomes very large, showing that a lot of signal is missing in the
‘MEAN-ALL’ column. For a more visual comparison of the two profiles,
you can use the command below to open both tables in TOPCAT:

$ astscript-fits-view collapsed-all-rad.fits \
collapsed-obj-rad.fits

Once TOPCAT has opened take the following steps:
1. Select ‘collapsed-all-rad.fits’
2. In the "Graphics" menu, select "Plane Plot".
3. Click on the "Axes" side-bar (by default, at the bottom half of the
window), and click on "Y Log" to view the vertical axis in
logarithmic scale.
4. In the "Layers" menu, select "Add Position Control". You will see
that at the bottom half, a new scatter plot information is
displayed.
5. Click on the scroll-down menu in front of "Table" and select ‘2:
collapsed-obj-rad.fits’. Afterwards, you will see the optimized
pseudo-narrow-band image radial profile as blue points.

---------- Footnotes ----------

(1) To derive the axis ratio and position angle automatically, you
can take the following steps. Note that we are not using NoiseChisel
because this crop has been intentionally selected to contain signal, so
there is no raw noise inside of it.
$ aststatistics collapsed-all.fits --sky --tilesize=5,5
$ astarithmetic collapsed-all.fits -h1 collapsed-all_sky.fits -hSKY_STD / 5 gt \
-ocollapsed-lab.fits
$ astmkcatalog collapsed-lab.fits -h1 --valuesfile=collapsed-all.fits \
--position-angle --axis-ratio
$ asttable collapsed-all_arith_cat.fits -Y
-78.817 0.694

File: gnuastro.info, Node: Color images with full dynamic range, Next: Zero point of an image, Prev: Detecting lines and extracting spectra in 3D data, Up: Tutorials

2.6 Color images with full dynamic range
========================================

Color images are fundamental tools to visualize astronomical datasets,
allowing to visualize valuable physical information within them. A
color image is a composite representation derived from different
channels. Each channel usually corresponding to different filters (each
showing wavelength intervals of the object's spectrum). In general,
most common color image formats (like JPEG, PNG or PDF) are defined from
a combination of Red-Green-Blue (RGB) channels (to cover the optical
range with normal cameras). These three filters are hard-wired in your
monitor and in most normal camera (for example smartphone or DSLR)
pixels. For more on the concept and usage of colors, see *note Color::
and *note Colormaps for single-channel pixels::.

However, normal images (that you take with your smartphone during the
day for example) have a very limited dynamic range (difference between
brightest and fainest part of an image). For example in an image you
take from a farm, the brightness pixel (the sky) cannot be more than 255
times the faintest/darkest shadow in the image (because normal cameras
produce unsigned 8 bit integers; containing $2^8=256$ levels; see *note
Numeric data types::).

However, astronomical sources span a much wider dynamic range such
that their central parts can be tens of millions of times brighter than
their larger outer regions. Our astronomical images in the FITS format
are therefore usually 32-bit floating points to preserve this
information. Therefore a simple linear scaling of 32-bit astronomical
data to the 8-bit range will put most of the pixels on the darkest level
and barely show anything! This presents a major challenge in
visualizing our astronomical images on a monitor, in print or for a
projector when showing slides.

In this tutorial, we review how to prepare your images and create
informative RGB images for your PDF reports. We start with aligning the
images to the same pixel grid (which is usually necessary!) and using
the low-level engine (Gnuastro's *note ConvertType:: program) directly
to create an RGB image. Afterwards, we will use a higher-level
installed script (*note Color images with gray faint regions::). This
is a high-level wrapper over ConvertType that does some pre-processing
and stretches the pixel values to enhance their 8-bit representation
before calling ConvertType.

* Menu:

* Color channels in same pixel grid:: Warping all inputs to the same pixel grid.
* Color image using linear transformation:: A linear color mapping won't show much!
* Color for bright regions and grayscale for faint:: Show the full dynamic range.
* Manually setting color-black-gray regions:: Physically motivated regions.
* Weights contrast markers and other customizations:: Nice ways to enhance visual appearance.

File: gnuastro.info, Node: Color channels in same pixel grid, Next: Color image using linear transformation, Prev: Color images with full dynamic range, Up: Color images with full dynamic range

2.6.1 Color channels in same pixel grid
---------------------------------------

In order to use different images as color channels, it is important that
the images be properly aligned and on the same pixel grid. When your
inputs are high-level products of the same survey, this is usually the
case. However, in many other situations the images you plan to use as
different color channels lie on different sky positions, even if they
may have the same number of pixels. In this section we will show how to
solve this problem.

For an example dataset, let's use the same SDSS field that we used in
*note Detecting large extended targets::: the field covering the outer
parts of the M51 group. With the commands below, we'll make an ‘inputs’
directory and download and prepare the three g, r and i band images of
SDSS over the same field there:

$ mkdir in
$ sdssurl=https://dr12.sdss.org/sas/dr12/boss/photoObj/frames
$ for f in g r i; do \
wget $sdssurl/301/3716/6/frame-$f-003716-6-0117.fits.bz2 \
-O$f.fits.bz2; \
bunzip2 $f.fits.bz2; \
astfits $f.fits --copy=0 -oin/$f-sdss.fits; \
rm $f.fits; \
done

Let's have a look at the three three images with the first command,
and get their number of pixels with the second:

## Open the images locked by image coordinates
$ astscript-fits-view in/*-sdss.fits

## Check the number of pixels along each axis of all images.
$ astfits in/*-sdss.fits --keyvalue=NAXIS1,NAXIS2
in/g-sdss.fits 2048 1489
in/i-sdss.fits 2048 1489
in/r-sdss.fits 2048 1489

From the first command, the images look like they cover the same
astronomical object (M51) in the same region of the sky, and with the
second, we see that they have the same number of pixels. But this
general visual inspection does not guarantee that the astronomical
objects within the pixel grid cover exactly the same positions (within a
pixel!) on the sky. Let's open the images again, but this time asking
DS9 to only show one at a time, and to "blink" between them:

$ astscript-fits-view in/*-sdss.fits \
--ds9extra="-single -zoom to fit -blink"

If you pay attention, you will see that the objects within each image
are at slightly different locations. If you don't immediately see it,
try zooming in to any star within the image and let DS9 continue
blinking. You will see that the star jumps a few pixels between each
blink.

In essence, the images are not aligned on the same pixel grid,
therefore, the same source does not share identical image coordinates
across these three images. As a consequence, it is necessary to align
the images before making the color image, otherwise this misalignment
will generate multiply-peaked point-sources (stars and centers of
galaxies) and artificial color gradients in the more diffuse parts. To
align the images to the same pixel grid, we will employ Gnuastro's *note
Warp:: program. In particular, its features to *note Align pixels with
WCS considering distortions::.

Let's take the middle (r band) filter as the reference to define our
grid. With the first command after building the ‘aligned/’ directory,
let's align the r filter to the celestial coordinates (so the M51
group's position angle doesn't depend on the orientation of the
telescope when it took this image). With for the other two filters, we
will use Warp's ‘--gridfile’ option to ensure that ensure that their
pixel grid and WCS exactly match the r band image, but the pixel values
come from the other two filters. Finally, in the last command, we'll
visualize the three aligned images.

## Put all three channels in the same pixel grid.
$ mkdir aligned
$ astwarp in/r-sdss.fits --output=aligned/r-sdss.fits
$ astwarp in/g-sdss.fits --output=aligned/g-sdss.fits \
--gridfile=aligned/r-sdss.fits
$ astwarp in/i-sdss.fits --output=aligned/i-sdss.fits \
--gridfile=aligned/r-sdss.fits

## Open the images locked by image coordinates
$ astscript-fits-view aligned/*-sdss.fits \
--ds9extra="-single -zoom to fit -blink"

As the images blink between each other, zoom in to some of the
smaller stars and you will see that they no longer jump from one blink
to the next. These images are now precisely pixel-aligned. We are now
equipped with the essential data to proceed with the color image
generation in *note Color image using linear transformation::.

File: gnuastro.info, Node: Color image using linear transformation, Next: Color for bright regions and grayscale for faint, Prev: Color channels in same pixel grid, Up: Color images with full dynamic range

2.6.2 Color image using linear transformation
---------------------------------------------

Previously (in *note Color channels in same pixel grid::), we downloaded
three SDSS filters of M51 and described how you can put them all in the
same pixel grid. In this section, we will explore the raw and low-level
process of generating color images using the input images (without
modifying the pixel value distributions). We will use Gnuastro's
ConvertType program (with executable name ‘astconvertt’).

Let's create our first color image using the aligned SDSS images
mentioned in the previous section. The order in which you provide the
images matters, so ensure that you sort the filters from redder to bluer
(iSDSS and gSDSS are respectively the reddest and bluest of the three
filters used here).

$ astconvertt aligned/i-sdss.fits aligned/r-sdss.fits \
aligned/g-sdss.fits -g1 --output=m51.pdf

*Other color formats:* In the example above, we are using PDF because
this is usually the best format to later also insert marks that are
commonly necessary in scientific publications (see *note Marking objects
for publication::. But you can also generate JPEG and TIFF outputs
simply by using a different suffix for your output file (for example
‘--output=m51.jpg’ or ‘--output=m51.tiff’).

Open the image with your PDF viewer and have a look. Do you see
something? Initially, it appears predominantly black. However, upon
closer inspection, you will discern very tiny points where some color is
visible. These points correspond to the brightest part of the brightest
sources in this field! The reason you saw much more structure when
looking at the image in DS9 previously in *note Color channels in same
pixel grid:: was that ‘astscript-fits-view’ used DS9's ‘-zscale’ option
to scale the values in a non-linear way! Let's have another look at the
images with the linear ‘minmax’ scaling of DS9:

$ astscript-fits-view aligned/*-sdss.fits \
--ds9extra="-scale minmax -lock scalelimits"

You see that it looks very similar to the PDF we generated above:
almost fully black! This phenomenon exemplifies the challenge discussed
at the start of this tutorial in *note Color images with full dynamic
range::). Given the vast number of pixels close to the sky background
level compared to the relatively few very bright pixels, visualizing the
entire dynamic range simultaneously is tricky.

To address this challenge, the low-level ConvertType program allows
you to selectively choose the pixel value ranges to be displayed in the
color image. This can be accomplished using the ‘--fluxlow’ and
‘--fluxhigh’ options of ConvertType. Pixel values below ‘--fluxlow’ are
mapped to the minimum value (displayed as black in the default
colormap), and pixel values above ‘--fluxhigh’ are mapped to the maximum
value (displayed as white)) The choice of these values depends on the
pixel value distribution of the images.

But before that, we have to account for an important differences
between the filters: the brightness of the background also has different
values in different filters (the sky has colors!) So before making more
progress, generally, first you have to subtract the sky from all three
images you want to feed to the color channels. In a previous tutorial
(*note Detecting large extended targets::) we used these same images as
a basis to show how you can do perfect sky subtraction in the presence
of large extended objects like M51. Here we are just doing a
visualization and bringing pixels to 8-bit, so we don't need that level
of precision reached there (we won't be doing photometry!). Therefore,
let's just keep the ‘--tilesize=100,100’ of NoiseChisel.

$ mkdir no-sky
$ for f in i r g; do \
astnoisechisel aligned/$f-sdss.fits --tilesize=100,100 \
--output=no-sky/$f-sdss.fits; \
done

*Accounting for zero points:* An important step that we have not
implemented in this section is to unify the zero points of the three
filters. In the case of SDSS (and some other surveys), the images have
already been brought to the same zero point, but that is not generally
the case. So before subtracting sky (and estimating the standard
deviation) you should also unify the zero points of your images (for
example through Arithmetic's ‘counts-to-nanomaggy’ or ‘counts-to-jy’
described in *note Unit conversion operators::). If you don't already
have the zero point of your images, see the dedicated tutorial: *note
Zero point of an image::.

Now that we know the noise fluctuates around zero in all three
images, we can start to define the values for the ‘--fluxlow’ and
‘--fluxhigh’. But the sky standard deviation comes from the sky
brightness in different filters and is therefore different! Let's have
a look by taking the median value of the ‘SKY_STD’ extension of
NoiseChisel's output:

$ aststatistics no-sky/i-sdss.fits -hSKY_STD --median
2.748338e-02

$ aststatistics no-sky/r-sdss.fits -hSKY_STD --median
1.678463e-02

$ aststatistics no-sky/g-sdss.fits -hSKY_STD --median
9.687680e-03

You see that the sky standard deviation of the reddest filter (i) is
almost three times the bluest filter (g)! This is usually the case in
any scenario (redder emission usually requires much less energy and gets
absorbed less, so the background is usually brighter in the reddest
filters). As a result, we should define our limits based on the noise
of the reddest filter. Let's set the minimum flux to 0 and the maximum
flux to ~50 times the noise of the i-band image ($0.027\times50=1.35$).

$ astconvertt no-sky/i-sdss.fits no-sky/r-sdss.fits no-sky/g-sdss.fits \
-g1 --fluxlow=0.0 --fluxhigh=1.35 --output=m51.pdf

After opening the new color image, you will observe that a spiral arm
of M51 and M51B (or NGC5195, which is interacting with M51), become
visible. However, the majority of the image remains black. Feel free
to experiment with different values for ‘--fluxhigh’ to set the maximum
value closer to the noise-level and see the more diffuse structures.
For instance, try with ‘--fluxhigh=0.27’ the brightest pixels will have
a signal-to-noise ratio of 10, or even ‘--fluxhigh=0.135’ for a
signal-to-noise ratio of 5. But you will notice that, the brighter
areas of the galaxy become "saturated": you don't see the structure of
brighter parts of the galaxy any more. As you bring down the maximum
threshold, the saturated areas also increase in size: loosing some
useful information on the bright side!

Let's go to the extreme and decrease the threshold to close the
noise-level (for example ‘--fluxhigh=0.027’ to have a signal-to-noise
ratio of 1)! You will see that the noise now becomes colored! You
generally don't want this because the difference in filter values of one
pixel are only physically meaningful when they have a high
signal-to-noise ratio. For lower signal-to-noise ratios, we should
avoid color.

Ideally, we want to see both the brighter parts of the central
galaxy, as well as the fainter diffuse parts together! But with the
simple linear transformation here, that is not possible! You need some
pre-processing (before calling ConvertType) to scale the images. For
example, you can experiment with taking the logarithm or the square root
of the images (using *note Arithmetic::) before creating the color
image.

These non-linear functions transform pixel values, mapping them to a
new range. After applying such transformations, you can use the
transformed images as inputs to ‘astconvertt’ to generate color images
(similar to how we subtracted the sky; which is a linear operation). In
addition to that, it is possible to use a different color schema for
showing the different brightness ranges as it is explained in the next
section. In the next section (*note Color for bright regions and
grayscale for faint::), we'll review one high-level installed script
which will simplify all these pre-processings and help you produce
images with more information in them.

File: gnuastro.info, Node: Color for bright regions and grayscale for faint, Next: Manually setting color-black-gray regions, Prev: Color image using linear transformation, Up: Color images with full dynamic range

2.6.3 Color for bright regions and grayscale for faint
------------------------------------------------------

In the previous sections we aligned three SDSS images of M51 group *note
Color channels in same pixel grid::, and created a linearly-scaled color
image (only using ‘astconvertt’ program) in *note Color image using
linear transformation::. But we saw that showing the brighter and
fainter parts of the galaxy in a single image is impossible in the
linear scale! In this section, we will use Gnuastro's
‘astscript-color-faint-gray’ installed script to address this problem
and create images which visualize a major fraction of the contents of
our astronomical data.

This script aims to solve the problems mentioned in the previous
section. See Infante-Sainz et al. 2024
(https://arxiv.org/abs/2401.03814), which first introduced this script,
for examples of the final images we will be producing in this tutorial.
This script uses a non-linear transformation to modify the bright input
values before combining them to produce the color image. Furthermore,
for the faint regions of the image, it will use grayscale and avoid
color over all (as we saw, colored noised is not too nice to look at!).
The faint regions are also inverted: so the brightest pixel in the faint
(black-and-white or grayscale) region is black and the faintest pixels
will be white. Black therefore creates a smooth transition from the
colored bright pixels: the faintest colored pixel is also black. Since
the background is white and the diffuse parts are black, the final
product will also show nice in print or show on a projector (the
background is not black, but white!).

The SDSS image we used in the previous sections doesn't show the full
glory of the M51 group! Therefore, in this section, we will use the
wider images from the J-PLUS survey (https://www.j-plus.es).
Fortunately J-PLUS includes the SDSS filters, so we can use the same
iSDSS, rSDSSS, and gSDSS filters of J-PLUS. As a consequence, similar to
the previous section, the R, G, and B channels are respectively mapped
to the iSDSS, rSDSS and gSDSS filters of J-PLUS.

The J-PLUS identification numbers for the images containing the M51
galaxy group are in these three filters are respectively: 92797, 92801,
92803. The J-PLUS images are already sky subtracted and aligned into
the same pixel grid (so we will not need the ‘astwarp’ and
‘astnoisechisel’ steps before). However, zero point magnitudes of the
J-PLUS images are different: 23.43, 23.74, 23.74. Also, the field of
view of the J-PLUS Camera is very large and we only need a small region
to see the M51 galaxy group. Therefore, we will crop the regions around
the M51 group with a width of 0.35 degree wide (or 21 arcmin) and put
the crops in the same ‘aligned/’ directory we made before (which
contains the inputs to the colored images). With all the above
information, let's download, crop, and have a look at the images to
check that everything is fine. Finally, let's run
‘astscript-color-faint-gray’ on the three cropped images.

## Download
$ url=https://archive.cefca.es/catalogues/vo/siap/jplus-dr3/get_fits?id=
$ wget "$url"92797 -Oin/i-jplus.fits.fz
$ wget "$url"92801 -Oin/r-jplus.fits.fz
$ wget "$url"92803 -Oin/g-jplus.fits.fz

## Crop
$ widthdeg=0.35
$ ra=202.4741207
$ dec=47.2171879
$ for f in i r g; do \
astcrop in/$f-jplus.fits.fz --center=$ra,$dec \
--width=$widthdeg --output=aligned/$f-jplus.fits; \
done

## Visual inspection of the images used for the color image
$ astscript-fits-view aligned/*-jplus.fits

## Create colored image.
$ R=aligned/i-jplus.fits
$ G=aligned/r-jplus.fits
$ B=aligned/g-jplus.fits
$ astscript-color-faint-gray $R $G $B -g1 --output=m51.pdf

After opening the PDF, you will notice that it is a color image with
a gray background, making the M51 group and background galaxies visible
together. However, the images does not look nice and there is
significant room for improvement! You will notice that at the end of
its operation, the script printed some numerical values for four options
in a table, to show automatically estimated parameter values. To
enhance the output, let's go through and explain these step by step.

The first important point to take into account is the photometric
calibration. If the images are photometrically calibrated, then it is
necessary to use the calibration to put the images in the same physical
units and create "real" colors. The script is able to do it through the
zero point magnitudes with the option ‘--zeropoint’ (or ‘-z’). With
this option, the images are internally transformed to have the same
pixel units and then create the color image. Since the magnitude zero
points are 23.43, 23.74, 23.74 for the i, r, and g images, let's use
them in the definition

$ astscript-color-faint-gray $R $G $B -g1 --output=m51.pdf \
-z23.43 -z23.74 -z23.74

Open the image and have a look. This image does not differ too much
from the one generated by default (not using the zero point magnitudes).
This is because the zero point values used here are similar for the
three images. But in other datasets the calibration could make a big
difference!

Let's consider another vital parameter: the minimum value to be
displayed (‘--minimum’ or ‘-m’). Pixel values below this number will
not be shown on the color image. In general, if the sky background has
been subtracted (see *note Color image using linear transformation::),
you can use the same value (0) for all three. However, it is possible
to consider different minimum values for the inputs (in this case use as
many ‘-m’ as input images). In this particular case, a minimum value of
zero for all images is suitable. To keep the command simple, we'll add
the zero point, minimum and HDU of each image in the variable that also
had its filename.

$ R="aligned/i-jplus.fits -h1 --zeropoint=23.43 --minimum=0.0"
$ G="aligned/r-jplus.fits -h1 --zeropoint=23.74 --minimum=0.0"
$ B="aligned/g-jplus.fits -h1 --zeropoint=23.74 --minimum=0.0"
$ astscript-color-faint-gray $R $G $B --output=m51.pdf

In contrast to the previous image, the new PDF (with a minimum value
of zero) exhibits a better background visualization because it is
avoiding negative pixels to be included in the scaling (they are white).

Now let's review briefly how the script modifies the pixel value
distribution in order to show the entire dynamical range in an
appropriate way. The script combines the three images into a single one
by using a the mean operator, as a consequence, the combined image is
the average of the three R, G, and B images. This averaged image is
used for performing the asinh transformation of Lupton et al. 2004
(https://ui.adsabs.harvard.edu/abs/2004PASP..116..133L) that is
controlled by two parameters: ‘--qbright’ ($q$) and ‘--stretch’ ($s$).

The asinh transformation consists in transforming the combined image
($I$) according to the expression: $f(I) =
asinh(q\times{}s\times{}I)/q$. When $q\rightarrow0$, the expression
becomes linear with a slope of the "stretch" ($s$) parameter: $f(I) =
s\times{}I$. In practice, we can use this characteristic to first set a
low value for ‘--qbright’ and see the brighter parts in color, while
changing the parameter ‘--stretch’ to show linearly the fainter regions
(outskirts of the galaxies for example. The image obtained previously
was computed by the default parameters (‘--qthresh=1.0’ and
‘--stretch=1.0’). So, let's set a lower value for ‘--qbright’ and check
the result.

$ astscript-color-faint-gray $R $G $B --output=m51-qlow.pdf \
--qbright=0.01

Comparing ‘m51.pdf’ and ‘m51-qlow.pdf’, you will see that a large
area of the previously colored colored pixels have become black. Only
the very brightest pixels (core of the galaxies and stars) are shown in
color. Now, let's bring out the fainter regions around the brightest
pixels linearly by increasing ‘--stretch’. This allows you to reveal
fainter regions, such as outer parts of galaxies, spiral arms, stellar
streams, and similar structures. Please, try different values to see
the effect of changing this parameter. Here, we will use the value of
‘--stretch=100’.

$ astscript-color-faint-gray $R $G $B --output=m51-qlow-shigh.pdf \
--qbright=0.01 --stretch=100

Do you see how the spiral arms and the outskirts of the galaxies have
become visible as ‘--stretch’ is increased? After some trials, you will
have the necessary feeling to see how it works. Please, play with these
two parameters until you obtain the desired results. Depending on the
absolute pixel values of the input images and the photometric
calibration, these two parameters will be different. So, when using
this script on your own data, take your time to study and analyze which
parameters are good for showing the entire dynamical range. For this
tutorial, we will keep it simple and use the previous parameters. Let's
define a new variable to keep the parameters already discussed so we
have short command-line examples.

$ params="--qbright=0.01 --stretch=100"
$ astscript-color-faint-gray $R $G $B $params --output=m51.pdf
$ rm m51-qlow.pdf m51-qlow-shigh.pdf

Having a separate color-map for the fainter parts is generally a good
thing, but for some reason you may not want it! To disable this
feature, you can use the ‘--coloronly’ option:

$ astscript-color-faint-gray $R $G $B $params --coloronly \
--output=m51-coloronly.pdf

Open the image and note that now the coloring has gone all the way
into the noise (producing a black background). In contrast with the
gray background images before, the fainter/smaller stars/galaxies and
the low surface brightness features are not visible anymore! These
regions show the interaction of two galaxies; as well as all the other
background galaxies and foreground stars. These structures were
entirely hidden in the "only-color" images. Consequently, the gray
background color scheme is particularly useful for visualizing the most
features of your data and you will rarely need to use the ‘--coloronly’
option. We will therefore not use this option anymore in this tutorial;
and let's clean the temporary file made before:

$ rm m51-coloronly.pdf

Now that we have the basic parameters are set, let's consider other
parameters that allow to fine tune the three ranges of values: color for
the brightest pixel values, black for intermediate pixel values, and
gray for the faintest pixel values:
• ‘--colorval’ defines the boundary between the color and black
regions (the lowest pixel value that is colored).
• ‘--grayval’ defines the boundary between the black and gray regions
(the highest gray value).
Looking at the last lines that the script prints, we see that the
default value estimated for ‘--colorval’ and ‘--grayval’ are roughly
1.4. What do they mean? To answer this question it is necessary to
have a look at the image that is used to separate those different
regions. By default, this image is computed internally by the script
and removed at the end. To have a look at it, you need to use the
option ‘--keeptmp’ to keep the temporary files. Let's put the temporary
files into the ‘tmp’ directory with the option ‘--tmpdir=tmp --keeptmp’.
The first will use the name ‘tmp’ for the temporary directory and with
the second, we ask the script to not delete (keep) it after all
operations are done.

$ astscript-color-faint-gray $R $G $B $params --output=m51.pdf \
--tmpdir=tmp --keeptmp

The image that defines the thresholds is
‘./tmp/colorgray_threshold.fits’. By default, this image is the
asinh-transformed image with the pixel values between 0 (faint) and 100
(bright). If you obtain the statistics of this image, you will see that
the median value is exactly the value that the script is giving as the
‘--colorval’.

$ aststatistics ./tmp/colorgray_threshold.fits

In other words, all pixels between 100 and this value (1.4) on the
threshold image will be shown in color. To see its effect, let's
increase this parameter to ‘--colorval=25’. By doing this, we expect
that only bright pixels (those between 100 and 25 in the threshold
image) will be in color.

$ astscript-color-faint-gray $R $G $B $params --colorval=25 \
--output=m51-colorval.pdf

Open ‘m51-colorval.pdf’ and check that it is true! Only the central
part of the objects (very bright pixels, those between 100 and 25 on the
threshold image) are shown in color. Fainter pixels (below 25 on the
threshold image) are shown in black and gray. However, in many
situations it is good to be able to show the outskirts of galaxies and
low surface brightness features in pure black, while showing the
background in gray. To do that, we can use another threshold that
separates the black and gray pixels: ‘--grayval’.

Similar to ‘--colorval’, the ‘--grayval’ option defines the
separation between the pure black and the gray pixels from the threshold
image. For example, by setting ‘--grayval=5’, those pixels below 5 in
the threshold image will be shown in gray, brighter pixels will be shown
in black until the value 25. Pixels brighter than 25 are shown in
color.

$ astscript-color-faint-gray $R $G $B $params --output=m51-check.pdf \
--colorval=25 --grayval=5

Open the image and check that the regions shown in color are smaller
(as before), and that now there is a region around those color pixels
that are in pure black. After the black pixels toward the fainter ones,
they are shown in gray. As explained above, in the gray region, the
brightest are black and the faintest are white. It is recommended to
experiment with different values around the estimated one to have a
feeling on how it changes the image. To have even better idea of those
regions, please run the following example to keep temporary files and
check the labeled image it has produced:

$ astscript-color-faint-gray $R $G $B $params --output=m51-check.pdf \
--colorval=25 --grayval=5 \
--tmpdir=tmp --keeptmp

$ astscript-fits-view tmp/total_mask-2color-1black-0gray.fits

In this segmentation image, pixels equal to 2 will be shown in color,
pixels equal to 1 will be shown as pure black, and pixels equal to zero
are shown in gray. By default, the script sets the same value for both
thresholds. That means that there is not many pure black pixels. By
adjusting the ‘--colorval’ and ‘--grayval’ parameters, you can obtain an
optimal result to show the bright and faint parts of your data within
one printable image. The values used here are somewhat extreme to
illustrate the logic of the procedure, but we encourage you to
experiment with values close to the estimated by default in order to
have a smooth transition between the three regions (color, black, and
gray). The script can provide additional information about the pixel
value distributions used to estimate the parameters by using the
‘--checkparams’ option.

To conclude this section of the tutorial, let's clean up the
temporary test files:

$ rm m51-check.pdf m51-colorval.pdf

File: gnuastro.info, Node: Manually setting color-black-gray regions, Next: Weights contrast markers and other customizations, Prev: Color for bright regions and grayscale for faint, Up: Color images with full dynamic range

2.6.4 Manually setting color-black-gray regions
-----------------------------------------------

In *note Color for bright regions and grayscale for faint::, we created
a non-linear colored image. We used the ‘--colorval’ and ‘--grayval’
options to specify which regions to show in gray (faintest values),
black (intermediate values) and color (brightest values). We also saw
that the script uses a labeled image with three possible values for each
pixel to identify how that pixel should be colored.

A useful feature of this script is the possibility of providing this
labeled image as an input directly. This expands the possibilities of
generating color images in a more quantitative way. In this section,
we'll use this feature to use a more physically motivated criteria to
select these three regions (the surface brightness in the reddest band).

First, let's generate a surface brightness image from the R channel.
That is, the value of each pixel will be in the units of surface
brightness (mag/arcsec$^2$). To do that, we need obtain the pixel area
in arcsec and use the zero point value of the image. Then, the
‘counts-to-sb’ operator of ‘astarithmetic’ is used. For more on the
conversion of NaN surface brightness values and the value to ‘R_sbl’
(which is roughly the surface brightness limit of this image), see *note
FITS images in a publication::.

$ sb_sbl=26
$ sb_zp=23.43
$ sb_img=aligned/i-jplus.fits
$ pixarea=$(astfits $sb_img --pixelareaarcsec2 --quiet)

# Compute the SB image (set NaNs to SB of 26!)
$ astarithmetic $sb_img $sb_zp $pixarea counts-to-sb set-sb \
sb sb isblank sb $sb_sbl gt or $sb_sbl where \
--output=sb.fits

# Have a look at the image
$ astscript-fits-view sb.fits --ds9scale=minmax \
--ds9extra="-invert"

Remember that because ‘sb.fits’ is a surface brightness image where
lower values are brighter and higher values are fainter. Let's build
the labeled image that defines the regions (‘regions.fits’) step-by-step
with the following criteria in surface brightness (SB)

$\rm{SB}<23$
These are the brightest pixels, we want these in color. In the
regions labeled image, these should get a value of 2.
$23<\rm{SB}<25$
These are the intermediate pixel values, to see the fainter parts
better, we want these in pure black (no change in color in this
range). In the regions labeled image, these should get a value of
1.
$\rm{SB}>25$
These are the faintest pixel values, we want these in a gray color
map (pixels with an SB of 25 will be black and as they become
fainter, they will become lighter shades of gray). In the regions
labeled image, these should get a value of 0.

# SB thresholds (low and high)
$ sb_faint=25
$ sb_bright=23

# Select the three ranges of pixels.
$ astarithmetic sb.fits set-sb \
sb $sb_bright lt set-color \
sb $sb_bright ge sb $sb_faint lt and set-black \
color 2 u8 x black + \
--output=regions.fits

# Check the images
$ astscript-fits-view regions.fits

We can now use this labeled image with the ‘--regions’ option for
obtaining the final image with the desired regions (the ‘R’, ‘G’, ‘B’
and ‘params’ shell variables were set previously in *note Color for
bright regions and grayscale for faint::):

$ astscript-color-faint-gray $R $G $B $params --output=m51-sb.pdf \
--regions=regions.fits

Open ‘m51-sb.pdf’ and have a look. Do you see how the different
regions (SB intervals) have been colored differently? They come from
the SB levels we defined, and because it is using absolute thresholds in
physical units of surface brightness, the visualization is not only a
nice looking color image, but can be used in scientific analysis.

This is really interesting because now it is possible to use color
images for detecting low surface brightness features at the same time
they provide quantitative measurements. Of course, here we have defined
this region label image just using two surface brightness values, but it
is possible to define any other labeled region image that you may need
for your particular purpose.

File: gnuastro.info, Node: Weights contrast markers and other customizations, Prev: Manually setting color-black-gray regions, Up: Color images with full dynamic range

2.6.5 Weights, contrast, markers and other customizations
---------------------------------------------------------

Previously (in *note Manually setting color-black-gray regions::) we
used an absolute (in units of surface brightness) thresholding for
selecting which regions to show by color, black and gray. To keep the
previous configurations and avoid long commands, let's add the previous
options to the ‘params’ shell variable. To help in readability, we will
repeat the other shell variables from previous sections also:

$ R="aligned/i-jplus.fits -h1 --zeropoint=23.43 --minimum=0.0"
$ G="aligned/r-jplus.fits -h1 --zeropoint=23.74 --minimum=0.0"
$ B="aligned/g-jplus.fits -h1 --zeropoint=23.74 --minimum=0.0"
$ params="--regions=regions.fits --qbright=0.01 --stretch=100"
$ astscript-color-faint-gray $R $G $B $params --output=m51.pdf

To modify the color balance of the output image, you can weigh the
three channels differently with the ‘--weight’ or ‘-w’ option. For
example, by using ‘-w1 -w1 -w2’, you give two times more weight to the
blue channel than to the red and green channels:

$ astscript-color-faint-gray $R $G $B $params -w1 -w1 -w2 \
--output=m51-weighted.pdf

The colored pixels of the output are much bluer now and the
distinction between the two merging galaxies is more clear. However,
keep in mind that altering the different filters can lead to incorrect
subsequent analyses by the readers/viewers of this work (for example
they will falsely think that the galaxy is blue, and not red!). If the
reduction and photometric calibration are correct, and the images
represent what you consider as the red, green, and blue channels, then
the output color image should be suitable without weights.

In certain situations, the combination of channels may not have a
traditional color interpretation. For instance, combining an X-ray
channel with an optical filter and a far-infrared image can complicate
the interpretation in terms of human understanding of color. But the
physical interpretation remains valid as the different channels (colors
in the output) represent different physical phenomena of astronomical
sources. Another easier example is the use of narrow-band filters such
as the H-alpha of J-PLUS survey. This is shown in the Bottom-right
panel of Figure 1 by Infante-Sainz et al. 2024
(https://arxiv.org/abs/2401.03814), in this case the G channel has been
substituted by the image corresponding to the H-alpha filter to show the
star formation regions. Therefore, please use the weights with caution,
as it can significantly affect the output and misinform your
readers/viewers.

If you do apply weights be sure to report the weights in the caption
of the image (beside the filters that were used for each channel). With
great power there must also come great responsibility!

Two additional transformations are available to modify the appearance
of the output color image. The linear transformation combines bias
adjustment and contrast enhancement through the ‘--bias’ and
‘--contrast’ options. In most cases, only the contrast adjustment is
necessary to improve the quality of the color image. To illustrate the
impact of adjusting image contrast, we will generate an image with
higher contrast and compare with the previous one.

$ astscript-color-faint-gray $R $G $B $params --contrast=2 \
--output=m51-contrast.pdf

When you compare this (‘m51-contrast.pdf’) with the previous output
(‘m51.pdf’), you see that the colored parts are now much more clear!
Use this option also with caution because it may happen that the bright
parts become saturated.

Another option available for transforming the image appearance is the
gamma correction, a non-linear transformation that can be useful in
specific cases. You can experiment with different gamma values to
observe the impact on the resulting image. Lower gamma values will
enhance faint structures, while higher values will emphasize brighter
regions. Let's have a look by giving two very different values to it
with the simple loop below:

$ for g in 0.4 2.0; do \
astscript-color-faint-gray $R $G $B $params --contrast=2 \
--gamma=$g --output=m51-gamma-$g.pdf; \
done

Comparing the last three files (‘m51-contrast.pdf’,
‘m51-gamma-0.4.pdf’ and ‘m51-gamma-2.0.pdf’), you will clearly see the
effect of the ‘--gamma’.

Instead of using a combination of the three input images for the gray
background, you can introduce a fourth image that will be used for
generating the gray background. This image is referred to as the "K"
channel and may be useful when a particular filter is deeper, has unique
characteristics, or you have built by some custom processing to show the
diffuse features better. In this case, this image will be used for
defining the ‘--colorval’ and ‘--grayval’ thresholds, but the rationale
remains the same as explained earlier.

Two additional options are available to smooth different regions by
convolving with a Gaussian kernel: ‘--colorkernelfwhm’ for smoothing
color regions and ‘--graykernelfwhm’ for convolving gray regions. The
value specified for these options represents the full width at half
maximum of the Gaussian kernel.

Finally, another commonly useful feature is ‘--markoptions’: it
allows you to mark and label the final output image with vector graphics
over the color image. The arguments passed through this option are
directly passed to ConvertType for the generation of the output image.
This feature was already used in *note Marking objects for publication::
of the *note General program usage tutorial::; see there for a more
complete introduction.

Let's create four marks/labels just to illustrate the procedure
within ‘astscript-color-faint-gray’. First we need to create a table
that contains the parameters for creating the marks (coordinates, shape,
size, colors, etc.). In order to have an example that could be easily
salable to more marks, with elaborated options let's create it by parts:
the header with the column names, and the parameters. With the
following commands, we'll create the header that contains the column
metadata.

echo "# Column 1: ra [pix, f32] RA coordinate" > markers.txt
echo "# Column 2: dec [pix, f32] Dec coordinate" >> markers.txt
echo "# Column 3: shape [none, u8] Marker shape" >> markers.txt
echo "# Column 4: size [pix, f32] Marker Size" >> markers.txt
echo "# Column 5: aratio [none, f32] Axis ratio" >> markers.txt
echo "# Column 6: angle [deg, f32] Position angle" >> markers.txt
echo "# Column 7: color [none, u8] Marker color" >> markers.txt

Next is to create the parameters that define the markers. In this
case, with the lines below we create four markers (cross, ellipse,
square, and line) at different positions, with different shapes, and
colors. These lines are appended to the header file created previously.
echo "400.00 400.00 3 60.000 0.50 0.000 8" >> markers.txt
echo "1800.0 400.00 4 120.00 0.30 45.00 58" >> markers.txt
echo "400.00 1800.0 6 180.00 1.00 0.000 85" >> markers.txt
echo "1800.0 1800.0 8 240.00 1.00 -45.0 25" >> markers.txt

Now that we have the table containing the definition of the markers,
we use the ‘--markoptions’ option of this script. This option will pass
what ever is given to it directly to ConvertType, so you can use all the
options in *note Drawing with vector graphics::. For this basic
example, let's give it the following options:

markoptions="--mode=img \
--sizeinarcsec \
--markshape=shape \
--markrotate=angle \
--markcolor=color \
--marks=markers.txt \
--markcoords=ra,dec \
--marksize=size,aratio"

The last step consists in executing the script with the option that
provides all the markers options.

$ astscript-color-faint-gray $R $G $B $params --contrast=2 \
--markoptions="$markoptions" \
--output=m51-marked.pdf

Open the ‘m51-marked.pdf’ and check that the four markers have been
printed on the image. With this quick example we just show the
possibility of drawing markers on images very easily. This task can be
automated, for example by plotting markers from a given catalog at
specific positions, and so on. Note that there are many other options
for customize your markers/drawings over an output of ConvertType, see
*note Drawing with vector graphics:: and *note Marking objects for
publication::.

Congratulations! By following the tutorial up to this point, we have
been able to reproduce three images of Infante-Sainz et al. 2024
(https://arxiv.org/abs/2401.03814). You can see the commands that were
used to generate them within the reproducible source of that paper at
<https://codeberg.org/gnuastro/papers/src/branch/color-faint-gray>.
Remember that this paper is exactly reproducible with Maneage, so you
can explore and build the entire paper by yourself. For more on
Maneage, see Akhlaghi et al. 2021
(https://ui.adsabs.harvard.edu/abs/2021CSE....23c..82A).

This tutorial provided a general overview of the various options to
construct a color image from three different FITS images using the
‘astscript-color-faint-gray’ script. Keep in mind that the optimal
parameters for generating the best color image depend on your specific
goals and the quality of your input images. We encourage you to follow
this tutorial with the provided J-PLUS images and later with your own
dataset. See *note Color images with gray faint regions:: for more
information, and please consider citing Infante-Sainz et al. 2024
(https://arxiv.org/abs/2401.03814) if you use this script in your work
(the full BibTeX entry of this paper will be given to you with the
‘--cite’ option).

File: gnuastro.info, Node: Zero point of an image, Next: Pointing pattern design, Prev: Color images with full dynamic range, Up: Tutorials

2.7 Zero point of an image
==========================

The "zero point" of an image is astronomical jargon for the calibration
factor of its pixel values; allowing us to convert the raw pixel values
to physical units. It is therefore a critical step during data
reduction. For more on the definition and importance of the zero point
magnitude, see *note Brightness flux magnitude:: and *note Zero point
estimation::.

In this tutorial, we will use Gnuastro's ‘astscript-zeropoint’, to
estimate the zero point of a single exposure image from the J-PLUS
survey (https://www.j-plus.es), while using an SDSS
(http://www.sdss.org) image as reference (recall that all SDSS images
have been calibrated to have a fixed zero point of 22.5). In this case,
both images that we are using were taken with the SDSS _r_ filter. See
Eskandarlou et al. 2023 (https://arxiv.org/abs/2312.04263).

*Same filters and SVO filter database:* It is very important that both
your images are taken with the same filter. When looking at filter
names, don't forget that different filter systems sometimes have the
same names for one filter, such as the name "R"; which is used in both
the Johnson and SDSS filter systems. Hence if you confront an image in
the "R" or "r" filter, double check to see exactly which filter system
it corresponds to. If you know which observatory your data came from,
you can use the SVO database (http://svo2.cab.inta-csic.es/theory/fps)
to confirm the similarity of the transmission curves of the filters of
your input and reference images. SVO contains the filter data for many
of the observatories world-wide.

* Menu:

* Zero point tutorial with reference image:: Using a reference image.
* Zero point tutorial with reference catalog:: Using a reference catalog.

File: gnuastro.info, Node: Zero point tutorial with reference image, Next: Zero point tutorial with reference catalog, Prev: Zero point of an image, Up: Zero point of an image

2.7.1 Zero point tutorial with reference image
----------------------------------------------

First, let’s create a directory named ‘tutorial-zeropoint’ to keep
things clean and work in that. Then, with the commands below, you can
download an image from J-PLUS and SDSS. To speed up the analysis, the
image is cropped to have a smaller region around its center.

$ mkdir tutorial-zeropoint
$ cd tutorial-zeropoint
$ jplusdr2=http://archive.cefca.es/catalogues/vo/siap/jplus-dr2/reduced
$ wget $jplusdr2/get_fits?id=771463 -O jplus.fits.fz
$ astcrop jplus.fits.fz --center=107.7263,40.1754 \
--width=0.6 --output=jplus-crop.fits

Although we cropped the J-PLUS image, it is still very large in
comparison with the SDSS image (the J-PLUS field of view is almost
$1.5\times1.5$ deg$^2$, while the field of view of SDSS in each filter
is almost $0.3\times0.5$ deg$^2$). Therefore, let's download two SDSS
images (and then decompress them) in the region of the cropped J-PLUS
image to have a more accurate result compared to a single SDSS
footprint: generally, your zero point estimation will have less scatter
with more overlap between your reference image(s) and your input image.

$ sdssbase=https://dr12.sdss.org/sas/dr12/boss/photoObj/frames
$ wget $sdssbase/301/6509/5/frame-r-006509-5-0115.fits.bz2 \
-O sdss1.fits.bz2
$ wget $sdssbase/301/6573/5/frame-r-006573-5-0174.fits.bz2 \
-O sdss2.fits.bz2
$ bunzip2 sdss1.fits.bz2
$ bunzip2 sdss2.fits.bz2

To have a feeling of the data, let's open the three images with
‘astscript-fits-view’ using the command below. Wait a few seconds to
see the three images "blinking" one after another. The largest one is
the J-PLUS crop and the two smaller ones that partially cover it in
different regions are from SDSS.

$ astscript-fits-view sdss1.fits sdss2.fits jplus-crop.fits \
--ds9extra="-lock frame wcs -single -zoom to fit -blink yes"

The test above showed that the three images are already
astrometrically calibrated (the coverage of the pixel positions on the
sky is correct in both). To confirm, you can zoom-in to a certain
object and confirm it on a pixel level. It is always good to do the
visual check above when you are confronted with new images (and may not
be confident about the accuracy of the astrometry). Do not forget that
the goal here is to find the calibration of pixel values; and that we
assume pixel positions are already calibrated (the image already has a
good astrometry).

The SDSS images are Sky subtracted, while this single-exposure J-PLUS
image still contains the counts related to the Sky emission within them.
In the J-PLUS survey, the sky-level in each pixel is kept in a separate
‘BACKGROUND_MODEL’ HDU of ‘jplus.fits.fz’; this allows you to use a
different sky if you like. The SDSS image FITS files also have multiple
extensions. To understand our inputs, let's have a fast look at the
basic info of each:

$ astfits sdss1.fits
Fits (GNU Astronomy Utilities) 0.22
Run on Fri Apr 14 11:24:03 2023
-----
HDU (extension) information: 'sdss1.fits'.
Column 1: Index (counting from 0, usable with '--hdu').
Column 2: Name ('EXTNAME' in FITS standard, usable with '--hdu').
('n/a': no name in HDU metadata)
Column 3: Image data type or 'table' format (ASCII or binary).
Column 4: Size of data in HDU.
Column 5: Units of data in HDU (only images).
('n/a': no unit in HDU metadata, or HDU is a table)
-----
0 n/a float32 2048x1489 nanomaggy
1 n/a float32 2048 n/a
2 n/a table_binary 1x3 n/a
3 n/a table_binary 1x31 n/a

$ astfits jplus.fits.fz
Fits (GNU Astronomy Utilities) 0.22
Run on Fri Apr 14 11:21:30 2023
-----
HDU (extension) information: 'jplus.fits.fz'.
Column 1: Index (counting from 0, usable with '--hdu').
Column 2: Name ('EXTNAME' in FITS standard, usable with '--hdu').
('n/a': no name in HDU metadata)
Column 3: Image data type or 'table' format (ASCII or binary).
Column 4: Size of data in HDU.
Column 5: Units of data in HDU (only images).
('n/a': no unit in HDU metadata, or HDU is a table)
-----
0 n/a no-data 0 n/a
1 IMAGE float32 9216x9232 adu
2 MASKED_PIXELS int16 9216x9232 n/a
3 BACKGROUND_MODEL float32 9216x9232 n/a
4 MASK_MODEL uint8 9216x9232 n/a

Therefore, in order to be able to compare the SDSS and J-PLUS images,
we should first subtract the sky from the J-PLUS image. To do that, we
can either subtract the ‘BACKGROUND_MODEL’ HDU from the ‘IMAGE’ HDU
using *note Arithmetic::, or we can use *note NoiseChisel:: to find a
good sky ourselves. As scientists we like to tweak and be creative, so
let's estimate it ourselves with the command below. Generally, you may
not have a pre-estimated Sky estimation like above, so you should be
prepared to subtract the sky yourself.

$ astnoisechisel jplus-crop.fits --output=jplus-nc.fits
$ astscript-fits-view jplus-nc.fits

Notice that there is a relatively bright star in the center-bottom of
the image. In the "Cube" window, click on the "Next" button to see the
‘DETECTIONS’ HDU. The large footprint of the bright star is obvious.
Press the "Next" button one more time to get to the ‘SKY’ HDU. You see
that in the center-bottom, the footprint of the large star is clearly
visible in the measured Sky level. This is not good! With Sky values
above 54 ADU in the center of the star (the white pixels). This
over-subtracted Sky level in part of the image will affect your
magnitude measurements and thus the zero point!

In *note General program usage tutorial::, we have a section on *note
NoiseChisel optimization for detection::, there is also a full tutorial
on this in *note Detecting large extended targets::. Therefore, we will
not go into the details of NoiseChisel optimization here. Given the
large images of J-PLUS, we will increase the tile-size to $100\times100$
pixels and the number of neighbors to identify outlying tiles to 50
(these are usually the first parameters you should start editing when
you are confronted with a new image). After the second command, check
the ‘SKY’ extension to confirm that there is no footprint of any bright
object there. You will still see a gradient, but note the minimum and
maximum values of the Sky level: their difference is more than 26 times
smaller than the noise standard deviation (so statistically speaking, it
is pretty flat!)

$ astnoisechisel jplus-crop.fits --output=jplus-nc.fits \
--tilesize=100,100 --outliernumngb=50
$ astscript-fits-view jplus-nc.fits

## Check that the gradient in the sky is statistically negligible.
$ aststatistics jplus-nc.fits -hSKY --minimum --maximum \
| awk '{print $2-$1}'
0.32809
$ aststatistics jplus-nc.fits -hSKY_STD --median
8.377977e+00

We are now ready to find the zero point! First, let's run the
‘astscript-zeropoint’ with ‘--help’ to see the option names (recall that
you can see more details of each option in *note Invoking
astscript-zeropoint::). For the first time, let's use the script in the
most simple state possible. We will keep only the essential options:
the names of the input and reference images (and their HDUs), the name
of the output, and also two apertures with radii of 3 arcsec to start
with:

$ astscript-zeropoint --help
$ astscript-zeropoint jplus-nc.fits --hdu=INPUT-NO-SKY \
--refimgs=sdss1.fits,sdss2.fits \
--output=jplus-zeropoint.fits \
--refimgszp=22.5,22.5 \
--refimgshdu=0,0 \
--aperarcsec=3

The output is a FITS table (because generally, you will give more
apertures and choose the best one based on a higher-level analysis).
Let's check the output's internal structure with Gnuastro's ‘astfits’
program.

$ astfits jplus-zeropoint.fits
-----
0 n/a no-data 0 n/a
1 ZEROPOINTS table_binary 1x3 n/a
2 APER-3 table_binary 321x2 n/a

You can see that there are two HDUs in this file. The HDU names give
a hint, so let's have a look at each extension with Gnuastro's
‘asttable’ program:

$ asttable jplus-zeropoint.fits --hdu=1 -i
--------
jplus-zeropoint.fits (hdu: 1)
------- ----- ---- -------
No.Name Units Type Comment
------- ----- ---- -------
1 APERTURE arcsec float32 n/a
2 ZEROPOINT mag float32 n/a
3 ZPSTD mag float32 n/a
--------
Number of rows: 1
--------

As you can see, in the first extension, for each of the apertures you
requested (‘APERTURE’), there is a zero point (‘ZEROPOINT’) and the
standard deviation of the measurements on the apertures (‘ZPSTD’). In
this case, we only requested one aperture, so it only has one row. Now,
let's have a look at the next extension:

$ asttable jplus-zeropoint.fits --hdu=2 -i
--------
jplus-zeropoint.fits (hdu: 2)
------- ----- ---- -------
No.Name Units Type Comment
------- ----- ---- -------
1 MAG-REF f32 float32 Magnitude of reference.
2 MAG-DIFF f32 float32 Magnitude diff with input.
--------
Number of rows: 321
--------

It contains a table of measurements for the aperture with the least
scatter. In this case, we only gave one aperture, so it is the same.
If you give multiple apertures, only the one with least scatter will be
present by default. In the ‘MAG-REF’ column you see the magnitudes
within each aperture on the reference (SDSS) image(s). The ‘MAG-DIFF’
column contains the difference of the input (J-PLUS) and reference
(SDSS) magnitudes for each aperture (see *note Zero point estimation::).
The two catalogs, created by the aperture photometry from the SDSS
images, are merged into one so that there are more stars to compare.
Therefore, no matter how many reference images you provide, there will
only be a single table here. If the two SDSS images overlapped, each
object in the overlap region would have two rows (one row for the
measurement from one SDSS image, and another from the measurement from
the other).

Now that we have obtained the zero point of the J-PLUS image, let's
go a little deeper into lower-level details of how this script operates.
This will help you better understand what happened and how to interpret
and improve the outputs when you are confronted with a new image and
strange outputs.

To keep intermediate results the ‘astscript-zeropoint’ script keeps
temporary files in a temporary directory and later deletes it (and all
the intermediate products). If you like to check the temporary files of
the intermediate steps, you can use ‘--keeptmp’ option to not remove
them.

Let's take a closer look into the contents of each HDU. First, we'll
use Gnuastro’s ‘asttable’ to see the measured zero point for this
aperture. We are using ‘-Y’ to have human-friendly (non-scientific!)
numbers (which are sufficient here) and ‘-O’ to also show the metadata
of each column at the start.

$ asttable jplus-zeropoint.fits -Y -O
# Column 1: APERTURE [arcsec,f32,] Aperture used.
# Column 2: ZEROPOINT [mag ,f32,] Zero point (sig-clip median).
# Column 3: ZPSTD [mag ,f32,] Zero point Standard deviation.
3.000 26.435 0.057

Now, let's have a look at the first 10 rows of the second (‘APER-3’)
extension. From the previous check we did above, we see that it
contains 321 rows!

$ asttable jplus-zeropoint.fits -Y -O --hdu=APER-3 --head=10
# Column 1: MAG-REF [f32,f32,] Magnitude of reference.
# Column 2: MAG-DIFF [f32,f32,] Magnitude diff with input.
16.461 30.035
16.243 28.209
15.427 26.427
20.064 26.459
17.334 26.425
20.518 26.504
17.100 26.400
16.919 26.428
17.654 26.373
15.392 26.429

But the table above is hard to interpret, so let's plot it. To do
this, we'll use the same ‘astscript-fits-view’ command above that we
used for images. It detects if the file has a image or table HDU and
will call DS9 or TOPCAT respectively. You can also use any other
plotter you like (TOPCAT is not part of Gnuastro), this script just
calls it.

$ astscript-fits-view jplus-zeropoint.fits --hdu=APER-3

After ‘TOPCAT’ opens, you can select the "Graphics" menu and then
"Plain plot". This will show a plot with the SDSS (reference image)
magnitude on the horizontal axis and the difference of magnitudes
between the the input and reference (the zero point) on the vertical
axis.

In an ideal world, the zero point should be independent of the
magnitude of the different stars that were used. Therefore, this plot
should be a horizontal line (with some scatter as we go to fainter
stars). But as you can see in the plot, in the real world, this
expected behavior is seen only for stars with magnitudes about 16 to 19
in the reference SDSS images. The stars that are brighter than 16 are
saturated in one (or both) surveys(1). Therefore, they do not have the
correct magnitude or mag-diff. You can check some of these stars
visually by using the blinking command above and zooming into some of
the brighter stars in the SDSS images.

On the other hand, it is natural that we cannot measure accurate
magnitudes for the fainter stars because the noise level (or "depth") of
each image is limited. As a result, the horizontal line becomes wider
(scattered) as we go to the right (fainter magnitudes on the horizontal
axis). So, let's limit the range of used magnitudes from the SDSS
catalog to calculate a more accurate zero point for the J-PLUS image.
For this reason, we have the ‘--magnituderange’ option in
‘astscript-zeropoint’.

*Necessity of sky subtraction:* To obtain this horizontal line, it is
very important that both your images have been sky subtracted. Please,
repeat the last ‘astscript-zeropoint’ command above only by changing the
input file to ‘jplus-crop.fits’. Then use Gnuastro’s
‘astscript-fits-view’ again to draw a plot with ‘TOPCAT’ (also same as
above). Instead of a horizontal line, you will see _a sloped line_ in
the magnitude range above! This happens because the sky level acts as a
source of constant signal in all apertures, so the magnitude difference
will not be independent of the star's magnitude, but dependent on it
(the measurement on a fainter star will be dominated by the sky level).

*Remember:* if you see a sloped line instead of a horizontal line,
the input or reference image(s) are not sky subtracted.

Another key parameter of this script is the aperture size
(‘--aperarcsec’) for the aperture photometry of images. On one hand, if
the selected aperture is too small, you will be at the mercy of the
differing PSFs between your input and reference image(s): part of the
light of the star will be lost in the image with the worse PSF. On the
other hand, with large aperture size, the light of neighboring objects
(stars/galaxies) can affect the photometry. We should select an
aperture radius of the same order than the one used in the reference
image, typically 2 to 3 times the PSF FWHM of the images. For now,
let's assume the values 2, 3, 4, 5, and 6 arcsec for the aperture sizes
parameter. The script will compare the result for several aperture
sizes and choose the one with least standard deviation value, ‘ZPSTD’
column of the ‘ZEROPOINTS’ HDU.

Let's re-run the script with the following changes:
• Using ‘--magnituderange’ to limit the stars used for estimating the
zero point.
• Giving more values for aperture size to find the best for these two
images as explained above.
• Call ‘--keepzpap’ option to keep the result of matching the
catalogs done with the selected apertures in the different
extensions of the output file.

$ astscript-zeropoint jplus-nc.fits --hdu=INPUT-NO-SKY \
--refimgs=sdss1.fits,sdss2.fits \
--output=jplus-zeropoint.fits \
--refimgszp=22.5,22.5 \
--aperarcsec=2,3,4,5,6 \
--magnituderange=16,18 \
--refimgshdu=0,0 \
--keepzpap

Now, check number of HDU extensions by ‘astfits’.

$ astfits jplus-zeropoint.fits
-----
0 n/a no-data 0 n/a
1 ZEROPOINTS table_binary 5x3 n/a
2 APER-2 table_binary 319x2 n/a
3 APER-3 table_binary 321x2 n/a
4 APER-4 table_binary 323x2 n/a
5 APER-5 table_binary 323x2 n/a
6 APER-6 table_binary 325x2 n/a

You can see that the output file now has a separate HDU for each
aperture (thanks to ‘--keepzpap’.) The ‘ZEROPOINTS’ hdu contains the
final zero point values for each aperture and their error. The best
zero point value belongs to the aperture that has the least scatter (has
the lowest standard deviation). The rest of extensions contain the zero
point value computed within each aperture (as discussed above).

Let's check the different tables by plotting all magnitude tables at
the same time with ‘TOPCAT’.

$ astscript-fits-view jplus-zeropoint.fits

After ‘TOPCAT’ has opened take the following steps:
1. From the "Graphics" menu, select "Plain plot". You will see the
last HDU's scatter plot open in a new window (for ‘APER-6’, with
red points). The Bottom-left panel has the logo of a red-blue
scatter plot that has written ‘6:jplus-zeropoint.fits’ in front of
it (showing that this is the 6th HDU of this file). In the
bottom-right panel, you see the names of the columns that are being
displayed.
2. In the "Layers" menu, Click on "Add Position Control". On the
bottom-left panel, you will notice that a new blue-red scatter plot
has appeared but it just says ‘<no table>’. In the bottom-right
panel, in front of "Table:", select any other extension. This will
plot the same two columns of that extension as blue points.
Zoom-in to the region of the horizontal line to see/compare the
different scatters.

Change the HDU given to "Table:" and see the distribution of zero
points for the different apertures.

The manual/visual operation above is critical if this is your first
time with a new dataset (it shows all kinds of systematic biases (like
the Sky issue above)! But once you know your data has no systematic
biases, choosing between the different apertures is not easy visually!
Let's have a look at the table the ‘ZEROPOINTS’ HDU (we don't need to
explicitly call this HDU since it is the first one):

$ asttable jplus-zeropoint.fits -O -Y
# Column 1: APERTURE [arcsec,f32,] Aperture used.
# Column 2: ZEROPOINT [mag ,f32,] Zero point (sig-clip median).
# Column 3: ZPSTD [mag ,f32,] Zero point Standard deviation.
2.000 26.405 0.028
3.000 26.436 0.030
4.000 26.448 0.035
5.000 26.458 0.042
6.000 26.466 0.056

The most accurate zero point is the one where ‘ZPSTD’ is the
smallest. In this case, minimum of ‘ZPSTD’ is with radii of 2 and 3
arcseconds. Run the ‘astscript-fits-view’ command above again to open
TOPCAT. Let's focus on the magnitude plots in these two apertures and
determine a more accurate range of magnitude. The more reliable option
is the range between 16.4 (where we have no saturated stars) and 18.5
mag (fainter than this, the scatter becomes too strong). Finally, let's
set some more apertures between 2 and 3 arcseconds radius:

$ astscript-zeropoint jplus-nc.fits --hdu=INPUT-NO-SKY \
--refimgs=sdss1.fits,sdss2.fits \
--output=jplus-zeropoint.fits \
--magnituderange=16.4,18.5 \
--refimgszp=22.5,22.5 \
--aperarcsec=2,2.5,3,3.5,4 \
--refimgshdu=0,0 \
--keepzpap

$ asttable jplus-zeropoint.fits -Y
2.000 26.405 0.037
2.500 26.425 0.033
3.000 26.436 0.034
3.500 26.442 0.039
4.000 26.449 0.044

The aperture with the least scatter is therefore the 2.5 arcsec
radius aperture, giving a zero point of 26.425 magnitudes for this
image. However, you can see that the scatter for the 3 arcsec aperture
is also acceptable. Actually, the ‘ZPSTD’ for of the 2.5 and 3 arcsec
apertures only have a difference of $3\%$ ($=
(0.034−0.0333)/0.033\times100$). So simply choosing the minimum is just
a first-order approximation (which is accurate within
$26.436−26.425=0.011$ magnitudes)

Note that in aperture photometry, the PSF plays an important role
(because the aperture is fixed but the two images can have very
different PSFs). The aperture with the least scatter should also
account for the differing PSFs. Overall, please, always check the
different and intermediate steps to make sure the parameters are the
good so the estimation of the zero point is correct.

If you are happy with the minimum, you don't have to search for the
minimum aperture or its corresponding zero point yourself. This script
has written it in ‘ZPVALUE’ keyword of the table. With the first
command, we also see the name of the file also, (you can use this on
many files for example). With the second command, we are only printing
the number by adding the ‘-q’ (or ‘--quiet’) option (this is useful in a
script where you want to write the value in a shell variable to use
later).

$ astfits jplus-zeropoint.fits --keyvalue=ZPVALUE
jplus-zeropoint.fits 2.642512e+01

$ astfits jplus-zeropoint.fits --keyvalue=ZPVALUE -q
2.642512e+01

Generally, this script will write the following FITS keywords (all
starting with ‘ZP’) for your future reference in its output:

$ astfits jplus-zeropoint.fits -h1 | grep ^ZP
ZPAPER = 2.5 / Best aperture.
ZPVALUE = 26.42512 / Best zero point.
ZPSTD = 0.03276644 / Best std. dev. of zeropoint.
ZPMAGMIN= 16.4 / Min mag for obtaining zeropoint.
ZPMAGMAX= 18.5 / Max mag for obtaining zeropoint.

Using the ‘--keyvalue’ option of the *note Fits:: program, you can
easily get multiple of the values in one run (where necessary):

$ astfits jplus-zeropoint.fits --hdu=1 --quiet \
--keyvalue=ZPAPER,ZPVALUE,ZPSTD
2.500000e+00 2.642512e+01 3.276644e-02

---------- Footnotes ----------

(1) To learn more about saturated pixels and recognition of the
saturated level of the image, please see *note Saturated pixels and
Segment's clumps::

File: gnuastro.info, Node: Zero point tutorial with reference catalog, Prev: Zero point tutorial with reference image, Up: Zero point of an image

2.7.2 Zero point tutorial with reference catalog
------------------------------------------------

In *note Zero point tutorial with reference image::, we explained how to
use the ‘astscript-zeropoint’ for estimating the zero point of one image
based on a reference image. Sometimes there is not a reference image
and we need to use a reference catalog. Fortunately,
‘astscript-zeropoint’ can also use the catalog instead of the image to
find the zero point.

To show this, let's download a catalog of SDSS in the area that
overlaps with the cropped J-PLUS image (used in the previous section).
For more on Gnuastro's Query program, please see *note Query::. The
columns of ID, RA, Dec and magnitude in the SDSS _r_ filter are called
by their name in the SDSS catalog.

$ astquery vizier \
--dataset=sdss12 \
--overlapwith=jplus-crop.fits \
--column=objID,RA_ICRS,DE_ICRS,rmag \
--output=sdss-catalog.fits

To visualize the position of the SDSS objects over the J-PLUS image,
let's use ‘astscript-ds9-region’ (for more details please see *note SAO
DS9 region files from table::) with the command below (it will
automatically open DS9 and load the regions it created):

$ astscript-ds9-region sdss-catalog.fits \
--column=RA_ICRS,DE_ICRS \
--color=red --width=3 --output=sdss.reg \
--command="ds9 jplus-nc.fits[INPUT-NO-SKY] \
-scale zscale"

Now, we are ready to estimate the zero point of the J-PLUS image
based on the SDSS catalog. To download the input image and understand
how to use the ‘astscript-zeropoint’, please see *note Zero point
tutorial with reference image::.

Many of the options (like the aperture size) and magnitude range are
the same so we will not discuss them further. You will notice that the
only substantive difference of the command below with the last command
in the previous section is that we are using ‘--refcat’ instead of
‘--refimgs’. There are also some cosmetic differences for example a new
output name, not using ‘--refimgszp’ since it is only necessary for
images) and the ‘--*column’ options which are used to identify the names
of the necessary columns of the input catalog:

$ astscript-zeropoint jplus-nc.fits --hdu=INPUT-NO-SKY \
--refcat=sdss-catalog.fits \
--refcatmag=rmag \
--refcatra=RA_ICRS \
--refcatdec=DE_ICRS \
--output=jplus-zeropoint-cat.fits \
--magnituderange=16.4,18.5 \
--aperarcsec=2,2.5,3,3.5,4 \
--keepzpap

Let's inspect the output with the command below.

$ asttable jplus-zeropoint-cat.fits -Y
2.000 26.337 0.034
2.500 26.386 0.036
3.000 26.417 0.041
3.500 26.439 0.043
4.000 26.455 0.050

As you see, the values and standard deviations are very similar to
the results we got previously in *note Zero point tutorial with
reference image::. The Standard deviations are generally a little
higher here because we didn't do the photometry ourselves, but they are
statistically similar.

Before we finish, let's open the two outputs (from a reference image
and reference catalog) with the command below. To confirm how they
compare, we are showing the result for ‘APER-3’ extension in both
(following the TOPCAT plotting recipe in *note Zero point tutorial with
reference image::).

$ astscript-fits-view jplus-zeropoint.fits jplus-zeropoint-cat.fits \
-hAPER-3