Angular coverage on the sky (GNU Astronomy Utilities)

Next: Cosmological coverage and visualizing tables, Previous: Dataset inspection and cropping, Up: General program usage tutorial [Contents][Index]

2.1.5 Angular coverage on the sky ¶

The cropped images in Dataset inspection and cropping are the deepest images we currently have of the sky. The first thing that comes to mind may be this: “How large is this field on the sky?”.

More accurate method: the steps mentioned in this section are primarily designed to help you get familiar with the FITS WCS standard and some shells scripting. The accuracy of this method will decrease as your image becomes large (on the scale of degrees). For an accurate method, see Area of non-blank pixels on sky.

You can get a fast and crude answer with Gnuastro’s Fits program, using this command:

$ astfits flat-ir/xdf-f160w.fits --skycoverage

It will print the sky coverage in two formats (all numbers are in units of degrees for this image): 1) the image’s central RA and Dec and full width around that center, 2) the range of RA and Dec covered by this image. You can use these values in various online query systems. You can also use this option to automatically calculate the area covered by this image. With the --quiet option, the printed output of --skycoverage will not contain human-readable text, making it easier for automatic (computer) processing:

$ astfits flat-ir/xdf-f160w.fits --skycoverage --quiet

The second row is the coverage range along RA and Dec (compare with the outputs before using --quiet). We can thus simply subtract the second from the first column and multiply it with the difference of the fourth and third columns to calculate the image area. We will also multiply each by 60 to have the area in arc-minutes squared.

$ astfits flat-ir/xdf-f160w.fits --skycoverage --quiet \
        | awk 'NR==2{print ($2-$1)*60*($4-$3)*60}'

The returned value is \(9.06711\) arcmin\(^2\). However, this method ignores the fact that many of the image pixels are blank! In other words, the image does cover this area, but there is no data in more than half of the pixels. So let’s calculate the area coverage over-which we actually have data.

The FITS world coordinate system (WCS) metadata standard contains the key to answering this question. Run the following command to see all the FITS keywords (metadata) for one of the images (almost identical with the other images because they are scaled to the same region of Sky):

$ astfits flat-ir/xdf-f160w.fits -h1

Look into the keywords grouped under the ‘World Coordinate System (WCS)’ title. These keywords define how the image relates to the outside world. In particular, the CDELT* keywords (or CDELT1 and CDELT2 in this 2D image) contain the “Coordinate DELTa” (or change in coordinate units) with a change in one pixel. But what is the units of each “world” coordinate? The CUNIT* keywords (for “Coordinate UNIT”) have the answer. In this case, both CUNIT1 and CUNIT1 have a value of deg, so both “world” coordinates are in units of degrees. We can thus conclude that the value of CDELT* is in units of degrees-per-pixel³¹.

With the commands below, we will use CDELT (along with the number of non-blank pixels) to find the answer of our initial question: “how much of the sky does this image cover?”. The lines starting with ## are just comments for you to read and understand each command. Do not type them on the terminal (no problem if you do, they will just not have any effect). The commands are intentionally repetitive in some places to better understand each step and also to demonstrate the beauty of command-line features like history, variables, pipes and loops (which you will commonly use as you become more proficient on the command-line).

Use shell history: Do not forget to make effective use of your shell’s history: you do not have to re-type previous command to add something to them (like the examples below). This is especially convenient when you just want to make a small change to your previous command. Press the “up” key on your keyboard (possibly multiple times) to see your previous command(s) and modify them accordingly.

Your locale does not use ‘.’ as decimal separator: on systems that do not use an English language environment, the dates, numbers, etc., can be printed in different formats (for example, ‘0.5’ can be written as ‘0,5’: with a comma). With the LC_NUMERIC line at the start of the script below, we are ensuring a unified format in the output of seq. For more, please see Numeric locale.

## Make sure that the decimal separator is a point in any environment.
$ export LC_NUMERIC=C

## See the general statistics of non-blank pixel values.
$ aststatistics flat-ir/xdf-f160w.fits

## We only want the number of non-blank pixels (add '--number').
$ aststatistics flat-ir/xdf-f160w.fits --number

## Keep the result of the command above in the shell variable `n'.
$ n=$(aststatistics flat-ir/xdf-f160w.fits --number)

## See what is stored the shell variable `n'.
$ echo $n

## Show all the FITS keywords of this image.
$ astfits flat-ir/xdf-f160w.fits -h1

## The resolution (in degrees/pixel) is in the `CDELT' keywords.
## Only show lines that contain these characters, by feeding
## the output of the previous command to the `grep' program.
$ astfits flat-ir/xdf-f160w.fits -h1 | grep CDELT

## Since the resolution of both dimensions is (approximately) equal,
## we will only read the value of one (CDELT1) with '--keyvalue'.
$ astfits flat-ir/xdf-f160w.fits -h1 --keyvalue=CDELT1

## We do not need the file name in the output (add '--quiet').
$ astfits flat-ir/xdf-f160w.fits -h1 --keyvalue=CDELT1 --quiet

## Save it as the shell variable `r'.
$ r=$(astfits flat-ir/xdf-f160w.fits -h1 --keyvalue=CDELT1 --quiet)

## Print the values of `n' and `r'.
$ echo $n $r

## Use the number of pixels (first number passed to AWK) and
## length of each pixel's edge (second number passed to AWK)
## to estimate the area of the field in arc-minutes squared.
$ echo $n $r | awk '{print $1 * ($2*60)^2}'

The output of the last command (area of this field) is 4.03817 (or approximately 4.04) arc-minutes squared. Just for comparison, this is roughly 175 times smaller than the average moon’s angular area (with a diameter of 30 arc-minutes or half a degree).

Some FITS writers do not use the CDELT convention, making it hard to use the steps above. In such cases, you can extract the pixel scale with the --pixelscale option of Gnuastro’s Fits program like the command below. Similar to the --skycoverage option above, you can also use the --quiet option to allow easy usage of the values in scripts.

$ astfits flat-ir/xdf-f160w.fits --pixelscale

AWK for table/value processing: As you saw above AWK is a powerful and simple tool for text processing. You will see it often in shell scripts. GNU AWK (the most common implementation) comes with a free and wonderful book in the same format as this book which will allow you to master it nicely. Just like this manual, you can also access GNU AWK’s manual on the command-line whenever necessary without taking your hands off the keyboard. Just run info awk.

Footnotes

(31)

With the FITS CDELT convention, rotation (PC or CD keywords) and scales (CDELT) are separated. In the FITS standard the CDELT keywords are optional. When CDELT keywords are not present, the PC matrix is assumed to contain both the coordinate rotation and scales. Note that not all FITS writers use the CDELT convention. So you might not find the CDELT keywords in the WCS metadata of some FITS files. However, all Gnuastro programs (which use the default FITS keyword writing format of WCSLIB) write their output WCS with the CDELT convention, even if the input does not have it. If your dataset does not use the CDELT convention, you can feed it to any (simple) Gnuastro program (for example, Arithmetic) and the output will have the CDELT keyword. See Section 8 of the FITS standard for more