# Licensed under a 3-clause BSD style license - see LICENSE.rstimporthtmlimportloggingimportwarningsfromenumimportEnum,autoimportnumpyasnpimportscipy.interpolateimportastropy.unitsasufromastropy.coordinatesimport(ICRS,AltAz,BaseCoordinateFrame,CartesianRepresentation,SkyCoord,UnitSphericalRepresentation,)fromastropy.ioimportfitsfromastropy.tableimportTablefromastropy.unitsimportQuantityfromastropy.utilsimportlazypropertyfromgammapy.utils.deprecationimportGammapyDeprecationWarningfromgammapy.utils.fitsimportearth_location_from_dict,earth_location_to_dictfromgammapy.utils.scriptsimportmake_pathfromgammapy.utils.timeimporttime_ref_from_dict,time_ref_to_dict,time_to_fits_headerlog=logging.getLogger(__name__)__all__=["FixedPointingInfo","PointingInfo","PointingMode"]def_check_coord_frame(coord_or_frame,expected_frame,name):"""Check if a skycoord or frame is given in expected_frame."""is_coord=isinstance(coord_or_frame,SkyCoord)is_frame=isinstance(coord_or_frame,BaseCoordinateFrame)ifnot(is_frameoris_coord):raiseTypeError(f"{name} must be a 'astropy.coordinates.SkyCoord'""or {expected_frame} instance")ifis_coord:frame=coord_or_frame.frameelse:frame=coord_or_frameifnotisinstance(frame,expected_frame):raiseValueError(f"{name} is in wrong frame, expected {expected_frame}, got {frame}")
[docs]classPointingMode(Enum):""" Describes the behavior of the pointing during the observation. See :ref:`gadf:iact-events-obs-mode`. For ground-based instruments, the most common options will be: * POINTING: The telescope observes a fixed position in the ICRS frame * DRIFT: The telescope observes a fixed position in the alt-az frame Gammapy only supports fixed pointing positions over the whole observation (either in equatorial or horizontal coordinates). OGIP also defines RASTER, SLEW and SCAN. These cannot be treated using a fixed pointing position in either frame, so they would require the pointing table, which is at the moment not supported by gammapy. Data releases based on gadf v0.2 do not have consistent OBS_MODE keyword e.g. the H.E.S.S. data releases uses the not-defined value "WOBBLE". For all gadf data, we assume OBS_MODE to be the same as "POINTING", unless it is set to "DRIFT", making the assumption that one observation only contains a single fixed position. """POINTING=auto()DRIFT=auto()@staticmethoddeffrom_gadf_string(val):"""Parse a string from the GADF header into a PointingMode."""# OBS_MODE is not well-defined and not mandatory in GADF 0.2# We always assume that the observations are pointing observations# unless the OBS_MODE is set to DRIFTifval.upper()=="DRIFT":returnPointingMode.DRIFTelse:returnPointingMode.POINTING
[docs]classFixedPointingInfo:"""IACT array pointing info. Data format specification: :ref:`gadf:iact-pnt` Parameters ---------- meta : `~astropy.table.Table.meta` Meta header info from Table on pointing. Passing this is deprecated, provide ``mode`` and ``fixed_icrs`` or ``fixed_altaz`` instead or use `FixedPointingInfo.from_fits_header` instead. mode : `PointingMode` How the telescope was pointing during the observation. fixed_icrs : `~astropy.coordinates.SkyCoord`, optional The coordinates of the observation in ICRS as a `~astropy.coordinates.SkyCoord` object. Default is None. Required if mode is `PointingMode.POINTING`. fixed_altaz : `~astropy.coordinates.SkyCoord`, optional The coordinates of the observation in alt-az as a `~astropy.coordinates.SkyCoord` object. Default is None. Required if mode is `PointingMode.DRIFT`. Examples -------- >>> from gammapy.data import FixedPointingInfo, PointingMode >>> from astropy.coordinates import SkyCoord >>> import astropy.units as u >>> fixed_icrs = SkyCoord(83.633 * u.deg, 22.014 * u.deg, frame="icrs") >>> pointing_info = FixedPointingInfo(fixed_icrs=fixed_icrs) >>> print(pointing_info) FixedPointingInfo: <BLANKLINE> mode: PointingMode.POINTING coordinates: <SkyCoord (ICRS): (ra, dec) in deg (83.633, 22.014)> """def__init__(self,meta=None,*,mode=None,fixed_icrs=None,fixed_altaz=None,# these have nothing really to do with pointing_info# needed for backwards compatibility but should be removed and accessed# from the observation, not the pointing info.location=None,time_start=None,time_stop=None,time_ref=None,legacy_altaz=None,# store altaz given for mode=POINTING separately so it can be removed easily):self._meta=None# TODO: for backwards compatibility, remove in 2.0# and make other keywards requiredifmetaisnotNone:warnings.warn("Initializing a FixedPointingInfo using a `meta` dict is deprecated",GammapyDeprecationWarning,)self._meta=metaself.__dict__.update(self.from_fits_header(meta).__dict__)returnifmodeisnotNone:warnings.warn("Passing mode is deprecated and the argument will be removed in Gammapy 1.3."" pointing mode is deduced from whether fixed_icrs or fixed_altaz is given",GammapyDeprecationWarning,)self._location=locationself._time_start=time_startself._time_stop=time_stopself._time_ref=time_refself._legacy_altaz=legacy_altazorAltAz(np.nan*u.deg,np.nan*u.deg)iffixed_icrsisnotNoneandfixed_altazisnotNone:raiseValueError("fixed_icrs and fixed_altaz are mutually exclusive")iffixed_icrsisnotNone:_check_coord_frame(fixed_icrs,ICRS,"fixed_icrs")ifnp.isnan(fixed_icrs.ra.value)ornp.isnan(fixed_icrs.dec.value):warnings.warn("In future, fixed_icrs must have non-nan values",GammapyDeprecationWarning,)self._mode=PointingMode.POINTINGself._fixed_icrs=fixed_icrsself._fixed_altaz=Noneelse:_check_coord_frame(fixed_altaz,AltAz,"fixed_altaz")self._mode=PointingMode.DRIFTself._fixed_icrs=Noneself._fixed_altaz=fixed_altazdef_repr_html_(self):try:returnself.to_html()exceptAttributeError:returnf"<pre>{html.escape(str(self))}</pre>"
[docs]@classmethoddeffrom_fits_header(cls,header):""" Parse `~gammapy.data.FixedPointingInfo` from the given FITS header. Parameters ---------- header : `astropy.fits.Header` Header to parse, e.g. from a GADF EVENTS HDU. Currently, only the GADF format is supported. Returns ------- pointing : `~gammapy.data.FixedPointingInfo` The FixedPointingInfo instance filled from the given header. """obs_mode=header.get("OBS_MODE","POINTING")mode=PointingMode.from_gadf_string(obs_mode)try:location=earth_location_from_dict(header)exceptKeyError:location=None# we allow missing RA_PNT / DEC_PNT in POINTING for some reason...# FIXME: actually enforce this to be present instead of using nanra=u.Quantity(header.get("RA_PNT",np.nan),u.deg)dec=u.Quantity(header.get("DEC_PNT",np.nan),u.deg)alt=header.get("ALT_PNT")az=header.get("AZ_PNT")legacy_altaz=Nonefixed_icrs=Nonepointing_altaz=None# we can be more strict with DRIFT, as support was only added recentlyifmodeisPointingMode.DRIFT:ifaltisNoneorazisNone:raiseIOError("Keywords ALT_PNT and AZ_PNT are required for OBSMODE=DRIFT")pointing_altaz=AltAz(alt=alt*u.deg,az=az*u.deg)else:fixed_icrs=SkyCoord(ra,dec)ifnp.isnan(ra.value)ornp.isnan(dec.value):warnings.warn("RA_PNT / DEC_PNT will be required in a future version of"" gammapy for pointing-mode POINTING",GammapyDeprecationWarning,)# store given altaz also for POINTING for backwards compatibility,# FIXME: remove in 2.0ifaltisnotNoneandazisnotNone:legacy_altaz=AltAz(alt=alt*u.deg,az=az*u.deg)time_start=header.get("TSTART")time_stop=header.get("TSTOP")time_ref=Noneiftime_startisnotNoneortime_stopisnotNone:time_ref=time_ref_from_dict(header)time_unit=u.Unit(header.get("TIMEUNIT","s"),format="fits")iftime_startisnotNone:time_start=time_ref+u.Quantity(time_start,time_unit)iftime_stopisnotNone:time_stop=time_ref+u.Quantity(time_stop,time_unit)returncls(location=location,fixed_icrs=fixed_icrs,fixed_altaz=pointing_altaz,time_start=time_start,time_stop=time_stop,time_ref=time_ref,legacy_altaz=legacy_altaz,)
[docs]defto_fits_header(self,format="gadf",version="0.3",time_ref=None):""" Convert this FixedPointingInfo object into a fits header for the given format. Parameters ---------- format : {"gadf"} Format, currently only "gadf" is supported. Default is "gadf". version : str, optional Version of the ``format``, this function currently supports gadf versions 0.2 and 0.3. Default is "0.3". time_ref : `astropy.time.Time`, optional Reference time for storing the time related information in fits format. Default is None. Returns ------- header : `astropy.fits.Header` Header with fixed pointing information filled for the requested format. """ifformat!="gadf":raiseValueError(f'Only the "gadf" format supported, got {format}')ifversionnotin{"0.2","0.3"}:raiseValueError(f"Unsupported version {version} for format {format}")ifself.mode==PointingMode.DRIFTandversion=="0.2":raiseValueError("mode=DRIFT is only supported by GADF 0.3")header=fits.Header()ifself.modeisPointingMode.POINTING:header["OBS_MODE"]="POINTING"header["RA_PNT"]=self.fixed_icrs.ra.deg,u.deg.to_string("fits")header["DEC_PNT"]=self.fixed_icrs.dec.deg,u.deg.to_string("fits")elifself.modeisPointingMode.DRIFT:header["OBS_MODE"]="DRIFT"header["AZ_PNT"]=self.fixed_altaz.az.deg,u.deg.to_string("fits")header["ALT_PNT"]=self.fixed_altaz.alt.deg,u.deg.to_string("fits")# FIXME: remove in 2.0ifself._legacy_altazisnotNoneandnotnp.isnan(self._legacy_altaz.alt.value):header["AZ_PNT"]=self._legacy_altaz.az.deg,u.deg.to_string("fits")header["ALT_PNT"]=self._legacy_altaz.alt.deg,u.deg.to_string("fits")ifself._time_startisnotNone:header["TSTART"]=time_to_fits_header(self._time_start,epoch=time_ref)ifself._time_stopisnotNone:header["TSTOP"]=time_to_fits_header(self._time_stop,epoch=time_ref)ifself._time_startisnotNoneorself._time_stopisnotNone:header.update(time_ref_to_dict(time_ref))ifself._locationisnotNone:header.update(earth_location_to_dict(self._location))returnheader
[docs]@classmethoddefread(cls,filename,hdu="EVENTS"):"""Read pointing information table from file to obtain the metadata. Parameters ---------- filename : str Filename. hdu : int or str, optional HDU number or name. Default is "EVENTS". Returns ------- pointing_info : `PointingInfo` Pointing information. """filename=make_path(filename)header=fits.getheader(filename,extname=hdu)returncls.from_fits_header(header)
@propertydefmode(self):"""See `PointingMode`, if not present, assume POINTING."""returnself._mode@propertydeffixed_altaz(self):"""The fixed coordinates of the observation in alt-az as a `~astropy.coordinates.SkyCoord` object. None if not a DRIFT observation. """returnself._fixed_altaz@propertydeffixed_icrs(self):""" The fixed coordinates of the observation in ICRS as a `~astropy.coordinates.SkyCoord` object. None if not a POINTING observation. """returnself._fixed_icrs
[docs]defget_icrs(self,obstime=None,location=None)->SkyCoord:""" Get the pointing position in ICRS frame for a given time. If the observation was performed tracking a fixed position in ICRS, the icrs pointing is returned with the given obstime attached. If the observation was performed in drift mode, the fixed alt-az coordinates are transformed to ICRS using the observation location and the given time. Parameters ---------- obstime : `astropy.time.Time`, optional Time for which to get the pointing position in ICRS frame. Default is None. location : `astropy.coordinates.EarthLocation`, optional Observatory location, only needed for drift observations to transform from horizontal coordinates to ICRS. Default is None. Returns ------- icrs : `astropy.coordinates.SkyCoord` Pointing position in ICRS frame. """ifself.mode==PointingMode.POINTING:returnSkyCoord(self._fixed_icrs.data,location=location,obstime=obstime)ifself.mode==PointingMode.DRIFT:ifobstimeisNone:obstime=self.obstimereturnself.get_altaz(obstime,location=location).icrsraiseValueError(f"Unsupported pointing mode: {self.mode}.")
[docs]defget_altaz(self,obstime=None,location=None)->SkyCoord:""" Get the pointing position in alt-az frame for a given time. If the observation was performed tracking a fixed position in ICRS, the icrs pointing is transformed at the given time using the location of the observation. If the observation was performed in drift mode, the fixed alt-az coordinate is returned with `obstime` attached. Parameters ---------- obstime : `astropy.time.Time`, optional Time for which to get the pointing position in alt-az frame. Default is None. location : `astropy.coordinates.EarthLocation`, optional Observatory location, only needed for pointing observations to transform from ICRS to horizontal coordinates. Default is None. Returns ------- altaz : `astropy.coordinates.SkyCoord` Pointing position in alt-az frame. """location=locationiflocationisnotNoneelseself._locationframe=AltAz(location=location,obstime=obstime)ifself.mode==PointingMode.POINTING:returnself.fixed_icrs.transform_to(frame)ifself.mode==PointingMode.DRIFT:# see https://github.com/astropy/astropy/issues/12965alt=self.fixed_altaz.altaz=self.fixed_altaz.azreturnSkyCoord(alt=u.Quantity(np.full(obstime.shape,alt.deg),u.deg,copy=False),az=u.Quantity(np.full(obstime.shape,az.deg),u.deg,copy=False),frame=frame,)raiseValueError(f"Unsupported pointing mode: {self.mode}.")
[docs]classPointingInfo:"""IACT array pointing info. Data format specification: :ref:`gadf:iact-pnt`. Parameters ---------- table : `~astropy.table.Table` Table (with meta header information) on pointing. Examples -------- >>> from gammapy.data import PointingInfo >>> pointing_info = PointingInfo.read('$GAMMAPY_DATA/tests/pointing_table.fits.gz') >>> print(pointing_info) Pointing info: <BLANKLINE> Location: GeodeticLocation(lon=<Longitude 16.50022222 deg>, lat=<Latitude -23.27177778 deg>, height=<Quantity 1835. m>) MJDREFI, MJDREFF, TIMESYS = (51910, 0.000742870370370241, 'TT') Time ref: 2001-01-01T00:01:04.184 Time ref: 51910.00074287037 MJD (TT) Duration: 1586.0000000000018 sec = 0.44055555555555603 hours Table length: 100 <BLANKLINE> START: Time: 2004-01-21T19:50:02.184 Time: 53025.826414166666 MJD (TT) RADEC: 83.6333 24.5144 deg ALTAZ: 11.4575 41.3409 deg <BLANKLINE> <BLANKLINE> END: Time: 2004-01-21T20:16:28.184 Time: 53025.844770648146 MJD (TT) RADEC: 83.6333 24.5144 deg ALTAZ: 3.44573 42.1319 deg <BLANKLINE> <BLANKLINE> Note: In order to reproduce the example you need the tests datasets folder. You may download it with the command ``gammapy download datasets --tests --out $GAMMAPY_DATA`` """def__init__(self,table):self.table=table
[docs]@classmethoddefread(cls,filename,hdu="POINTING"):"""Read `PointingInfo` table from file. Parameters ---------- filename : str Filename. hdu : int or str, optional HDU number or name. Default is "POINTING". Returns ------- pointing_info : `PointingInfo` Pointing information. """filename=make_path(filename)table=Table.read(filename,hdu=hdu)returncls(table=table)
def__str__(self):ss="Pointing info:\n\n"ss+=f"Location: {self.location.geodetic}\n"m=self.table.metass+="MJDREFI, MJDREFF, TIMESYS = {}\n".format((m["MJDREFI"],m["MJDREFF"],m["TIMESYS"]))ss+=f"Time ref: {self.time_ref.fits}\n"ss+=f"Time ref: {self.time_ref.mjd} MJD (TT)\n"sec=self.duration.to("second").valuehour=self.duration.to("hour").valuess+=f"Duration: {sec} sec = {hour} hours\n"ss+="Table length: {}\n".format(len(self.table))ss+="\nSTART:\n"+self._str_for_index(0)+"\n"ss+="\nEND:\n"+self._str_for_index(-1)+"\n"returnssdef_str_for_index(self,idx):"""Information for one point in the pointing table."""ss="Time: {}\n".format(self.time[idx].fits)ss+="Time: {} MJD (TT)\n".format(self.time[idx].mjd)ss+="RADEC: {} deg\n".format(self.radec[idx].to_string())ss+="ALTAZ: {} deg\n".format(self.altaz[idx].to_string())returnss@lazypropertydeflocation(self):"""Observatory location as an `~astropy.coordinates.EarthLocation` object."""returnearth_location_from_dict(self.table.meta)@lazypropertydeftime_ref(self):"""Time reference as a `~astropy.time.Time` object."""returntime_ref_from_dict(self.table.meta)@lazypropertydefduration(self):"""Pointing table duration as a `~astropy.time.TimeDelta` object. The time difference between the first and last entry. """returnself.time[-1]-self.time[0]@lazypropertydeftime(self):"""Time array as a `~astropy.time.Time` object."""met=Quantity(self.table["TIME"].astype("float64"),"second")time=self.time_ref+metreturntime.tt@lazypropertydefradec(self):"""RA / DEC position from table as a `~astropy.coordinates.SkyCoord`."""lon=self.table["RA_PNT"]lat=self.table["DEC_PNT"]returnSkyCoord(lon,lat,unit="deg",frame="icrs")@lazypropertydefaltaz_frame(self):"""ALT / AZ frame as a `~astropy.coordinates.AltAz` object."""returnAltAz(obstime=self.time,location=self.location)@lazypropertydefaltaz(self):"""ALT / AZ position computed from RA / DEC as a`~astropy.coordinates.SkyCoord`."""returnself.radec.transform_to(self.altaz_frame)@lazypropertydefaltaz_from_table(self):"""ALT / AZ position from table as a `~astropy.coordinates.SkyCoord`."""lon=self.table["AZ_PNT"]lat=self.table["ALT_PNT"]returnSkyCoord(lon,lat,unit="deg",frame=self.altaz_frame)@staticmethoddef_interpolate_cartesian(mjd_support,coord_support,mjd):xyz=coord_support.cartesianx_new=scipy.interpolate.interp1d(mjd_support,xyz.x)(mjd)y_new=scipy.interpolate.interp1d(mjd_support,xyz.y)(mjd)z_new=scipy.interpolate.interp1d(mjd_support,xyz.z)(mjd)returnCartesianRepresentation(x_new,y_new,z_new).represent_as(UnitSphericalRepresentation)
[docs]defaltaz_interpolate(self,time):"""Interpolate pointing for a given time."""altaz_frame=AltAz(obstime=time,location=self.location)returnSkyCoord(self._interpolate_cartesian(self.time.mjd,self.altaz,time.mjd),frame=altaz_frame,)
[docs]defget_icrs(self,obstime):""" Get the pointing position in ICRS frame for a given time. Parameters ---------- obstime : `astropy.time.Time` Time for which to get the pointing position in ICRS frame. Returns ------- icrs : `astropy.coordinates.SkyCoord` Pointing position in ICRS frame. """returnSkyCoord(self._interpolate_cartesian(self.time.mjd,self.radec,obstime.mjd),obstime=obstime,frame="icrs",)
[docs]defget_altaz(self,obstime):""" Get the pointing position in alt-az frame for a given time. If the observation was performed tracking a fixed position in ICRS, the icrs pointing is transformed at the given time using the location of the observation. If the observation was performed in drift mode, the fixed alt-az coordinate is returned with `obstime` attached. Parameters ---------- obstime : `astropy.time.Time` Time for which to get the pointing position in alt-az frame. Returns ------- altaz : `astropy.coordinates.SkyCoord` Pointing position in alt-az frame. """# give precedence to ALT_PNT / AZ_PNT if presentif"ALT_PNT"inself.tableand"AZ_PNT"inself.table:altaz=self.altaz_from_tableframe=AltAz(obstime=obstime,location=self.location)returnSkyCoord(self._interpolate_cartesian(self.time.mjd,altaz,obstime.mjd),frame=frame,)# fallback to transformation from required ICRS if notreturnself.altaz_interpolate(time=obstime)