# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpfromastropyimportunitsasufromastropy.coordinatesimportSkyCoordfromastropy.ioimportfitsfromastropy.nddataimportNoOverlapErrorfromregionsimportPointSkyRegionfromgammapy.mapsimportHpxNDMap,Map,MapAxis,RegionNDMapfromgammapy.maps.hpx.ioimportHPX_FITS_CONVENTIONS,HpxConvfromgammapy.utils.scriptsimportmake_pathfromgammapy.utils.timeimporttime_ref_from_dictfrom.importLightCurveTemplateTemporalModel,Models,SkyModel,TemplateSpatialModel__all__=["read_hermes_cube"]def_template_model_from_cta_sdc(filename):"""To create a `LightCurveTemplateTemporalModel` from the energy-dependent temporal model files of the cta-sdc1. This format is subject to change. """filename=str(make_path(filename))withfits.open(filename)ashdul:frame=hdul[0].header.get("frame","icrs")position=SkyCoord(hdul[0].header["LONG"]*u.deg,hdul[0].header["LAT"]*u.deg,frame=frame)energy_hdu=hdul["ENERGIES"]energy_axis=MapAxis.from_nodes(nodes=energy_hdu.data,unit=energy_hdu.header["TUNIT1"],name="energy",interp="log",)time_hdu=hdul["TIMES"]time_header=time_hdu.headertime_header.setdefault("MJDREFF",0.5)time_header.setdefault("MJDREFI",55555)time_header.setdefault("scale","tt")time_min=time_hdu.data["Initial Time"]time_max=time_hdu.data["Final Time"]edges=np.append(time_min,time_max[-1])*u.Unit(time_header["TUNIT1"])data=hdul["SPECTRA"]time_ref=time_ref_from_dict(time_header)time_axis=MapAxis.from_edges(edges=edges,name="time",interp="log")reg_map=RegionNDMap.create(region=PointSkyRegion(center=position),axes=[energy_axis,time_axis],data=np.array(list(data.data)*u.Unit(data.header["UNITS"])),)returnLightCurveTemplateTemporalModel(reg_map,t_ref=time_ref,filename=filename)
[docs]defread_hermes_cube(filename):"""Read 3d templates produced with hermes."""# add hermes conventions to the list used by gammapyhermes_conv=HpxConv(convname="hermes-template",colstring="TFLOAT",hduname="xtension",frame="COORDTYPE",quantity_type="differential",)HPX_FITS_CONVENTIONS["hermes-template"]=hermes_convmaps=[]energy_nodes=[]withfits.open(filename)ashdulist:# cannot read directly in 3d with Map.read because BANDS HDU is missing# https://gamma-astro-data-formats.readthedocs.io/en/v0.2/skymaps/index.html#bands-hdu# so we have to loop over hdus and create the energy axisforhduinhdulist[1:]:template=HpxNDMap.from_hdu(hdu,format="hermes-template")# fix missing/incompatible infostemplate._unit=u.Unit(hdu.header["TUNIT1"])# .from_hdu expect "BUNIT"iftemplate.geom.frame=="G":template._geom._frame="galactic"maps.append(template)energy_nodes.append(hdu.header["ENERGY"])# SI unit (see header comment)# create energy axis and set unitenergy_nodes*=u.Jouleenergy_nodes=energy_nodes.to("GeV")axis=MapAxis(energy_nodes,interp="log",name="energy_true",node_type="center",unit="GeV")returnMap.from_stack(maps,axis=axis)
defcutout_template_models(models,cutout_kwargs,datasets_names=None):"""apply cutout to template models"""models_cut=Models()forminmodels:ifisinstance(m.spatial_model,TemplateSpatialModel):try:map_=m.spatial_model.map.cutout(**cutout_kwargs)except(NoOverlapError,ValueError):continuetemplate_cut=TemplateSpatialModel(map_,normalize=m.spatial_model.normalize,)model_cut=SkyModel(spatial_model=template_cut,spectral_model=m.spectral_model,datasets_names=datasets_names,)models_cut.append(model_cut)else:models_cut.append(m)returnmodels_cut