Source code for gammapy.data.utils

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import astropy.units as u
from astropy.table import Table
from gammapy.utils.cluster import standard_scaler

__all__ = ["get_irfs_features"]


def get_irfs_features(
    observations,
    energy_true,
    position=None,
    fixed_offset=None,
    names=None,
    containment_fraction=0.68,
    apply_standard_scaler=False,
):
    """Get features from IRFs properties at a given position. Used for observations clustering.

    Parameters
    ----------
    observations : `~gammapy.data.Observations`
        Container holding a list of `~gammapy.data.Observation`.
    energy_true : `~astropy.units.Quantity`
        Energy true at which to compute the containment radius.
    position : `~astropy.coordinates.SkyCoord`, optional
        Sky position. Default is None.
    fixed_offset : `~astropy.coordinates.Angle`, optional
        Offset calculated from the pointing position. Default is None.
        If neither the `position` nor the `fixed_offset` is specified,
        it uses the position of the center of the map by default.
    names : {"edisp-bias", "edisp-res", "psf-radius"}
        IRFs properties to be considered.
        Default is None. If None, all the features are computed.
    containment_fraction : float, optional
        Containment_fraction to compute the `psf-radius`.
        Default is 68%.
    apply_standard_scaler : bool, optional
        Compute standardize features by removing the mean and scaling to unit variance.
        Default is False.

    Returns
    -------
    features : `~astropy.table.Table`
        Features table.
    """
    from gammapy.irf import EDispKernelMap, PSFMap

    if names is None:
        names = ["edisp-bias", "edisp-res", "psf-radius"]

    if position and fixed_offset:
        raise ValueError(
            "`position` and `fixed_offset` arguments are mutually exclusive"
        )

    rows = []

    for obs in observations:
        psf_kwargs = dict(fraction=containment_fraction, energy_true=energy_true)
        if isinstance(obs.psf, PSFMap) and isinstance(obs.edisp, EDispKernelMap):
            if position is None:
                position = obs.psf.psf_map.geom.center_skydir
            edisp_kernel = obs.edisp.get_edisp_kernel(position=position)
            psf_kwargs["position"] = position
        else:
            if fixed_offset is None:
                if position is None:
                    offset = 0 * u.deg
                else:
                    offset_max = np.minimum(
                        obs.psf.axes["offset"].center[-1],
                        obs.edisp.axes["offset"].center[-1],
                    )
                    offset = np.minimum(
                        position.separation(obs.get_pointing_icrs(obs.tmid)), offset_max
                    )
            else:
                offset = fixed_offset
            edisp_kernel = obs.edisp.to_edisp_kernel(offset)
            psf_kwargs["offset"] = offset

        data = {}
        for name in names:
            if name == "edisp-bias":
                data[name] = edisp_kernel.get_bias(energy_true)[0]
            if name == "edisp-res":
                data[name] = edisp_kernel.get_resolution(energy_true)[0]
            if name == "psf-radius":
                data[name] = obs.psf.containment_radius(**psf_kwargs).to("deg")
            data["obs_id"] = obs.obs_id

        rows.append(data)

    features = Table(rows)

    if apply_standard_scaler:
        features = standard_scaler(features)

    features.meta = dict(
        energy_true=energy_true,
        fixed_offset=fixed_offset,
        containment_fraction=containment_fraction,
        apply_standard_scaler=apply_standard_scaler,
    )

    if position:
        features.meta["lon"] = position.galactic.l
        features.meta["lat"] = position.galactic.b
        features.meta["frame"] = "galactic"

    return features