Search https://cxc.harvard.edu/sherpa/ Contact the CXC HelpDesk
Last modified: 17 Dec 2024

URL: https://cxc.cfa.harvard.edu/sherpa/threads/expmap/
/ sherpa / threads / expmap /
Return to Threads Page: Top | All | Fitting

Using an Exposure Map in Fitting Image Data

Sherpa Threads (CIAO 4.17 Sherpa)

Overview

Synopsis:

This thread shows how to use an exposure map when fitting 2-D spatial data. The exposure map file is input to Sherpa as a file-based exposure map model via the load_table_model function.

Last Update: 17 Dec 2024 - reviewed for CIAO 4.17, updated returned defined model string formatting.

Contents

Getting Started

Please follow the Sherpa Getting Started thread.

Reading and Plotting 2-D FITS Data

We are using 2-D spatial data from the FITS data file img.fits. This data set is input to Sherpa with the load_image command: 

sherpa> load_image("img.fits")
sherpa> show_data()
Data Set: 1
Filter: 
name      = img.fits
x0        = Float64[6400]
x1        = Float64[6400]
y         = Float64[6400]
shape     = (80, 80)
staterror = None
syserror  = None
sky       = physical
 crval    = [3944.,3920.]
 crpix    = [0.5,0.5]
 cdelt    = [5.,5.]
eqpos     = world
 crval    = [40.0117,59.9967]
 crpix    = [4096.5,4096.5]
 cdelt    = [-0.0001, 0.0001]
 crota    = 0
 epoch    = 2000
 equinox  = 2000
coord     = logical

The data set may be viewed as a contour plot (contour_data) or an image (image_data). Here we show the contour plot method (Figure 1): 

sherpa> contour_data()

Figure 1: Contour plot of the data

[Contour plot of image data]
[Print media version: Contour plot of image data]

Figure 1: Contour plot of the data

Contour plot of image data set created with the contour_data function.

Setting the Exposure Map

We define a file-based exposure map model by loading an exposure map file with the load_table_model command: 

sherpa> load_table_model("emap", "expmap.fits")

To display the status of the model emap, use the print() command; notice that Sherpa identifies the exposure map model as "tablemodel.emap": 

sherpa> print(emap)
tablemodel.emap
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   emap.ampl    thawed            1 -3.40282e+38  3.40282e+38

Defining and Fitting the Source

One can now define a model to be used as a source model. After viewing Figure 1, the beta2d model is found to be a promising candidate for the source. 

sherpa> set_source(beta2d.b1 * emap) 
sherpa> show_model()
Model: 1
beta2d.b1 * tablemodel.emap
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   b1.r0        thawed           10  1.17549e-38  3.40282e+38           
   b1.xpos      thawed            0 -3.40282e+38  3.40282e+38           
   b1.ypos      thawed            0 -3.40282e+38  3.40282e+38           
   b1.ellip     frozen            0            0        0.999           
   b1.theta     frozen            0     -6.28319      6.28319    radians
   b1.ampl      thawed            1 -3.40282e+38  3.40282e+38           
   b1.alpha     thawed            1          -10           10           
   emap.ampl    thawed            1 -3.40282e+38  3.40282e+38           

sherpa> b1.r0 = 30
sherpa> b1.xpos = 40
sherpa> b1.ypos = 40
sherpa> b1.ellip = 0.3
sherpa> b1.theta = 5
sherpa> b1.ampl = 3.0
sherpa> b1.alpha = 1.5

sherpa> thaw(b1)
sherpa> freeze(emap.ampl)

sherpa> show_model()
Model: 1
beta2d.b1 * tablemodel.emap
   Param        Type          Value          Min          Max      Units
   -----        ----          -----          ---          ---      -----
   b1.r0        thawed           30  1.17549e-38  3.40282e+38           
   b1.xpos      thawed           40 -3.40282e+38  3.40282e+38           
   b1.ypos      thawed           40 -3.40282e+38  3.40282e+38           
   b1.ellip     thawed          0.3            0        0.999           
   b1.theta     thawed            5     -6.28319      6.28319    radians
   b1.ampl      thawed            3 -3.40282e+38  3.40282e+38           
   b1.alpha     thawed          1.5          -10           10           
   emap.ampl    frozen            1 -3.40282e+38  3.40282e+38

Next, we fit the model to the data. Since the data is Poisson, we switch to the Cash statistic and use the Nelder-Mead optimiser: 

sherpa> set_method('neldermead')
sherpa> set_stat('cash')
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = cash
Initial fit statistic = 398229
Final fit statistic   = -723845 at function evaluation 1499
Data points           = 6400
Degrees of freedom    = 6393
Change in statistic   = 1.12207e+06
   b1.r0          12.0514     
   b1.xpos        39.3166     
   b1.ypos        40.9343     
   b1.ellip       0.0193133   
   b1.theta       -1.44989    
   b1.ampl        1.33832     
   b1.alpha       1.58933

To display the fit and residuals of the plot, we first try the contour command (Figure 2): 

sherpa> contour('fit', 'resid')

sherpa> for ax in plt.gcf().axes:
   ...:     ax.set_aspect('equal')
   ...:

Figure 2: Contour plot: fit and residuals

[The top plot shows contours around the source (also drawn as contours), and the bottom plot residuals (as contours, which renders the plot hard to read due to the Poisson nature of the data).]
[Print media version: The top plot shows contours around the source (also drawn as contours), and the bottom plot residuals (as contours, which renders the plot hard to read due to the Poisson nature of the data).]

Figure 2: Contour plot: fit and residuals

The contour plot does not show this data well.

The data can be better viewed in DS9 with the image_fit and image_resid functions. First, image_fit, which displays the data image, model image, and fit image (Figure 3): 

sherpa> image_fit()

and the latter the (data - model) fit residuals image (Figure 4): 

sherpa> image_resid()

Figure 3: Data, model, and fit images

[DS9 display of data image, model image, and fit image]
[Print media version: DS9 display of data image, model image, and fit image]

Figure 3: Data, model, and fit images

DS9 display of data image, model image, and fit image.

Figure 4: Fit residuals image

[DS9 display of fit residuals image]
[Print media version: DS9 display of fit residuals image]

Figure 4: Fit residuals image

DS9 display of the (data-model) fit residuals image.

Saving a Sherpa Session

To save the Sherpa session in order to return to the analysis at a later point: 

sherpa> save("expmap.save")

sherpa> save_all("expmap.ascii")

where the save function records all the information about the current session to the binary file expmap.save, and the save_all function records the session settings to an editable ASCII file.

To restore the session that was saved to the binary file expmap.save or ASCII file expmap.ascii: 

sherpa> restore("expmap.save")

sherpa> %run -i expmap.ascii

Scripting It

The file fit.py is a Python script which performs the primary commands used above; it can be executed by typing %run -i fit.py on the Sherpa command line.

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, et cetera.)

Summary

This thread is complete, so we can exit the Sherpa session: 

sherpa> quit

History

14 Jan 2005 reviewed for CIAO 3.2: no changes
21 Dec 2005 reviewed for CIAO 3.3: no changes
01 Dec 2006 reviewed for CIAO 3.4: no changes
07 Dec 2008 updated for Sherpa 4.1
29 Apr 2009 new script command is available with CIAO 4.1.2
07 Jan 2010 updated for CIAO 4.2
13 Jul 2010 updated for CIAO 4.2 Sherpa v2: removal of S-Lang version of thread.
30 Jan 2012 reviewed for CIAO 4.4 (no changes)
13 Dec 2012 reviewed for CIAO 4.5 (no changes)
05 Dec 2013 reviewed for CIAO 4.6: no changes
18 Mar 2015 reviewed for CIAO 4.7, updated model parameter boundaries.
10 Dec 2015 reviewed for CIAO 4.8, updated screen output
01 Nov 2016 reviewed for CIAO 4.9, no content change
01 Jun 2018 reviewed for CIAO 4.10, no content change
12 Dec 2018 reviewed for CIAO 4.11, no content change
13 Dec 2019 Updated for CIAO 4.12: use of Matplotlib rather than ChIPS and switched to a Poisson-based statistic for the fit.
28 Mar 2022 reviewed for CIAO 4.14, no content change
08 Dec 2022 reviewed for CIAO 4.15, no content change
04 Dec 2023 reviewed for CIAO 4.16, no content change
17 Dec 2024 reviewed for CIAO 4.17, updated returned defined model string formatting.
Return to Threads Page: Top | All | Fitting
/ sherpa / threads / expmap /