# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpfromastropyimportunitsasufromastropy.coordinatesimportSkyCoordfromastropy.ioimportfitsfromastropy.nddataimportNoOverlapErrorfromastropy.timeimportTimefromregionsimportPointSkyRegionfromgammapy.mapsimportHpxNDMap,Map,MapAxis,RegionNDMapfromgammapy.maps.hpx.ioimportHPX_FITS_CONVENTIONS,HpxConvfromgammapy.utils.scriptsimportmake_pathfromgammapy.utils.timeimporttime_ref_from_dict,time_ref_to_dictfrom.importLightCurveTemplateTemporalModel,Models,SkyModel,TemplateSpatialModel__all__=["read_hermes_cube"]def_template_model_from_cta_sdc(filename,t_ref=None):"""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.headerift_refisNone:t_ref=Time(55555.5,format="mjd",scale="tt")time_header.update(time_ref_to_dict(t_ref,t_ref.scale))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,scale=t_ref.scale)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()ifmodelsisNone:returnmodels_cutforminmodels: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