Fitting FITS Image Data
![[CXC Logo]](/ciao/imgs/cxc-logo.gif){kind=logo}
Sherpa Threads (CIAO 4.4 Sherpa v1)
Overview
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.
Last Update: 30 Jan 2012 - reviewed for CIAO 4.4 (no changes)
Contents
- Create and Read Fits Image Data
- Filtering Image Data
- Defining a Source Model
- Fitting the model
- Adding an extra component
- Checking the Sherpa session
- Scripting It
- Summary
- History
-
Images
- Figure 1: Image of the data
- Figure 2: Filter region defined on the image
- Figure 3: Filtered data
- Figure 4: FITS image file containing filter information
- Figure 5: Residuals of the best-fit model using: g1+bgnd
- Figure 6: Residuals of the best-fit model using: g1+g2+bgnd
- Figure 7: The g2 model contribution to the best-fit model g1+g2+bgnd
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.
![[ds9 display of the source image]](data_orig.png)
Figure 1: Image of the data
Note that a logarithmic scaling has been chosen to bring out the faint emission of the source.
Setting the coordinate system
The default coordinate system for fitting images is that of logical (also called image) which uses the FITS convention 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 cannot subtract a background data set but either have to fit it - as we do below - or ignore the background signal.
sherpa> set_stat("cash")
The CSTAT statistic - equivalent to the XSPEC C statistic - may also be used in fitting image data. The advantage of CSTAT over 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 Simplex optimizer instead (the NelderMead method is the same as Simplex).
sherpa> set_method("simplex")
The fit should also be tried with the moncar optimization method to ensure that the global best-fit has been found, rather than a local one.
The current Sherpa status may be reviewed by using the show_all command:
sherpa> show_all() Data Set: 1 Filter: name = image2.fits x0 = Float64[56376] x1 = Float64[56376] y = Float64[56376] shape = (216, 261) staterror = None syserror = None sky = physical crval = [ 3798.5 4019.5] crpix = [ 0.5 0.5] cdelt = [ 2. 2.] eqpos = world crval = [ 278.386 -10.5899] crpix = [ 4096.5 4096.5] cdelt = [-0.0001 0.0001] crota = 0 epoch = 2000 equinox = 2000 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.
![[A circle surrounds the source emission and the details of this region are also shown.]](fit_region.png)
Figure 2: Filter region defined on the image
The region displayed on the image is a circle with center at (4072.46,4249.34) and a radius of 108 (using the physical coordinate system). This region will be used to fit the source.
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 saved 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:
One may also specify a FITS or ASCII 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.
![[The region outside the selected circle is no-longer displayed.]](data_filtered.png)
Figure 3: Filtered data
This image shows the data that is to be fit (again using a logarithmic scaling to display the full dynamic range of the data). The excluded pixels have been set to NaN which ds9 displays as the background color of white (using its default preference settings).
Load filter information from file
As of the release of Sherpa version 2 in CIAO 4.2, it is possible to filter 2-D data sets with filter information read from FITS image files using the load_filter function. The filter image input to load_filter should contain 1s and 0s to indicate which pixels should be noticed/ignored, and should match the source image in shape and number of pixels. When the new 'ignore' argument of load_filter is set to 'False' (default), the pixels in the data image which correspond to those marked by 1s in the filter image will be noticed, and the reverse when 'ignore' is 'True'.
The filter information contained in a FITS image filter file would be applied to the source image data set using the load_filter function as follows:
sherpa> load_filter("load_filter_image.fits", ignore=False)
An example of a FITS image file which can be used to filter the source image might look like that shown in Figure 4, where the pixels to be noticed are those which have a value of 1 (because 'ignore=False'), and those which should be excluded from the analysis are the pixels with a value of 0.
![[The pixels in the central point source region, and in the background surrounding the apparent extent of the diffuse source emission, are set to 0 to be excluded from the analysis.]](load_filter_image.png)
Figure 4: FITS image file containing filter information
The filter contained in the example FITS image filter file 'load_filter_image.fits' excludes from the analysis the pixels in the central point source region, and those in the background region surrounding the apparent extent of the diffuse source emission.
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_source(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. For simple, single-component models, the Sherpa guess() function may be used to automatically set sensible starting values and ranges for the model parameters. To have Sherpa automatically query for the initial parameter values when a model is established, set 'paramprompt(True)' (it is 'False' by default).
Note that the ellipticity and its angle (theta) are already 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
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> show_model() Model: 1 gauss2d.g1 Param Type Value Min Max Units ----- ---- ----- --- --- ----- g1.fwhm thawed 20 1.17549e-38 4300 g1.xpos thawed 4065.5 -3.40282e+38 4300 g1.ypos thawed 4250.5 -3.40282e+38 4300 g1.ellip frozen 0 0 0.999 g1.theta frozen 0 0 6.28319 radians g1.ampl thawed 20 1 1000
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_source(g1+const2d.bgnd) sherpa> bgnd.c0 = 0.2 sherpa> bgnd.c0.max = 100
The model is now:
sherpa> show_model(); Model: 1 (gauss2d.g1 + const2d.bgnd) Param Type Value Min Max Units ----- ---- ----- --- --- ----- g1.fwhm thawed 20 1.17549e-38 4300 g1.xpos thawed 4065.5 -3.40282e+38 4300 g1.ypos thawed 4250.5 -3.40282e+38 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() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = 20661.5 Final fit statistic = -48907.8 at function evaluation 510 Data points = 9171 Degrees of freedom = 9166 Change in statistic = 69569.3 g1.fwhm 57.9479 g1.xpos 4070.4 g1.ypos 4251.11 g1.ampl 23.3561 bgnd.c0 0.266362
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. This can be calculated at the Sherpa prompt since we can easily access parameter values:
sherpa> g1.fwhm.val * 0.492
(If you just say g1.fwhm * 0.492 then you will not get the correct answer; the .val suffix for the parameter name is required to read the value)
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()
![[ds9 now shows the residual (data-model) image and graphs along the X and Y axes show profiles of this image]](initial-resid.png)
Figure 5: Residuals of the best-fit model using: g1+bgnd
The residual image of the best-fitting two-dimensional gaussian plus constant model. The source is obviously more complex than this description: an excess in the core is visible as well as spatially-correlated residuals (e.g. regions of positive or negative pixels) at larger scales.
The horizontal and vertical graphs - accessed from the ds9 "View" menu - plot the data along the crosshair axes. The color scale has also been changed from that used in Figure 3.
The resulting image is shown in Figure 5. 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 insufficient to describe all the structure in the source; however it will suffice for this example.
The save_model and save_resid functions may be used to save the 2-D fit results to FITS image files:
sherpa> save_model("model.fits")
sherpa> save_resid("resid.fits")
Here, save_model saves the model array associated with 2-D data set id=1 to the FITS image file 'model.fits', and the save_resid function saves the simple data-minus-model counts residuals array to the FITS image file 'resid.fits'.
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 conf() command to estimate errors on these parameters, since this relies on the change in the statistic value rather than the absolute value:
sherpa> conf() g1.ypos lower bound: -0.178832 g1.fwhm lower bound: -0.28347 g1.fwhm upper bound: 0.28347 g1.xpos lower bound: -0.1807 g1.xpos upper bound: 0.1807 g1.ypos upper bound: 0.178832 g1.ampl lower bound: -0.243301 g1.ampl upper bound: 0.245803 bgnd.c0 lower bound: -0.00931359 bgnd.c0 upper bound: 0.00946505 Dataset = 1 Confidence Method = confidence Iterative Fit Method = None Fitting Method = neldermead Statistic = cash confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- g1.fwhm 57.9477 -0.28347 0.28347 g1.xpos 4070.4 -0.1807 0.1807 g1.ypos 4251.11 -0.178832 0.178832 g1.ampl 23.3562 -0.243301 0.245803 bgnd.c0 0.266365 -0.00931359 0.00946505
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_source(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> fit() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = -50788.3 Final fit statistic = -53918.2 at function evaluation 544 Data points = 9171 Degrees of freedom = 9166 Change in statistic = 3129.92 g2.fwhm 6.08633 g2.xpos 4070.78 g2.ypos 4249.34 g2.ampl 228.72 bgnd.c0 0.260449
Since the g1 model was frozen, those values did not change in the fit:
sherpa> show_model() Model: 1 ((gauss2d.g2 + gauss2d.g1) + const2d.bgnd) Param Type Value Min Max Units ----- ---- ----- --- --- ----- g2.fwhm thawed 6.08633 1.17549e-38 4300 g2.xpos thawed 4070.78 -3.40282e+38 4300 g2.ypos thawed 4249.34 -3.40282e+38 4300 g2.ellip frozen 0 0 0.999 g2.theta frozen 0 0 6.28319 radians g2.ampl thawed 228.72 1 1000 g1.fwhm frozen 57.9479 1.17549e-38 4300 g1.xpos frozen 4070.4 -3.40282e+38 4300 g1.ypos frozen 4251.11 -3.40282e+38 4300 g1.ellip frozen 0 0 0.999 g1.theta frozen 0 0 6.28319 radians g1.ampl frozen 23.3561 1 1000 bgnd.c0 thawed 0.260449 0 100
To see if the central component is elliptical - as suggested by the earlier residual image (Figure 5) - we "thaw" the ellipticity and position-angle parameters of the second gaussian. They are set to non-zero values to help the fit.
sherpa> thaw(g2.ellip) sherpa> thaw(g2.theta) sherpa> g2.ellip = 0.1 sherpa> g2.theta = 1 sherpa> fit() Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = -53955.8 Final fit statistic = -54058.3 at function evaluation 616 Data points = 9171 Degrees of freedom = 9164 Change in statistic = 102.484 g2.fwhm 7.42265 g2.xpos 4070.73 g2.ypos 4249.34 g2.ellip 0.315914 g2.theta 0.803453 g2.ampl 230.245 bgnd.c0 0.2603
The residuals of this fit are shown in Figure 6.
![[The residuals are smaller in amplitude,]](g2-resid.png)
Figure 6: Residuals of the best-fit model using: g1+g2+bgnd
Using two gaussians together with a constant provides a better model of the data than Figure 5. However, we note that the residual image still shows correlated areas. You could add further model components to the fit to account for these regions, or decide to use a different set of functions to try and describe the emission.
To see how reliable the ellipticity and theta values are we calculate the one-sigma errors using the conf() method. In this case we explicitly list the parameters we are interested in (g2.ellip and g2.theta) to avoid calculating the limits for the other parameters.
sherpa> conf(g2.ellip, g2.theta) g2.theta lower bound: -0.0399627 g2.theta upper bound: 0.0399627 g2.ellip lower bound: -0.0222117 g2.ellip upper bound: 0.0215095 Dataset = 1 Confidence Method = confidence Fitting Method = neldermead Statistic = cash confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- g2.ellip 0.315914 -0.0222117 0.0215095 g2.theta 0.803453 -0.0399627 0.0399627
We can use the image_source_component command to image one, or a combination of, individual unconvolved source model components in order to visualize the contribution to the full model being used to fit the image data. (The image_model_component command is available for visualizing PSF-convolved model components, and get_model_component_image/get_source_component_image for accessing the data arrays and preferences which define the corresponding image.)
Here we image the extra g2 gaussian component to quickly visualize its contribution to the fit of the emission.
sherpa> image_source_component(g2)
The g2 model component image is shown in Figure 7.
![[Example of image_source_component functionality]](g2-modcomp.png)
Figure 7: The g2 model contribution to the best-fit model g1+g2+bgnd
An image of the g2 model contribution to the fit of the emission using best-fit model g1+g2+bgnd.
Checking the Sherpa session
The final overall status of this Sherpa session may be viewed as follows:
sherpa> show_all() Data Set: 1 Filter: Circle(4072.46,4249.34,108) name = image2.fits x0 = Float64[56376] x1 = Float64[56376] y = Float64[56376] shape = (216, 261) staterror = None syserror = None sky = physical crval = [ 3798.5 4019.5] crpix = [ 0.5 0.5] cdelt = [ 2. 2.] eqpos = world crval = [ 278.386 -10.5899] crpix = [ 4096.5 4096.5] cdelt = [-0.0001 0.0001] crota = 0 epoch = 2000 equinox = 2000 coord = physical Model: 1 ((gauss2d.g2 + gauss2d.g1) + const2d.bgnd) Param Type Value Min Max Units ----- ---- ----- --- --- ----- g2.fwhm thawed 7.42265 1.17549e-38 4300 g2.xpos thawed 4070.73 -3.40282e+38 4300 g2.ypos thawed 4249.34 -3.40282e+38 4300 g2.ellip thawed 0.315914 0 0.999 g2.theta thawed 0.803453 0 6.28319 radians g2.ampl thawed 230.245 1 1000 g1.fwhm frozen 57.9479 1.17549e-38 4300 g1.xpos frozen 4070.4 -3.40282e+38 4300 g1.ypos frozen 4251.11 -3.40282e+38 4300 g1.ellip frozen 0 0 0.999 g1.theta frozen 0 0 6.28319 radians g1.ampl frozen 23.3561 1 1000 bgnd.c0 thawed 0.2603 0 100 Optimization Method: NelderMead name = simplex ftol = 1.19209289551e-07 maxfev = None initsimplex = 0 finalsimplex = 9 step = None iquad = 1 verbose = 0 Statistic: Cash Maximum likelihood function Fit:Dataset = 1 Method = neldermead Statistic = cash Initial fit statistic = -53955.8 Final fit statistic = -54058.3 at function evaluation 616 Data points = 9171 Degrees of freedom = 9164 Change in statistic = 102.484 g2.fwhm 7.42265 g2.xpos 4070.73 g2.ypos 4249.34 g2.ellip 0.315914 g2.theta 0.803453 g2.ampl 230.245 bgnd.c0 0.2603 Confidence:Dataset = 1 Confidence Method = confidence Fitting Method = neldermead Statistic = cash confidence 1-sigma (68.2689%) bounds: Param Best-Fit Lower Bound Upper Bound ----- -------- ----------- ----------- g2.ellip 0.315914 -0.0222117 0.0215095 g2.theta 0.803453 -0.0399627 0.0399627
Scripting It
The commands used to run this thread may be saved in a text file such as fit.py, which can then be executed as a script with execfile("fit.py"). Alternatively, the Sherpa session can be saved to a binary file with the save command (restored with the restore command), or to an editable ASCII file with save_all.
The Sherpa script command may be used to save everything typed on the command line in a Sherpa session:
sherpa> script(filename="sherpa.log", clobber=False)
(Note that restoring a Sherpa session from such a file could be problematic since it may include syntax errors, unwanted fitting trials, etcetera.)
The CXC is committed to helping Sherpa users transition to new syntax as smoothly as possible. If you have existing Sherpa scripts or save files, submit them to us via the CXC Helpdesk and we will provide the CIAO/Sherpa 4.4 syntax to you.
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 |
| 16 Jul 2008 | minor update to use simplex rather than the neldermead optimisation method as they are the same |
| 11 Dec 2008 | updated for CIAO 4.1 |
| 05 Jan 2010 | updated for CIAO 4.2 |
| 07 Jun 2010 | updated to include save_model and save_resid functions for saving 2-D fit results |
| 12 Jun 2010 | updated for CIAO 4.2 Sherpa v2: load_filter now accepts filter files in FITS image format for filtering 2-D data sets. S-Lang version of thread removed |
| 15 Dec 2010 | updated for CIAO 4.3: the new image_source_component function is available for visualizing the contribution of individual unconvolved model components to a fit |
| 30 Jan 2012 | reviewed for CIAO 4.4 (no changes) |
