Merging Central
Combining Datasets for Imaging & Spatial Analysis
This provides an overview of the considerations involved with combining observations for spatial/imaging analysis. For the most part, the most common tasks is automated using through the merge_obs/reproject_obs/flux_obs family of CIAO scripts. A broad discussion that began pre-dating these scripts and expanded upon their introduction is found in the Why combine data? 'why' page.
Astrometric Corrections (optional)
While updating the astrometry is generally unnecessary, if your science requires the best possible source positions, then you may want to apply astrometric offsets to individual observations to get the correct registration. Chandra's absolute pointing accuracy is generally better than 0.4 arcsec but can be improved in an absolute sense by cross-matching Chandra sources with other higher precision catalogs or in a relative sense by cross matching one Chandra observation with another observation. The CIAO tools reproject_aspect, wcs_match, and wcs_update can be used to compute the fine astrometic shifts between two source lists and apply the offsets to various Chandra files, but is also a major weakness, since there are more sophisticated algorithms available in the literature than what are presently implemented in CIAO.
-
The basics of astrometric corrections with Chandra data is covered in the Correcting Absolute Astrometry thread.
-
For a discussion on caveats and nuances involved with the updating astrometry with CIAO, see the Revising Astrometry for Imaging & Spatial Analysis page.
Reprojection:
When spatially combining data sets with similar sky coordinate definitions, the effects of the differing tangent planes will be difficult to detect; however, an [extreme] example of how things can go awry without reprojecting to a common sky coordinates grid is by looking at data sets near the celestial poles (e.g. observations within 10 deg of the south celestial pole). The combined data sets without reprojecting the observations to a common tangent plane, a single observation will define the WCS used to transform each observation's sky coordinates to celestial coordinates in the combined data sets, which is clearly incorrect as seen in Figure 1 .
TL;DR it's a fair assumption that the sky coordinates system and its respective WCS transform are unique for each observation.
Co-adding images within CIAO is typically performed in terms of the real-valued "sky" (also referred as "physical") coordinates by the tools.
Recall that the "sky" coordinate system records the event position on a fictitious tangent plane to a nominal, fixed celestial pointing direction. That is, the "sky" plane is the projection of the celestial sphere onto a 2D plane where the plane intersects orthogonally with the 3D unit sphere in the pointing direction.
For a typical observation the tangent point direction is defined by its nominal right ascension and declination which is generally the same as the telescope's aimpoint location—the time-averaged RA and Dec of the telescope's optical-axis due to spacecraft dithering—albeit, the nominal tangent point and aimpoint can differ if the Science Instrument Module is offset for an observation.
The nominal RA and Dec for the tangent point for an observation is then set to have the sky coordinate position to approximately be (4096.5,4096.5) for ACIS observations, (16384.5,16384.5) for HRC-I observations, and (32768.5,32768.5) for HRC-S observations. The default aimpoint is different each proposal cycle because of thermal relaxation of spacecraft structures.
A practical consequence of the sky coordinate system being specific to a given observation is that the WCS transform between sky and celestial coordinates will also be specific to that observation. This means that to correctly stack observations for imaging analysis, the observations should all share a common sky coordinate system and WCS transform, which is achieved through reprojecting the data sets to a common spherical projection onto a plane.
In this example, we use 15 ACIS observations with aimpoint declinations, \({-80}^{\circ} \geq \delta > {-90}^{\circ}\), so they should be unevenly distributed about the pole at various declinations. The most simplistic approach to combining files is to use dmmerge on a set of FITS tables, e.g. event files with a common set of columns.
We can perform the query for the observations in ChaSeR or with find_chandra_obsid, which can also automatically download the datasets, which we then proceed to reprocess with a common CalDB version using chandra_repro, while neglecting any consideration about the chosen ACIS data mode for the purposes of this example.
unix% find_chandra_obsid arg=0 dec=-89.9999 radius=600 instrument=acis detail=basic download=none verbose=1 # obsid sepn inst grat time obsdate piname target 10365 420.4 ACIS-S NONE 5.6 2009-11-18 Marshall "PKS 1303-827" 3477 378.0 ACIS-S NONE 19.8 2002-03-31 Fox GRB020321 8266 384.4 ACIS-I NONE 8.0 2007-06-24 Murray "RXC J1539.5-8335" 8272 465.8 ACIS-I NONE 7.9 2007-03-22 Murray "ACO S 405" 19534 515.9 ACIS-S NONE 26.7 2018-02-06 Lansbury NuSTARJ151253-8124.3 20959 515.9 ACIS-S NONE 20.8 2019-01-04 Lansbury NuSTARJ151253-8124.3 21421 202.1 ACIS-S NONE 21.9 2019-02-18 Elvis ESO005-G004 10143 261.2 ACIS-S NONE 2.0 2008-12-10 Fox 1RXSJ200924.1-853911 16283 466.7 ACIS-I NONE 7.4 2014-09-17 Murray "ACO S 405" 15124 289.4 ACIS-I NONE 22.8 2013-12-10 Jones "ACO 4023" 14925 585.1 ACIS-S NONE 38.5 2013-08-27 Wolter "AM 0642-801" 15156 584.0 ACIS-I NONE 8.0 2013-06-08 Reiprich A2837 18296 536.6 ACIS-I NONE 12.6 2016-09-02 Rossetti "PSZ2 G295.25-21.55" 22293 571.9 ACIS-S NONE 19.7 2019-12-11 Wolk "HD 39091" 22326 267.6 ACIS-I NONE 37.4 2020-12-18 Safdi "RE J0317-853" unix% find_chandra_obsid arg=0 dec=-89.9999 radius=600 instrument=acis download=all... lots of screen output ...unix% chandra_repro "*" outdir=""... lots of screen output ...
Table of ACIS ObsIDs within 10 degrees of the southern celestial pole.
ObsID | RA | Declination |
---|---|---|
15156 | 00:52:44.90 | -80:15:57.60 |
22326 | 03:17:15.85 | -85:32:25.56 |
8272 | 03:51:28.00 | -82:14:11.00 |
16283 | 03:51:31.70 | -82:13:20.00 |
22293 | 05:37:09.89 | -80:28:08.83 |
21421 | 06:05:41.60 | -86:37:55.00 |
14925 | 06:38:33.60 | -80:14:52.00 |
18296 | 09:19:02.50 | -81:03:25.20 |
10365 | 13:08:38.00 | -82:59:34.80 |
19534 | 15:12:54.00 | -81:24:06.00 |
20959 | 15:12:54.00 | -81:24:06.00 |
8266 | 15:39:25.20 | -83:35:34.00 |
3477 | 16:11:02.40 | -83:42:00.00 |
10143 | 20:09:13.00 | -85:38:46.80 |
15124 | 23:39:44.00 | -85:10:33.00 |
dmmerge-ing these event files will result in will result in using the sky to WCS transform defined in the first input file to the tool.
An example of directly using the reprocessed data sets and merging event files that have not been projected to a common sky plane.
unix% ls */repro/*evt* > evt.lis unix% dmmerge infile="@evt.lis[cols time,expno,ccd_id,node_id,chip,\ ? tdet,det,sky,pha,energy,fltgrade,grade,status]" \ ? outfile="merge_no-reproj.evt" clobber+
[Version: full-size]
Figure 1: Stacked Observations without Reprojecting
So you'll need to redefine the sky coordinates of all these observations by reprojecting the events to a common tangent plan, which will then be displayed as the set of obervations unevenly scattered around the celestial pole at varying hour angles and declinations.
Stacking the observations to a common sky tangent point centered near the sourthern celestial pole returns the expected mosaic about the pole.
unix% merge_obs infiles="*/repro/*evt*.fits" \ ? outroot="merge_reproj_Dec-90/" bands=broad binsize="16" \ ? refcoord="0:0:0 -89:59:59.999999" clobber+
[Version: full-size]
Figure 2: Stacked, Reprojected Observations
Source Detection with Merged PSF maps:
A primary reason users typically [spatially] stack observations is to detect faint sources. The CIAO celldetect and wavdetect tools can make use of information about the PSF size at a given point to either determine the size of a detect cell at that point or determine the optimum wavelet size that matches the size of the sources being detected (for point sources, this scale is comparable to the size of the PSF), for the respective tool using their psffile parameter with PSF maps.
A discussion may be found in the Running wavdetect or celldetect on Merged Data: Choosing psffile thread and a brief overview of the strengths and weaknesses between the two available detection method (vtpdetect, a Voronoi tessellation and percolation source detection technique does not make use of PSF maps) in CIAO are summarized below.
celldetect
Pros
- Fast and Robust
- Works well for point sources
- Detailed PSF shape not needed, only approximate size
Cons
- Extended sources can be difficult, without careful selection of cell sizes
- Can get confused in crowded fields
- Exposure maps needed to eliminate edge sources
- Not very sensitive, unless background maps are used, but these can be difficult to make
wavdetect
Pros
- Works well in crowded fields
- Works well for point sources superimposed on extended emission
- Only approximate PSF shape is needed
- Edge of field effects not a problem
Cons
- Slow, especially if many wavelet scales used
- Memory-intensive and very large image sizes can pose a problem
In wavdetect, providing an exposure map will help with eliminating the vast majority of false detection along the chip edges. The most notable feature of the above thread is the discussion regarding the various methods to combine PSF maps.
The pixel values of the PSF map represents the radius of a circular PSF at that pixel for the given energy and ecf parameter values passed to mkpsfmap in the units that are also specified in the tool assuming a flat detector. However, realize that the circular assumption used will completely fall apart at large off-axis angles, just based on the changes to the PSF morphology and size.
Scattered throughout the existing documentation related to generating PSF maps (with mkpsfmap and dependent scripts such as fluximage/merge_obs/flux_obs) there are instances where a PSF 'encircled counts fraction' of ecf=0.393 is used to generate PSF maps. An ECF of 39.3% corresponds to the 1\(\sigma\) integrated volume of a 2D Gaussian and in the examples, it is chosen so that point sources with only a few counts which will be concentrated in the center of the PSF to be better detected.
The discussion will be focused on wavdetect herein since it is the most widely used CIAO source detection tool. While the tool cannot compute point spread function sizes if wavdetect is not provided with a PSF map, it will still run. The source characteristics will not be reliable and if the size is wrong, then the wrong wavelet scales may be used to determine the source parameters, and thus the source properties may be wrongly estimated. However, the source detection process is unaffected. That is, when a PSF map is provided, the PSF sizes are used by the tool to determine the detected source significance and whether or not it is included in the source list, thus limiting the number of spurious detections.
wavdetect is basically a wrapper script around the wtransform and wrecon tools. Source detection is performed with wtransform, which provides a list of source candidates and wrecon tests these candidates for significance, providing wavdetect with a source list file of significant detections.
In the detection step, wtransform does not use the PSF information but uses a Ricker (aka "Mexican Hat") wavelet function at different size-scales to correlate the input image with the wavelets. Pixels with a high-correlation are considered sources. The tool then uses an estimated background and the selected detection threshold (sigthresh) to determine whether or not the sources from the correlated pixels can then be considered candidate source detections.
In the next step the PSF map is important. Since wrecon will take the candidate sources and smooth them at various wavelet scales. The source is only statistically significant—and kept as a source—when smoothed with a wavelet scale greater than the PSF size and is reconstructed using the wavelet scale closest to the PSF size. The source size and position are also better estimated during this process. Without a PSF map, the wavelets are scaled with the the image pixels, and only reliable in the case where the PSF size does not vary in the region being examined.
As an aside—and point of clarification that is frequently asked—pertaining to the ellsigma parameter found in both wavdetect and celldetect is that it has no relationship to mkpsfmap's ecf value used to generate PSF maps and it does not affect the detections themselves. ellsigma is an optional parameter with a default value of 3 standard deviations of the distribution of counts within the source cell (the observed counts), with the distribution assumed to be Gaussian. For example, if "ellsigma" is 1, then the elliptical source region will contain 39.3% of the total observed source counts but is only used as a scaling factor on the major- and minor-axes of the ellipses for each detection to enhance the visualization of the sources regions when the source list is overlaid onto the data and does not affect the source detection process or other determined source properties.
Ultimately, providing a PSF will affect the wavdetect/celldetect results, giving a shorter source list with better estimated source properties. While the PSF size information does not affect the actual detection process, since it is not used, but the PSF map definitely affects the final source list, as it is used for the significance test to reject false detections.
CAVEATS
caveats: merging exposure maps/exposure-corrected images... the merge_obs approach was the simplest