/
irf_map.py
148 lines (118 loc) · 4.57 KB
/
irf_map.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
from copy import deepcopy
import numpy as np
from astropy.io import fits
from gammapy.maps import Map
class IRFMap:
"""IRF map base class"""
def __init__(self, irf_map, exposure_map):
self._irf_map = irf_map
self.exposure_map = exposure_map
@classmethod
def from_hdulist(
cls,
hdulist,
hdu=None,
hdu_bands=None,
exposure_hdu=None,
exposure_hdu_bands=None,
):
"""Create from `~astropy.io.fits.HDUList`.
Parameters
----------
hdulist : `~astropy.fits.HDUList`
HDU list.
hdu : str
Name or index of the HDU with the IRF map.
hdu_bands : str
Name or index of the HDU with the IRF map BANDS table.
exposure_hdu : str
Name or index of the HDU with the exposure map data.
exposure_hdu_bands : str
Name or index of the HDU with the exposure map BANDS table.
Returns
-------
irf_map : `IRFMap`
IRF map.
"""
if hdu is None:
hdu = cls._hdu_name
irf_map = Map.from_hdulist(hdulist, hdu=hdu, hdu_bands=hdu_bands)
if exposure_hdu is None:
exposure_hdu = cls._hdu_name + "_exposure"
if exposure_hdu in hdulist:
exposure_map = Map.from_hdulist(
hdulist, hdu=exposure_hdu, hdu_bands=exposure_hdu_bands
)
else:
exposure_map = None
return cls(irf_map, exposure_map)
@classmethod
def read(cls, filename):
"""Read an IRF_map from file and create corresponding object"""
with fits.open(filename, memmap=False) as hdulist:
return cls.from_hdulist(hdulist)
def to_hdulist(self):
"""Convert to `~astropy.io.fits.HDUList`.
Returns
-------
hdu_list : `~astropy.io.fits.HDUList`
HDU list.
"""
hdulist = self._irf_map.to_hdulist(hdu=self._hdu_name)
exposure_hdu = self._hdu_name + "_exposure"
if self.exposure_map is not None:
new_hdulist = self.exposure_map.to_hdulist(hdu=exposure_hdu)
hdulist.extend(new_hdulist[1:])
return hdulist
def write(self, filename, overwrite=False, **kwargs):
"""Write IRF map to fits"""
hdulist = self.to_hdulist(**kwargs)
hdulist.writeto(filename, overwrite=overwrite)
def stack(self, other, weights=None):
"""Stack IRF map with another one in place.
Parameters
----------
other : `~gammapy.irf.IRFMap`
Energy dispersion map to be stacked with this one.
weights : `~gammapy.maps.Map`
Map with stacking weights.
"""
if self.exposure_map is None or other.exposure_map is None:
raise ValueError("Missing exposure map for PSFMap.stack")
cutout_info = other._irf_map.geom.cutout_info
if cutout_info is not None:
slices = cutout_info["parent-slices"]
parent_slices = Ellipsis, slices[0], slices[1]
else:
parent_slices = None
self._irf_map.data[parent_slices] *= self.exposure_map.data[parent_slices]
self._irf_map.stack(other._irf_map * other.exposure_map.data, weights=weights)
# stack exposure map
if weights and "energy" in weights.geom.axes_names:
weights = weights.slice_by_idx({"energy": slice(0, 1)})
self.exposure_map.stack(other.exposure_map, weights=weights)
with np.errstate(invalid="ignore"):
self._irf_map.data[parent_slices] /= self.exposure_map.data[parent_slices]
self._irf_map.data = np.nan_to_num(self._irf_map.data)
def copy(self):
"""Copy IRF map"""
return deepcopy(self)
def cutout(self, position, width, mode="trim"):
"""Cutout IRF map.
Parameters
----------
position : `~astropy.coordinates.SkyCoord`
Center position of the cutout region.
width : tuple of `~astropy.coordinates.Angle`
Angular sizes of the region in (lon, lat) in that specific order.
If only one value is passed, a square region is extracted.
mode : {'trim', 'partial', 'strict'}
Mode option for Cutout2D, for details see `~astropy.nddata.utils.Cutout2D`.
Returns
-------
cutout : `IRFMap`
Cutout IRF map.
"""
irf_map = self._irf_map.cutout(position, width, mode)
exposure_map = self.exposure_map.cutout(position, width, mode)
return self.__class__(irf_map, exposure_map=exposure_map)