ChIPS - contouring images

Introduction

The Starting ChIPS thread describes how to start ChIPS.

The data used in this thread is available in the chips_data.tar.gz file.

unix% ls -1 chips/contour/
a4059_chandra.fits
a4059_nvss.fits
chandra.con
contour-overlay.state
nvss.con
projection.dat

The plots used in this thread are based on the data in the images a4059_chandra.fits and a4059_nvss.fits which contain images of the galaxy cluster Abell 4059 from Chandra and the NVSS. Both datasets were taken from the archive, with no re-processing applied. The Chandra image has been binned heavily (by a factor of 16) to increase the signal-to-noise ratio for each pixel. The projection.dat file contains three columns (x axis, y axis, statistic) from a projection run in Sherpa (they are from a spectral analysis of a galaxy cluster, but not Abell 4059).

The contour-overlay.state file is a ChIPS state file, and represents the finished version of Figure 7. This file allows you to view, and edit, the visualization without having to enter in all the preceeding commands. The chandra.con and nvss.con files are contour files created by ds9, and were used to create Figure 8.

Handling images with more-than two dimensions

The a4059_nvss.fits image is stored as a 4D array, as shown below:

unix% dmlist a4059_nvss.fits cols
 
--------------------------------------------------------------------------------
Columns for Image Block PRIMARY
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   PRIMARY[121,121,1,1] JY/BEAM      Real4(121x121x1x1) -Inf:+Inf            
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block PRIMARY
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    POS(X) = (#1) 
                 (Y)   (#2)
   2   3      Z                    = #3 
   3   4      #AXIS4               = #4 
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block PRIMARY
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+359.1696) +SIN[(-0.0042)* (POS(X)-(+61.0))]
                   (DEC)   (-34.6717 )      (+0.0042)  (   (Y) (+61.0)) 
   2   3      STOKES               = Z 
   3   4      FREQ                 = +1.4E+09  +100000000.0 * (#AXIS4  -1.0)

If this file is used as is, the presence of the extra dimensions will confuse make_figure() and add_contour(), resulting in the following error message:

chips ERROR: Data for contours must be regularly gridded

We will make use of the Data Model's bin specifier to remove these extra dimensions when we read in the file. Since the dmlist output above shows the sky coordinates are called POS we can use the syntax "[bin pos]", as shown below:

unix% dmlist "a4059_nvss.fits[bin pos]" cols
 
--------------------------------------------------------------------------------
Columns for Image Block PRIMARY_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   PRIMARY_IMAGE[121,121]              Real4(121x121) -Inf:+Inf            
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block PRIMARY_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    POS(X) = (#1) 
                 (Y)   (#2)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block PRIMARY_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+359.1696) +SIN[(-0.0042)* (POS(X)-(+61.0))]
                   (DEC)   (-34.6717 )      (+0.0042)  (   (Y) (+61.0)) 

Creating a contour plot of an image

In this section we are going to create a contour plot of the a4059_chandra.fits image, and edit the properties of the axes to better display the data. We finish the section by adding a labelled region, used to indicate the size scale of the image.

chips> make_figure("a4059_chandra.fits");
chips> info;
Window [win1]
  Frame [frm1]
    Plot [plot1]   (0.15,0.15)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Contour [ctr1]

chips> print (get_contour.levels);
[500.0, 1000.0, 1500.0, 2000.0, 2500.0]

which uses the make_figure() command to create a contour plot (Figure 1). We will show in later sections how the add_contour() command can also be used.

As we did not specify the contour levels, the command calculated values automatically that span the data range of the image. The levels being displayed can be found by inspecting the levels field of the object returned by get_contour(); in this case there are 5 contours, starting at 500 counts per pixel and extending up to 2500 counts per pixel.

The axis limits are automatically chosen to cover the region of the image that contains contours; there are a few high-count pixels to the North of the cluster - at a position close to (359.22,-34.705) - that cause the Y axis to be extended this far.

[Version: full-size, PNG, postscript]

Figure 1: Basic contour plot of a4059_chandra.fits

The make_figure command has created a contour plot of the image data in a4059_chandra.fits. The contour levels were calculated from the image, and the command chose to display the image using the WCS information for the axis ranges (both these defaults can be manually over-ridden).

The axis display is not ideal in this case, in that the labels are in decimal degrees, there are no labels, and the X axis is not displayed with the correct orientation. These will be addressed when creating Figure 2.

We now customize the plot to correct the orientation of the X axis and improve the labelling of the axes. The resulting figure is shown in Figure 2:

chips> reverse_axes(X_AXIS);

chips> set_xaxis({"tickformat","ra"});
chips> set_yaxis({"tickformat","dec"});

chips> set_plot_xlabel("RA");
chips> set_plot_ylabel("Declination");

chips> yax = ChipsAxis;
chips> yax.ticklabel.angle = 90;
chips> yax.ticklabel.offset = 20;
chips> yax.ticklabel.halign = 0.5;
chips> set_yaxis (yax);

[Version: full-size, PNG, postscript]

Figure 2: Creating a publication-quality plot

The axes from Figure 1 have been adjusted to produce a publication-quality plot: the X axis has been reversed, labels added, the numeric labels have been changed from degrees to sexagesimal format, and the ones on the Y axis have been rotated by 90 degrees to avoid overlapping the label and edge of the window.

We now add a region to the plot to act as a scale bar. We choose a width of 30 arcseconds, which corresponds to a maximum radius of 15 * sqrt(2) arcseconds for a square (the remaining factor of 3600 in the add_region call is to convert to degrees).

chips> add_region(4, 359.28, -34.73, 15*sqrt(2)/3600);

A label is added to indicate the width of the square. The move_label() command is used to move the label closer to the square (setting the third argument to 1 states that the label should be shifted relative to its old location), and the valign attribute is set to 0.5 to center the label vertically. The result is shown in Figure 3.

chips> add_label(359.27, -34.73, "30\"", {"color", "green"});
chips> move_label(0.004, 0, 1);
chips> print (get_label.valign);
0
chips> set_label({"valign", 0.5, "size", 14});

[Version: full-size, PNG, postscript]

Figure 3: Adding a region and label

The green square has a width of 30 arcseconds.

The pick() and get_pick() commands can be used to identify the positions of features in the contour: pick() is useful for a simple interactive session since it prints the results to the screen, whereas get_pick() returns the results as variables, so is of more use for scripted analysis. In the following we show how you can use both commands; first we set an undo tag so that we can easily return to the same state as Figure 3.

chips> set_undo_tag("pick");

Once pick() has been called you can move the mouse around within the ChIPS plot window and - depending on your window manager - see the current coordinates of the mouse in the window title. The values shown use the current coordinate system; in this case this is the data coordinate system of the image and so decimal degrees are used (the values are not formatted in sexagesimal format). For this example a point near the peak of the cluster emission was chosen. Since pick() was called with no arguments only one position was displayed before control was returned; if pick(2) had been used then two positions would have been selected, and pick(0) would have continued until the Escape key was pressed.

chips> pick;
(359.254,-34.7596)

We now use get_pick() to return an array of points; here we selected two corners of a rectangle surrounding the core emission of the cluster. The return value is then used to select new limits for the plot, with the result shown in Figure 4.

chips> c = get_pick(2);
chips> print(c);
[array([ 359.25998797,  359.24346934]), array([-34.754112  , -34.76977369]), array([2, 2], dtype=uint32)]
chips> x = c[0];
chips> y = c[1];
chips> limits(X_AXIS, x[0], x[1]);
chips> limits(Y_AXIS, y[0], y[1]);

The pick_limits() command automates this procedure - i.e. it replaces the need to call get_pick() and then limits() - and can also restrict the limits change to just the X or Y axis.

[Version: full-size, PNG, postscript]

Figure 4: Manually zooming into the contour

The contour plot shows the core emission from the cluster. The order of the axis limits passed to the limits() call does not matter in this case: both

chips> limits(X_AXIS, x[0], x[1]);

and

chips> limits(X_AXIS, x[1], x[0]);

would produce the same plot (i.e. the orientation of the axis is not determined by the order of the values sent to limits()).

We then revert these changes with the undo() command, using the previously-created tag of pick; the outcome of this command should be the same as Figure 3.

chips> undo("pick");

Overlaying a contour

In this section we shall overlay contours from the file a4059_nvss.fits onto Figure 3.

Our first attempt uses add_contour() to add a second contour plot. The resulting plot, shown in Figure 5, is not useful because the second image covers a larger area of the sky, and the contours chosen for it are not informative.

[Version: full-size, PNG, postscript]

Figure 5: A first attempt at overlaying radio contours

The axis limits have increased greatly since the radio image covers a much larger region of the sky than the X-ray image and the default contour levels go to far down into the noise. For reference, the levels are:

chips> print (get_contour.levels);
[0.0, 0.10000000149011612, 0.20000000298023224, 0.30000000447034836, 0.40000000596046448, 0.5000000074505806, 0.60000000894069672]

We could call undo() and repeat the call, using more sensible contour levels, but here we show how you can change the levels displayed by a contour, using the set_contour() call (the levels were chosen by trial and error based on the range of pixel values in the image):

chips> set_contour({"levels", [0.05, 0.1, 0.3, 0.5]});

We could manually change the axis limits, but instead we decide to use the limits derived from the first contour, using the object, id form of the limits() command. First we use the info() command to find out the name of the first contour (i.e. ctr1), and then use this along with the term chips_contour in the limits call. The resulting ranges do not quite match the limits used in Figure 3 since the limits call does not pad the values (e.g. the axis.pad setting) as the make_figure() and add_contour() calls do. We therefore choose to change the Y axis range so as not to include the small contours to the North of the cluster and to center the emission in the plot. The result is shown in Figure 6.

chips> info;
Window [win1]
  Frame [frm1]
    Plot [plot1]   (0.15,0.15)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Contour [ctr1]
      Region [reg1]
      Label [lbl1]
      Contour [ctr2]

chips> limits(chips_contour, "ctr1");
chips> print (get_plot_yrange);
[-34.778116881434684, -34.705309436694833]
chips> limits (Y_AXIS, -34.79, -34.72);

[Version: full-size, PNG, postscript]

Figure 6: Re-sizing the plot

The plot now shows contours of the X-ray and radio emission from Abell 4059, although it is not clear which is which.

To clarify the plot contents we decide to change the color of the radio contours to red and add labels to indicate the contour types. The labels are positioned using the "plot-normalized" coordinate system, where (0,0) and (1,1) refer to the bottom-left and top-right corners of the plot (i.e. where the X and Y axes meet).

chips> set_contour({"color", "red"});

chips> add_label(0.95, 0.1, "Chandra", {"coordsys", "plot_norm", "halign", 1});
chips> add_label(0.95, 0.15, "NVSS", {"coordsys", "plot_norm", "halign", 1, "color", "red"});

chips> save_state("contour-overlay.state");

The save_state() call creates a binary file which can be loaded into another ChIPS session (even on a different machine) by the load_state() command. The result will be the same plot, which can then be edited, changed, or printed. The ChIPS data tar file contains a copy of the state file for you to try out.

The final figure is shown in Figure 7.

[Version: full-size, PNG, postscript]

Figure 7: Cleaning up the plot

As the labels have been added using plot-normalized coordinates, they will not move if the axis limits are changed: for example after

chips> limits(Y_AXIS, -35, -34);

the labels remain in the bottom-right of the plot.

The contour plot can be compared to that created by ds9, shown in Figure 8.

[Version: full-size]

Figure 8: A ds9 version of the contour plot

The white contours are for the X-ray emission, the red contours are the radio emission, and the background image is the Chandra image.

The contour levels were manually set to those used to the values used by ChIPS, and the smoothing scale was set to 1. The contours used are available in the ChIPS tar file as chandra.con and nvss.con.

Controlling the contour levels

In this section we shall show how you can change the levels that are used to contour an image, as well as highlight other ways that contours can be created.

First we read in the image data using the read_file() command, which returns a Crate object representing the most-interesting block in this file (in this case the image). We use this object in the call to add_contour, and specify the levels coordinate system to use.

chips> img = read_file("a4059_chandra.fits");

chips> clear;
chips> levels = [400, 600, 800, 1200, 1800];
chips> add_contour(img, levels, {"wcs", "sky"});
chips> set_plot_xlabel("X");
chips> set_plot_ylabel("Y");

In order to center the cluster emission in the plot, and make the axes have the same length, we find the X and Y ranges of the contour, and then use limits() to set the width of each axis to be 500 sky pixels. The result is Figure 9.

chips> xr = get_contour_xrange;
chips> xr[1] - xr[0];
426.854675293

chips> yr = get_contour_yrange;
chips> yr[1] - yr[0];
728.061660767

chips> limits(X_AXIS, 3800, 4300);
chips> limits(Y_AXIS, 4000, 4500);

[Version: full-size, PNG, postscript]

Figure 9: Manually creating a contour plot

This contour was created by manually specifying the contour levels, over-riding the default choice of coordinates (so that the sky coordinates were used rather than the WCS ones), and the data used for the contour was taken from a Crate rather than a file name.

An image array could also have been used, as shown below, but then there would have been no coordinate transform, so the axes would have displayed the logical (i.e. pixel) coordinate system of the array:

chips> ivals = get_imagevals;
chips> ivals;
Integer_Type[64,65]
chips> add_contour(ivals, 65, 64);

Note that the dimensions of the image array need to be given in this form of the call, and that you give the size of the second dimension (the X axis) first.

The properties of a contour can be read with the get_contour() command and set with set_contour(). In the following we change the mode field to illustrate different ways that contour levels can be calculated. First we show the contour object:

chips> print(get_contour);
algorithm = 1
color = default
depth = 100
id = None
interval = 10.0
levels = [400.0, 600.0, 800.0, 1200.0, 1800.0]
mode = arbitrary
numlevels = 5
stem = None
style = solid
thickness = 1.0
wcs = sky

The mode is set to arbitrary since we manually set the levels when we created the plot. We start by changing the mode to nice, which attempts to use levels at human-readable values:

chips> set_contour({"mode", "nice"});
chips> get_contour.levels;
[500.0, 1000.0, 1500.0, 2000.0, 2500.0]

The levels field contains the currently-displayed levels, and can be accessed directly, as shown above. The resulting plot is shown in Figure 10.

[Version: full-size, PNG, postscript]

Figure 10: Using the nice mode for contour levels

The nice mode is the default mode for contour plots, so the levels here are the same as those shown in Figure 1. The axes are different since here we explicitly used the SKY coordinate system, rather than the celestial system that is picked by default.

We now change to the interval mode, and set the interval between contours to 200, which creates Figure 11.

chips> set_contour({"mode", "interval", "interval", 200});
chips> get_contour.levels;
[200.0, 400.0, 600.0, 800.0, 1000.0, 1200.0, 1400.0, 1600.0, 1800.0, 2000.0, 2200.0, 2400.0]

[Version: full-size, PNG, postscript]

Figure 11: Using the interval mode for contour levels

The interval mode creates levels that are separated by the contour's interval value (here set to 200), starting at the minimum value of the image (in this case 0).

chips> min(get_imagevals(img));
0

The count mode tries to create a fixed number of levels, given by the numlevels attribute of the contour. The actual number created depends on the data range, as the levels created are meant to be human-friendly in the same way that the nice mode levels are. This is why the following ends up creating a plot with only 2 contours (Figure 12):

chips> set_contour({"mode", "count", "numlevels", 4});
chips> get_contour.levels;
[1000.0, 2000.0]

[Version: full-size, PNG, postscript]

Figure 12: Using the count mode for contour levels

Using the count mode has resulted in a contour with two levels.

Finally we return to the arbitrary mode by specifying the levels to use (this automatically sets the mode field to arbitrary).

chips> set_contour({"levels", [500, 1000, 1400, 1800, 2200]});
chips> get_contour;
algorithm = 1
color = default
depth = 100
id = None
interval = 200.0
levels = [500.0, 1000.0, 1400.0, 1800.0, 2200.0]
mode = arbitrary
numlevels = 4
stem = None
style = solid
thickness = 1.0
wcs = sky

The plots in the next section - for example Figure 13 - use these levels.

Controlling the appearance of a contour plot

We start by by changing the line color and style of the contour plot from the previous section, using a dotted orange line. The result is Figure 13.

chips> set_contour({"color", "orange", "style", "dot"});

[Version: full-size, PNG, postscript]

Figure 13: Changing the line style for contours

The line styles available are the same as for other objects. Due to the way the patterns are constructed, the best results are generally found with the "solid" and "dot" styles.

The results can also depend on the thickness of the contour (here we have used the default value of 1). The on-screen display can look poor when the thickness is set to 0.5 on some graphics cards, but the hard-copy versions created by print_window() are displayed correctly in this case.

We finish by showing how you can create a contour of data in (x,y,z) form - in this case the contents of the text file projection.dat. We start by reading in the data from the file (as the first line of the file is a comment line with three words on it, these words are used as the column names):

chips> !head -5 projection.dat
# kTx Abundance Cash
4.0 1.0e-01 80.36740
4.550 1.0e-01 71.87300
5.100 1.0e-01 66.36260
5.650 1.0e-01 62.58350

chips> cr = read_file("projection.dat");
chips> print (get_col_names(cr));
['kTx' 'Abundance' 'Cash']
chips> kt = copy_colvals(cr, "ktx");
chips> abun = copy_colvals(cr, "abundance");
chips> sval = copy_colvals(cr, "cash");

We use split() to create a second plot, arranged below the first plot; the third argument gives the spacing between the plots. The set_preferences() call makes all future labels default to a size of 14, which saves having to add "size", 14 to the attribute list in the three add_label() calls below.

To identify the contour levels - we use the values returned by Sherpa when the projection was created - we create three contours of the same data, using color and labels to distinguish them. By default, ChIPS will create a name for each new object it makes - such as "plot1" and "cont1" - but this can be over-ridden by setting the id attribute when calling one of the add_<object>() calls. We use this to create meaningful names for the contour/label pairs (of 1sig, 2sig, and 3sig).

chips> split(2, 1, 0.1);
chips> set_preferences({"label.size", 14});

chips> add_contour(kt, abun, sval, [53.6454], {"id", "1sig"});
chips> add_label (7, 1, "1\sigma"R, {"id", "1sig"});

chips> add_contour(kt, abun, sval, [57.5304], {"id", "2sig", "color", "cyan"});
chips> add_label (9.5, 1.3, "2\sigma"R, {"id", "2sig", "color", "cyan"});

chips> add_contour(kt, abun, sval, [63.1794], {"id", "3sig", "color", "brown"});
chips> add_label (11.5, 1.5, "3\sigma"R, {"id", "3sig", "color", "brown"});

chips> set_plot_xlabel("kT_x");
chips> set_plot_ylabel("Abundance");

The info() command shows the full list of ChIPS objects that have been created in this visualization. We then change to the first plot, using the current_plot command, and add a title to it. If we had not changed plot then the title would have appeared above the second plot and overlap the "X" axis label from the first plot.

chips> info;
Window [win1]
  Frame [frm1]
    Plot [plot1]   (0.15,0.57)  .. (0.90,0.90)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Contour [ctr1]
    Plot [plot2]   (0.15,0.15)  .. (0.90,0.48)
      Border bottom [bx1]  top [bx2]  left [by1]  right [by2]
      X Axis [ax1]
      Y Axis [ay1]
      Contour [1sig]
      Label [1sig]
      Contour [2sig]
      Label [2sig]
      Contour [3sig]
      Label [3sig]

chips> current_plot("plot1");
chips> set_plot_title("Multiple contours");

chips> print_window("mplot", {"fittopage", 1});
chips> print_window("mplot", {"fittopage", 1, "format", "pdf"});

The print_window() calls create postscript and PDF versions of the plot; by setting the fittopage attribute to 1 the plot is expanded to fill the page (the keepaspect attribute has been left at its default value of 1 so that the aspect ratio of the plot is maintained by this expansion). The final plot is shown in Figure 14.

[Version: full-size, PNG, postscript, encapsulated postscript, PDF]

Figure 14: Multiple contour plots

The two plots each contain contoured data of an image; in the top plot there is one contour with multiple levels of a FITS image, and in the second there are three contours (one per level) of data read in from an ASCII file.

If we wanted to change the properties of the contours in the second plot, we could use the names we gave as id values to say something like the following (remembering to change the current plot first):

chips> current_plot("plot2");
chips> set_contour("2sig", {"color", green"});
chips> set_label("2sig", {"color", "green"});