Last modified: 13 Jan 2022

URL: https://cxc.cfa.harvard.edu/ciao/threads/dm_cube/

Using Data Cubes

CIAO 4.16 Science Threads


Overview

Synopsis:

The CIAO Data Model allows you to filter and manipulate 3-dimensional images, known as "data cubes." You can select 3-dimensional subsets, or slice out 2-dimensional pieces. If two dimensions of the cube represent a pair of position axes, you can apply a region filter to those axes.

Purpose:

To create and filter data cubes.

Last Update: 13 Jan 2022 - Reviewed for CIAO 4.14. Minor tweaks only.


Contents


Get Started

Download the sample data: 1463 (ACIS-S, Jupiter)

unix% download_chandra_obsid 1463 evt2

The thread also uses a VLA radio image from the NRAO Archive (project AZ0128).


Creating a Position-time Data Cube from an Event File

While we choose to create an (x,y,time) data cube, users may bin on any three columns that make sense in the analysis. For instance, you may want to create a PHA or energy axis to see how the spectral characteristics of a source change over time.

A common analysis involving data cubes is to create a file with two position axes and one time axis. This example shows how one might choose the binning parameters to create such a cube by inspecting various possibilities via dmlist before writing out the file.

Here we are using an observation of jupiter - acisf01463N006_evt2.fits. We need to determine a suitable time range and step size, and select a spatial range for the filtering.

What if we were to simply bin all three axes by a factor of one? dmcopy can be used to examine the effects by creating a virtual output file:

unix% dmlist "acisf01463N006_evt2.fits[bin x,y,time]" cols
# DMLIST (CIAO): WARNING: Creating large image: 1745 MB. Current max set at 500 MB.
Increase maximum using [opt mem=n] or increase blocking to reduce size.
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[8192,8192,22378]              Int2(8192x8192x22378) -                    
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (#1) [pixel]
                 (y)   (#2)
   2   3      time                 = +59969682.285558 [s] +1.0 * (#3  -0.50)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6834)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6789 )           (+0.000136667)  (   (y) (+4096.50)) 

The result is an 8192 x 8192 x 22378 cube, which would about 3 terabytes in size.

The time axes is 22378 pixels because the default pixel size is one second and this is a 22 ks observation. We can shorten this axis to only 22 steps in time by binning in units of 1000 seconds. Additionally, the spatial size can be reduced by binning x and y by a factor of 8:

unix% dmlist "acisf01463N006_evt2.fits[bin x=::8,y=::8,time=::1000]" cols
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[1024,1024,23]              Int2(1024x1024x23) -                    
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (+0.50)[pixel] +(+8.0)* ((#1)-(+0.50))
                 (y)   (+0.50)         (+8.0)  ((#2) (+0.50))
   2   3      time                 = +59969682.285558 [s] +1000.0 * (#3  -0.50)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6834)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6789 )           (+0.000136667)  (   (y) (+4096.50)) 

The output is now a 1024 x 1024 x 23 cube, which is more reasonable.

An alternate approach is to keep the full spatial resolution, but use a small region of the file:

unix% dmlist "acisf01463N006_evt2.fits[bin x=3900:4400,y=4100:4600,time=::1000]" cols
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[500,500,23]              Int2(500x500x23) -                    
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (+3900.0)[pixel] +(+1.0)* ((#1)-(+0.50))
                 (y)   (+4100.0)         (+1.0)  ((#2) (+0.50))
   2   3      time                 = +59969682.285558 [s] +1000.0 * (#3  -0.50)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6834)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6789 )           (+0.000136667)  (   (y) (+4096.50)) 

which shows that we will get a 500 x 500 x 23 image.

Finally, we decide to use this binning with dmcopy to create the file:

unix% dmcopy "acisf01463N006_evt2.fits[bin x=3900:4400,y=4100:4600,time=::1000]" jupiter_cube.fits

unix% dmlist jupiter_cube.fits blocks
 
--------------------------------------------------------------------------------
Dataset: jupiter_cube.fits
--------------------------------------------------------------------------------
 
     Block Name                          Type         Dimensions
--------------------------------------------------------------------------------
Block    1: EVENTS_IMAGE                   Image      Int2(500x500x23)
Block    2: GTI7                           Table         2 cols x 4        rows
Block    3: GTI0                           Table         2 cols x 2        rows
Block    4: GTI1                           Table         2 cols x 2        rows
Block    5: GTI2                           Table         2 cols x 3        rows
Block    6: GTI3                           Table         2 cols x 5        rows

As expected, the output is a 500 x 500 x 23 image.

How to display the file

SAOImage ds9, the default imager distributed with CIAO, has the capability to display data cubes.

unix% ds9 jupiter_cube.fits &

When the file is opened, ds9 automatically detects that it is a cube and launches the data cube dialog box, as shown in Figure 1. (If the data cube dialog box doesn't launch, open it from the "Frame" menu.) The spatial axes are displayed, while the third axis - time, in this case - is accessible via the dialog box. When "Play" is chosen, ds9 cycles through the bins of the time axis, essentially creating a movie of the object. The speed of the frame changes is controlled from the "Interval" menu of the dialog box.

Figure 1: Viewing a data cube

[Thumbnail image: One interval of the data cube is visible at a time when displayed in ds9. The data cube control window is also open.]

[Version: full-size]

[Print media version: One interval of the data cube is visible at a time when displayed in ds9. The data cube control window is also open.]

Figure 1: Viewing a data cube

This data cube has 23 intervals. Choosing "Play" from the data cube window steps through the intervals.

The other buttons in the dialog box (e.g. "Prev" and "Next") allow the user to move back and forth manually as well.

ds9 also has the option to render 3D datasets, in 3D.

unix% ds9 -3d jupiter_cube.fits &

Figure 2: 3-D rendering menu

[Thumbnail image: 3D render menu]

[Version: full-size]

[Print media version: 3D render menu]

Figure 2: 3-D rendering menu

A new 3D menu is now available that allows the user to rotate the image up-and-down (Elevation) as well as left-to-right (Azimuth). The current 3D slice is highlighted in cyan.

The data cube can be rotated using the sliders or by typing directly in the text box at the end of each slider as is shown in Figure 3

Figure 3: 3-D rendering example

[Thumbnail image: ds9 now supports real-time 3D rendering of data cubes. In this sequence of screen shots, the image has been rotated from -60 to +75 degrees.]

[Version: full-size]

[Print media version: ds9 now supports real-time 3D rendering of data cubes. In this sequence of screen shots, the image has been rotated from -60 to +75 degrees.]

Figure 3: 3-D rendering example

ds9 now supports real-time 3D rendering of data cubes. In this sequence of screen shots, the image has been rotated from -60 (top left) to +75 degrees (bottom right).


Manipulating the Data Cube

The data cube can be filtered using DM syntax in the same way as 2D files. All of these examples use the data cube created in the previous section.

Range filtering

A filter is applied to the time column of the cube to select a range of 3500 s:

unix% dmcopy "jupiter_cube.fits[time=59969680:59973180]" range_cube.fits

unix% dmlist range_cube.fits cols
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[500,500,4]              Int2(500x500x4) -                    
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (+3900.0) +(+1.0)* ((#1)-(+0.50))
                 (y)   (+4100.0)  (+1.0)  ((#2) (+0.50))
   2   3      time                 = +59969682.285558 [s] +1000.0 * (#3  -0.50)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6834)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6789 )           (+0.000136667)  (   (y) (+4096.50)) 

The output file is 500 x 500 x 4, since the filter spanned four of the 1000 s time bins. Displaying the file in ds9 (Figure 4) looks similar to the unfiltered file; there are just fewer steps available in the data cube dialog box.

Figure 4: Time-filtered data cube

[Thumbnail image: One interval of the data cube is visible at a time when displayed in ds9. The data cube control window is also open.]

[Version: full-size]

[Print media version: One interval of the data cube is visible at a time when displayed in ds9. The data cube control window is also open.]

Figure 4: Time-filtered data cube

This data cube has 4 intervals.


Region filtering

One may decide to apply a region filter to restrict the spatial axes of the file further. The region file used here, shown on the data in Figure 5 is:

unix% cat circle.reg 
# Region file format: CIAO version 1.0
circle(4143.5,4266.5,137.16489)

Note that if you are working with an object that moves in time, such as the solar system object we're using, make sure the region is large enough that the object won't drift out of the field of view. Define the region, then use the data cube dialog box to step forward and confirm that the object will still lie within the region.

Figure 5: Region filter overlaid on the data

[Thumbnail image: A green circle is defined on the data with the center to the west (right) of Jupiter.]

[Version: full-size]

[Print media version: A green circle is defined on the data with the center to the west (right) of Jupiter.]

Figure 5: Region filter overlaid on the data

The region is large enough that Jupiter won't drift out of it over the time interval.

dmcopy is used to apply the region filter to the unfiltered data cube:

unix% dmcopy "jupiter_cube.fits[sky=region(circle.reg)][opt null=-99]" region_cube.fits

Note that if you don't want to use a region file, the equivalent syntax for defining the region on the command line is:

unix% dmcopy "jupiter_cube.fits[sky=circle(4143.5,4266.5,137.16489)][opt null=-99]" region_cube_2.fits

The resulting file has the same time axis as the unfiltered cube, but smaller spatial axes:

unix% dmlist region_cube.fits cols
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range            Null
   1   EVENTS_IMAGE[275,275,23]              Int2(275x275x23) -                    -99        
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (+4006.0) +(+1.0)* ((#1)-(+0.50))
                 (y)   (+4129.0)  (+1.0)  ((#2) (+0.50))
   2   3      time                 = +59969682.285558 [s] +1000.0 * (#3  -0.50)
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6834)[deg] +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6789 )           (+0.000136667)  (   (y) (+4096.50)) 

It may be helpful to think of this file as a cylinder: a stack of (x,y) circles with a height of the time axis. Figure 6 shows the file in ds9.

Figure 6: Region-filtered data cube

[Thumbnail image: One interval of the spatially-filtered cube looks like a circle of data in ds9. The data cube control window is also open.]

[Version: full-size]

[Print media version: One interval of the spatially-filtered cube looks like a circle of data in ds9. The data cube control window is also open.]

Figure 6: Region-filtered data cube

The spatial axes have been filtered by the region file, but the time axis is the same (23 intervals).


Slicing a position-time data cube

It is also possible to slice time planes out of the file:

unix% dmcopy "jupiter_cube.fits[#3=10][bin x,y]" plane10.fits

This filters the cube by selecting pixels where the logical axis 3 ("#3") coordinate is equal to 10.

The output file in an (x,y) image representing that slice of time:

unix% dmlist plane10.fits cols
 
--------------------------------------------------------------------------------
Columns for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   EVENTS_IMAGE[500,500]              Int2(500x500)  -                    
 
--------------------------------------------------------------------------------
Physical Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    sky(x) = (+3900.0) +(+1.0)* ((#1)-(+0.50))
                 (y)   (+4100.0)  (+1.0)  ((#2) (+0.50))
 
--------------------------------------------------------------------------------
World Coordinate Axis Transforms for Image Block EVENTS_IMAGE
--------------------------------------------------------------------------------
 
Group# Axis# 
   1   1,2    EQPOS(RA ) = (+24.6832) +TAN[(-0.000136667)* (sky(x)-(+4096.50))]
                   (DEC)   (+8.6784 )      (+0.000136667)  (   (y) (+4096.50)) 

This is a two-dimensional image, shown in Figure 7.

Figure 7: A time plane from the data cube

[Thumbnail image: Displaying the time-filtered event file looks the same as viewing a single interval of the cube.]

[Version: full-size]

[Print media version: Displaying the time-filtered event file looks the same as viewing a single interval of the cube.]

Figure 7: A time plane from the data cube

The time axis has been filtered to one interval and the spatial axes were binned by a factor of 1 into a 500x500 image.


Removing Extra Axes: a 4D image that's really a 2D image

Some data, particularly radio images, have additional coordinate axes that are only one pixel wide, which are used to convey extra metadata. While useful, these extra axes may confuse applications that are designed for 2-dimensional data. dmcopy can be used to strip away the extra axes (at the cost of losing those metadata).

For instance, data obtained with the VLA looks like:

unix% dmlist vla_radio.img cols
 
-----------------------------------------------------------------------
Columns for Image Block PRIMARY
-----------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   PRIMARY[256,256,1,1] JY/BEAM      Real4(256x256x1x1) -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 ) = (+265.6222) +SIN[(-0.000277778)* (POS(X)-(+128.0))]
                   (DEC)   (-28.9884 )      (+0.000277778)  (   (Y) (+129.0)) 
   2   3      FREQ                 = +4.33399E+10  +100000000.0 * (Z  -1.0)
   3   4      STOKES               = AXIS4 

This is a 256 x 256 x 1 x 1 image whose physical axes are known to CIAO as X, Y, Z, AXIS4 and whose world coordinate axes are called RA, DEC, FREQ and STOKES. If we want to make this into a simple RA, DEC image:

unix% dmcopy "vla_radio.img[bin x,y]" vla_ra_dec.img

The file vla_ra_dec.img contains only the axes which were included in the binning specification; the additional information from the input file is discarded. Since we binned by a factor of one, the output image is also 256 x 256.

unix% dmlist vla_ra_dec.img cols
 
--------------------------------------------------------------------------------
Columns for Image Block PRIMARY_IMAGE
--------------------------------------------------------------------------------
 
ColNo  Name                 Unit        Type             Range
   1   PRIMARY_IMAGE[256,256]              Real4(256x256) -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 ) = (+265.6222) +SIN[(-0.000277778)* (POS(X)-(+128.0))]
                   (DEC)   (-28.9884 )      (+0.000277778)  (   (Y) (+129.0)) 

Figure 8 shows the VLA image in ds9.

Figure 8: VLA image

[Thumbnail image: The VLA image is displayed in (RA,Dec)]

[Version: full-size]

[Print media version: The VLA image is displayed in (RA,Dec)]

Figure 8: VLA image

The image was binned by a factor of one, so the output is still 256 x 256.


Final Note: Responses Files

While it is possible with CIAO tools to create 3d and higher-dimensional data cubes, it is not possible to create arbitrary dimension response files. For example when the 3D spatial-time image is created in this example:

unix% dmcopy "acisf01463N006_evt2.fits[bin x=3900:4400,y=4100:4600,time=::1000]" jupiter_cube.fits

It is not possible (or at least not easy) to create a 3D exposure map to accompany this image that gives the amount of good-time per time-slice. Having valid response information may not be critical for some types of analysis or for some observations but it may be something to consider when dealing with non-traditional data products.


History

27 Jan 2006 new for CIAO 3.3: original version
01 Dec 2006 updated for CIAO 3.4: CIAO version
08 Jan 2008 updated for CIAO 4.0: changed filename in examples from jupiter.fits to acisf01463N002_evt2.fits
14 Jan 2009 updated for CIAO 4.1: images are inline
05 Feb 2010 reviewed for CIAO 4.2: no changes
13 Jan 2011 reviewed for CIAO 4.3: no changes
11 Jan 2012 reviewed for CIAO: updated filename to acisf01463N006_evt2.fits (version N006); minor changes in screen output due to the data reprocessing
13 Dec 2012 Review for CIAO 4.5; added example of 3D data cube rendering.
25 Nov 2013 Review for CIAO 4.6. No Changes.
07 Apr 2014 Added a final note about response files.
17 Dec 2014 Reviewed for CIAO 4.7; no changes.
01 Feb 2016 Updated ds9 links.
03 Jan 2017 Reviewed for CIAO 4.9; no changes.
13 Jan 2022 Reviewed for CIAO 4.14. Minor tweaks only.