This is a fixed-text formatted version of a Jupyter notebook.

Fitting 2D images with Sherpa

Introduction

Sherpa is the X-ray satellite Chandra modeling and fitting application. It enables the user to construct complex models from simple definitions and fit those models to data, using a variety of statistics and optimization methods. The issues of constraining the source position and morphology are common in X- and Gamma-ray astronomy. This notebook will show you how to apply Sherpa to CTA data.

Here we will set up Sherpa to fit the counts map and loading the ancillary images for subsequent use. A relevant test statistic for data with Poisson fluctuations is the one proposed by Cash (1979). The simplex (or Nelder-Mead) fitting algorithm is a good compromise between efficiency and robustness. The source fit is best performed in pixel coordinates.

This tutorial has 2 important parts 1. Generating the Maps 2. The actual fitting with sherpa.

Since sherpa deals only with 2-dim images, the first part of this tutorial shows how to prepare gammapy maps to make classical images.

Necessary imports

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
import astropy.units as u
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from gammapy.extern.pathlib import Path
from gammapy.data import DataStore
from gammapy.irf import EnergyDispersion, make_mean_psf
from gammapy.maps import WcsGeom, MapAxis, Map
from gammapy.cube import MapMaker, PSFKernel

Generate the required Maps

We first generate the required maps using 3 simulated runs on the Galactic center, exactly as in the analysis_3d tutorial.

It is always advisable to make the maps on fine energy bins, and then sum them over to get an image.

In [2]:
# Define which data to use
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
obs_ids = [110380, 111140, 111159]
observations = data_store.get_observations(obs_ids)
In [3]:
energy_axis = MapAxis.from_edges(
    np.logspace(-1, 1.0, 10), unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=(0, 0),
    binsz=0.02,
    width=(10, 8),
    coordsys="GAL",
    proj="CAR",
    axes=[energy_axis],
)
In [4]:
%%time
maker = MapMaker(geom, offset_max=4.0 * u.deg)
maps = maker.run(observations)
CPU times: user 9.8 s, sys: 1.9 s, total: 11.7 s
Wall time: 11.7 s

Making a PSF Map

Make a PSF map and weigh it with the exposure at the source position to get a 2D PSF

In [5]:
# mean PSF
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")
table_psf = make_mean_psf(observations, src_pos)

# PSF kernel used for the model convolution
psf_kernel = PSFKernel.from_table_psf(table_psf, geom, max_radius="0.3 deg")

# get the exposure at the source position
exposure_at_pos = maps["exposure"].get_by_coord(
    {
        "lon": src_pos.l.value,
        "lat": src_pos.b.value,
        "energy": energy_axis.center,
    }
)

# now compute the 2D PSF
psf2D = psf_kernel.make_image(exposures=exposure_at_pos)

Make 2D images from 3D ones

Since sherpa image fitting works only with 2-dim images, we convert the generated maps to 2D images using make_images() and save them as fits files. The exposure is weighed with the spectrum before averaging (assumed to be a power law by default).

In [6]:
maps = maker.make_images()
In [7]:
Path("analysis_3d").mkdir(exist_ok=True)

maps["counts"].write("analysis_3d/counts_2D.fits", overwrite=True)
maps["background"].write("analysis_3d/background_2D.fits", overwrite=True)
maps["exposure"].write("analysis_3d/exposure_2D.fits", overwrite=True)
fits.writeto("analysis_3d/psf_2D.fits", psf2D.data, overwrite=True)

Read the maps and store them in a sherpa model

We now have the prepared files which sherpa can read. This part of the notebook shows how to do image analysis using sherpa

In [8]:
import sherpa.astro.ui as sh

sh.set_stat("cash")
sh.set_method("simplex")

sh.load_image("analysis_3d/counts_2D.fits")
sh.set_coord("logical")

sh.load_table_model("expo", "analysis_3d/exposure_2D.fits")
sh.load_table_model("bkg", "analysis_3d/background_2D.fits")
sh.load_psf("psf", "analysis_3d/psf_2D.fits")
WARNING: imaging routines will not be available,
failed to import sherpa.image.ds9_backend due to
'RuntimeErr: DS9Win unusable: Could not find ds9 on your PATH'
WARNING: failed to import sherpa.astro.xspec; XSPEC models will not be available

In principle one might first want to fit the background amplitude. However the background estimation method already yields the correct normalization, so we freeze the background amplitude to unity instead of adjusting it. The (smoothed) residuals from this background model are then computed and shown.

In [9]:
sh.set_full_model(bkg)
bkg.ampl = 1
sh.freeze(bkg)
In [10]:
resid = Map.read("analysis_3d/counts_2D.fits")
resid.data = sh.get_data_image().y - sh.get_model_image().y
resid_smooth = resid.smooth(width=4)
resid_smooth.plot(add_cbar=True);
../_images/notebooks_image_fitting_with_sherpa_17_0.png

Find and fit the brightest source

We then find the position of the maximum in the (smoothed) residuals map, and fit a (symmetrical) Gaussian source with that initial position:

In [11]:
yp, xp = np.unravel_index(
    np.nanargmax(resid_smooth.data), resid_smooth.data.shape
)
ampl = resid_smooth.get_by_pix((xp, yp))[0]

# creates g0 as a gauss2d instance
sh.set_full_model(bkg + psf(sh.gauss2d.g0) * expo)
g0.xpos, g0.ypos = xp, yp
sh.freeze(g0.xpos, g0.ypos)  # fix the position in the initial fitting step

# fix exposure amplitude so that typical exposure is of order unity
expo.ampl = 1e-9
sh.freeze(expo)
sh.thaw(g0.fwhm, g0.ampl)  # in case frozen in a previous iteration

g0.fwhm = 10  # give some reasonable initial values
g0.ampl = ampl
In [12]:
%%time
sh.fit()
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<timed eval> in <module>

<string> in fit(id, *otherids, **kwargs)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in fit(self, id, *otherids, **kwargs)
   9449         """
   9450         kwargs['bkg_only'] = False
-> 9451         self._fit(id, *otherids, **kwargs)
   9452
   9453     def fit_bkg(self, id=None, *otherids, **kwargs):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _fit(self, id, *otherids, **kwargs)
   9546             ids, f = self._get_bkg_fit(id, otherids)
   9547         else:
-> 9548             ids, f = self._get_fit(id, otherids)
   9549
   9550         if 'filter_nan' in kwargs and kwargs.pop('filter_nan'):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_fit(self, id, otherids, estmethod)
   7467     def _get_fit(self, id, otherids=(), estmethod=None):
   7468
-> 7469         fit_to_ids, datasets, models = self._prepare_fit(id, otherids)
   7470
   7471         self._add_extra_data_and_models(fit_to_ids, datasets, models)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _prepare_fit(self, id, otherids)
   7442             mod = None
   7443             if i in self._models or i in self._sources:
-> 7444                 mod = self._get_model(i)
   7445
   7446             # The issue with putting a try/catch here is that if an exception

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_model(self, id)
   5507         data = self.get_data(id)
   5508         model, is_source = self._get_model_status(id)
-> 5509         return self._add_convolution_models(id, data, model, is_source)
   5510
   5511     def _get_source(self, id=None):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   8394         model = \
   8395             sherpa.ui.utils.Session._add_convolution_models(self, id, data,
-> 8396                                                             model, is_source)
   8397         id = self._fix_id(id)
   8398         if (isinstance(data, sherpa.astro.data.DataPHA) and is_source):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in <listcomp>(.0)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/instrument.py in fold(self, data)
   1097
   1098     def fold(self, data):
-> 1099         _PSFModel.fold(self, data)
   1100
   1101         # Set WCS coordinates of kernel data set to match source data set.

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in fold(self, data)
    555         # FIXME how will we know the native dimensionality of the
    556         # raveled model without the values?
--> 557         self._check_pixel_size(data)
    558
    559         kargs = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in _check_pixel_size(self, data)
    718         if hasattr(self.kernel, "sky"):
    719             # This corresponds to the case when the kernel is actually a psf image, not just a model.
--> 720             psf_pixel_size = self.kernel.sky.cdelt
    721             data_pixel_size = data.sky.cdelt
    722

AttributeError: 'NoneType' object has no attribute 'cdelt'

Fit all parameters of this Gaussian component, fix them and re-compute the residuals map.

In [13]:
sh.thaw(g0.xpos, g0.ypos)
sh.fit()
sh.freeze(g0)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-13-d0c708c71933> in <module>
      1 sh.thaw(g0.xpos, g0.ypos)
----> 2 sh.fit()
      3 sh.freeze(g0)

<string> in fit(id, *otherids, **kwargs)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in fit(self, id, *otherids, **kwargs)
   9449         """
   9450         kwargs['bkg_only'] = False
-> 9451         self._fit(id, *otherids, **kwargs)
   9452
   9453     def fit_bkg(self, id=None, *otherids, **kwargs):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _fit(self, id, *otherids, **kwargs)
   9546             ids, f = self._get_bkg_fit(id, otherids)
   9547         else:
-> 9548             ids, f = self._get_fit(id, otherids)
   9549
   9550         if 'filter_nan' in kwargs and kwargs.pop('filter_nan'):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_fit(self, id, otherids, estmethod)
   7467     def _get_fit(self, id, otherids=(), estmethod=None):
   7468
-> 7469         fit_to_ids, datasets, models = self._prepare_fit(id, otherids)
   7470
   7471         self._add_extra_data_and_models(fit_to_ids, datasets, models)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _prepare_fit(self, id, otherids)
   7442             mod = None
   7443             if i in self._models or i in self._sources:
-> 7444                 mod = self._get_model(i)
   7445
   7446             # The issue with putting a try/catch here is that if an exception

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_model(self, id)
   5507         data = self.get_data(id)
   5508         model, is_source = self._get_model_status(id)
-> 5509         return self._add_convolution_models(id, data, model, is_source)
   5510
   5511     def _get_source(self, id=None):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   8394         model = \
   8395             sherpa.ui.utils.Session._add_convolution_models(self, id, data,
-> 8396                                                             model, is_source)
   8397         id = self._fix_id(id)
   8398         if (isinstance(data, sherpa.astro.data.DataPHA) and is_source):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in <listcomp>(.0)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/instrument.py in fold(self, data)
   1097
   1098     def fold(self, data):
-> 1099         _PSFModel.fold(self, data)
   1100
   1101         # Set WCS coordinates of kernel data set to match source data set.

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in fold(self, data)
    555         # FIXME how will we know the native dimensionality of the
    556         # raveled model without the values?
--> 557         self._check_pixel_size(data)
    558
    559         kargs = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in _check_pixel_size(self, data)
    718         if hasattr(self.kernel, "sky"):
    719             # This corresponds to the case when the kernel is actually a psf image, not just a model.
--> 720             psf_pixel_size = self.kernel.sky.cdelt
    721             data_pixel_size = data.sky.cdelt
    722

AttributeError: 'NoneType' object has no attribute 'cdelt'
In [14]:
resid.data = sh.get_data_image().y - sh.get_model_image().y
resid_smooth = resid.smooth(width=3)
resid_smooth.plot();
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-14-8dacc4a753bc> in <module>
----> 1 resid.data = sh.get_data_image().y - sh.get_model_image().y
      2 resid_smooth = resid.smooth(width=3)
      3 resid_smooth.plot();

<string> in get_model_image(id)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in get_model_image(self, id)
  14355
  14356         """
> 14357         self._prepare_imageobj(id, self._modelimage)
  14358         return self._modelimage
  14359

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _prepare_imageobj(self, id, imageobj, model)
  14264             imageobj.prepare_image(self.get_data(id), self.get_source(id))
  14265         else:
> 14266             imageobj.prepare_image(self.get_data(id), self.get_model(id))
  14267         return imageobj
  14268

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in get_model(self, id)
   5611
   5612         """
-> 5613         return self._get_model(id)
   5614
   5615     def _runparamprompt(self, pars):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_model(self, id)
   5507         data = self.get_data(id)
   5508         model, is_source = self._get_model_status(id)
-> 5509         return self._add_convolution_models(id, data, model, is_source)
   5510
   5511     def _get_source(self, id=None):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   8394         model = \
   8395             sherpa.ui.utils.Session._add_convolution_models(self, id, data,
-> 8396                                                             model, is_source)
   8397         id = self._fix_id(id)
   8398         if (isinstance(data, sherpa.astro.data.DataPHA) and is_source):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in <listcomp>(.0)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/instrument.py in fold(self, data)
   1097
   1098     def fold(self, data):
-> 1099         _PSFModel.fold(self, data)
   1100
   1101         # Set WCS coordinates of kernel data set to match source data set.

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in fold(self, data)
    555         # FIXME how will we know the native dimensionality of the
    556         # raveled model without the values?
--> 557         self._check_pixel_size(data)
    558
    559         kargs = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in _check_pixel_size(self, data)
    718         if hasattr(self.kernel, "sky"):
    719             # This corresponds to the case when the kernel is actually a psf image, not just a model.
--> 720             psf_pixel_size = self.kernel.sky.cdelt
    721             data_pixel_size = data.sky.cdelt
    722

AttributeError: 'NoneType' object has no attribute 'cdelt'

Iteratively find and fit additional sources

Instantiate additional Gaussian components, and use them to iteratively fit sources, repeating the steps performed above for component g0. (The residuals map is shown after each additional source included in the model.) This takes some time…

In [15]:
# initialize components with fixed, zero amplitude
for i in range(1, 10):
    model = sh.create_model_component("gauss2d", "g" + str(i))
    model.ampl = 0
    sh.freeze(model)

gs = [g0, g1, g2]
sh.set_full_model(bkg + psf(g0 + g1 + g2) * expo)
In [16]:
%%time
for i in range(1, len(gs)):
    yp, xp = np.unravel_index(
        np.nanargmax(resid_smooth.data), resid_smooth.data.shape
    )
    ampl = resid_smooth.get_by_pix((xp, yp))[0]
    gs[i].xpos, gs[i].ypos = xp, yp
    gs[i].fwhm = 10
    gs[i].ampl = ampl

    sh.thaw(gs[i].fwhm)
    sh.thaw(gs[i].ampl)
    sh.fit()

    sh.thaw(gs[i].xpos)
    sh.thaw(gs[i].ypos)
    sh.fit()
    sh.freeze(gs[i])

    resid.data = sh.get_data_image().y - sh.get_model_image().y
    resid_smooth = resid.smooth(width=6)
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<timed exec> in <module>

<string> in fit(id, *otherids, **kwargs)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in fit(self, id, *otherids, **kwargs)
   9449         """
   9450         kwargs['bkg_only'] = False
-> 9451         self._fit(id, *otherids, **kwargs)
   9452
   9453     def fit_bkg(self, id=None, *otherids, **kwargs):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _fit(self, id, *otherids, **kwargs)
   9546             ids, f = self._get_bkg_fit(id, otherids)
   9547         else:
-> 9548             ids, f = self._get_fit(id, otherids)
   9549
   9550         if 'filter_nan' in kwargs and kwargs.pop('filter_nan'):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_fit(self, id, otherids, estmethod)
   7467     def _get_fit(self, id, otherids=(), estmethod=None):
   7468
-> 7469         fit_to_ids, datasets, models = self._prepare_fit(id, otherids)
   7470
   7471         self._add_extra_data_and_models(fit_to_ids, datasets, models)

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _prepare_fit(self, id, otherids)
   7442             mod = None
   7443             if i in self._models or i in self._sources:
-> 7444                 mod = self._get_model(i)
   7445
   7446             # The issue with putting a try/catch here is that if an exception

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_model(self, id)
   5507         data = self.get_data(id)
   5508         model, is_source = self._get_model_status(id)
-> 5509         return self._add_convolution_models(id, data, model, is_source)
   5510
   5511     def _get_source(self, id=None):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   8394         model = \
   8395             sherpa.ui.utils.Session._add_convolution_models(self, id, data,
-> 8396                                                             model, is_source)
   8397         id = self._fix_id(id)
   8398         if (isinstance(data, sherpa.astro.data.DataPHA) and is_source):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in <listcomp>(.0)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/instrument.py in fold(self, data)
   1097
   1098     def fold(self, data):
-> 1099         _PSFModel.fold(self, data)
   1100
   1101         # Set WCS coordinates of kernel data set to match source data set.

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in fold(self, data)
    555         # FIXME how will we know the native dimensionality of the
    556         # raveled model without the values?
--> 557         self._check_pixel_size(data)
    558
    559         kargs = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in _check_pixel_size(self, data)
    718         if hasattr(self.kernel, "sky"):
    719             # This corresponds to the case when the kernel is actually a psf image, not just a model.
--> 720             psf_pixel_size = self.kernel.sky.cdelt
    721             data_pixel_size = data.sky.cdelt
    722

AttributeError: 'NoneType' object has no attribute 'cdelt'
In [17]:
resid_smooth.plot(add_cbar=True);
../_images/notebooks_image_fitting_with_sherpa_27_0.png

Generating output table and Test Statistics estimation

When adding a new source, one needs to check the significance of this new source. A frequently used method is the Test Statistics (TS). This is done by comparing the change of statistics when the source is included compared to the null hypothesis (no source ; in practice here we fix the amplitude to zero).

\(TS = Cstat(source) - Cstat(no source)\)

The criterion for a significant source detection is typically that it should improve the test statistic by at least 25 or 30. We have added only 3 sources to save time, but you should keep doing this till del(stat) is less than the required number.

In [18]:
from astropy.stats import gaussian_fwhm_to_sigma
from astropy.table import Table

rows = []
for g in gs:
    ampl = g.ampl.val
    g.ampl = 0
    stati = sh.get_stat_info()[0].statval
    g.ampl = ampl
    statf = sh.get_stat_info()[0].statval
    delstat = stati - statf

    geom = resid.geom
    coord = geom.pix_to_coord((g.xpos.val, g.ypos.val))
    pix_scale = geom.pixel_scales.mean().deg
    sigma = g.fwhm.val * pix_scale * gaussian_fwhm_to_sigma
    rows.append(
        dict(delstat=delstat, glon=coord[0], glat=coord[1], sigma=sigma)
    )

table = Table(rows=rows, names=rows[0])
for name in table.colnames:
    table[name].format = ".5g"
table
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-18-bdd7ade33e59> in <module>
      6     ampl = g.ampl.val
      7     g.ampl = 0
----> 8     stati = sh.get_stat_info()[0].statval
      9     g.ampl = ampl
     10     statf = sh.get_stat_info()[0].statval

<string> in get_stat_info()

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in get_stat_info(self)
   7614
   7615         """
-> 7616         return self._get_stat_info()
   7617
   7618     def get_fit_results(self):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _get_stat_info(self)
   9560     def _get_stat_info(self):
   9561
-> 9562         ids, datasets, models = self._prepare_fit(None)
   9563
   9564         extra_ids = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _prepare_fit(self, id, otherids)
   7442             mod = None
   7443             if i in self._models or i in self._sources:
-> 7444                 mod = self._get_model(i)
   7445
   7446             # The issue with putting a try/catch here is that if an exception

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _get_model(self, id)
   5507         data = self.get_data(id)
   5508         model, is_source = self._get_model_status(id)
-> 5509         return self._add_convolution_models(id, data, model, is_source)
   5510
   5511     def _get_source(self, id=None):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   8394         model = \
   8395             sherpa.ui.utils.Session._add_convolution_models(self, id, data,
-> 8396                                                             model, is_source)
   8397         id = self._fix_id(id)
   8398         if (isinstance(data, sherpa.astro.data.DataPHA) and is_source):

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in _add_convolution_models(self, id, data, model, is_source)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/ui/utils.py in <listcomp>(.0)
   5500             model = self._add_psf(id, data, model)
   5501         else:
-> 5502             [psf.fold(data) for psf in self._psf_models if psf in model]
   5503
   5504         return model

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/astro/instrument.py in fold(self, data)
   1097
   1098     def fold(self, data):
-> 1099         _PSFModel.fold(self, data)
   1100
   1101         # Set WCS coordinates of kernel data set to match source data set.

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in fold(self, data)
    555         # FIXME how will we know the native dimensionality of the
    556         # raveled model without the values?
--> 557         self._check_pixel_size(data)
    558
    559         kargs = {}

~/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/sherpa/instrument.py in _check_pixel_size(self, data)
    718         if hasattr(self.kernel, "sky"):
    719             # This corresponds to the case when the kernel is actually a psf image, not just a model.
--> 720             psf_pixel_size = self.kernel.sky.cdelt
    721             data_pixel_size = data.sky.cdelt
    722

AttributeError: 'NoneType' object has no attribute 'cdelt'

Exercises

  1. Keep adding sources till there are no more significat ones in the field. How many Gaussians do you need?
  2. Use other morphologies for the sources (eg: disk, shell) rather than only Gaussian.
  3. Compare the TS between different models

More about sherpa

These are good resources to learn more about Sherpa:

You could read over the examples there, and try to apply a similar analysis to this dataset here to practice.

If you want a deeper understanding of how Sherpa works, then these proceedings are good introductions: