Fitting FITS Image Data

Sherpa Threads (CIAO 4.0 Beta)

[Python Syntax]

Overview

Last Update: 15 May 2008 - rewritten for Sherpa Beta

Synopsis:

This thread models image data of the supernova remnant G21.5-0.9 using a source model expression that involves multiple model components. The Create and Read Fits Image Data section shows how the input file was created.

Proceed to the HTML or hardcopy (PDF: A4 | letter) version of the thread.

Create and Read Fits Image Data

Sample ObsID used: 1838 (ACIS-S, G21.5-0.9)

File types needed: evt2

First, an image of the region of interest is created:

unix% punlearn dmcopy
unix% dmcopy \
      "acisf01838N001_evt2.fits[sky=box(4059.25,4235.75,521.5,431.5)][bin x=::2,y=::2]" \
      image2.fits

Because we plan to model the extended structure, not detailed features in the image, we can use coarser binning and reduce the time taken for the fit. A binning factor of 2 is used.

The choice of binning should be determined by the quality of the data and the models you wish to fit. Increasing the binning factor when creating the image is likely to lose small-scale information.

After starting a Sherpa session, the dataset is input using the load_image() command. The image is then displayed by the image_data() command:

sherpa> load_image("image2.fits")
sherpa> image_data()

which creates Figure 1 . The ds9 imager is used by default; see the Using SAOImage ds9 CIAO thread for more information.

Setting the coordinate system

The default coordinate system for fitting images is that of logical (also called image) which uses the FITS convension for numbering pixels: the first pixel is centered at (1,1) with the bottom-left corner being at (0.5,0.5) and the top-right corner being at (1.5,1.5).

Although a simple system, it can be difficult when converting the best-fit parameters, such as source location, back to more-meaningful coordinate systems. For this reason we choose to use the physical coordinate system for fitting. This is done with the set_coord() function:

sherpa> set_coord("physical")

Choosing the appropriate statistics

When fitting Chandra or XMM-Newton imaging data it is almost certain that many pixels will only contain a small number of counts. For this reason, we use the Cash statistic to fit the source model. This means that we can not subtract a background dataset but either have to fit it - as we do below - or ignore the background signal.

sherpa> set_stat("cash")

The CSTAT statistic may also be used in fitting image data. CSTAT is equivalent to the XSPEC implementation of the Cash statistic. The advantage of CSTAT over Sherpa's implementation of CASH is that one can assign an approximate goodness-of-fit measure to a given value of the CSTAT statistic.

Setting the fit method

The default fitting method (Levenberg-Marquardt) is not well suited for use with the Cash statistic, so we choose to use the NelderMead optimizer instead.

sherpa> set_method("neldermead")

The fit should be tried with a variety of optimization methods to ensure that the global best-fit has been found, rather than a local one.

The current Sherpa status may be reviewed:

sherpa> show_all()
Optimization Method: NelderMead
Statistic:           Cash

Data Set: 1
name      = image2.fits
x0        = Float64[56376]
x1        = Float64[56376]
y         = Int16[56376]
shape     = (216, 261)
staterror = None
syserror  = None
sky       = <pytransform.WCSTANTransform; proxy of C++ WCSTANTransform instance at 0x815ec70>
eqpos     = <pytransform.WCSTANTransform; proxy of C++ WCSTANTransform instance at 0x8335ae8>
coord     = physical

Filtering Image Data

Here we illustrate several ways of filtering image data within Sherpa.

Define the region to filter

After the data has been displayed, the filter region may be determined from the imager. Alternatively, an existing region file may be used.

Use the "Region -> Shape" menu in ds9 to choose a region shape, then left-click on the display to draw the desired shape. the Region file format should be set to "CIAO" and the coordinate system to "Physical". For further instructions on how to create regions in ds9, see the Using CIAO Region Files thread.

After the desired region size is set, as shown in Figure 2 , use the "Region -> List Regions" menu to get the region definition:

# Region file format: CIAO version 1.0
circle(4072.46,4249.34,108)

The region may also be save to a file; see the ds9 thread for instructions.

Define the filter on the command line

The filter is now set on the command line with the region definition:

sherpa> notice2d("circle(4072.46,4249.34,108)")

One may also specify a file which contains the region definition:

sherpa> notice2d("filter.reg")

Using the image_data() command again displays the filtered data, as shown in Figure 3 .

Defining a Source Model

We wish to fit the image using a source model expression that involves multiple model components. By the end of the thread we will use two two-dimensional Gaussian functions together with a constant, but for now we begin with only one Gaussian. Since this is an example these choices are not physically motivated but just an empirical choice to describe the emission.

sherpa> set_model(gauss2d.g1)

It may be necessary to modify the starting parameter values and parameter limits before performing the fit. For complex models, physical knowledge of the problem can guide the initial settings of parameters and their ranges and improve the final convergence of the fit; note that there is no guess function in the Beta version of Sherpa.

After setting the values, the ellipticity and its angle (theta) are frozen so they do not vary during the fit.

sherpa> g1.ampl = 20
sherpa> g1.fwhm = 20
sherpa> g1.xpos = 4065.5
sherpa> g1.ypos = 4250.5
sherpa> freeze(g1.ellip)
sherpa> freeze(g1.theta)

The upper limit of the fwhm, xpos, and ypos parameters is constrained to roughly the size of the image in physical coordinates (300). ampl is restricted to the range 1:1000.

sherpa> g1.fwhm.max = 4300
sherpa> g1.xpos.max = 4300
sherpa> g1.ypos.max = 4300
sherpa> g1.ampl.min = 1
sherpa> g1.ampl.max = 1000

The current source model definition may be displayed:

sherpa> print(get_model())
g1
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   g1.fwhm      thawed           20       1e-120          4300
   g1.xpos      thawed       4065.5      -1e+120          4300
   g1.ypos      thawed       4250.5      -1e+120          4300
   g1.ellip     frozen            0            0        0.999
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      thawed           20            1         1000
s

The image_fit() command displays a montage of the data, current model, and residuals that can be used to see how well the model describes the data.

We decide to include a background component before fitting, and so add in a constant model to the source expression. The value of the background model (0.2) was estimated from a radial profile. The upper limit is set to 100 as a reasonable constraint.

sherpa> set_model(gauss2d.g1+const2d.bgnd)
sherpa> bgnd.c0 = 0.2
sherpa> bgnd.c0.max = 100

The model is now:

sherpa> print(get_model())
(g1 + bgnd)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   g1.fwhm      thawed           20       1e-120          4300
   g1.xpos      thawed       4065.5      -1e+120          4300
   g1.ypos      thawed       4250.5      -1e+120          4300
   g1.ellip     frozen            0            0        0.999
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      thawed           20            1         1000
   bgnd.c0      thawed          0.2            0          100

Fitting the model

We are now ready to fit the data:

sherpa> fit()
Statistic value = -48907.8 at function evaluation 449
Data points = 9171
Degrees of freedom = 9166
   g1.fwhm        57.9477
   g1.xpos        4070.4
   g1.ypos        4251.11
   g1.ampl        23.3562
   bgnd.c0        0.266363

The best-fit values for the gaussian FWHM and amplitude and the background level are close to the values we guessed previously. Since we are using the PHYSICAL coordinate system the xpos, ypos, and fwhm parameters are given in the physical system. For this ACIS image this is the SKY coordinate system, which means the FWHM corresponds to 57.9 * 0.492 = 28.5 arcseconds.

The fit can be visualised using image_fit(), which displays the data, model, and residual image. Here we are interested in just the residuals:

sherpa> image_resid()

The resulting image is shown in Figure 4 . Although the gaussian does a reasonable job at describing the radial profile of the large-scale emission, the residuals do appear to be spatially correlated. This suggests that the model used here is not sufficient to describe all the structure in the source; however it will suffice for this example.

Calculating errors on the fit parameters

Note that the statistic value is negative because we are using the Cash statistic - which is a maximum-likelihood statistic - rather than a chi-square statistic. This means that we cannot get an "absolute" measure of the fit - i.e. how well the model actually fits the data. However, we can still use the projection() command to estimate errors on these parameters:

sherpa> projection()
projection 1-sigma bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g1.fwhm           57.9477    -0.283035     0.283835
   g1.xpos            4070.4    -0.180546     0.180499
   g1.ypos           4251.11    -0.178588     0.178737
   g1.ampl           23.3562    -0.243317     0.245693
   bgnd.c0          0.266363  -0.00930696   0.00946291

Adding an extra component

The residual image and profile of the fit show excess emission at the center of the source. We will add another gaussian component to try and account for this emission. The upper parameter limits used for g1 are also used for g2.

sherpa> freeze(g1)

sherpa> set_model(gauss2d.g2+g1+bgnd)
sherpa> g2.fwhm = 10
sherpa> g2.ampl = 100
sherpa> g2.xpos = 4065.5
sherpa> g2.ypos = 4250.5
sherpa> g2.fwhm.max=4300
sherpa> g2.xpos.max=4300
sherpa> g2.ypos.max=4300
sherpa> g2.ampl.min=1
sherpa> g2.ampl.max=1000
sherpa> freeze(g2.ellip)
sherpa> freeze(g2.theta)

sherpa> fit()
Statistic value = -53918.2 at function evaluation 657
Data points = 9171
Degrees of freedom = 9166
   g2.fwhm        6.08638
   g2.xpos        4070.78
   g2.ypos        4249.34
   g2.ampl        228.719
   bgnd.c0        0.260454

Since the g1 model was frozen, those values did not change in the fit:

sherpa> print(get_model())
((gauss2d[g2] + gauss2d[g1]) + const2d[bgnd])
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   g2.fwhm      thawed      6.08638       1e-120         4300
   g2.xpos      thawed      4070.78      -1e+120         4300
   g2.ypos      thawed      4249.34      -1e+120         4300
   g2.ellip     frozen            0            0        0.999
   g2.theta     frozen            0            0      6.28319    radians
   g2.ampl      thawed      228.719            1         1000
   g1.fwhm      frozen      57.9477       1e-120         4300
   g1.xpos      frozen       4070.4      -1e+120         4300
   g1.ypos      frozen      4251.11      -1e+120         4300
   g1.ellip     frozen            0            0        0.999
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      frozen      23.3562            1         1000
   bgnd.c0      thawed     0.260454            0          100

To see if the central component is elliptical - as suggested by the earlier residual image - we "thaw" the ellipticity and position-angle parameters of the second gaussian. They are set to non-zero values to move them away from their minimum values (which helps the fit).

sherpa> thaw(g2.ellip)
sherpa> thaw(g2.theta)
sherpa> g2.ellip = 0.1
sherpa> g2.theta = 1

sherpa> fit()
Statistic value = -54058.2 at function evaluation 649
Data points = 9171
Degrees of freedom = 9164
   g2.fwhm        7.42257
   g2.xpos        4070.73
   g2.ypos        4249.34
   g2.ellip       0.315899
   g2.theta       0.803431
   g2.ampl        230.247
   bgnd.c0        0.260308

The residuals of this fit are shown in Figure 5 .

To see how reliable the ellipticity and theta values are we calculate the one-sigma errors using theprojection() method:

sherpa> projection(g2.ellip, g2.theta)
projection 1-sigma bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g2.ellip         0.315899    -0.022182    0.0215196
   g2.theta         0.803431   -0.0398697    0.0401327

Checking the Sherpa session

The final overall status of this Sherpa session may be viewed as follows:

sherpa> show_all()
Optimization Method: NelderMead
Statistic:           Cash

Data Set: 1
name      = image2.fits
x0        = Float64[56376]
x1        = Float64[56376]
y         = Int16[56376]
shape     = (216, 261)
staterror = None
syserror  = None
sky       = <pytransform.WCSTANTransform; proxy of C++ WCSTANTransform instance at 0x8433558>
eqpos     = <pytransform.WCSTANTransform; proxy of C++ WCSTANTransform instance at 0x841e3a0>
coord     = physical

Model: 1
((g2 + g1) + bgnd)
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   g2.fwhm      thawed      7.42257       1e-120         4300
   g2.xpos      thawed      4070.73      -1e+120         4300
   g2.ypos      thawed      4249.34      -1e+120         4300
   g2.ellip     thawed     0.315899            0        0.999
   g2.theta     thawed     0.803431            0      6.28319    radians
   g2.ampl      thawed      230.247            1         1000
   g1.fwhm      frozen      57.9477       1e-120         4300
   g1.xpos      frozen       4070.4      -1e+120         4300
   g1.ypos      frozen      4251.11      -1e+120         4300
   g1.ellip     frozen            0            0        0.999
   g1.theta     frozen            0            0      6.28319    radians
   g1.ampl      frozen      23.3562            1         1000

Fit: 1
Statistic value = -54058.2 at function evaluation 649
Data points = 9171
Degrees of freedom = 9164
   g2.fwhm        7.42257
   g2.xpos        4070.73
   g2.ypos        4249.34
   g2.ellip       0.315899
   g2.theta       0.803431
   g2.ampl        230.247
   bgnd.c0        0.260308

Projection: 1
projection 1-sigma bounds:
   Param            Best-Fit  Lower Bound  Upper Bound
   -----            --------  -----------  -----------
   g2.ellip         0.315899    -0.022182    0.0215196
   g2.theta         0.803431   -0.0398697    0.0401327

Summary of Commands

The following commands are used to run this thread. They may be saved in a text file and executed as a script: execfile("script.py").

load_image("image2.fits")
image_data()

set_coord("physical")
set_stat("cash")
set_method("neldermead")
show_all()

notice2d("circle(4072.46,4249.34,108)")

set_model(gauss2d.g1)
g1.ampl = 20
g1.fwhm = 20
g1.xpos = 4065.5
g1.ypos = 4250.5
freeze(g1.ellip)
freeze(g1.theta)
g1.fwhm.max = 4300
g1.xpos.max = 4300
g1.ypos.max = 4300
g1.ampl.min = 1
g1.ampl.max = 1000
print(get_model())

set_model(gauss2d.g1+const2d.bgnd)
bgnd.c0 = 0.2
bgnd.c0.max = 100
print(get_model())

fit()
image_resid()

projection()

freeze(g1)
set_model(gauss2d.g2+g1+bgnd)
g2.fwhm = 10
g2.ampl = 100
g2.xpos = 4065.5
g2.ypos = 4250.5
g2.fwhm.max = 4300
g2.xpos.max = 4300
g2.ypos.max = 4300
g2.ampl.min = 1
g2.ampl.max = 1000
freeze(g2.ellip)
freeze(g2.theta)

fit()
print(get_model())

thaw(g2.ellip)
thaw(g2.theta)
g2.ellip = 0.1
g2.theta = 1

fit()

projection(g2.ellip, g2.theta)
show_all()

Summary

In this thread we have shown you how you can fit a two-dimensional model to your image data. As with fitting one-dimensional data, care must be taken to avoid reaching a local, rather than global, minimum. It is suggested that you fit your data using different optimization methods - for instance the fit above could be re-done with the Monte Carlo ("moncar") method - and to try different initial parameter values when fitting.

In the example above the fit was not very complex: possible additions would have been to link the centers of the two gaussians to be the same; simultaneously fit the same model to more than one dataset (e.g. for multiple observations of a source or when analysing the data from the three imaging cameras in XMM-Newton); include an exposure map in the fit to account for instrumental features such as chip gaps and bad columns; or convolve the model by the PSF during fitting.

History

15 May 2008

rewritten for Sherpa Beta