# Licensed under a 3-clause BSD style license - see LICENSE.rstimporthtmlimportnumpyasnpimportastropy.unitsasuimportmatplotlib.pyplotaspltfromgammapy.mapsimportMapfromgammapy.modeling.modelsimportPowerLawSpectralModel__all__=["PSFKernel"]
[docs]classPSFKernel:"""PSF kernel for `~gammapy.maps.Map`. This is a container class to store a PSF kernel that can be used to convolve `~gammapy.maps.WcsNDMap` objects. It is usually computed from an `~gammapy.irf.PSFMap`. Parameters ---------- psf_kernel_map : `~gammapy.maps.Map` PSF kernel stored in a Map. Examples -------- :: from gammapy.maps import Map, WcsGeom, MapAxis from gammapy.irf import PSFMap from astropy import units as u # Define energy axis energy_axis_true = MapAxis.from_energy_bounds( "0.1 TeV", "10 TeV", nbin=4, name="energy_true ) # Create WcsGeom and map geom = WcsGeom.create(binsz=0.02 * u.deg, width=2.0 * u.deg, axes=[energy_axis_true]) some_map = Map.from_geom(geom) # Fill map at three positions some_map.fill_by_coord([[0, 0, 0], [0, 0, 0], [0.3, 1, 3]]) psf = PSFMap.from_gauss( energy_axis_true=energy_axis_true, sigma=[0.3, 0.2, 0.1] * u.deg ) kernel = psf.get_psf_kernel(geom=geom, max_radius=1*u.deg) # Do the convolution some_map_convolved = some_map.convolve(kernel) some_map_convolved.plot_grid(); """def__init__(self,psf_kernel_map,normalize=True):self._psf_kernel_map=psf_kernel_mapifnormalize:self.normalize()def_repr_html_(self):try:returnself.to_html()exceptAttributeError:returnf"<pre>{html.escape(str(self))}</pre>"
[docs]defnormalize(self):"""Force normalisation of the kernel."""data=self.psf_kernel_map.dataifself.psf_kernel_map.geom.is_image:axis=(0,1)else:axis=(1,2)withnp.errstate(divide="ignore",invalid="ignore"):data=np.nan_to_num(data/data.sum(axis=axis,keepdims=True))self.psf_kernel_map.data=data
@propertydefdata(self):"""Access the PSFKernel numpy array."""returnself._psf_kernel_map.data@propertydefpsf_kernel_map(self):"""The map object holding the kernel as a `~gammapy.maps.Map`."""returnself._psf_kernel_map
[docs]@classmethoddefread(cls,*args,**kwargs):"""Read kernel Map from file."""psf_kernel_map=Map.read(*args,**kwargs)returncls(psf_kernel_map)
[docs]@classmethoddeffrom_spatial_model(cls,model,geom,max_radius=None,factor=4):"""Create PSF kernel from spatial model. Parameters ---------- geom : `~gammapy.maps.WcsGeom` Map geometry. model : `~gammapy.modeling.models.SpatiaModel` Gaussian width. max_radius : `~astropy.coordinates.Angle`, optional Desired kernel map size. Default is None. factor : int, optional Oversample factor to compute the PSF. Default is 4. Returns ------- kernel : `~gammapy.irf.PSFKernel` The kernel Map with reduced geometry according to the max_radius. """ifmax_radiusisNone:max_radius=model.evaluation_radiusgeom=geom.to_odd_npix(max_radius=max_radius)model.position=geom.center_skydirgeom=geom.upsample(factor=factor)map=model.integrate_geom(geom,oversampling_factor=1)returncls(psf_kernel_map=map.downsample(factor=factor))
[docs]@classmethoddeffrom_gauss(cls,geom,sigma,max_radius=None,factor=4):"""Create Gaussian PSF. This is used for testing and examples. The map geometry parameters (pixel size, energy bins) are taken from ``geom``. The Gaussian width ``sigma`` is a scalar. TODO : support array input if it should vary along the energy axis. Parameters ---------- geom : `~gammapy.maps.WcsGeom` Map geometry. sigma : `~astropy.coordinates.Angle` Gaussian width. max_radius : `~astropy.coordinates.Angle`, optional Desired kernel map size. Default is None. factor : int, optional Oversample factor to compute the PSF. Default is 4. Returns ------- kernel : `~gammapy.irf.PSFKernel` The kernel Map with reduced geometry according to the max_radius. """fromgammapy.modeling.modelsimportGaussianSpatialModelgauss=GaussianSpatialModel(sigma=sigma)returncls.from_spatial_model(model=gauss,geom=geom,max_radius=max_radius,factor=factor)
[docs]defwrite(self,*args,**kwargs):"""Write the Map object which contains the PSF kernel to file."""self.psf_kernel_map.write(*args,**kwargs)
[docs]defto_image(self,spectrum=None,exposure=None,keepdims=True):"""Transform 3D PSFKernel into a 2D PSFKernel. Parameters ---------- spectrum : `~gammapy.modeling.models.SpectralModel`, optional Spectral model to compute the weights. Default is power-law with spectral index of 2. exposure : `~astropy.units.Quantity` or `~numpy.ndarray`, optional 1D array containing exposure in each true energy bin. It must have the same size as the PSFKernel energy axis. Default is uniform exposure over energy. keepdims : bool, optional If true, the resulting PSFKernel will keep an energy axis with one bin. Default is True. Returns ------- weighted_kernel : `~gammapy.irf.PSFKernel` The weighted kernel summed over energy. """map=self.psf_kernel_mapifspectrumisNone:spectrum=PowerLawSpectralModel(index=2.0)ifexposureisNone:exposure=np.ones(map.geom.axes[0].center.shape)exposure=u.Quantity(exposure)ifexposure.shape!=map.geom.axes[0].center.shape:raiseValueError("Incorrect exposure_array shape")# Compute weights vectorfornameinmap.geom.axes.names:if"energy"inname:energy_name=nameenergy_edges=map.geom.axes[energy_name].edgesweights=spectrum.integral(energy_min=energy_edges[:-1],energy_max=energy_edges[1:],intervals=True)weights*=exposureweights/=weights.sum()spectrum_weighted_kernel=map.copy()spectrum_weighted_kernel.quantity*=weights[:,np.newaxis,np.newaxis]returnself.__class__(spectrum_weighted_kernel.sum_over_axes(keepdims=keepdims))
[docs]defslice_by_idx(self,slices):"""Slice by index."""kernel=self.psf_kernel_map.slice_by_idx(slices=slices)returnself.__class__(psf_kernel_map=kernel)
[docs]defplot_kernel(self,ax=None,energy=None,add_cbar=False,**kwargs):"""Plot PSF kernel. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. energy : `~astropy.units.Quantity`, optional If None, the PSF kernel is summed over the energy axis. Otherwise, the kernel corresponding to the energy bin including the given energy is shown. add_cbar : bool, optional Add a colorbar. Default is False. kwargs: dict Keyword arguments passed to `~matplotlib.axes.Axes.imshow`. Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxifenergyisnotNone:kernel_map=self.psf_kernel_mapenergy_center=kernel_map.geom.axes["energy_true"].center.to(energy.unit)idx=np.argmin(np.abs(energy_center.value-energy.value))kernel_map.get_image_by_idx([idx]).plot(ax=ax,add_cbar=add_cbar,**kwargs)else:self.to_image().psf_kernel_map.plot(ax=ax,add_cbar=add_cbar,**kwargs)returnax
[docs]defpeek(self,figsize=(15,5)):"""Quick-look summary plots. Parameters ---------- figsize : tuple, optional Size of the figure. Default is (15, 5). """fig,axes=plt.subplots(nrows=1,ncols=2,figsize=figsize)axes[0].set_title("Energy-integrated PSF kernel")self.plot_kernel(ax=axes[0],add_cbar=True)axes[1].set_title("PSF kernel at 1 TeV")self.plot_kernel(ax=axes[1],add_cbar=True,energy=1*u.TeV)plt.tight_layout()