gammapy.cube.make.MapDatasetMaker.make_exposure_irf() - Code Metrics - Inspection of "Rename and refactor `MapMakerObs`" - gammapy/gammapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Pull Request — master (#2450)

by Axel

created 2019-10-14 08:02 UTC

MapDatasetMaker.make_exposure_irf() A

↳ Parent: gammapy.cube.make

Complexity

Conditions

Size

Total Lines	22
Code Lines	9

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
eloc	9
dl	0
loc	22
rs	9.95
c	0
b	0
f	0
cc	1
nop	2

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
from functools import lru_cache
import numpy as np
from astropy.coordinates import Angle
from astropy.nddata.utils import NoOverlapError
from astropy.utils import lazyproperty
from gammapy.irf import EnergyDependentMultiGaussPSF
from gammapy.maps import Map, WcsGeom
from gammapy.modeling.models import BackgroundModel
from .background import make_map_background_irf
from .counts import fill_map_counts
from .edisp_map import make_edisp_map
from .exposure import _map_spectrum_weight, make_map_exposure_true_energy
from .fit import MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, BINSZ_IRF, MapDataset
from .psf_map import make_psf_map

__all__ = ["MapDatasetMaker", "MapMakerRing"]

log = logging.getLogger(__name__)


class MapDatasetMaker:
    """Make maps for a single IACT observation.

    Parameters
    ----------
    geom : `~gammapy.maps.WcsGeom`
        Reference image geometry in reco energy, used for counts and background maps
    offset_max : `~astropy.coordinates.Angle`
        Maximum offset angle
    geom_true : `~gammapy.maps.WcsGeom`
        Reference image geometry in true energy, used for IRF maps. It can have a coarser
        spatial bins than the counts geom.
        If none, the same as geom is assumed
    migra_axis : `~gammapy.maps.MapAxis`
        Migration axis for edisp map
    rad_axis : `~gammapy.maps.MapAxis`
        Radial axis for psf map.
    cutout : bool
         Whether to cutout the observation.
    cutout_mode : {'trim', 'partial', 'strict'}
        Mode option for cutting out the observation,
        for details see `~astropy.nddata.utils.Cutout2D`.
    """

    def __init__(
        self,
        geom,
        offset_max,
        geom_true=None,
        background_oversampling=None,
        migra_axis=None,
        rad_axis=None,
        cutout_mode="trim",
        cutout=True
    ):
        self.geom = geom
        self.geom_true = geom_true if geom_true else geom.to_binsz(BINSZ_IRF)
        self.offset_max = Angle(offset_max)
        self.background_oversampling = background_oversampling
        self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT
        self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT
        self.cutout_mode = cutout_mode
        self.cutout_width = 2 * self.offset_max
        self.cutout = cutout

    def _cutout_geom(self, geom, observation):
        if self.cutout:
            cutout_kwargs = {
                "position": observation.pointing_radec,
                "width": self.cutout_width,
                "mode": self.cutout_mode,
            }
            return geom.cutout(**cutout_kwargs)
        else:
            return geom

    @lazyproperty
    def geom_exposure(self):
        """Exposure map geom (`Geom`)"""
        energy_axis = self.geom_true.get_axis_by_name("energy")
        geom_exposure = self.geom.to_image().to_cube([energy_axis])
        return geom_exposure

    @lazyproperty
    def geom_psf(self):
        """PSFMap geom (`Geom`)"""
        energy_axis = self.geom_true.get_axis_by_name("ENERGY")
        geom_psf = self.geom_true.to_image().to_cube([self.rad_axis, energy_axis])
        return geom_psf

    @lazyproperty
    def geom_edisp(self):
        """EdispMap geom (`Geom`)"""
        energy_axis = self.geom_true.get_axis_by_name("ENERGY")
        geom_edisp = self.geom_true.to_image().to_cube([self.migra_axis, energy_axis])
        return geom_edisp

    def make_counts(self, observation):
        """Make counts map.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        counts : `Map`
            Counts map.
        """
        geom = self._cutout_geom(self.geom, observation)
        counts = Map.from_geom(geom)
        fill_map_counts(counts, observation.events)
        return counts

    def make_exposure(self, observation):
        """Make exposure map.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        exposure : `Map`
            Exposure map.
        """
        geom = self._cutout_geom(self.geom_exposure, observation)
        exposure = make_map_exposure_true_energy(
            pointing=observation.pointing_radec,
            livetime=observation.observation_live_time_duration,
            aeff=observation.aeff,
            geom=geom,
        )
        return exposure

    @lru_cache(maxsize=1)
    def make_exposure_irf(self, observation):
        """Make exposure map with irf geometry.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        exposure : `Map`
            Exposure map.
        """
        geom = self._cutout_geom(self.geom_true, observation)
        exposure = make_map_exposure_true_energy(
            pointing=observation.pointing_radec,
            livetime=observation.observation_live_time_duration,
            aeff=observation.aeff,
            geom=geom,
        )
        return exposure

    def make_background(self, observation):
        """Make background map.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        background : `Map`
            Background map.
        """
        geom = self._cutout_geom(self.geom, observation)

        bkg_coordsys = observation.bkg.meta.get("FOVALIGN", "ALTAZ")
        if bkg_coordsys == "ALTAZ":
            pointing = observation.fixed_pointing_info
        elif bkg_coordsys == "RADEC":
            pointing = observation.pointing_radec
        else:
            raise ValueError(
                f"Invalid background coordinate system: {bkg_coordsys!r}\n"
                "Options: ALTAZ, RADEC"
            )
        background = make_map_background_irf(
            pointing=pointing,
            ontime=observation.observation_time_duration,
            bkg=observation.bkg,
            geom=geom,
            oversampling=self.background_oversampling,
        )
        return background

    def make_edisp(self, observation):
        """Make edisp map.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        edisp : `EdispMap`
            Edisp map.
        """
        geom = self._cutout_geom(self.geom_edisp, observation)

        exposure = self.make_exposure_irf(observation)

        edisp = make_edisp_map(
            edisp=observation.edisp,
            pointing=observation.pointing_radec,
            geom=geom,
            max_offset=self.offset_max,
            exposure_map=exposure,
        )
        return edisp

    def make_psf(self, observation):
        """Make psf map.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        psf : `PSFMap`
            Psf map.
        """
        psf = observation.psf
        geom = self._cutout_geom(self.geom_psf, observation)

        if isinstance(psf, EnergyDependentMultiGaussPSF):
            psf = psf.to_psf3d(rad=self.rad_axis.center)

        exposure = self.make_exposure_irf(observation)

        psf = make_psf_map(
            psf=psf,
            pointing=observation.pointing_radec,
            geom=geom,
            max_offset=self.offset_max,
            exposure_map=exposure,
        )
        return psf

    @lru_cache(maxsize=1)
    def make_mask_safe(self, observation):
        """Make offset mask.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        mask : `~numpy.ndarray`
            Mask
        """
        geom = self._cutout_geom(self.geom, observation)
        offset = geom.separation(observation.pointing_radec)
        return offset >= self.offset_max

    @lru_cache(maxsize=1)
    def make_mask_safe_irf(self, observation):
        """Make offset mask with irf geometry.

        Parameters
        ----------
        observation : `DataStoreObservation`
            Observation container.

        Returns
        -------
        mask : `~numpy.ndarray`
            Mask
        """
        geom = self._cutout_geom(self.geom_true, observation)
        offset = geom.separation(observation.pointing_radec)
        return offset >= self.offset_max

    def run(self, observation, selection=None):
        """Make map dataset.

        Parameters
        ----------
        observation : `~gammapy.data.DataStoreObservation`
            Observation
        selection : list
            List of str, selecting which maps to make.
            Available: 'counts', 'exposure', 'background', 'psf', 'edisp'
            By default, all maps are made.

        Returns
        -------
        dataset : `MapDataset`
            Map dataset.

        """
        selection = _check_selection(selection)

        kwargs = {}
        kwargs["name"] = "obs_{}".format(observation.obs_id)
        kwargs["gti"] = observation.gti

        mask_safe = self.make_mask_safe(observation)

        # expand mask safe into 3rd dimension
        nbin = self.geom.get_axis_by_name("energy").nbin
        kwargs["mask_safe"] = ~mask_safe & np.ones(nbin, dtype=bool)[:, np.newaxis, np.newaxis]

        mask_safe_irf = self.make_mask_safe_irf(observation)

        if "counts" in selection:
            counts = self.make_counts(observation)
            # TODO: remove masking out the values here and instead handle the safe mask only when
            #  fitting and / or stacking datasets?
            counts.data[..., mask_safe] = 0
            kwargs["counts"] = counts

        if "exposure" in selection:
            exposure = self.make_exposure(observation)
            exposure.data[..., mask_safe] = 0
            kwargs["exposure"] = exposure

        if "background" in selection:
            background_map = self.make_background(observation)
            background_map.data[..., mask_safe] = 0
            kwargs["background_model"] = BackgroundModel(background_map)

        if "psf" in selection:
            psf = self.make_psf(observation)
            psf.exposure_map.data[..., mask_safe_irf] = 0
            kwargs["psf"] = psf

        if "edisp" in selection:
            edisp = self.make_edisp(observation)
            psf.exposure_map.data[..., mask_safe_irf] = 0

            kwargs["edisp"] = edisp

        return MapDataset(**kwargs)


def _check_selection(selection):
    """Handle default and validation of selection"""
    available = ["counts", "exposure", "background", "psf", "edisp"]
    if selection is None:
        selection = available

    if not isinstance(selection, list):
        raise TypeError("Selection must be a list of str")

    for name in selection:
        if name not in available:
            raise ValueError(f"Selection not available: {name!r}")

    return selection


class MapMakerRing:
    """Make maps from IACT observations.

    The main motivation for this class in addition to the `MapMaker`
    is to have the common image background estimation methods,
    like `~gammapy.cube.RingBackgroundEstimator`,
    that work using on and off maps.

    To ensure adequate statistics, only observations that are fully
    contained within the reference geometry will be analysed

    Parameters
    ----------
    geom : `~gammapy.maps.WcsGeom`
        Reference image geometry
    offset_max : `~astropy.coordinates.Angle`
        Maximum offset angle
    exclusion_mask : `~gammapy.maps.Map`
        Exclusion mask
    background_estimator : `~gammapy.cube.RingBackgroundEstimator`
        or `~gammapy.cube.AdaptiveRingBackgroundEstimator`
        Ring background estimator or something with an equivalent API.

    Examples
    --------
    Here is an example how to ise the MapMakerRing with H.E.S.S. DL3 data::

        import numpy as np
        import astropy.units as u
        from astropy.coordinates import SkyCoord
        from regions import CircleSkyRegion
        from gammapy.maps import Map, WcsGeom, MapAxis
        from gammapy.cube import MapMakerRing, RingBackgroundEstimator
        from gammapy.data import DataStore

        # Create observation list
        data_store = DataStore.from_file(
            "$GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz"
        )
        data_sel = data_store.obs_table["TARGET_NAME"] == "MSH 15-52"
        obs_table = data_store.obs_table[data_sel]
        observations = data_store.get_observations(obs_table["OBS_ID"])

        # Define the geom
        pos = SkyCoord(228.32, -59.08, unit="deg")
        energy_axis = MapAxis.from_edges(np.logspace(0, 5.0, 5), unit="TeV", name="energy")
        geom = WcsGeom.create(skydir=pos, binsz=0.02, width=(5, 5), axes=[energy_axis])

        # Make a region mask
        regions = CircleSkyRegion(center=pos, radius=0.3 * u.deg)
        mask = Map.from_geom(geom)
        mask.data = mask.geom.region_mask([regions], inside=False)

        # Run map maker with ring background estimation
        ring_bkg = RingBackgroundEstimator(r_in="0.5 deg", width="0.3 deg")
        maker = MapMakerRing(
            geom=geom, offset_max="2 deg", exclusion_mask=mask, background_estimator=ring_bkg
        )
        images = maker.run_images(observations)
    """

    def __init__(
        self, geom, offset_max, exclusion_mask=None, background_estimator=None
    ):

        self.geom = geom
        self.offset_max = Angle(offset_max)
        self.exclusion_mask = exclusion_mask
        self.background_estimator = background_estimator

    def _get_empty_maps(self, selection):
        # Initialise zero-filled maps
        maps = {}
        for name in selection:
            if name == "exposure":
                maps[name] = Map.from_geom(self.geom, unit="m2 s")
            else:
                maps[name] = Map.from_geom(self.geom, unit="")
        return maps

    def _get_obs_maker(self, obs):
        # Compute cutout geometry and slices to stack results back later

        # Make maps for this observation
        return MapDatasetMaker(
            geom=self.geom,
            offset_max=self.offset_max,
        )

    @staticmethod
    def _maps_sum_over_axes(maps, spectrum, keepdims):
        """Compute weighted sum over map axes.

        Parameters
        ----------
        spectrum : `~gammapy.modeling.models.SpectralModel`
            Spectral model to compute the weights.
            Default is power-law with spectral index of 2.
        keepdims : bool, optional
            If this is set to True, the energy axes is kept with a single bin.
            If False, the energy axes is removed
        """
        images = {}
        for name, map in maps.items():
            if name == "exposure":
                map = _map_spectrum_weight(map, spectrum)
            images[name] = map.sum_over_axes(keepdims=keepdims)
        # TODO: PSF (and edisp) map sum_over_axis

        return images


    def _run(self, observations, sum_over_axis=False, spectrum=None, keepdims=False):
        selection = ["on", "exposure_on", "off", "exposure_off", "exposure"]
        maps = self._get_empty_maps(selection)
        if sum_over_axis:
            maps = self._maps_sum_over_axes(maps, spectrum, keepdims)

        for obs in observations:
            try:
                obs_maker = self._get_obs_maker(obs)
            except NoOverlapError:
                log.info(f"Skipping obs_id: {obs.obs_id} (no map overlap)")
                continue

            dataset = obs_maker.run(obs, selection=["counts", "exposure", "background"])
            maps_obs = {}
            maps_obs["counts"] = dataset.counts
            maps_obs["exposure"] = dataset.exposure
            maps_obs["background"] = dataset.background_model.map
            maps_obs["exclusion"] = self.exclusion_mask.cutout(
                position=obs.pointing_radec, width=2 * self.offset_max, mode="trim"
                )

            if sum_over_axis:
                maps_obs = self._maps_sum_over_axes(maps_obs, spectrum, keepdims)
                maps_obs["exclusion"] = maps_obs["exclusion"].sum_over_axes(
                    keepdims=keepdims
                )
                maps_obs["exclusion"].data = (
                    maps_obs["exclusion"].data / self.geom.axes[0].nbin
                )

            maps_obs_bkg = self.background_estimator.run(maps_obs)
            maps_obs.update(maps_obs_bkg)
            maps_obs["exposure_on"] = maps_obs.pop("background")
            maps_obs["on"] = maps_obs.pop("counts")

            # Now paste the returned maps on the ref geom
            for name in selection:
                data = maps_obs[name].quantity.to_value(maps[name].unit)
                maps[name].fill_by_coord(maps_obs[name].geom.get_coord(), data)

        self._maps = maps
        return maps

    def run_images(self, observations, spectrum=None, keepdims=False):
        """Run image making.

        The maps are summed over on the energy axis for a classical image analysis.

        Parameters
        ----------
        observations : `~gammapy.data.Observations`
            Observations to process
        spectrum : `~gammapy.modeling.models.SpectralModel`, optional
            Spectral model to compute the weights.
            Default is power-law with spectral index of 2.
        keepdims : bool, optional
            If this is set to True, the energy axes is kept with a single bin.
            If False, the energy axes is removed

        Returns
        -------
        maps : dict of `~gammapy.maps.Map`
            Dictionary containing the following maps:

            * ``"on"``: counts map
            * ``"exposure_on"``: on exposure map, which is just the
              template background map from the IRF
            * ``"exposure_off"``: off exposure map convolved with the ring
            * ``"off"``: off map
        """
        return self._run(
            observations, sum_over_axis=True, spectrum=spectrum, keepdims=keepdims
        )

    def run(self, observations):
        """Run map making.

        Parameters
        ----------
        observations : `~gammapy.data.Observations`
            Observations to process

        Returns
        -------
        maps : dict of `~gammapy.maps.Map`
            Dictionary containing the following maps:

            * ``"on"``: counts map
            * ``"exposure_on"``: on exposure map, which is just the
              template background map from the IRF
            * ``"exposure_off"``: off exposure map convolved with the ring
            * ``"off"``: off map
        """
        return self._run(observations, sum_over_axis=False)


1			# Licensed under a 3-clause BSD style license - see LICENSE.rst
2			import logging
3			from functools import lru_cache
4			import numpy as np
5			from astropy.coordinates import Angle
6			from astropy.nddata.utils import NoOverlapError
7			from astropy.utils import lazyproperty
8			from gammapy.irf import EnergyDependentMultiGaussPSF
9			from gammapy.maps import Map, WcsGeom
10			from gammapy.modeling.models import BackgroundModel
11			from .background import make_map_background_irf
12			from .counts import fill_map_counts
13			from .edisp_map import make_edisp_map
14			from .exposure import _map_spectrum_weight, make_map_exposure_true_energy
15			from .fit import MIGRA_AXIS_DEFAULT, RAD_AXIS_DEFAULT, BINSZ_IRF, MapDataset
16			from .psf_map import make_psf_map
17
18			__all__ = ["MapDatasetMaker", "MapMakerRing"]
19
20			log = logging.getLogger(__name__)
21
22
23			class MapDatasetMaker:
24			"""Make maps for a single IACT observation.
25
26			Parameters
27			----------
28			geom : `~gammapy.maps.WcsGeom`
29			Reference image geometry in reco energy, used for counts and background maps
30			offset_max : `~astropy.coordinates.Angle`
31			Maximum offset angle
32			geom_true : `~gammapy.maps.WcsGeom`
33			Reference image geometry in true energy, used for IRF maps. It can have a coarser
34			spatial bins than the counts geom.
35			If none, the same as geom is assumed
36			migra_axis : `~gammapy.maps.MapAxis`
37			Migration axis for edisp map
38			rad_axis : `~gammapy.maps.MapAxis`
39			Radial axis for psf map.
40			cutout : bool
41			Whether to cutout the observation.
42			cutout_mode : {'trim', 'partial', 'strict'}
43			Mode option for cutting out the observation,
44			for details see `~astropy.nddata.utils.Cutout2D`.
45			"""
46
47			def __init__(
48			self,
49			geom,
50			offset_max,
51			geom_true=None,
52			background_oversampling=None,
53			migra_axis=None,
54			rad_axis=None,
55			cutout_mode="trim",
56			cutout=True
57			):
58			self.geom = geom
59			self.geom_true = geom_true if geom_true else geom.to_binsz(BINSZ_IRF)
60			self.offset_max = Angle(offset_max)
61			self.background_oversampling = background_oversampling
62			self.migra_axis = migra_axis if migra_axis else MIGRA_AXIS_DEFAULT
63			self.rad_axis = rad_axis if rad_axis else RAD_AXIS_DEFAULT
64			self.cutout_mode = cutout_mode
65			self.cutout_width = 2 * self.offset_max
66			self.cutout = cutout
67
68			def _cutout_geom(self, geom, observation):
69			if self.cutout:
70			cutout_kwargs = {
71			"position": observation.pointing_radec,
72			"width": self.cutout_width,
73			"mode": self.cutout_mode,
74			}
75			return geom.cutout(**cutout_kwargs)
76			else:
77			return geom
78
79			@lazyproperty
80			def geom_exposure(self):
81			"""Exposure map geom (`Geom`)"""
82			energy_axis = self.geom_true.get_axis_by_name("energy")
83			geom_exposure = self.geom.to_image().to_cube([energy_axis])
84			return geom_exposure
85
86			@lazyproperty
87			def geom_psf(self):
88			"""PSFMap geom (`Geom`)"""
89			energy_axis = self.geom_true.get_axis_by_name("ENERGY")
90			geom_psf = self.geom_true.to_image().to_cube([self.rad_axis, energy_axis])
91			return geom_psf
92
93			@lazyproperty
94			def geom_edisp(self):
95			"""EdispMap geom (`Geom`)"""
96			energy_axis = self.geom_true.get_axis_by_name("ENERGY")
97			geom_edisp = self.geom_true.to_image().to_cube([self.migra_axis, energy_axis])
98			return geom_edisp
99
100			def make_counts(self, observation):
101			"""Make counts map.
102
103			Parameters
104			----------
105			observation : `DataStoreObservation`
106			Observation container.
107
108			Returns
109			-------
110			counts : `Map`
111			Counts map.
112			"""
113			geom = self._cutout_geom(self.geom, observation)
114			counts = Map.from_geom(geom)
115			fill_map_counts(counts, observation.events)
116			return counts
117
118			def make_exposure(self, observation):
119			"""Make exposure map.
120
121			Parameters
122			----------
123			observation : `DataStoreObservation`
124			Observation container.
125
126			Returns
127			-------
128			exposure : `Map`
129			Exposure map.
130			"""
131			geom = self._cutout_geom(self.geom_exposure, observation)
132			exposure = make_map_exposure_true_energy(
133			pointing=observation.pointing_radec,
134			livetime=observation.observation_live_time_duration,
135			aeff=observation.aeff,
136			geom=geom,
137			)
138			return exposure
139
140			@lru_cache(maxsize=1)
141			def make_exposure_irf(self, observation):
142			"""Make exposure map with irf geometry.
143
144			Parameters
145			----------
146			observation : `DataStoreObservation`
147			Observation container.
148
149			Returns
150			-------
151			exposure : `Map`
152			Exposure map.
153			"""
154			geom = self._cutout_geom(self.geom_true, observation)
155			exposure = make_map_exposure_true_energy(
156			pointing=observation.pointing_radec,
157			livetime=observation.observation_live_time_duration,
158			aeff=observation.aeff,
159			geom=geom,
160			)
161			return exposure
162
163			def make_background(self, observation):
164			"""Make background map.
165
166			Parameters
167			----------
168			observation : `DataStoreObservation`
169			Observation container.
170
171			Returns
172			-------
173			background : `Map`
174			Background map.
175			"""
176			geom = self._cutout_geom(self.geom, observation)
177
178			bkg_coordsys = observation.bkg.meta.get("FOVALIGN", "ALTAZ")
179			if bkg_coordsys == "ALTAZ":
180			pointing = observation.fixed_pointing_info
181			elif bkg_coordsys == "RADEC":
182			pointing = observation.pointing_radec
183			else:
184			raise ValueError(
185			f"Invalid background coordinate system: {bkg_coordsys!r}\n"
186			"Options: ALTAZ, RADEC"
187			)
188			background = make_map_background_irf(
189			pointing=pointing,
190			ontime=observation.observation_time_duration,
191			bkg=observation.bkg,
192			geom=geom,
193			oversampling=self.background_oversampling,
194			)
195			return background
196
197			def make_edisp(self, observation):
198			"""Make edisp map.
199
200			Parameters
201			----------
202			observation : `DataStoreObservation`
203			Observation container.
204
205			Returns
206			-------
207			edisp : `EdispMap`
208			Edisp map.
209			"""
210			geom = self._cutout_geom(self.geom_edisp, observation)
211
212			exposure = self.make_exposure_irf(observation)
213
214			edisp = make_edisp_map(
215			edisp=observation.edisp,
216			pointing=observation.pointing_radec,
217			geom=geom,
218			max_offset=self.offset_max,
219			exposure_map=exposure,
220			)
221			return edisp
222
223			def make_psf(self, observation):
224			"""Make psf map.
225
226			Parameters
227			----------
228			observation : `DataStoreObservation`
229			Observation container.
230
231			Returns
232			-------
233			psf : `PSFMap`
234			Psf map.
235			"""
236			psf = observation.psf
237			geom = self._cutout_geom(self.geom_psf, observation)
238
239			if isinstance(psf, EnergyDependentMultiGaussPSF):
240			psf = psf.to_psf3d(rad=self.rad_axis.center)
241
242			exposure = self.make_exposure_irf(observation)
243
244			psf = make_psf_map(
245			psf=psf,
246			pointing=observation.pointing_radec,
247			geom=geom,
248			max_offset=self.offset_max,
249			exposure_map=exposure,
250			)
251			return psf
252
253			@lru_cache(maxsize=1)
254			def make_mask_safe(self, observation):
255			"""Make offset mask.
256
257			Parameters
258			----------
259			observation : `DataStoreObservation`
260			Observation container.
261
262			Returns
263			-------
264			mask : `~numpy.ndarray`
265			Mask
266			"""
267			geom = self._cutout_geom(self.geom, observation)
268			offset = geom.separation(observation.pointing_radec)
269			return offset >= self.offset_max
270
271			@lru_cache(maxsize=1)
272			def make_mask_safe_irf(self, observation):
273			"""Make offset mask with irf geometry.
274
275			Parameters
276			----------
277			observation : `DataStoreObservation`
278			Observation container.
279
280			Returns
281			-------
282			mask : `~numpy.ndarray`
283			Mask
284			"""
285			geom = self._cutout_geom(self.geom_true, observation)
286			offset = geom.separation(observation.pointing_radec)
287			return offset >= self.offset_max
288
289			def run(self, observation, selection=None):
290			"""Make map dataset.
291
292			Parameters
293			----------
294			observation : `~gammapy.data.DataStoreObservation`
295			Observation
296			selection : list
297			List of str, selecting which maps to make.
298			Available: 'counts', 'exposure', 'background', 'psf', 'edisp'
299			By default, all maps are made.
300
301			Returns
302			-------
303			dataset : `MapDataset`
304			Map dataset.
305
306			"""
307			selection = _check_selection(selection)
308
309			kwargs = {}
310			kwargs["name"] = "obs_{}".format(observation.obs_id)
311			kwargs["gti"] = observation.gti
312
313			mask_safe = self.make_mask_safe(observation)
314
315			# expand mask safe into 3rd dimension
316			nbin = self.geom.get_axis_by_name("energy").nbin
317			kwargs["mask_safe"] = ~mask_safe & np.ones(nbin, dtype=bool)[:, np.newaxis, np.newaxis]
318
319			mask_safe_irf = self.make_mask_safe_irf(observation)
320
321			if "counts" in selection:
322			counts = self.make_counts(observation)
323			# TODO: remove masking out the values here and instead handle the safe mask only when
324			# fitting and / or stacking datasets?
325			counts.data[..., mask_safe] = 0
326			kwargs["counts"] = counts
327
328			if "exposure" in selection:
329			exposure = self.make_exposure(observation)
330			exposure.data[..., mask_safe] = 0
331			kwargs["exposure"] = exposure
332
333			if "background" in selection:
334			background_map = self.make_background(observation)
335			background_map.data[..., mask_safe] = 0
336			kwargs["background_model"] = BackgroundModel(background_map)
337
338			if "psf" in selection:
339			psf = self.make_psf(observation)
340			psf.exposure_map.data[..., mask_safe_irf] = 0
341			kwargs["psf"] = psf
342
343			if "edisp" in selection:
344			edisp = self.make_edisp(observation)
345			psf.exposure_map.data[..., mask_safe_irf] = 0
			0 ignored issues – show introduced 2019-10-10 20:20 UTC by Report Bug Copy Issue Report The variable `psf` does not seem to be defined in case `"psf" in selection` on line `338` is `False`. Are you sure this can never be the case? Loading history...
346			kwargs["edisp"] = edisp
347
348			return MapDataset(**kwargs)
349
350
351			def _check_selection(selection):
352			"""Handle default and validation of selection"""
353			available = ["counts", "exposure", "background", "psf", "edisp"]
354			if selection is None:
355			selection = available
356
357			if not isinstance(selection, list):
358			raise TypeError("Selection must be a list of str")
359
360			for name in selection:
361			if name not in available:
362			raise ValueError(f"Selection not available: {name!r}")
363
364			return selection
365
366
367			class MapMakerRing:
368			"""Make maps from IACT observations.
369
370			The main motivation for this class in addition to the `MapMaker`
371			is to have the common image background estimation methods,
372			like `~gammapy.cube.RingBackgroundEstimator`,
373			that work using on and off maps.
374
375			To ensure adequate statistics, only observations that are fully
376			contained within the reference geometry will be analysed
377
378			Parameters
379			----------
380			geom : `~gammapy.maps.WcsGeom`
381			Reference image geometry
382			offset_max : `~astropy.coordinates.Angle`
383			Maximum offset angle
384			exclusion_mask : `~gammapy.maps.Map`
385			Exclusion mask
386			background_estimator : `~gammapy.cube.RingBackgroundEstimator`
387			or `~gammapy.cube.AdaptiveRingBackgroundEstimator`
388			Ring background estimator or something with an equivalent API.
389
390			Examples
391			--------
392			Here is an example how to ise the MapMakerRing with H.E.S.S. DL3 data::
393
394			import numpy as np
395			import astropy.units as u
396			from astropy.coordinates import SkyCoord
397			from regions import CircleSkyRegion
398			from gammapy.maps import Map, WcsGeom, MapAxis
399			from gammapy.cube import MapMakerRing, RingBackgroundEstimator
400			from gammapy.data import DataStore
401
402			# Create observation list
403			data_store = DataStore.from_file(
404			"$GAMMAPY_DATA/hess-dl3-dr1/hess-dl3-dr3-with-background.fits.gz"
405			)
406			data_sel = data_store.obs_table["TARGET_NAME"] == "MSH 15-52"
407			obs_table = data_store.obs_table[data_sel]
408			observations = data_store.get_observations(obs_table["OBS_ID"])
409
410			# Define the geom
411			pos = SkyCoord(228.32, -59.08, unit="deg")
412			energy_axis = MapAxis.from_edges(np.logspace(0, 5.0, 5), unit="TeV", name="energy")
413			geom = WcsGeom.create(skydir=pos, binsz=0.02, width=(5, 5), axes=[energy_axis])
414
415			# Make a region mask
416			regions = CircleSkyRegion(center=pos, radius=0.3 * u.deg)
417			mask = Map.from_geom(geom)
418			mask.data = mask.geom.region_mask([regions], inside=False)
419
420			# Run map maker with ring background estimation
421			ring_bkg = RingBackgroundEstimator(r_in="0.5 deg", width="0.3 deg")
422			maker = MapMakerRing(
423			geom=geom, offset_max="2 deg", exclusion_mask=mask, background_estimator=ring_bkg
424			)
425			images = maker.run_images(observations)
426			"""
427
428			def __init__(
429			self, geom, offset_max, exclusion_mask=None, background_estimator=None
430			):
431
432			self.geom = geom
433			self.offset_max = Angle(offset_max)
434			self.exclusion_mask = exclusion_mask
435			self.background_estimator = background_estimator
436
437			def _get_empty_maps(self, selection):
438			# Initialise zero-filled maps
439			maps = {}
440			for name in selection:
441			if name == "exposure":
442			maps[name] = Map.from_geom(self.geom, unit="m2 s")
443			else:
444			maps[name] = Map.from_geom(self.geom, unit="")
445			return maps
446
447			def _get_obs_maker(self, obs):
448			# Compute cutout geometry and slices to stack results back later
449
450			# Make maps for this observation
451			return MapDatasetMaker(
452			geom=self.geom,
453			offset_max=self.offset_max,
454			)
455
456			@staticmethod
457			def _maps_sum_over_axes(maps, spectrum, keepdims):
458			"""Compute weighted sum over map axes.
459
460			Parameters
461			----------
462			spectrum : `~gammapy.modeling.models.SpectralModel`
463			Spectral model to compute the weights.
464			Default is power-law with spectral index of 2.
465			keepdims : bool, optional
466			If this is set to True, the energy axes is kept with a single bin.
467			If False, the energy axes is removed
468			"""
469			images = {}
470			for name, map in maps.items():
471			if name == "exposure":
472			map = _map_spectrum_weight(map, spectrum)
473			images[name] = map.sum_over_axes(keepdims=keepdims)
474			# TODO: PSF (and edisp) map sum_over_axis
475
476			return images
477
478
479			def _run(self, observations, sum_over_axis=False, spectrum=None, keepdims=False):
480			selection = ["on", "exposure_on", "off", "exposure_off", "exposure"]
481			maps = self._get_empty_maps(selection)
482			if sum_over_axis:
483			maps = self._maps_sum_over_axes(maps, spectrum, keepdims)
484
485			for obs in observations:
486			try:
487			obs_maker = self._get_obs_maker(obs)
488			except NoOverlapError:
489			log.info(f"Skipping obs_id: {obs.obs_id} (no map overlap)")
490			continue
491
492			dataset = obs_maker.run(obs, selection=["counts", "exposure", "background"])
493			maps_obs = {}
494			maps_obs["counts"] = dataset.counts
495			maps_obs["exposure"] = dataset.exposure
496			maps_obs["background"] = dataset.background_model.map
497			maps_obs["exclusion"] = self.exclusion_mask.cutout(
498			position=obs.pointing_radec, width=2 * self.offset_max, mode="trim"
499			)
500
501			if sum_over_axis:
502			maps_obs = self._maps_sum_over_axes(maps_obs, spectrum, keepdims)
503			maps_obs["exclusion"] = maps_obs["exclusion"].sum_over_axes(
504			keepdims=keepdims
505			)
506			maps_obs["exclusion"].data = (
507			maps_obs["exclusion"].data / self.geom.axes[0].nbin
508			)
509
510			maps_obs_bkg = self.background_estimator.run(maps_obs)
511			maps_obs.update(maps_obs_bkg)
512			maps_obs["exposure_on"] = maps_obs.pop("background")
513			maps_obs["on"] = maps_obs.pop("counts")
514
515			# Now paste the returned maps on the ref geom
516			for name in selection:
517			data = maps_obs[name].quantity.to_value(maps[name].unit)
518			maps[name].fill_by_coord(maps_obs[name].geom.get_coord(), data)
519
520			self._maps = maps
521			return maps
522
523			def run_images(self, observations, spectrum=None, keepdims=False):
524			"""Run image making.
525
526			The maps are summed over on the energy axis for a classical image analysis.
527
528			Parameters
529			----------
530			observations : `~gammapy.data.Observations`
531			Observations to process
532			spectrum : `~gammapy.modeling.models.SpectralModel`, optional
533			Spectral model to compute the weights.
534			Default is power-law with spectral index of 2.
535			keepdims : bool, optional
536			If this is set to True, the energy axes is kept with a single bin.
537			If False, the energy axes is removed
538
539			Returns
540			-------
541			maps : dict of `~gammapy.maps.Map`
542			Dictionary containing the following maps:
543
544			* ``"on"``: counts map
545			* ``"exposure_on"``: on exposure map, which is just the
546			template background map from the IRF
547			* ``"exposure_off"``: off exposure map convolved with the ring
548			* ``"off"``: off map
549			"""
550			return self._run(
551			observations, sum_over_axis=True, spectrum=spectrum, keepdims=keepdims
552			)
553
554			def run(self, observations):
555			"""Run map making.
556
557			Parameters
558			----------
559			observations : `~gammapy.data.Observations`
560			Observations to process
561
562			Returns
563			-------
564			maps : dict of `~gammapy.maps.Map`
565			Dictionary containing the following maps:
566
567			* ``"on"``: counts map
568			* ``"exposure_on"``: on exposure map, which is just the
569			template background map from the IRF
570			* ``"exposure_off"``: off exposure map convolved with the ring
571			* ``"off"``: off map
572			"""
573			return self._run(observations, sum_over_axis=False)
574

gammapy / gammapy

Pull Request — master (#2450)

MapDatasetMaker.make_exposure_irf() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like