Source code for gammapy.image.models.core

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
import numpy as np
import astropy.units as u
from astropy.coordinates.angle_utilities import angular_separation
from astropy.coordinates import Angle, Longitude, Latitude, SkyCoord
from ...utils.fitting import Parameter, Model
from ...maps import Map
from scipy.integrate import quad
from scipy.special import erf


__all__ = [
    "SkySpatialModel",
    "SkyPointSource",
    "SkyGaussian",
    "SkyDisk",
    "SkyEllipse",
    "SkyShell",
    "SkyDiffuseConstant",
    "SkyDiffuseMap",
]


log = logging.getLogger(__name__)

EDGE_WIDTH_95 = 2.326174307353347


def smooth_edge(x, width):
    value = (x / width).to_value("")
    return 0.5 * (1 - erf(value * EDGE_WIDTH_95))


[docs]class SkySpatialModel(Model): """Sky spatial model base class."""
[docs] def __call__(self, lon, lat): """Call evaluate method""" kwargs = dict() for par in self.parameters.parameters: kwargs[par.name] = par.quantity return self.evaluate(lon, lat, **kwargs)
@property def position(self): """Spatial model center position""" try: lon = self.lon_0.quantity lat = self.lat_0.quantity return SkyCoord(lon, lat, frame=self.frame) except IndexError: raise ValueError("Model does not have a defined center position")
[docs]class SkyPointSource(SkySpatialModel): r"""Point Source. .. math:: \phi(lon, lat) = \delta{(lon - lon_0, lat - lat_0)} Parameters ---------- lon_0 : `~astropy.coordinates.Longitude` :math:`lon_0` lat_0 : `~astropy.coordinates.Latitude` :math:`lat_0` frame : {"galactic", "icrs"} Coordinate frame of `lon_0` and `lat_0`. """ __slots__ = ["frame", "lon_0", "lat_0"] def __init__(self, lon_0, lat_0, frame="galactic"): self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 ) self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) super().__init__([self.lon_0, self.lat_0]) @property def evaluation_radius(self): """Returns the effective radius of the sky region where the model evaluates to non-zero. For a Gaussian source, we fix it to :math:`0`. Returns ------- radius : `~astropy.units.Quantity` Radius """ return 0 * u.deg
[docs] @staticmethod def evaluate(lon, lat, lon_0, lat_0): """Evaluate the model (static function).""" wrapval = lon_0 + 180 * u.deg lon = Angle(lon).wrap_at(wrapval) _, grad_lon = np.gradient(lon) grad_lat, _ = np.gradient(lat) lon_diff = np.abs((lon - lon_0) / grad_lon) lat_diff = np.abs((lat - lat_0) / grad_lat) lon_val = np.select([lon_diff < 1], [1 - lon_diff], 0) / np.abs(grad_lon) lat_val = np.select([lat_diff < 1], [1 - lat_diff], 0) / np.abs(grad_lat) return lon_val * lat_val
[docs]class SkyGaussian(SkySpatialModel): r"""Two-dimensional symmetric Gaussian model .. math:: \phi(\text{lon}, \text{lat}) = N \times \text{exp}\left\{-\frac{1}{2} \frac{1-\text{cos}\theta}{1-\text{cos}\sigma}\right\}\,, where :math:`\theta` is the angular separation between the center of the Gaussian and the evaluation point. This angle is calculated on the celestial sphere using the function `angular.separation` defined in `astropy.coordinates.angle_utilities`. The Gaussian is normalized to 1 on the sphere: .. math:: N = \frac{1}{4\pi a\left[1-\text{exp}(-1/a)\right]}\,,\,\,\,\, a = 1-\text{cos}\sigma\,. The normalization factor is in units of :math:`\text{sr}^{-1}`. In the limit of small :math:`\theta` and :math:`\sigma`, this definition reduces to the usual form: .. math:: \phi(\text{lon}, \text{lat}) = \frac{1}{2\pi\sigma^2} \exp{\left(-\frac{1}{2} \frac{\theta^2}{\sigma^2}\right)} Parameters ---------- lon_0 : `~astropy.coordinates.Longitude` :math:`\text{lon}_0` lat_0 : `~astropy.coordinates.Latitude` :math:`\text{lat}_0` sigma : `~astropy.coordinates.Angle` :math:`\sigma` frame : {"galactic", "icrs"} Coordinate frame of `lon_0` and `lat_0`. """ __slots__ = ["frame", "lon_0", "lat_0", "sigma"] def __init__(self, lon_0, lat_0, sigma, frame="galactic"): self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 ) self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) self.sigma = Parameter("sigma", Angle(sigma), min=0) super().__init__([self.lon_0, self.lat_0, self.sigma]) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. For a Gaussian source, we fix it to :math:`5\sigma`. Returns ------- radius : `~astropy.coordinates.Angle` Radius in angular units """ return 5 * self.parameters["sigma"].quantity
[docs] @staticmethod def evaluate(lon, lat, lon_0, lat_0, sigma): """Evaluate the model (static function).""" sep = angular_separation(lon, lat, lon_0, lat_0) a = 1.0 - np.cos(sigma) norm = 1 / (4 * np.pi * a * (1.0 - np.exp(-1.0 / a))) exponent = -0.5 * ((1 - np.cos(sep)) / a) return u.Quantity(norm.value * np.exp(exponent).value, "sr-1", copy=False)
[docs]class SkyDisk(SkySpatialModel): r"""Constant radial disk model. .. math:: \phi(lon, lat) = \frac{1}{2 \pi (1 - \cos{r_0}) } \cdot \begin{cases} 1 & \text{for } \theta \leq r_0 \\ 0 & \text{for } \theta > r_0 \end{cases} where :math:`\theta` is the sky separation. To improve fit convergence of the model, the sharp edges is smoothed using `~scipy.special.erf`. Parameters ---------- lon_0 : `~astropy.coordinates.Longitude` :math:`lon_0` lat_0 : `~astropy.coordinates.Latitude` :math:`lat_0` r_0 : `~astropy.coordinates.Angle` :math:`r_0` edge : `~astropy.coordinates.Angle` Width of the edge. The width is defined as the range within the smooth edges of the model drops from 95% to 5% of its amplitude (see illustration below). frame : {"galactic", "icrs"} Coordinate frame of `lon_0` and `lat_0`. Examples -------- Here is an illustration of the definition of the edge parameter: .. plot:: :include-source: import matplotlib.pyplot as plt from astropy import units as u from gammapy.image.models import SkyDisk lons = np.linspace(0, 0.3, 500) * u.deg r_0 = 0.2 * u.deg edge = 0.1 * u.deg disk = SkyDisk(lon_0="0 deg", lat_0="0 deg", r_0=r_0, edge=edge) profile = disk(lons, 0 * u.deg) plt.plot(lons, profile / profile.max(), alpha=0.5) plt.xlabel("Radius (deg)") plt.ylabel("Profile (A.U.)") edge_min, edge_max = (r_0 - edge / 2.).value, (r_0 + edge / 2.).value plt.vlines([edge_min, edge_max], 0, 1, linestyles=["--"]) plt.annotate("", xy=(edge_min, 0.5), xytext=(edge_min + edge.value, 0.5), arrowprops=dict(arrowstyle="<->", lw=2)) plt.text(0.2, 0.52, "Edge width", ha="center", size=12) plt.hlines([0.95], edge_min - 0.02, edge_min + 0.02, linestyles=["-"]) plt.text(edge_min + 0.02, 0.95, "95%", size=12, va="center") plt.hlines([0.05], edge_max - 0.02, edge_max + 0.02, linestyles=["-"]) plt.text(edge_max - 0.02, 0.05, "5%", size=12, va="center", ha="right") plt.show() """ __slots__ = ["frame", "lon_0", "lat_0", "r_0"] def __init__(self, lon_0, lat_0, r_0, edge="0.01 deg", frame="galactic"): self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 ) self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) self.r_0 = Parameter("r_0", Angle(r_0)) self.edge = Parameter("edge", Angle(edge), min=0.01, frozen=True) super().__init__([self.lon_0, self.lat_0, self.r_0, self.edge]) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. For a Disk source, we fix it to :math:`r_0`. Returns ------- radius : `~astropy.coordinates.Angle` Radius in angular units """ return self.parameters["r_0"].quantity
[docs] @staticmethod def evaluate(lon, lat, lon_0, lat_0, r_0, edge): """Evaluate the model (static function).""" sep = angular_separation(lon, lat, lon_0, lat_0) # Surface area of a spherical cap, see https://en.wikipedia.org/wiki/Spherical_cap norm = 1.0 / (2 * np.pi * (1 - np.cos(r_0))) in_disk = smooth_edge(sep - r_0, edge) return u.Quantity(norm.value * in_disk, "sr-1", copy=False)
[docs]class SkyEllipse(SkySpatialModel): r"""Constant elliptical model. .. math:: \phi(\text{lon}, \text{lat}) = \begin{cases} N & \text{for } \,\,\,dist(F_1,P)+dist(F_2,P)\leq 2 a \\ 0 & \text{otherwise }\,, \end{cases} where :math:`F_1` and :math:`F_2` represent the foci of the ellipse, :math:`P` is a generic point of coordinates :math:`(\text{lon}, \text{lat})`, :math:`a` is the major semiaxis of the ellipse and N is the model's normalization, in units of :math:`\text{sr}^{-1}`. The model is defined on the celestial sphere, with a normalization defined by: .. math:: \int_{4\pi}\phi(\text{lon}, \text{lat}) \,d\Omega = 1\,. Parameters ---------- lon_0 : `~astropy.coordinates.Longitude` :math:`\text{lon}_0`: `lon` coordinate for the center of the ellipse. lat_0 : `~astropy.coordinates.Latitude` :math:`\text{lat}_0`: `lat` coordinate for the center of the ellipse. semi_major : `~astropy.coordinates.Angle` :math:`a`: length of the major semiaxis, in angular units. e : `float` Eccentricity of the ellipse (:math:`0< e< 1`). theta : `~astropy.coordinates.Angle` :math:`\theta`: Rotation angle of the major semiaxis. The rotation angle increases clockwise (i.e., East of North) from the positive `lon` axis. edge : `~astropy.coordinates.Angle` Width of the edge. The width is defined as the range within the smooth edges of the model drops from 95% to 5% of its amplitude (see illustration for `SkyDisk`). frame : {"galactic", "icrs"} Coordinate frame of `lon_0` and `lat_0`. Examples -------- .. plot:: :include-source: import numpy as np import matplotlib.pyplot as plt import astropy.units as u from gammapy.image.models.core import SkyEllipse from gammapy.maps import Map, WcsGeom model = SkyEllipse("2 deg", "2 deg", "1 deg", 0.8, "30 deg") m_geom = WcsGeom.create(binsz=0.01, width=(3, 3), skydir=(2, 2), coordsys="GAL", proj="AIT") coords = m_geom.get_coord() lon = coords.lon * u.deg lat = coords.lat * u.deg vals = model(lon, lat) skymap = Map.from_geom(m_geom, data=vals.value) _, ax, _ = skymap.smooth("0.05 deg").plot() transform = ax.get_transform('galactic') ax.scatter(2, 2, transform=transform, s=20, edgecolor='red', facecolor='red') ax.text(2.0, 1.8, r"$(l_0, b_0)$", transform=transform, ha="center") ax.plot([2, 2 + np.cos(np.pi / 6)], [2, 2 + np.sin(np.pi / 6)], color="r", transform=transform) ax.hlines(y=2, color='r', linestyle='--', transform=transform, xmin=0, xmax=5) ax.text(2.5, 2.06, r"$\theta$", transform=transform); plt.show() """ __slots__ = ["frame", "lon_0", "lat_0", "semi_major", "e", "theta", "_offset_by"] def __init__( self, lon_0, lat_0, semi_major, e, theta, edge="0.01 deg", frame="galactic" ): try: from astropy.coordinates.angle_utilities import offset_by self._offset_by = offset_by except ImportError: raise ImportError("The SkyEllipse model requires astropy>=3.1") self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 ) self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) self.semi_major = Parameter("semi_major", Angle(semi_major)) self.e = Parameter("e", e, min=0, max=1) self.theta = Parameter("theta", Angle(theta)) self.edge = Parameter("edge", Angle(edge), frozen=True, min=0.01) super().__init__( [self.lon_0, self.lat_0, self.semi_major, self.e, self.theta, self.edge] ) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. For an elliptical source, we fix it to the length of the semi-major axis. Returns ------- radius : `~astropy.coordinates.Angle` Radius in angular units """ radius = self.parameters["semi_major"].quantity return radius
[docs] @staticmethod def compute_norm(semi_major, e): """Compute the normalization factor.""" semi_minor = semi_major * np.sqrt(1 - e ** 2) def integral_fcn(x, a, b): A = 1 / np.sin(a) ** 2 B = 1 / np.sin(b) ** 2 C = A - B cs2 = np.cos(x) ** 2 return 1 - np.sqrt(1 - 1 / (B + C * cs2)) return ( 2 * quad(lambda x: integral_fcn(x, semi_major, semi_minor), 0, np.pi)[0] ) ** -1
[docs] def evaluate(self, lon, lat, lon_0, lat_0, semi_major, e, theta, edge): """Evaluate the model (static function).""" # find the foci of the ellipse c = semi_major * e lon_1, lat_1 = self._offset_by(lon_0, lat_0, 90 * u.deg - theta, c) lon_2, lat_2 = self._offset_by(lon_0, lat_0, 270 * u.deg - theta, c) sep_1 = angular_separation(lon, lat, lon_1, lat_1) sep_2 = angular_separation(lon, lat, lon_2, lat_2) in_ellipse = smooth_edge(sep_1 + sep_2 - 2 * semi_major, 2 * edge) norm = SkyEllipse.compute_norm(semi_major, e) return u.Quantity(norm * in_ellipse, "sr-1", copy=False)
[docs]class SkyShell(SkySpatialModel): r"""Shell model. .. math:: \phi(lon, lat) = \frac{3}{2 \pi (r_{out}^3 - r_{in}^3)} \cdot \begin{cases} \sqrt{r_{out}^2 - \theta^2} - \sqrt{r_{in}^2 - \theta^2} & \text{for } \theta \lt r_{in} \\ \sqrt{r_{out}^2 - \theta^2} & \text{for } r_{in} \leq \theta \lt r_{out} \\ 0 & \text{for } \theta > r_{out} \end{cases} where :math:`\theta` is the sky separation and :math:`r_{\text{out}} = r_{\text{in}}` + width Note that the normalization is a small angle approximation, although that approximation is still very good even for 10 deg radius shells. Parameters ---------- lon_0 : `~astropy.coordinates.Longitude` :math:`lon_0` lat_0 : `~astropy.coordinates.Latitude` :math:`lat_0` radius : `~astropy.coordinates.Angle` Inner radius, :math:`r_{in}` width : `~astropy.coordinates.Angle` Shell width frame : {"galactic", "icrs"} Coordinate frame of `lon_0` and `lat_0`. """ __slots__ = ["frame", "lon_0", "lat_0", "radius", "width"] def __init__(self, lon_0, lat_0, radius, width, frame="galactic"): self.frame = frame self.lon_0 = Parameter( "lon_0", Longitude(lon_0).wrap_at("180d"), min=-180, max=180 ) self.lat_0 = Parameter("lat_0", Latitude(lat_0), min=-90, max=90) self.radius = Parameter("radius", Angle(radius)) self.width = Parameter("width", Angle(width)) super().__init__([self.lon_0, self.lat_0, self.radius, self.width]) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. Given by :math:`r_\text{out}`. Returns ------- radius : `~astropy.coordinates.Angle` Radius in angular units """ return self.parameters["radius"].quantity + self.parameters["width"].quantity
[docs] @staticmethod def evaluate(lon, lat, lon_0, lat_0, radius, width): """Evaluate the model (static function).""" sep = angular_separation(lon, lat, lon_0, lat_0) radius_out = radius + width norm = 3 / (2 * np.pi * (radius_out ** 3 - radius ** 3)) with np.errstate(invalid="ignore"): # np.where and np.select do not work with quantities, so we use the # workaround with indexing value = np.sqrt(radius_out ** 2 - sep ** 2) mask = [sep < radius] value[mask] = (value - np.sqrt(radius ** 2 - sep ** 2))[mask] value[sep > radius_out] = 0 return norm * value
[docs]class SkyDiffuseConstant(SkySpatialModel): """Spatially constant (isotropic) spatial model. Parameters ---------- value : `~astropy.units.Quantity` Value """ __slots__ = ["value"] frame = None def __init__(self, value=1): self.value = Parameter("value", value) super().__init__([self.value]) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. For a constant diffuse model, we fix it to None. Returns ------- radius : `~astropy.coordinates.Angle` None """ return None
[docs] @staticmethod def evaluate(lon, lat, value): return value
[docs]class SkyDiffuseMap(SkySpatialModel): """Spatial sky map template model (2D). This is for a 2D image. Use `~gammapy.cube.models.SkyDiffuseCube` for 3D cubes with an energy axis. Parameters ---------- map : `~gammapy.maps.Map` Map template norm : float Norm parameter (multiplied with map values) meta : dict, optional Meta information, meta['filename'] will be used for serialization normalize : bool Normalize the input map so that it integrates to unity. interp_kwargs : dict Interpolation keyword arguments passed to `gammapy.maps.Map.interp_by_coord`. Default arguments are {'interp': 'linear', 'fill_value': 0}. """ __slots__ = ["map", "norm", "meta", "_interp_kwargs"] def __init__(self, map, norm=1, meta=None, normalize=True, interp_kwargs=None): if (map.data < 0).any(): log.warning("Diffuse map has negative values. Check and fix this!") self.map = map if normalize: self.normalize() self.norm = Parameter("norm", norm) self.meta = dict() if meta is None else meta interp_kwargs = {} if interp_kwargs is None else interp_kwargs interp_kwargs.setdefault("interp", "linear") interp_kwargs.setdefault("fill_value", 0) self._interp_kwargs = interp_kwargs super().__init__([self.norm]) @property def evaluation_radius(self): r"""Returns the effective radius of the sky region where the model evaluates to non-zero. For a DiffuseMap source, we fix it to half of the maximal dimension of the map. Returns ------- radius : `~astropy.coordinates.Angle` Radius in angular units. """ return np.max(self.map.geom.width) / 2.0
[docs] def normalize(self): """Normalize the diffuse map model so that it integrates to unity.""" data = self.map.data / self.map.data.sum() data /= self.map.geom.solid_angle().to_value("sr") self.map = self.map.copy(data=data, unit="sr-1")
[docs] @classmethod def read(cls, filename, normalize=True, **kwargs): """Read spatial template model from FITS image. The default unit used if none is found in the file is ``sr-1``. Parameters ---------- filename : str FITS image filename. normalize : bool Normalize the input map so that it integrates to unity. kwargs : dict Keyword arguments passed to `Map.read()`. """ m = Map.read(filename, **kwargs) if m.unit == "": m.unit = "sr-1" return cls(m, normalize=normalize)
[docs] def evaluate(self, lon, lat, norm): """Evaluate model.""" coord = {"lon": lon.to_value("deg"), "lat": lat.to_value("deg")} val = self.map.interp_by_coord(coord, **self._interp_kwargs) return u.Quantity(norm.value * val, self.map.unit, copy=False)
@property def position(self): """`~astropy.coordinates.SkyCoord`""" return self.map.geom.center_skydir @property def frame(self): return self.position.frame.name