2.1.20 FITS images in a publication ¶

In the previous section (Reddest clumps, cutouts and parallelization), we visually inspected the positions of the reddest objects using DS9. That is very good for an interactive inspection of the objects: you can zoom-in and out, you can do measurements, etc. Once the experimentation phase of your project is complete, you want to show these objects over the whole image in a report, paper or slides.

One solution is to use DS9 itself! For example, run the astscript-fits-view command of the previous section to open DS9 with the regions over-plotted. Click on the “File” menu and select “Save Image”. In the side-menu that opens, you have multiple formats to select from. Usually for publications, we want to show the regions and text (in the colorbar) in vector graphics, so it is best to export to EPS. Once you have made the EPS, you can then convert it to PDF with the epspdf command.

Another solution is to use Gnuastro’s ConvertType program. The main difference is that DS9 is a Graphic User Interface (GUI) program, so it takes relatively long (about a second) to load, and it requires many dependencies. This will slow-down automatic conversion of many files, and will make your code hard to move to another operating system. DS9 does have a command-line interface that you can use to automate the creation of each file, however, it has a very peculiar command-line interface and formats (like the “region” files). However, in ConvertType, there is no graphic interface, so it has very few dependencies, it is fast, and finally, it takes normal tables (in plain-text or FITS) as input. So in this concluding step of the analysis, let’s build a nice publication-ready plot, showing the positions of the reddest objects in the image for our paper.

In Reddest clumps, cutouts and parallelization, we already used ConvertType to make JPEG postage stamps. Here, we will use it to make a PDF image of the whole deep region. To start, let’s simply run ConvertType on the F160W image:

$ astconvertt flat-ir/xdf-f160w.fits -oxdf.pdf

Open the output in a PDF viewer. You see that it is almost fully black! Let’s see why this happens! First, with the two commands below, let’s calculate the maximum value, and the standard deviation of the sky in this image (using NoiseChisel’s output, which we found at the end of NoiseChisel optimization for detection). Note that NoiseChisel writes the median sky standard deviation before interpolation in the MEDSTD keyword of the SKY_STD HDU. This is more robust than the median of the Sky standard deviation image (which has gone through interpolation).

$ max=$(aststatistics nc/xdf-f160w.fits -hINPUT-NO-SKY --maximum)
$ skystd=$(astfits nc/xdf-f160w.fits -hSKY_STD --keyvalue=MEDSTD -q)

$ echo $max $skystd
58.8292 0.000410282

$ echo $max $skystd | awk '{print $1/$2}'
143387

In the last command above, we divided the maximum by the sky standard deviation. You see that the maximum value is more than \(140000\) times larger than the noise level! On the other hand common monitors or printers, usually have a maximum dynamic range of 8-bits, only allowing for \(2^8=256\) layers. This is therefore the maximum number of “layers” you can have in a common display formats like JPEG, PDF or PNG! Dividing the result above by 256, we get a layer spacing of

$ echo $max $skystd | awk '{print $1/$2/256}'
560.106

In other words, the first layer (which is black) will contain all the pixel values below \(\sim560\)! So all pixels with a signal-to-noise ratio lower than \(\sim560\) will have a black color since they fall in the first layer of an 8-bit PDF (or JPEG) image. This happens because by default we are assuming a linear mapping from floating point to 8-bit integers.

To fix this, we should move to a different mapping. A good, physically motivated, mapping is Surface Brightness (which is in log-scale, see Brightness, Flux, Magnitude and Surface brightness). Fortunately this is very easy to do with Gnuastro’s Arithmetic program, as shown in the commands below (using the known zero point⁴¹, and after calculating the pixel area in units of arcsec\(^2\)):

$ zeropoint=25.94
$ pixarcsec2=$(astfits nc/xdf-f160w.fits --pixelareaarcsec2)
$ astarithmetic nc/xdf-f160w.fits $zeropoint $pixarcsec2 counts-to-sb \
                --output=xdf-f160w-sb.fits

With the two commands below, first, let’s look at the dynamic range of the image now (dividing the maximum by the minimum), and then let’s open the image and have a look at it:

$ aststatistics xdf-f160w-sb.fits --minimum --maximum
$ astscript-fits-view xdf-f160w-sb.fits

The good news is that the dynamic range has now decreased to about 2! In other words, we can distribute the 256 layers of an 8-bit display over a much smaller range of values, and therefore better visualize the data. However, there are two important points to consider from the output of the first command and a visual inspection of the second.

• The largest pixel value (faintest surface brightness level) in the image is \(\sim43\)! This is far too low to be realistic, and is just due to noise. As discussed in Measuring the dataset limits, the \(3\sigma\) surface brightness limit of this image, over 100 arcsec\(^2\) is roughly 32.66 mag/arcsec\(^2\).
• You see many NaN pixels in between the galaxies! These are due to the fact that the magnitude is defined on a logarithmic scale and the logarithm of a negative number is not defined.

In other words, we should replace all NaN pixels, and pixels with a surface brightness value fainter than the image surface brightness limit to this limit. With the first command below, we will first extract the surface brightness limit from the catalog headers that we calculated before, and then call Arithmetic to use this limit.

$ sblimit=$(astfits cat/xdf-f160w.fits --keyvalue=SBLMAG -q)
$ astarithmetic nc/xdf-f160w.fits $zeropoint $pixarcsec2 \
                counts-to-sb set-sb \
                sb sb $sblimit gt sb isblank or $sblimit where \
                --output=xdf-f160w-sb.fits

Let’s convert this image into a PDF with the command below:

$ astconvertt xdf-f160w-sb.fits --output=xdf-f160w-sb.pdf

It is much better now and we can visualize many features of the FITS file (from the central structures of the galaxies and stars, to a little into the noise and their low surface brightness features. However, the image generally looks a little too gray! This is because of that bright star in the bottom half of the image! Stars are very sharp! So let’s manually tell ConvertType to set any pixel with a value less than (brighter than) 20 to black (and not use the minimum). We do this with the --fluxlow option:

$ astconvertt xdf-f160w-sb.fits --output=xdf-f160w-sb.pdf --fluxlow=20

We are still missing some of the diffuse flux in this PDF. This is because of those negative pixels that were set to NaN. To better show these structures, we should warp the image to larger pixels. So let’s warp it to a pixel grid where the new pixels are \(4\times4\) larger than the original pixels. But be careful that warping should be done on the original image, not on the surface brightness image. We should re-calculate the surface brightness image after the warping is one. This is because \(log(a+b)\ne log(a)+log(b)\). Recall that surface brightness calculation involves a logarithm, and warping involves addition of pixel values.

$ astwarp nc/xdf-f160w.fits --scale=1/4 --centeroncorner \
          --output=xdf-f160w-warped.fits

$ pixarcsec2=$(astfits xdf-f160w-warped.fits --pixelareaarcsec2)

$ astarithmetic xdf-f160w-warped.fits $zeropoint $pixarcsec2 \
                counts-to-sb set-sb \
                sb sb $sblimit gt sb isblank or $sblimit where \
                --output=xdf-f160w-sb.fits

$ astconvertt xdf-f160w-sb.fits --output=xdf-f160w-sb.pdf --fluxlow=20

Above, we needed to re-calculate the pixel area of the warpped image, but we did not need to re-calculate the surface brightness limit! The reason is that the surface brightness limit is independent of the pixel area (in its derivation, the pixel area has been accounted for). As a side-effect of the warping, the number of pixels in the image also dramatically decreased, therefore the volume of the output PDF (in bytes) is also smaller, making your paper/report easier to upload/download or send by email. This visual resolution is still more than enough for including on top of a column in your paper!

To summarize, and to keep the image for the next section in a separate directory, here are the necessary commands:

$ zeropoint=25.94
$ mkdir report-image
$ cd report-image
$ sblimit=$(astfits cat/xdf-f160w.fits --keyvalue=SBLMAG -q)
$ astwarp nc/xdf-f160w.fits --scale=1/4 --centeroncorner \
          --output=warped.fits
$ pixarcsec2=$(astfits warped.fits --pixelareaarcsec2)
$ astarithmetic warped.fits $zeropoint $pixarcsec2 \
                counts-to-sb set-sb \
                sb sb $sblimit gt sb isblank or $sblimit where \
                --output=sb.fits
$ astconvertt sb.fits --output=sb.pdf --fluxlow=20

Finally, let’s remove all the temporary files we built in the top-level tutorial directory:

$ rm *.fits *.pdf