# Licensed under a 3-clause BSD style license - see LICENSE.rstimportloggingimportnumpyasnpimportastropy.unitsasufromastropy.coordinatesimportangular_separationfromastropy.visualizationimportquantity_supportimportmatplotlib.pyplotaspltfrommatplotlib.colorsimportLogNormfromgammapy.mapsimportMapAxes,MapAxisfromgammapy.maps.axesimportUNIT_STRING_FORMATfromgammapy.visualization.utilsimportadd_colorbarfrom.coreimportIRFfrom.ioimportgadf_is_pointlike__all__=["Background3D","Background2D","BackgroundIRF"]log=logging.getLogger(__name__)
[docs]classBackgroundIRF(IRF):"""Background IRF base class. Parameters ---------- axes : list of `MapAxis` or `MapAxes` object data : `~np.ndarray` Data array. unit : str or `~astropy.units.Unit` Data unit usually ``s^-1 MeV^-1 sr^-1``. meta : dict Metadata dictionary. """default_interp_kwargs=dict(bounds_error=False,fill_value=0.0,values_scale="log")"""Default Interpolation kwargs to extrapolate."""
[docs]@classmethoddeffrom_table(cls,table,format="gadf-dl3"):"""Read from `~astropy.table.Table`. Parameters ---------- table : `~astropy.table.Table` Table with background data. format : {"gadf-dl3"} Format specification. Default is "gadf-dl3". Returns ------- bkg : `Background2D` or `Background2D` Background IRF class. """# TODO: some of the existing background files have missing HDUCLAS keywords# which are required to define the correct Gammapy axis namesif"HDUCLAS2"notintable.meta:log.warning("Missing 'HDUCLAS2' keyword assuming 'BKG'")table=table.copy()table.meta["HDUCLAS2"]="BKG"axes=MapAxes.from_table(table,format=format)[cls.required_axes]# TODO: spec says key should be "BKG", but there are files around# (e.g. CTA 1DC) that use "BGD". For now we support bothif"BKG"intable.colnames:bkg_name="BKG"elif"BGD"intable.colnames:bkg_name="BGD"else:raiseValueError("Invalid column names. Need 'BKG' or 'BGD'.")data=table[bkg_name].quantity[0].Tifdata.unit==""orisinstance(data.unit,u.UnrecognizedUnit):data=u.Quantity(data.value,"s-1 MeV-1 sr-1",copy=False)log.warning("Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)")# TODO: The present HESS and CTA background fits files# have a reverse order (lon, lat, E) than recommended in GADF(E, lat, lon)# For now, we support both.ifaxes.shape==axes.shape[::-1]:log.error("Ambiguous axes order in Background fits files!")ifnp.shape(data)!=axes.shape:log.debug("Transposing background table on read")data=data.transpose()returncls(axes=axes,data=data.value,meta=table.meta,unit=data.unit,is_pointlike=gadf_is_pointlike(table.meta),fov_alignment=table.meta.get("FOVALIGN","RADEC"),)
[docs]classBackground3D(BackgroundIRF):"""Background 3D. Data format specification: :ref:`gadf:bkg_3d`. Parameters ---------- axes : list of `MapAxis` or `MapAxes` object Required data axes: ["energy", "fov_lon", "fov_lat"] in the given order. data : `~np.ndarray` Data array. unit : str or `~astropy.units.Unit` Data unit usually ``s^-1 MeV^-1 sr^-1``. fov_alignment : `~gammapy.irf.FoVAlignment` The orientation of the field of view coordinate system. meta : dict Metadata dictionary. Examples -------- Here's an example you can use to learn about this class: >>> from gammapy.irf import Background3D >>> filename = '$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits' >>> bkg_3d = Background3D.read(filename, hdu='BACKGROUND') >>> print(bkg_3d) Background3D ------------ <BLANKLINE> axes : ['energy', 'fov_lon', 'fov_lat'] shape : (21, 36, 36) ndim : 3 unit : 1 / (MeV s sr) dtype : >f4 <BLANKLINE> """tag="bkg_3d"required_axes=["energy","fov_lon","fov_lat"]default_unit=u.s**-1*u.MeV**-1*u.sr**-1
[docs]defto_2d(self):"""Convert to `Background2D`. This takes the values at Y = 0 and X >= 0. """# TODO: this is incorrect as it misses the Jacobian?idx_lon=self.axes["fov_lon"].coord_to_idx(0*u.deg)[0]idx_lat=self.axes["fov_lat"].coord_to_idx(0*u.deg)[0]data=self.quantity[:,idx_lon:,idx_lat].copy()offset=self.axes["fov_lon"].edges[idx_lon:]offset_axis=MapAxis.from_edges(offset,name="offset")returnBackground2D(axes=[self.axes["energy"],offset_axis],data=data.value,unit=data.unit)
[docs]defpeek(self,figsize=(10,8)):"""Quick-look summary plots. Parameters ---------- figsize : tuple, optional Size of the figure. Default is (10, 8). """returnself.to_2d().peek(figsize)
[docs]defplot_at_energy(self,energy=1*u.TeV,add_cbar=True,ncols=3,figsize=None,axes_loc=None,kwargs_colorbar=None,**kwargs,):"""Plot the background rate in FoV coordinates at a given energy. Parameters ---------- energy : `~astropy.units.Quantity`, optional List of energies. Default is 1 TeV. add_cbar : bool, optional Add color bar. Default is True. ncols : int, optional Number of columns to plot. Default is 3. figsize : tuple, optional Figure size. Default is None. axes_loc : dict, optional Keyword arguments passed to `~mpl_toolkits.axes_grid1.axes_divider.AxesDivider.append_axes`. kwargs_colorbar : dict, optional Keyword arguments passed to `~matplotlib.pyplot.colorbar`. **kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.pcolormesh`. """kwargs_colorbar=kwargs_colorbaror{}n=len(energy)cols=min(ncols,n)rows=1+(n-1)//colswidth=12cfraction=0.0ifadd_cbar:cfraction=0.15iffigsizeisNone:figsize=(width,rows*width//(cols*(1+cfraction)))fig,axes=plt.subplots(ncols=cols,nrows=rows,figsize=figsize,gridspec_kw={"hspace":0.2,"wspace":0.3},)x=self.axes["fov_lat"].edgesy=self.axes["fov_lon"].edgesX,Y=np.meshgrid(x,y)fori,eeinenumerate(energy):iflen(energy)==1:ax=axeselse:ax=axes.flat[i]bkg=self.evaluate(energy=ee)bkg_unit=bkg.unitbkg=bkg.valuewithquantity_support():caxes=ax.pcolormesh(X,Y,bkg.squeeze(),**kwargs)self.axes["fov_lat"].format_plot_xaxis(ax)self.axes["fov_lon"].format_plot_yaxis(ax)ax.set_title(str(ee))ifadd_cbar:label=f"Background [{bkg_unit.to_string(UNIT_STRING_FORMAT)}]"kwargs_colorbar.setdefault("label",label)cbar=add_colorbar(caxes,ax=ax,axes_loc=axes_loc,**kwargs_colorbar)cbar.formatter.set_powerlimits((0,0))row,col=np.unravel_index(i,shape=(rows,cols))ifcol>0:ax.set_ylabel("")ifrow<rows-1:ax.set_xlabel("")ax.set_aspect("equal","box")
[docs]classBackground2D(BackgroundIRF):"""Background 2D. Data format specification: :ref:`gadf:bkg_2d` Parameters ---------- axes : list of `MapAxis` or `MapAxes` object Required data axes: ["energy", "offset"] in the given order. data : `~np.ndarray` Data array. unit : str or `~astropy.units.Unit` Data unit usually ``s^-1 MeV^-1 sr^-1``. meta : dict Metadata dictionary. """tag="bkg_2d"required_axes=["energy","offset"]default_unit=u.s**-1*u.MeV**-1*u.sr**-1
[docs]defto_3d(self):"""Convert to Background3D."""offsets=self.axes["offset"].edgesedges_neg=np.negative(offsets)[::-1]edges_neg=edges_neg[edges_neg<=0]edges=np.concatenate((edges_neg,offsets[offsets>0]))fov_lat=MapAxis.from_edges(edges=edges,name="fov_lat")fov_lon=MapAxis.from_edges(edges=edges,name="fov_lon")axes=MapAxes([self.axes["energy"],fov_lon,fov_lat])coords=axes.get_coord()offset=angular_separation(0*u.rad,0*u.rad,coords["fov_lon"],coords["fov_lat"])data=self.evaluate(offset=offset,energy=coords["energy"])returnBackground3D(axes=axes,data=data,)
[docs]defplot_at_energy(self,energy=1*u.TeV,add_cbar=True,ncols=3,figsize=None,**kwargs):"""Plot the background rate in FoV coordinates at a given energy. Parameters ---------- energy : `~astropy.units.Quantity`, optional List of energy. Default is 1 TeV. add_cbar : bool, optional Add color bar. Default is True. ncols : int, optional Number of columns to plot. Default is 3. figsize : tuple, optional Figure size. Default is None. **kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.pcolormesh`. """bkg_3d=self.to_3d()bkg_3d.plot_at_energy(energy=energy,add_cbar=add_cbar,ncols=ncols,figsize=figsize,**kwargs)
[docs]defplot(self,ax=None,add_cbar=True,axes_loc=None,kwargs_colorbar=None,**kwargs):"""Plot energy offset dependence of the background model. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. add_cbar : bool, optional Add a colorbar to the plot. Default is True. axes_loc : dict, optional Keyword arguments passed to `~mpl_toolkits.axes_grid1.axes_divider.AxesDivider.append_axes`. kwargs_colorbar : dict, optional Keyword arguments passed to `~matplotlib.pyplot.colorbar`. kwargs : dict Keyword arguments passed to `~matplotlib.pyplot.pcolormesh`. Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxenergy_axis,offset_axis=self.axes["energy"],self.axes["offset"]data=self.quantity.valuekwargs.setdefault("cmap","GnBu")kwargs.setdefault("edgecolors","face")kwargs.setdefault("norm",LogNorm())kwargs_colorbar=kwargs_colorbaror{}withquantity_support():caxes=ax.pcolormesh(energy_axis.edges,offset_axis.edges,data.T,**kwargs)energy_axis.format_plot_xaxis(ax=ax)offset_axis.format_plot_yaxis(ax=ax)ifadd_cbar:label=(f"Background rate [{self.quantity.unit.to_string(UNIT_STRING_FORMAT)}]")kwargs_colorbar.setdefault("label",label)add_colorbar(caxes,ax=ax,axes_loc=axes_loc,**kwargs_colorbar)
[docs]defplot_offset_dependence(self,ax=None,energy=None,**kwargs):"""Plot background rate versus offset for a given energy. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. energy : `~astropy.units.Quantity`, optional Energy. Default is None. Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxifenergyisNone:energy_axis=self.axes["energy"]e_min,e_max=np.log10(energy_axis.center.value[[0,-1]])energy=np.logspace(e_min,e_max,4)*energy_axis.unitoffset_axis=self.axes["offset"]foreeinenergy:bkg=self.evaluate(offset=offset_axis.center,energy=ee)ifnp.isnan(bkg).all():continuelabel=f"energy = {ee:.1f}"withquantity_support():ax.plot(offset_axis.center,bkg,label=label,**kwargs)offset_axis.format_plot_xaxis(ax=ax)ax.set_ylabel(f"Background rate [{ax.yaxis.units.to_string(UNIT_STRING_FORMAT)}]")ax.set_yscale("log")ax.legend(loc="upper right")returnax
[docs]defplot_energy_dependence(self,ax=None,offset=None,**kwargs):"""Plot background rate versus energy for a given offset. Parameters ---------- ax : `~matplotlib.axes.Axes`, optional Matplotlib axes. Default is None. offset : `~astropy.coordinates.Angle`, optional Offset. Default is None. kwargs : dict Forwarded to plt.plot(). Returns ------- ax : `~matplotlib.axes.Axes` Matplotlib axes. """ax=plt.gca()ifaxisNoneelseaxifoffsetisNone:offset_axis=self.axes["offset"]off_min,off_max=offset_axis.center.value[[0,-1]]offset=np.linspace(off_min,off_max,4)*offset_axis.unitenergy_axis=self.axes["energy"]foroffinoffset:bkg=self.evaluate(offset=off,energy=energy_axis.center)label=f"offset = {off:.2f}"withquantity_support():ax.plot(energy_axis.center,bkg,label=label,**kwargs)energy_axis.format_plot_xaxis(ax=ax)ax.set_yscale("log")ax.set_ylabel(f"Background rate [{ax.yaxis.units.to_string(UNIT_STRING_FORMAT)}]")ax.legend(loc="best")returnax