gammapy.maps.region.ndmap.RegionNDMap.plot_hist() - Code Metrics - Inspection of "Some more docs fixes" - gammapy/gammapy - Measure and Improve Code Quality continuously with Scrutinizer

Passed

Push — master ( fb4dcc...ae5621 )

by Axel

created 2021-11-20 03:10 UTC

gammapy.maps.region.ndmap.RegionNDMap.plot_hist() A

↳ Parent: gammapy.maps.region.ndmap

Complexity

Conditions

Size

Total Lines	41
Code Lines	17

Duplication

Lines	0
Ratio	0 %

Importance

Changes

Metric	Value
cc	5
eloc	17
nop	3
dl	0
loc	41
rs	9.0833
c	0
b	0
f	0

from itertools import product
import numpy as np
from scipy.ndimage.measurements import label as ndi_label
from astropy import units as u
from astropy.io import fits
from astropy.nddata import block_reduce
from astropy.table import Table
from astropy.visualization import quantity_support
from gammapy.utils.interpolation import ScaledRegularGridInterpolator, StatProfileScale
from gammapy.utils.scripts import make_path
from ..axes import MapAxes
from ..core import Map
from ..geom import pix_tuple_to_idx
from ..region import RegionGeom
from ..utils import INVALID_INDEX

__all__ = ["RegionNDMap"]


class RegionNDMap(Map):
    """N-dimensional region map.
    A `~RegionNDMap` owns a `~RegionGeom` instance as well as a data array
    containing the values associated to that region in the sky along the non-spatial
    axis, usually an energy axis. The spatial dimensions of a `~RegionNDMap`
    are reduced to a single spatial bin with an arbitrary shape,
    and any extra dimensions are described by an arbitrary number of non-spatial axes.

    Parameters
    ----------
    geom : `~gammapy.maps.RegionGeom`
        Region geometry object.
    data : `~numpy.ndarray`
        Data array. If none then an empty array will be allocated.
    dtype : str, optional
        Data type, default is float32
    meta : `dict`
        Dictionary to store meta data.
    unit : str or `~astropy.units.Unit`
        The map unit
    """

    def __init__(self, geom, data=None, dtype="float32", meta=None, unit=""):
        if data is None:
            data = np.zeros(geom.data_shape, dtype=dtype)

        if meta is None:
            meta = {}

        self._geom = geom
        self.data = data
        self.meta = meta
        self.unit = u.Unit(unit)

    def plot(self, ax=None, axis_name=None, **kwargs):
        """Plot the data contained in region map along the non-spatial axis.

        Parameters
        ----------
        ax : `~matplotlib.pyplot.Axis`
            Axis used for plotting
        axis_name : str
            Which axis to plot on the x axis. Extra axes will be plotted as
            additional lines.
        **kwargs : dict
            Keyword arguments passed to `~matplotlib.pyplot.errorbar`

        Returns
        -------
        ax : `~matplotlib.pyplot.Axis`
            Axis used for plotting
        """
        import matplotlib.pyplot as plt

        ax = ax or plt.gca()

        if axis_name is None:
            if self.geom.axes.is_unidimensional:
                axis_name = self.geom.axes.primary_axis.name
            else:
                raise ValueError(
                    "Plotting a region map with multiple extra axes requires "
                    "specifying the 'axis_name' keyword."
                )

        axis = self.geom.axes[axis_name]

        kwargs.setdefault("marker", "o")
        kwargs.setdefault("markersize", 4)
        kwargs.setdefault("ls", "None")
        kwargs.setdefault("xerr", axis.as_plot_xerr)

        yerr_nd, yerr = kwargs.pop("yerr", None), None
        uplims_nd, uplims = kwargs.pop("uplims", None), None
        label_default = kwargs.pop("label", None)

        labels = product(
            *[ax.as_plot_labels for ax in self.geom.axes if ax.name != axis.name]
        )

        for label_axis, (idx, quantity) in zip(
            labels, self.iter_by_axis(axis_name=axis.name)
        ):
            if isinstance(yerr_nd, tuple):
                yerr = yerr_nd[0][idx], yerr_nd[1][idx]
            elif isinstance(yerr_nd, np.ndarray):
                yerr = yerr_nd[idx]

            if uplims_nd is not None:
                uplims = uplims_nd[idx]

            label = " ".join(label_axis) if label_default is None else label_default

            with quantity_support():
                ax.errorbar(
                    x=axis.as_plot_center,
                    y=quantity,
                    yerr=yerr,
                    uplims=uplims,
                    label=label,
                    **kwargs,
                )

        axis.format_plot_xaxis(ax=ax)

        if "energy" in axis_name:
            ax.set_yscale("log", nonpositive="clip")

        if len(self.geom.axes) > 1:
            plt.legend()

        return ax

    def plot_hist(self, ax=None, **kwargs):
        """Plot as histogram.

        kwargs are forwarded to `~matplotlib.pyplot.hist`

        Parameters
        ----------
        ax : `~matplotlib.axis` (optional)
            Axis instance to be used for the plot
        **kwargs : dict
            Keyword arguments passed to `~matplotlib.pyplot.hist`

        Returns
        -------
        ax : `~matplotlib.pyplot.Axis`
            Axis used for plotting
        """
        import matplotlib.pyplot as plt

        ax = plt.gca() if ax is None else ax

        kwargs.setdefault("histtype", "step")
        kwargs.setdefault("lw", 1)

        if not self.geom.axes.is_unidimensional:
            raise ValueError("Plotting is only supported for unidimensional maps")

        axis = self.geom.axes[0]

        with quantity_support():
            weights = self.data[:, 0, 0]
            ax.hist(
                axis.as_plot_center, bins=axis.as_plot_edges, weights=weights, **kwargs
            )

        if not self.unit.is_unity():
            ax.set_ylabel(f"Data [{self.unit}]")

        axis.format_plot_xaxis(ax=ax)
        ax.set_yscale("log")
        return ax

    def plot_interactive(self):
        raise NotImplementedError(
            "Interactive plotting currently not support for RegionNDMap"
        )

    def plot_region(self, ax=None, **kwargs):
        """Plot region

        Parameters
        ----------
        ax : `~astropy.visualization.WCSAxes`
            Axes to plot on. If no axes are given,
            the region is shown using the minimal
            equivalent WCS geometry.
        **kwargs : dict
            Keyword arguments forwarded to `~regions.PixelRegion.as_artist`
        """
        ax = self.geom.plot_region(ax, **kwargs)
        return ax

    def plot_mask(self, ax=None, **kwargs):
        """Plot the mask as a shaded area in a xmin-xmax range

        Parameters
        ----------
        ax : `~matplotlib.axis`
            Axis instance to be used for the plot.
        **kwargs : dict
            Keyword arguments passed to `~matplotlib.pyplot.axvspan`

        Returns
        -------
        ax : `~matplotlib.pyplot.Axis`
            Axis used for plotting
        """
        import matplotlib.pyplot as plt

        if not self.is_mask:
            raise ValueError("This is not a mask and cannot be plotted")

        kwargs.setdefault("color", "k")
        kwargs.setdefault("alpha", 0.05)
        kwargs.setdefault("label", "mask")

        ax = plt.gca() if ax is None else ax

        edges = self.geom.axes["energy"].edges.reshape((-1, 1, 1))

        labels, nlabels = ndi_label(self.data)

        for idx in range(1, nlabels + 1):
            mask = labels == idx
            xmin = edges[:-1][mask].min().value
            xmax = edges[1:][mask].max().value
            ax.axvspan(xmin, xmax, **kwargs)

        return ax

    @classmethod
    def create(
        cls,
        region,
        axes=None,
        dtype="float32",
        meta=None,
        unit="",
        wcs=None,
        binsz_wcs="0.1deg",
        data=None,
    ):
        """Create an empty region map object.

        Parameters
        ----------
        region : str or `~regions.SkyRegion`
            Region specification
        axes : list of `MapAxis`
            Non spatial axes.
        dtype : str
            Data type, default is 'float32'
        unit : str or `~astropy.units.Unit`
            Data unit.
        meta : `dict`
            Dictionary to store meta data.
        wcs : `~astropy.wcs.WCS`
            WCS projection to use for local projections of the region
        data : `~numpy.ndarray`
            Data array

        Returns
        -------
        map : `RegionNDMap`
            Region map
        """
        geom = RegionGeom.create(region=region, axes=axes, wcs=wcs, binsz_wcs=binsz_wcs)
        return cls(geom=geom, dtype=dtype, unit=unit, meta=meta, data=data)

    def downsample(
        self, factor, preserve_counts=True, axis_name="energy", weights=None
    ):
        """Downsample the non-spatial dimension by a given factor.

        Parameters
        ----------
        factor : int
            Downsampling factor.
        preserve_counts : bool
            Preserve the integral over each bin.  This should be true
            if the map is an integral quantity (e.g. counts) and false if
            the map is a differential quantity (e.g. intensity).
        axis_name : str
            Which axis to downsample. Default is "energy".
        weights : `RegionNDMap`
            Contains the weights to apply to the axis to reduce. Default
            is just weighs of one.

        Returns
        -------
        map : `RegionNDMap`
            Downsampled region map.
        """
        if axis_name is None:
            return self.copy()

        geom = self.geom.downsample(factor=factor, axis_name=axis_name)

        block_size = [1] * self.data.ndim
        idx = self.geom.axes.index_data(axis_name)
        block_size[idx] = factor

        if weights is None:
            weights = 1
        else:
            weights = weights.data

        func = np.nansum if preserve_counts else np.nanmean
        if self.is_mask:
            func = np.all
        data = block_reduce(self.data * weights, tuple(block_size), func=func)

        return self._init_copy(geom=geom, data=data)

    def upsample(self, factor, preserve_counts=True, axis_name="energy"):
        """Upsample the non-spatial dimension by a given factor.

        Parameters
        ----------
        factor : int
            Upsampling factor.
        preserve_counts : bool
            Preserve the integral over each bin.  This should be true
            if the RegionNDMap is an integral quantity (e.g. counts) and false if
            the RegionNDMap is a differential quantity (e.g. intensity).
        axis_name : str
            Which axis to upsample. Default is "energy".

        Returns
        -------
        map : `RegionNDMap`
            Upsampled region map.
        """
        geom = self.geom.upsample(factor=factor, axis_name=axis_name)
        data = self.interp_by_coord(geom.get_coord())

        if preserve_counts:
            data /= factor

        return self._init_copy(geom=geom, data=data)

    def iter_by_axis(self, axis_name):
        """Iterate data by axis

        Parameters
        ----------
        axis_name : str
            Axis name

        Returns
        -------
        idx, data : tuple, `~astropy.units.Quantity`
            Data and index
        """
        idx_axis = self.geom.axes.index_data(axis_name)
        shape = list(self.data.shape)
        shape[idx_axis] = 1

        for idx in np.ndindex(*shape):
            idx = list(idx)
            idx[idx_axis] = slice(None)
            yield tuple(idx), self.quantity[tuple(idx)]

    def fill_by_idx(self, idx, weights=None):

        # TODO: too complex, simplify!
        idx = pix_tuple_to_idx(idx)

        msk = np.all(np.stack([t != INVALID_INDEX.int for t in idx]), axis=0)
        idx = [t[msk] for t in idx]

        if weights is not None:
            if isinstance(weights, u.Quantity):
                weights = weights.to_value(self.unit)
            weights = weights[msk]

        idx = np.ravel_multi_index(idx, self.data.T.shape)
        idx, idx_inv = np.unique(idx, return_inverse=True)
        weights = np.bincount(idx_inv, weights=weights).astype(self.data.dtype)
        self.data.T.flat[idx] += weights

    def get_by_idx(self, idxs):
        return self.data[idxs[::-1]]

    def interp_by_coord(self, coords, **kwargs):
        pix = self.geom.coord_to_pix(coords)
        return self.interp_by_pix(pix, **kwargs)

    def interp_by_pix(self, pix, **kwargs):
        grid_pix = [np.arange(n, dtype=float) for n in self.data.shape[::-1]]

        if np.any(np.isfinite(self.data)):
            data = self.data.copy().T
            data[~np.isfinite(data)] = 0.0
        else:
            data = self.data.T

        scale = kwargs.get("values_scale", "lin")

        if scale == "stat-profile":
            axis = 2 + self.geom.axes.index("norm")
            kwargs["values_scale"] = StatProfileScale(axis=axis)

        fn = ScaledRegularGridInterpolator(grid_pix, data, **kwargs)
        return fn(tuple(pix), clip=False)

    def set_by_idx(self, idx, value):
        self.data[idx[::-1]] = value

    @classmethod
    def read(cls, filename, format="gadf", ogip_column=None, hdu=None):
        """Read from file.

        Parameters
        ----------
        filename : `pathlib.Path` or str
            Filename.
        format : {"gadf", "ogip", "ogip-arf"}
            Which format to use.
        ogip_column : {None, "COUNTS", "QUALITY", "BACKSCAL"}
            If format 'ogip' is chosen which table hdu column to read.
        hdu : str
            Name or index of the HDU with the map data.

        Returns
        -------
        region_map : `RegionNDMap`
            Region nd map
        """
        filename = make_path(filename)
        with fits.open(filename, memmap=False) as hdulist:
            return cls.from_hdulist(
                hdulist, format=format, ogip_column=ogip_column, hdu=hdu
            )

    def write(self, filename, overwrite=False, format="gadf", hdu="SKYMAP"):
        """Write map to file

        Parameters
        ----------
        filename : `pathlib.Path` or str
            Filename.
        format : {"gadf", "ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"}
            Which format to use.
        overwrite : bool
            Overwrite existing files?
        """
        filename = make_path(filename)
        self.to_hdulist(format=format, hdu=hdu).writeto(filename, overwrite=overwrite)

    def to_hdulist(self, format="gadf", hdu="SKYMAP", hdu_bands=None, hdu_region=None):
        """Convert to `~astropy.io.fits.HDUList`.

        Parameters
        ----------
        format : {"gadf", "ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"}
            Format specification
        hdu : str
            Name of the HDU with the map data, used for "gadf" format.
        hdu_bands : str
            Name or index of the HDU with the BANDS table, used for "gadf" format.
        hdu_region : str
            Name or index of the HDU with the region table.

        Returns
        -------
        hdulist : `~astropy.fits.HDUList`
            HDU list
        """
        hdulist = fits.HDUList()
        table = self.to_table(format=format)

        if hdu_bands is None:
            hdu_bands = f"{hdu.upper()}_BANDS"
        if hdu_region is None:
            hdu_region = f"{hdu.upper()}_REGION"

        if format in ["ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"]:
            hdulist.append(fits.BinTableHDU(table))
        elif format == "gadf":
            table.meta.update(self.geom.axes.to_header())
            hdulist.append(fits.BinTableHDU(table, name=hdu))
        else:
            raise ValueError(f"Unsupported format '{format}'")

        if format in ["ogip", "ogip-sherpa", "gadf"]:
            hdulist_geom = self.geom.to_hdulist(
                format=format, hdu_bands=hdu_bands, hdu_region=hdu_region
            )
            hdulist.extend(hdulist_geom[1:])

        return hdulist

    @classmethod
    def from_table(cls, table, format="", colname=None):
        """Create region map from table

        Parameters
        ----------
        table : `~astropy.table.Table`
            Table with input data
        format : {"gadf-sed", "lightcurve"}
            Format to use
        colname : str
            Column name to take the data from.

        Returns
        -------
        region_map : `RegionNDMap`
            Region map
        """
        if format == "gadf-sed":
            if colname is None:
                raise ValueError("Column name required")

            axes = MapAxes.from_table(table=table, format=format)

            if colname == "stat_scan":
                names = ["norm", "energy"]
            # TODO: this is not officially supported by GADF...
            elif colname in ["counts", "npred", "npred_excess"]:
                names = ["dataset", "energy"]
            else:
                names = ["energy"]

            axes = axes[names]
            data = table[colname].data
            unit = table[colname].unit or ""
        elif format == "lightcurve":
            axes = MapAxes.from_table(table=table, format=format)

            if colname == "stat_scan":
                names = ["norm", "energy", "time"]
            # TODO: this is not officially supported by GADF...
            elif colname in ["counts", "npred", "npred_excess"]:
                names = ["dataset", "energy", "time"]
            else:
                names = ["energy", "time"]

            axes = axes[names]
            data = table[colname].data
            unit = table[colname].unit or ""
        else:
            raise ValueError(f"Format not supported {format}")

        geom = RegionGeom.create(region=None, axes=axes)
        return cls(geom=geom, data=data, unit=unit, meta=table.meta, dtype=data.dtype)

    @classmethod
    def from_hdulist(cls, hdulist, format="gadf", ogip_column=None, hdu=None, **kwargs):
        """Create from `~astropy.io.fits.HDUList`.

        Parameters
        ----------
        hdulist : `~astropy.io.fits.HDUList`
            HDU list.
        format : {"gadf", "ogip", "ogip-arf"}
            Format specification
        ogip_column : {"COUNTS", "QUALITY", "BACKSCAL"}
            If format 'ogip' is chosen which table hdu column to read.
        hdu : str
            Name or index of the HDU with the map data.

        Returns
        -------
        region_nd_map : `RegionNDMap`
            Region map.
        """
        defaults = {
            "ogip": {"hdu": "SPECTRUM", "column": "COUNTS"},
            "ogip-arf": {"hdu": "SPECRESP", "column": "SPECRESP"},
            "gadf": {"hdu": "SKYMAP", "column": "DATA"},
        }

        if hdu is None:
            hdu = defaults[format]["hdu"]

        if ogip_column is None:
            ogip_column = defaults[format]["column"]

        geom = RegionGeom.from_hdulist(hdulist, format=format, hdu=hdu)

        table = Table.read(hdulist[hdu])
        quantity = table[ogip_column].quantity

        if ogip_column == "QUALITY":
            data, unit = np.logical_not(quantity.value.astype(bool)), ""
        else:
            data, unit = quantity.value, quantity.unit

        return cls(geom=geom, data=data, meta=table.meta, unit=unit, dtype=data.dtype)

    def _pad_spatial(self, *args, **kwargs):
        raise NotImplementedError("Spatial padding is not supported by RegionNDMap")

    def crop(self):
        raise NotImplementedError("Crop is not supported by RegionNDMap")

    def stack(self, other, weights=None, nan_to_num=True):
        """Stack other region map into map.

        Parameters
        ----------
        other : `RegionNDMap`
            Other map to stack
        weights : `RegionNDMap`
            Array to be used as weights. The spatial geometry must be equivalent
            to `other` and additional axes must be broadcastable.
        nan_to_num: bool
            Non-finite values are replaced by zero if True (default).
        """
        data = other.quantity.to_value(self.unit).astype(self.data.dtype)

        # TODO: re-think stacking of regions. Is making the union reasonable?
        # self.geom.union(other.geom)
        if nan_to_num:
            data = data.copy()
            data[~np.isfinite(data)] = 0
        if weights is not None:
            if not other.geom.to_image() == weights.geom.to_image():
                raise ValueError("Incompatible geoms between map and weights")
            data = data * weights.data

        self.data += data

    def to_table(self, format="gadf"):
        """Convert to `~astropy.table.Table`.

        Data format specification: :ref:`gadf:ogip-pha`

        Parameters
        ----------
        format : {"gadf", "ogip", "ogip-arf", "ogip-arf-sherpa"}
            Format specification

        Returns
        -------
        table : `~astropy.table.Table`
            Table
        """
        data = np.nan_to_num(self.quantity[:, 0, 0])

        if format == "ogip":
            if len(self.geom.axes) > 1:
                raise ValueError(
                    f"Writing to format '{format}' only supports a "
                    f"single energy axis. Got {self.geom.axes.names}"
                )

            energy_axis = self.geom.axes[0]
            energy_axis.assert_name("energy")
            table = Table()
            table["CHANNEL"] = np.arange(energy_axis.nbin, dtype=np.int16)
            table["COUNTS"] = np.array(data, dtype=np.int32)

            # see https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html
            table.meta = {
                "EXTNAME": "SPECTRUM",
                "telescop": "unknown",
                "instrume": "unknown",
                "filter": "None",
                "exposure": 0,
                "corrfile": "",
                "corrscal": "",
                "ancrfile": "",
                "hduclass": "OGIP",
                "hduclas1": "SPECTRUM",
                "hduvers": "1.2.1",
                "poisserr": True,
                "chantype": "PHA",
                "detchans": energy_axis.nbin,
                "quality": 0,
                "backscal": 0,
                "grouping": 0,
                "areascal": 1,
            }

        elif format in ["ogip-arf", "ogip-arf-sherpa"]:
            if len(self.geom.axes) > 1:
                raise ValueError(
                    f"Writing to format '{format}' only supports a "
                    f"single energy axis. Got {self.geom.axes.names}"
                )

            energy_axis = self.geom.axes[0]
            table = energy_axis.to_table(format=format)
            table.meta = {
                "EXTNAME": "SPECRESP",
                "telescop": "unknown",
                "instrume": "unknown",
                "filter": "None",
                "hduclass": "OGIP",
                "hduclas1": "RESPONSE",
                "hduclas2": "SPECRESP",
                "hduvers": "1.1.0",
            }

            if format == "ogip-arf-sherpa":
                data = data.to("cm2")

            table["SPECRESP"] = data

        elif format == "gadf":
            table = Table()
            data = self.quantity.flatten()
            table["CHANNEL"] = np.arange(len(data), dtype=np.int16)
            table["DATA"] = data
        else:
            raise ValueError(f"Unsupported format: '{format}'")

        meta = {k: self.meta.get(k, v) for k, v in table.meta.items()}
        table.meta.update(meta)
        return table

    def get_spectrum(self, *args, **kwargs):
        """Return self"""
        return self

    def to_region_nd_map(self, *args, **kwargs):
        return self

    def cutout(self, *args, **kwargs):
        """Return self"""
        return self


1		from itertools import product
2		import numpy as np
3		from scipy.ndimage.measurements import label as ndi_label
4		from astropy import units as u
5		from astropy.io import fits
6		from astropy.nddata import block_reduce
7		from astropy.table import Table
8		from astropy.visualization import quantity_support
9		from gammapy.utils.interpolation import ScaledRegularGridInterpolator, StatProfileScale
10		from gammapy.utils.scripts import make_path
11		from ..axes import MapAxes
12		from ..core import Map
13		from ..geom import pix_tuple_to_idx
14		from ..region import RegionGeom
15		from ..utils import INVALID_INDEX
16
17		__all__ = ["RegionNDMap"]
18
19
20		class RegionNDMap(Map):
21		"""N-dimensional region map.
22		A `~RegionNDMap` owns a `~RegionGeom` instance as well as a data array
23		containing the values associated to that region in the sky along the non-spatial
24		axis, usually an energy axis. The spatial dimensions of a `~RegionNDMap`
25		are reduced to a single spatial bin with an arbitrary shape,
26		and any extra dimensions are described by an arbitrary number of non-spatial axes.
27
28		Parameters
29		----------
30		geom : `~gammapy.maps.RegionGeom`
31		Region geometry object.
32		data : `~numpy.ndarray`
33		Data array. If none then an empty array will be allocated.
34		dtype : str, optional
35		Data type, default is float32
36		meta : `dict`
37		Dictionary to store meta data.
38		unit : str or `~astropy.units.Unit`
39		The map unit
40		"""
41
42		def __init__(self, geom, data=None, dtype="float32", meta=None, unit=""):
43		if data is None:
44		data = np.zeros(geom.data_shape, dtype=dtype)
45
46		if meta is None:
47		meta = {}
48
49		self._geom = geom
50		self.data = data
51		self.meta = meta
52		self.unit = u.Unit(unit)
53
54		def plot(self, ax=None, axis_name=None, **kwargs):
55		"""Plot the data contained in region map along the non-spatial axis.
56
57		Parameters
58		----------
59		ax : `~matplotlib.pyplot.Axis`
60		Axis used for plotting
61		axis_name : str
62		Which axis to plot on the x axis. Extra axes will be plotted as
63		additional lines.
64		**kwargs : dict
65		Keyword arguments passed to `~matplotlib.pyplot.errorbar`
66
67		Returns
68		-------
69		ax : `~matplotlib.pyplot.Axis`
70		Axis used for plotting
71		"""
72		import matplotlib.pyplot as plt
73
74		ax = ax or plt.gca()
75
76		if axis_name is None:
77		if self.geom.axes.is_unidimensional:
78		axis_name = self.geom.axes.primary_axis.name
79		else:
80		raise ValueError(
81		"Plotting a region map with multiple extra axes requires "
82		"specifying the 'axis_name' keyword."
83		)
84
85		axis = self.geom.axes[axis_name]
86
87		kwargs.setdefault("marker", "o")
88		kwargs.setdefault("markersize", 4)
89		kwargs.setdefault("ls", "None")
90		kwargs.setdefault("xerr", axis.as_plot_xerr)
91
92		yerr_nd, yerr = kwargs.pop("yerr", None), None
93		uplims_nd, uplims = kwargs.pop("uplims", None), None
94		label_default = kwargs.pop("label", None)
95
96		labels = product(
97		*[ax.as_plot_labels for ax in self.geom.axes if ax.name != axis.name]
98		)
99
100		for label_axis, (idx, quantity) in zip(
101		labels, self.iter_by_axis(axis_name=axis.name)
102		):
103		if isinstance(yerr_nd, tuple):
104		yerr = yerr_nd[0][idx], yerr_nd[1][idx]
105		elif isinstance(yerr_nd, np.ndarray):
106		yerr = yerr_nd[idx]
107
108		if uplims_nd is not None:
109		uplims = uplims_nd[idx]
110
111		label = " ".join(label_axis) if label_default is None else label_default
112
113		with quantity_support():
114		ax.errorbar(
115		x=axis.as_plot_center,
116		y=quantity,
117		yerr=yerr,
118		uplims=uplims,
119		label=label,
120		**kwargs,
121		)
122
123		axis.format_plot_xaxis(ax=ax)
124
125		if "energy" in axis_name:
126		ax.set_yscale("log", nonpositive="clip")
127
128		if len(self.geom.axes) > 1:
129		plt.legend()
130
131		return ax
132
133		def plot_hist(self, ax=None, **kwargs):
134		"""Plot as histogram.
135
136		kwargs are forwarded to `~matplotlib.pyplot.hist`
137
138		Parameters
139		----------
140		ax : `~matplotlib.axis` (optional)
141		Axis instance to be used for the plot
142		**kwargs : dict
143		Keyword arguments passed to `~matplotlib.pyplot.hist`
144
145		Returns
146		-------
147		ax : `~matplotlib.pyplot.Axis`
148		Axis used for plotting
149		"""
150		import matplotlib.pyplot as plt
151
152		ax = plt.gca() if ax is None else ax
153
154		kwargs.setdefault("histtype", "step")
155		kwargs.setdefault("lw", 1)
156
157		if not self.geom.axes.is_unidimensional:
158		raise ValueError("Plotting is only supported for unidimensional maps")
159
160		axis = self.geom.axes[0]
161
162		with quantity_support():
163		weights = self.data[:, 0, 0]
164		ax.hist(
165		axis.as_plot_center, bins=axis.as_plot_edges, weights=weights, **kwargs
166		)
167
168		if not self.unit.is_unity():
169		ax.set_ylabel(f"Data [{self.unit}]")
170
171		axis.format_plot_xaxis(ax=ax)
172		ax.set_yscale("log")
173		return ax
174
175		def plot_interactive(self):
176		raise NotImplementedError(
177		"Interactive plotting currently not support for RegionNDMap"
178		)
179
180		def plot_region(self, ax=None, **kwargs):
181		"""Plot region
182
183		Parameters
184		----------
185		ax : `~astropy.visualization.WCSAxes`
186		Axes to plot on. If no axes are given,
187		the region is shown using the minimal
188		equivalent WCS geometry.
189		**kwargs : dict
190		Keyword arguments forwarded to `~regions.PixelRegion.as_artist`
191		"""
192		ax = self.geom.plot_region(ax, **kwargs)
193		return ax
194
195		def plot_mask(self, ax=None, **kwargs):
196		"""Plot the mask as a shaded area in a xmin-xmax range
197
198		Parameters
199		----------
200		ax : `~matplotlib.axis`
201		Axis instance to be used for the plot.
202		**kwargs : dict
203		Keyword arguments passed to `~matplotlib.pyplot.axvspan`
204
205		Returns
206		-------
207		ax : `~matplotlib.pyplot.Axis`
208		Axis used for plotting
209		"""
210		import matplotlib.pyplot as plt
211
212		if not self.is_mask:
213		raise ValueError("This is not a mask and cannot be plotted")
214
215		kwargs.setdefault("color", "k")
216		kwargs.setdefault("alpha", 0.05)
217		kwargs.setdefault("label", "mask")
218
219		ax = plt.gca() if ax is None else ax
220
221		edges = self.geom.axes["energy"].edges.reshape((-1, 1, 1))
222
223		labels, nlabels = ndi_label(self.data)
224
225		for idx in range(1, nlabels + 1):
226		mask = labels == idx
227		xmin = edges[:-1][mask].min().value
228		xmax = edges[1:][mask].max().value
229		ax.axvspan(xmin, xmax, **kwargs)
230
231		return ax
232
233		@classmethod
234		def create(
235		cls,
236		region,
237		axes=None,
238		dtype="float32",
239		meta=None,
240		unit="",
241		wcs=None,
242		binsz_wcs="0.1deg",
243		data=None,
244		):
245		"""Create an empty region map object.
246
247		Parameters
248		----------
249		region : str or `~regions.SkyRegion`
250		Region specification
251		axes : list of `MapAxis`
252		Non spatial axes.
253		dtype : str
254		Data type, default is 'float32'
255		unit : str or `~astropy.units.Unit`
256		Data unit.
257		meta : `dict`
258		Dictionary to store meta data.
259		wcs : `~astropy.wcs.WCS`
260		WCS projection to use for local projections of the region
261		data : `~numpy.ndarray`
262		Data array
263
264		Returns
265		-------
266		map : `RegionNDMap`
267		Region map
268		"""
269		geom = RegionGeom.create(region=region, axes=axes, wcs=wcs, binsz_wcs=binsz_wcs)
270		return cls(geom=geom, dtype=dtype, unit=unit, meta=meta, data=data)
271
272		def downsample(
273		self, factor, preserve_counts=True, axis_name="energy", weights=None
274		):
275		"""Downsample the non-spatial dimension by a given factor.
276
277		Parameters
278		----------
279		factor : int
280		Downsampling factor.
281		preserve_counts : bool
282		Preserve the integral over each bin. This should be true
283		if the map is an integral quantity (e.g. counts) and false if
284		the map is a differential quantity (e.g. intensity).
285		axis_name : str
286		Which axis to downsample. Default is "energy".
287		weights : `RegionNDMap`
288		Contains the weights to apply to the axis to reduce. Default
289		is just weighs of one.
290
291		Returns
292		-------
293		map : `RegionNDMap`
294		Downsampled region map.
295		"""
296		if axis_name is None:
297		return self.copy()
298
299		geom = self.geom.downsample(factor=factor, axis_name=axis_name)
300
301		block_size = [1] * self.data.ndim
302		idx = self.geom.axes.index_data(axis_name)
303		block_size[idx] = factor
304
305		if weights is None:
306		weights = 1
307		else:
308		weights = weights.data
309
310		func = np.nansum if preserve_counts else np.nanmean
311		if self.is_mask:
312		func = np.all
313		data = block_reduce(self.data * weights, tuple(block_size), func=func)
314
315		return self._init_copy(geom=geom, data=data)
316
317		def upsample(self, factor, preserve_counts=True, axis_name="energy"):
318		"""Upsample the non-spatial dimension by a given factor.
319
320		Parameters
321		----------
322		factor : int
323		Upsampling factor.
324		preserve_counts : bool
325		Preserve the integral over each bin. This should be true
326		if the RegionNDMap is an integral quantity (e.g. counts) and false if
327		the RegionNDMap is a differential quantity (e.g. intensity).
328		axis_name : str
329		Which axis to upsample. Default is "energy".
330
331		Returns
332		-------
333		map : `RegionNDMap`
334		Upsampled region map.
335		"""
336		geom = self.geom.upsample(factor=factor, axis_name=axis_name)
337		data = self.interp_by_coord(geom.get_coord())
338
339		if preserve_counts:
340		data /= factor
341
342		return self._init_copy(geom=geom, data=data)
343
344		def iter_by_axis(self, axis_name):
345		"""Iterate data by axis
346
347		Parameters
348		----------
349		axis_name : str
350		Axis name
351
352		Returns
353		-------
354		idx, data : tuple, `~astropy.units.Quantity`
355		Data and index
356		"""
357		idx_axis = self.geom.axes.index_data(axis_name)
358		shape = list(self.data.shape)
359		shape[idx_axis] = 1
360
361		for idx in np.ndindex(*shape):
362		idx = list(idx)
363		idx[idx_axis] = slice(None)
364		yield tuple(idx), self.quantity[tuple(idx)]
365
366	View Code Duplication	def fill_by_idx(self, idx, weights=None):
		0 ignored issues – show Duplication introduced 2020-02-13 15:21 UTC by Report Bug Copy Issue Report This code seems to be duplicated in your project. Loading history...
367		# TODO: too complex, simplify!
368		idx = pix_tuple_to_idx(idx)
369
370		msk = np.all(np.stack([t != INVALID_INDEX.int for t in idx]), axis=0)
371		idx = [t[msk] for t in idx]
372
373		if weights is not None:
374		if isinstance(weights, u.Quantity):
375		weights = weights.to_value(self.unit)
376		weights = weights[msk]
377
378		idx = np.ravel_multi_index(idx, self.data.T.shape)
379		idx, idx_inv = np.unique(idx, return_inverse=True)
380		weights = np.bincount(idx_inv, weights=weights).astype(self.data.dtype)
381		self.data.T.flat[idx] += weights
382
383		def get_by_idx(self, idxs):
384		return self.data[idxs[::-1]]
385
386		def interp_by_coord(self, coords, **kwargs):
387		pix = self.geom.coord_to_pix(coords)
388		return self.interp_by_pix(pix, **kwargs)
389
390		def interp_by_pix(self, pix, **kwargs):
391		grid_pix = [np.arange(n, dtype=float) for n in self.data.shape[::-1]]
392
393		if np.any(np.isfinite(self.data)):
394		data = self.data.copy().T
395		data[~np.isfinite(data)] = 0.0
396		else:
397		data = self.data.T
398
399		scale = kwargs.get("values_scale", "lin")
400
401		if scale == "stat-profile":
402		axis = 2 + self.geom.axes.index("norm")
403		kwargs["values_scale"] = StatProfileScale(axis=axis)
404
405		fn = ScaledRegularGridInterpolator(grid_pix, data, **kwargs)
406		return fn(tuple(pix), clip=False)
407
408		def set_by_idx(self, idx, value):
409		self.data[idx[::-1]] = value
410
411		@classmethod
412		def read(cls, filename, format="gadf", ogip_column=None, hdu=None):
413		"""Read from file.
414
415		Parameters
416		----------
417		filename : `pathlib.Path` or str
418		Filename.
419		format : {"gadf", "ogip", "ogip-arf"}
420		Which format to use.
421		ogip_column : {None, "COUNTS", "QUALITY", "BACKSCAL"}
422		If format 'ogip' is chosen which table hdu column to read.
423		hdu : str
424		Name or index of the HDU with the map data.
425
426		Returns
427		-------
428		region_map : `RegionNDMap`
429		Region nd map
430		"""
431		filename = make_path(filename)
432		with fits.open(filename, memmap=False) as hdulist:
433		return cls.from_hdulist(
434		hdulist, format=format, ogip_column=ogip_column, hdu=hdu
435		)
436
437		def write(self, filename, overwrite=False, format="gadf", hdu="SKYMAP"):
438		"""Write map to file
439
440		Parameters
441		----------
442		filename : `pathlib.Path` or str
443		Filename.
444		format : {"gadf", "ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"}
445		Which format to use.
446		overwrite : bool
447		Overwrite existing files?
448		"""
449		filename = make_path(filename)
450		self.to_hdulist(format=format, hdu=hdu).writeto(filename, overwrite=overwrite)
451
452		def to_hdulist(self, format="gadf", hdu="SKYMAP", hdu_bands=None, hdu_region=None):
453		"""Convert to `~astropy.io.fits.HDUList`.
454
455		Parameters
456		----------
457		format : {"gadf", "ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"}
458		Format specification
459		hdu : str
460		Name of the HDU with the map data, used for "gadf" format.
461		hdu_bands : str
462		Name or index of the HDU with the BANDS table, used for "gadf" format.
463		hdu_region : str
464		Name or index of the HDU with the region table.
465
466		Returns
467		-------
468		hdulist : `~astropy.fits.HDUList`
469		HDU list
470		"""
471		hdulist = fits.HDUList()
472		table = self.to_table(format=format)
473
474		if hdu_bands is None:
475		hdu_bands = f"{hdu.upper()}_BANDS"
476		if hdu_region is None:
477		hdu_region = f"{hdu.upper()}_REGION"
478
479		if format in ["ogip", "ogip-sherpa", "ogip-arf", "ogip-arf-sherpa"]:
480		hdulist.append(fits.BinTableHDU(table))
481		elif format == "gadf":
482		table.meta.update(self.geom.axes.to_header())
483		hdulist.append(fits.BinTableHDU(table, name=hdu))
484		else:
485		raise ValueError(f"Unsupported format '{format}'")
486
487		if format in ["ogip", "ogip-sherpa", "gadf"]:
488		hdulist_geom = self.geom.to_hdulist(
489		format=format, hdu_bands=hdu_bands, hdu_region=hdu_region
490		)
491		hdulist.extend(hdulist_geom[1:])
492
493		return hdulist
494
495		@classmethod
496		def from_table(cls, table, format="", colname=None):
497		"""Create region map from table
498
499		Parameters
500		----------
501		table : `~astropy.table.Table`
502		Table with input data
503		format : {"gadf-sed", "lightcurve"}
504		Format to use
505		colname : str
506		Column name to take the data from.
507
508		Returns
509		-------
510		region_map : `RegionNDMap`
511		Region map
512		"""
513		if format == "gadf-sed":
514		if colname is None:
515		raise ValueError("Column name required")
516
517		axes = MapAxes.from_table(table=table, format=format)
518
519		if colname == "stat_scan":
520		names = ["norm", "energy"]
521		# TODO: this is not officially supported by GADF...
522		elif colname in ["counts", "npred", "npred_excess"]:
523		names = ["dataset", "energy"]
524		else:
525		names = ["energy"]
526
527		axes = axes[names]
528		data = table[colname].data
529		unit = table[colname].unit or ""
530		elif format == "lightcurve":
531		axes = MapAxes.from_table(table=table, format=format)
532
533		if colname == "stat_scan":
534		names = ["norm", "energy", "time"]
535		# TODO: this is not officially supported by GADF...
536		elif colname in ["counts", "npred", "npred_excess"]:
537		names = ["dataset", "energy", "time"]
538		else:
539		names = ["energy", "time"]
540
541		axes = axes[names]
542		data = table[colname].data
543		unit = table[colname].unit or ""
544		else:
545		raise ValueError(f"Format not supported {format}")
546
547		geom = RegionGeom.create(region=None, axes=axes)
548		return cls(geom=geom, data=data, unit=unit, meta=table.meta, dtype=data.dtype)
549
550		@classmethod
551		def from_hdulist(cls, hdulist, format="gadf", ogip_column=None, hdu=None, **kwargs):
552		"""Create from `~astropy.io.fits.HDUList`.
553
554		Parameters
555		----------
556		hdulist : `~astropy.io.fits.HDUList`
557		HDU list.
558		format : {"gadf", "ogip", "ogip-arf"}
559		Format specification
560		ogip_column : {"COUNTS", "QUALITY", "BACKSCAL"}
561		If format 'ogip' is chosen which table hdu column to read.
562		hdu : str
563		Name or index of the HDU with the map data.
564
565		Returns
566		-------
567		region_nd_map : `RegionNDMap`
568		Region map.
569		"""
570		defaults = {
571		"ogip": {"hdu": "SPECTRUM", "column": "COUNTS"},
572		"ogip-arf": {"hdu": "SPECRESP", "column": "SPECRESP"},
573		"gadf": {"hdu": "SKYMAP", "column": "DATA"},
574		}
575
576		if hdu is None:
577		hdu = defaults[format]["hdu"]
578
579		if ogip_column is None:
580		ogip_column = defaults[format]["column"]
581
582		geom = RegionGeom.from_hdulist(hdulist, format=format, hdu=hdu)
583
584		table = Table.read(hdulist[hdu])
585		quantity = table[ogip_column].quantity
586
587		if ogip_column == "QUALITY":
588		data, unit = np.logical_not(quantity.value.astype(bool)), ""
589		else:
590		data, unit = quantity.value, quantity.unit
591
592		return cls(geom=geom, data=data, meta=table.meta, unit=unit, dtype=data.dtype)
593
594		def _pad_spatial(self, args, *kwargs):
595		raise NotImplementedError("Spatial padding is not supported by RegionNDMap")
596
597		def crop(self):
598		raise NotImplementedError("Crop is not supported by RegionNDMap")
599
600		def stack(self, other, weights=None, nan_to_num=True):
601		"""Stack other region map into map.
602
603		Parameters
604		----------
605		other : `RegionNDMap`
606		Other map to stack
607		weights : `RegionNDMap`
608		Array to be used as weights. The spatial geometry must be equivalent
609		to `other` and additional axes must be broadcastable.
610		nan_to_num: bool
611		Non-finite values are replaced by zero if True (default).
612		"""
613		data = other.quantity.to_value(self.unit).astype(self.data.dtype)
614
615		# TODO: re-think stacking of regions. Is making the union reasonable?
616		# self.geom.union(other.geom)
617		if nan_to_num:
618		data = data.copy()
619		data[~np.isfinite(data)] = 0
620		if weights is not None:
621		if not other.geom.to_image() == weights.geom.to_image():
622		raise ValueError("Incompatible geoms between map and weights")
623		data = data * weights.data
624
625		self.data += data
626
627		def to_table(self, format="gadf"):
628		"""Convert to `~astropy.table.Table`.
629
630		Data format specification: :ref:`gadf:ogip-pha`
631
632		Parameters
633		----------
634		format : {"gadf", "ogip", "ogip-arf", "ogip-arf-sherpa"}
635		Format specification
636
637		Returns
638		-------
639		table : `~astropy.table.Table`
640		Table
641		"""
642		data = np.nan_to_num(self.quantity[:, 0, 0])
643
644		if format == "ogip":
645		if len(self.geom.axes) > 1:
646		raise ValueError(
647		f"Writing to format '{format}' only supports a "
648		f"single energy axis. Got {self.geom.axes.names}"
649		)
650
651		energy_axis = self.geom.axes[0]
652		energy_axis.assert_name("energy")
653		table = Table()
654		table["CHANNEL"] = np.arange(energy_axis.nbin, dtype=np.int16)
655		table["COUNTS"] = np.array(data, dtype=np.int32)
656
657		# see https://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/spectra/ogip_92_007/node6.html
658		table.meta = {
659		"EXTNAME": "SPECTRUM",
660		"telescop": "unknown",
661		"instrume": "unknown",
662		"filter": "None",
663		"exposure": 0,
664		"corrfile": "",
665		"corrscal": "",
666		"ancrfile": "",
667		"hduclass": "OGIP",
668		"hduclas1": "SPECTRUM",
669		"hduvers": "1.2.1",
670		"poisserr": True,
671		"chantype": "PHA",
672		"detchans": energy_axis.nbin,
673		"quality": 0,
674		"backscal": 0,
675		"grouping": 0,
676		"areascal": 1,
677		}
678
679		elif format in ["ogip-arf", "ogip-arf-sherpa"]:
680		if len(self.geom.axes) > 1:
681		raise ValueError(
682		f"Writing to format '{format}' only supports a "
683		f"single energy axis. Got {self.geom.axes.names}"
684		)
685
686		energy_axis = self.geom.axes[0]
687		table = energy_axis.to_table(format=format)
688		table.meta = {
689		"EXTNAME": "SPECRESP",
690		"telescop": "unknown",
691		"instrume": "unknown",
692		"filter": "None",
693		"hduclass": "OGIP",
694		"hduclas1": "RESPONSE",
695		"hduclas2": "SPECRESP",
696		"hduvers": "1.1.0",
697		}
698
699		if format == "ogip-arf-sherpa":
700		data = data.to("cm2")
701
702		table["SPECRESP"] = data
703
704		elif format == "gadf":
705		table = Table()
706		data = self.quantity.flatten()
707		table["CHANNEL"] = np.arange(len(data), dtype=np.int16)
708		table["DATA"] = data
709		else:
710		raise ValueError(f"Unsupported format: '{format}'")
711
712		meta = {k: self.meta.get(k, v) for k, v in table.meta.items()}
713		table.meta.update(meta)
714		return table
715
716		def get_spectrum(self, args, *kwargs):
717		"""Return self"""
718		return self
719
720		def to_region_nd_map(self, args, *kwargs):
721		return self
722
723		def cutout(self, args, *kwargs):
724		"""Return self"""
725		return self
726

gammapy / gammapy

Push — master ( fb4dcc...ae5621 )

gammapy.maps.region.ndmap.RegionNDMap.plot_hist() A

Complexity

Size

Duplication

Importance

Duplication Side-by-Side

Filter issues like