# Licensed under a 3-clause BSD style license - see LICENSE.rst"""Gammacat open TeV source catalog.https://github.com/gammapy/gamma-cat"""importloggingimportnumpyasnpfromastropyimportunitsasufromastropy.tableimportTablefromgammapy.estimatorsimportFluxPointsfromgammapy.modeling.modelsimportModel,SkyModelfromgammapy.utils.scriptsimportmake_pathfrom.coreimportSourceCatalog,SourceCatalogObject,format_flux_points_table__all__=["SourceCatalogGammaCat","SourceCatalogObjectGammaCat"]log=logging.getLogger(__name__)
[docs]classSourceCatalogObjectGammaCat(SourceCatalogObject):"""One object from the gamma-cat source catalog. Catalog is represented by `~gammapy.catalog.SourceCatalogGammaCat`. """_source_name_key="common_name"def__str__(self):returnself.info()@propertydefflux_points(self):"""Flux points (`~gammapy.estimators.FluxPoints`)."""returnFluxPoints.from_table(table=self.flux_points_table,reference_model=self.sky_model(),format="gadf-sed",)
[docs]definfo(self,info="all"):"""Information string. Parameters ---------- info : {'all', 'basic', 'position, 'model'} Comma separated list of options """ifinfo=="all":info="basic,position,model"ss=""ops=info.split(",")if"basic"inops:ss+=self._info_basic()if"position"inops:ss+=self._info_position()if"model"inops:ss+=self._info_morph()ss+=self._info_spectral_fit()ss+=self._info_spectral_points()returnss
[docs]defspectral_model(self):"""Source spectral model (`~gammapy.modeling.models.SpectralModel`). Parameter errors are statistical errors only. """data=self.dataspec_type=data["spec_type"]ifspec_type=="pl":tag="PowerLawSpectralModel"pars={"amplitude":data["spec_pl_norm"],"index":data["spec_pl_index"],"reference":data["spec_pl_e_ref"],}errs={"amplitude":data["spec_pl_norm_err"],"index":data["spec_pl_index_err"],}elifspec_type=="pl2":e_max=data["spec_pl2_e_max"]DEFAULT_E_MAX=u.Quantity(1e5,"TeV")ifnp.isnan(e_max.value)ore_max.value==0:e_max=DEFAULT_E_MAXtag="PowerLaw2SpectralModel"pars={"amplitude":data["spec_pl2_flux"],"index":data["spec_pl2_index"],"emin":data["spec_pl2_e_min"],"emax":e_max,}errs={"amplitude":data["spec_pl2_flux_err"],"index":data["spec_pl2_index_err"],}elifspec_type=="ecpl":tag="ExpCutoffPowerLawSpectralModel"pars={"amplitude":data["spec_ecpl_norm"],"index":data["spec_ecpl_index"],"lambda_":1.0/data["spec_ecpl_e_cut"],"reference":data["spec_ecpl_e_ref"],}errs={"amplitude":data["spec_ecpl_norm_err"],"index":data["spec_ecpl_index_err"],"lambda_":data["spec_ecpl_e_cut_err"]/data["spec_ecpl_e_cut"]**2,}elifspec_type=="none":returnNoneelse:raiseValueError(f"Invalid spec_type: {spec_type}")model=Model.create(tag,"spectral",**pars)forname,valueinerrs.items():model.parameters[name].error=valuereturnmodel
[docs]defspatial_model(self):"""Source spatial model (`~gammapy.modeling.models.SpatialModel`). TODO: add parameter errors! """d=self.datamorph_type=d["morph_type"]pars={"lon_0":d["glon"],"lat_0":d["glat"],"frame":"galactic"}errs={"lat_0":self.data["pos_err"],"lon_0":self.data["pos_err"]/np.cos(self.data["glat"]),}ifmorph_type=="point":tag="PointSpatialModel"elifmorph_type=="gauss":# TODO: add infos back once model support elongation# pars['x_stddev'] = d['morph_sigma']# pars['y_stddev'] = d['morph_sigma']# if not np.isnan(d['morph_sigma2']):# pars['y_stddev'] = d['morph_sigma2']# if not np.isnan(d['morph_pa']):# # TODO: handle reference frame for rotation angle# pars['theta'] = Angle(d['morph_pa'], 'deg').radtag="GaussianSpatialModel"pars["sigma"]=d["morph_sigma"]elifmorph_type=="shell":tag="ShellSpatialModel"# TODO: probably we shouldn't guess a shell width here!pars["radius"]=0.8*d["morph_sigma"]pars["width"]=0.2*d["morph_sigma"]elifmorph_type=="none":returnNoneelse:raiseValueError(f"Invalid morph_type: {morph_type!r}")model=Model.create(tag,"spatial",**pars)forname,valueinerrs.items():model.parameters[name].error=valuereturnmodel
[docs]defsky_model(self):"""Source sky model (`~gammapy.modeling.models.SkyModel`)."""returnSkyModel(spatial_model=self.spatial_model(),spectral_model=self.spectral_model(),name=self.name,)
def_add_source_meta(self,table):"""Copy over some information to `table.meta`."""d=self.datam=table.metam["origin"]="Data from gamma-cat"m["source_id"]=d["source_id"]m["common_name"]=d["common_name"]m["reference_id"]=d["reference_id"]@propertydefflux_points_table(self):"""Differential flux points (`~gammapy.estimators.FluxPoints`)."""d=self.datatable=Table()table.meta["SED_TYPE"]="dnde"self._add_source_meta(table)valid=np.isfinite(d["sed_e_ref"].value)ifvalid.sum()==0:returnNonetable["e_ref"]=d["sed_e_ref"]table["e_min"]=d["sed_e_min"]table["e_max"]=d["sed_e_max"]table["dnde"]=d["sed_dnde"]table["dnde_err"]=d["sed_dnde_err"]table["dnde_errn"]=d["sed_dnde_errn"]table["dnde_errp"]=d["sed_dnde_errp"]table["dnde_ul"]=d["sed_dnde_ul"]# Only keep rows that actually contain informationtable=table[valid]forcolnameintable.colnames:ifnotnp.isfinite(table[colname]).any():table.remove_column(colname)returntable
[docs]classSourceCatalogGammaCat(SourceCatalog):"""Gammacat open TeV source catalog. See: https://github.com/gammapy/gamma-cat One source is represented by `~gammapy.catalog.SourceCatalogObjectGammaCat`. Parameters ---------- filename : str Path to the gamma-cat fits file. Examples -------- Load the catalog data: >>> import matplotlib.pyplot as plt >>> import astropy.units as u >>> from gammapy.catalog import SourceCatalogGammaCat >>> cat = SourceCatalogGammaCat() Access a source by name: >>> source = cat['Vela Junior'] Access source spectral data and plot it: >>> ax = plt.subplot() >>> energy_range = [1, 10] * u.TeV >>> source.spectral_model().plot(energy_range, ax=ax) # doctest: +SKIP >>> source.spectral_model().plot_error(energy_range, ax=ax) # doctest: +SKIP >>> source.flux_points.plot(ax=ax) # doctest: +SKIP """tag="gamma-cat"description="An open catalog of gamma-ray sources"source_object_class=SourceCatalogObjectGammaCatdef__init__(self,filename="$GAMMAPY_DATA/catalogs/gammacat/gammacat.fits.gz"):filename=make_path(filename)table=Table.read(filename,hdu=1)source_name_key="common_name"source_name_alias=("other_names","gamma_names")super().__init__(table=table,source_name_key=source_name_key,source_name_alias=source_name_alias,)