Last modified: 12 Dec 2018


Using an Exposure Map in Fitting Image Data

Sherpa Threads (CIAO 4.11 Sherpa v1)



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: 12 Dec 2018 - reviewed for CIAO 4.11, no content change


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
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, creating a PostScript file of the output as well:

sherpa> contour_data()
sherpa> print_window("contour_plot")

This creates Figure 1.

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)
   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. The BETA2D model is defined for the source with set_source, then the initial parameter values are specified:

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.ellip, b1.theta)
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. By default, Sherpa uses Chi2Gehrels as the fit statistic and LevMar as the fit optimization method. We decide to change the optimization method for this fit to Nelder-Mead:

sherpa> set_method("neldermead")
sherpa> fit()
Dataset               = 1
Method                = neldermead
Statistic             = chi2gehrels
Initial fit statistic = 4.88095e+06
Final fit statistic   = 3255.75 at function evaluation 2294
Data points           = 6400
Degrees of freedom    = 6393
Probability [Q-value] = 1
Reduced statistic     = 0.509268
Change in statistic   = 4.8777e+06
   b1.r0          12.4624     
   b1.xpos        39.5139     
   b1.ypos        40.8959     
   b1.ellip       0.0259202   
   b1.theta       -4.6965     
   b1.ampl        1.31312     
   b1.alpha       1.66641     

To display the fit and residuals of the plot, we use image_fit and image_resid; the former opens in DS9 a display of the the data image, model image, and fit image, and the latter the (data - model) fit residuals image.

sherpa> image_fit()

sherpa> image_resid()

These commands create Figure 2 and Figure 3.

Figure 2: 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 2: Data, model, and fit images

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

Figure 3: Fit residuals image

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

Figure 3: 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:

where the save function records all the information about the current session to the binary file, 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 or ASCII file expmap.ascii:

sherpa> restore("")

sherpa> exec(open("expmap.ascii").read())

Scripting It

The file is a Python script which performs the primary commands used above; it can be executed by typing exec(open("").read()) 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.)

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.11 syntax to you.


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

sherpa> quit


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