IRAF paso a paso: 4 - Working with images


Image formats

IRAF can read, write and work with a variety of formats. Multidimensional images are allowed (for instance with 3, 4 up to 7 axes). Here is a list of the formats it uses:
 

Format Description
FITS Flexible Image Transport System - The "Universal" format.
Starting from V2.11, IRAF can use FITS in virtually any tasks dealing with images.
The FITS format now allows multiextension format, i.e. a number of images in the same physical file (used for mosaic data, for instance).
OIF (.imh + .pix) Old IRAF Format
This was the "standard" IRAF format. Now FITS is the de-facto format used by IRAF.
File image.imh contains the image header, file image.pix contains the pixel data.
In image.imh is the information about in which directory image.pix is.
STF (.hhh + .hhd) Space Telescope Format
The standard format used for HST images; analogous to OIF.
image.hhh is a plain ascii file; image.hhd contains pixel data.
Also, a variety of extension of the type ".??h" and ".??d" exists for HST images, depending on the instrument (WFPC2, FOC, etc.) and the type of data (raw, data quality file, etc.).
Multigroup format is often used (ex. the images of the 4 WFPC2 chips are in the same physical file).
Mask (.pl) A format used to identify bad pixels. Values is 0 for good pixels, non-zero for bad (masked) pixels.

Unless there are specific reasons to favour another format, you are better off using FITS as the default format. Indeed, FITS files can also be used directly in other applications (for instance read directly from DS9, or opened with PyFITS).

Also, IRAF can operate on tabular data both in ASCII and TABLE formats. The STSDAS ttools contains a variety of tasks to work on tables.
STSDAS tables are in binary format, and their header contains all the info (name of columns, units, format, etc.) about the table content.
The table I/O software can read/write also plain ASCII tables.


IRAF Coordinate System

IRAF supports 3 different coordinate systems:

Coord. system Description
Logical  Pixel units, relative to the current image. Starts from (1,1) = leftmost bottom pixel.
Physical  Similar to logical, but reference system is the one of the original image. Example:
cl> imcopy image.imh[10:500,20:500] subimage.imh
pixel (1,1) in logical system of subimage.imh has physical coordinates (10,20).
World Coordinate system having astronomical meaning. Might be Wavelenght, (RA,DEC), etc.
The image header contains all the info on the type of world coordinate system (=WCS), and the parameters for the transformation between pixel units and world coordinates.

Package imcoords contains tasks that deal with the World Coordinate System of an image. Also, several tasks may use the WCS, ex. imexamine, implot, listpixels, etc.
As an example, two images of the same field, but taken by two different instruments with different scales and orientations, may be registered to a common reference system by using only the WCSs of the two images (if accurate enough).
Sometimes, a task (or a user) might get confused by the physical or world coordinates. They may be deleted with the task wcsreset.

Image dev$wpix contains WCS keywords in its header, so it can be used as a test image to explore the WCS.


Input/Output from tape

Probably one of the most problem-prone and dreaded operations with IRAF.

Normally, images on tapes are in FITS format. To read them from a tape, do the following things (assuming you are already with the cl> prompt, are in the right directory, and the "dataio" package is loaded):

  1. Mount the tape on the appropriate device. Take note of the device name and of the machine it is connected to.
  2. cl> allocate machine!device
    where "machine" is the name of the workstation that owns the tape unit, "device" the name of the device (ex: nido!mtexa)
    If unsure of the device's name, do cl> devices to have a list of the devices, or consult the SIC Web pages.
  3. Read the images with the command:
    cl> rfits machine!device 1-10,15-98,110 WHT_
    Files 1 through 10, 15 through 98, and file 110 will be read. Name of files will be: WHT_0001, WHT_0002 ..., WHT_0015 ..., WHT_0110. Format will be the one defined by variable imtype
    Usually, the default hidden parameters for rfits are OK.
  4. Once finished, deallocate the tape (IMPORTANT).
    cl> deallocate machine!device
    and remove the tape.

Writing to tape is very similar. Do the same allocation/deallocation operations described above. Now to write to the tape:

You are not compelled to store files on tape in FITS format. In most circumstances, a tar archive is faster and simpler.


Getting info on an image

All the information about an image is contained in its header (except the pixel data, of course). You can see and also modify the header content.


Display an image

First you have to select your display tool (typically DS9). Then, use the display task to load the image into the image display.
Since there are a number of things you can do with this task (not to speak of the image display itself), let's have a look at a few illustrative examples:

Display also allows to mosaic more images on the image window, allowing a direct comparison, for instance, of the same field observed through different filters, or same filter but different epochs. Here you can find an example of script to do this.
There is also a CL task to do this: stsdas.graphics.sdisplay.mosaic_display

Do not forget to set the variable stdimage to the appropriate size. A complete list of the sizes available is in the file dev$imtoolrc.


Working with image sections

It is generally possible to work with just a region of an image. To define what region to work on, use the syntax: [Xstart:Xend,Ystart:Yend].
Ex: dev$pix[100:300,150:350] goes from pixel X=100 to X=300, Y=150 to Y=350, so it's a 201 x 201 region.
A "*" means: all pixels. Ex. [*,100:200] is a horizontal slice of the image, of same length in X as the original, but reduced Y size.

Many tasks accept image regions (in input; rarely on output); for instance:

Moreover, the image section notation, coupled with the imcopy and imtranspose tasks, allows simple coordinate transformations.

Command Description
imcopy dev$pix[*,-*] new1 flip the rows (image is upside down)
imcopy dev$pix[-*,*] new2 flip the columns (image is inverted left-right)
imcopy dev$pix[-*,-*] new3 image is rotated 180 degrees
imtranspose dev$pix[*,-*] new4 image is rotated 90 degrees counterclockwise
imtranspose dev$pix[-*,*] new5 image is rotated 90 degrees clockwise

Create empty images

Sometimes, it may be necessary to create empty images, i.e. with all pixels set to 0 or some other value.
The easiest thing is to use an already existing images, multiplying it by 0, using the imarith task.
Alternatively, you can use other tasks to create an image from scratch, such as:


There are essentially two ways of printing out an image:

  1. Using the "print option" of the image display (ex. DS9)
  2. Using the export task in package "dataio". This task provides a number of output formats (EPS, GIF, etc.), and a variety of lookup tables, color transformations, etc.

Moreover, a special task color.rgbsun creates a "quasitruecolor" image.
It needs three images of the same field, taken in three filters at different wavelengths. One image is associated to "blue", one to "green", one to "red". The output is a Sun 24-bit RGB rasterfile.
If the correspondence between filters and the above colors is good (for instance, a B, V, R set), the resulting image will have colors similar to the ones that would be observed by the naked eye. For an example, see:
file to create color image of a WFPC2 field

There is also a task called import which creates an IRAF image from a GIF, RAS, or other format. Its help is even more obscure than export's one...

Other options at the host level permit to save a displayed window as a raster image; for instance, in Unix the xv tool.


Working with tabular data

The TABLES package contains many tasks to handle tabular data of different types. Many analysis tasks in STSDAS indeed produce more or less complex tables as output. STSDAS tables are generally binary files with data stored in rows and columns; also, each element may be itself a vector, therefore forming sort of a 3-D table. Their default extension is ".tab".

Most of the tables-related tasks are in the ttools package. While they can unleash all of their power on well formatted and structured tables, many of those tasks can also work on plain ASCII tables, that is simple text files in row and column format. Here we illustrate a few basic operations we can do on ASCII tables (the file ascii_table.dat may be used as an example).

We can get info on the table by using the commands:
cl> tinfo ascii_table.dat ttout=yes  and 
cl>
tlcol ascii_table.dat nlist=4

To perform an arithmetic operation, we use the task tcalc:
cl> tcalc ascii_table.dat norm "(c7-390.0)/c3"
As we have not given any id to the column tables, we use "c7" and "c3" to refer to their position. The output column is called "norm" and is placed as last column.

We can select a subset of rows using task tselect:
cl> tselect ascii_table.dat ascii_table_subset.dat 'c5=="2003-03-31"&&c4>19.0'
In this specific case, we select all flat-fields taken on March 31 in the evening.

cl> tcalc ascii_table_subset.dat tdiff '60*(c4-19:28)'
We create column tdiff which says how many minutes since the sunset the image was taken. Notice that we can mix decimal and sexagesimal formats.

cl> tproject ascii_table_subset.dat ascii_table_project.dat "c1,c2,tdiff,norm"
Select relevant columns from the input table.

Now we can use any plotting program to display the sky brightness, through the various filters, as a function of time since sunset. Or, we can use task sgraph, briefly described in The sgraph task

Exercise:
The file tablafinalV.dat refers to a set of flat-field images taken before sunset. Column 2 lists their mean value, column 3 the exposure time and column 4 the time the image was taken. Can you produce a plot with in X the time before sunset (for simplicity assume sunset = 7:30:00), and in Y the sky brightness (in arbitrary units: just normalize the observed counts by the exposure time)?

Another interesting task is tjoin. It permits to join two tables into a new table on the basis of one (or more) common columns. For instance, if you have a table with a list of stars, with, say, ID, positions and magnitude, and another table with a similar, but different, list of objects with ID and radial velocities, you can produce a table with all objects common to both tables (from their ID) together with positions, magnitude and radial velocity.

As en exercise, we can join tabla_coord.dat and tabla_param.dat, and produce a final table containing all objects common to both tables:
cl> tjoin tabla_coord.dat tabla_param.dat tabla_out.dat c1 c1

A similar task is tmatch.