feat: add skymaps plots, and skymaps from array instantiation

.pre-commit-config.yaml

@@ -13,7 +13,6 @@ repos:
     rev: 'v1.10.0'
     hooks:
      - id: python-no-eval
-     - id: rst-backticks
   # For python files
   - repo: https://github.com/psf/black
     # Code style

CHANGELOG.md

@@ -8,7 +8,8 @@
 * `cdshaelpix.nested.vertices` can now take depth as a `numpy.ndarray` instead of only
   accepting a single depth
 * new module `skymap` added.
-* read nested all-sky skymaps in the implicit scheme from fits files with `Skymap.from_fits`
+* read/write, and plot nested all-sky skymaps in the implicit scheme from fits files with
+  `Skymap.from_fits`, `Skymap.from_array`, `Skymap.quick_plot`, and `Skymap.to_fits`
 
 ### Fixed
 

docs/api.rst

@@ -87,4 +87,14 @@ While in the nested scheme, ``nside`` is a power of two
         vertices
         vertices_skycoord
 
+cdshealpix.skymap
+~~~~~~~~~~~~~~~~~
+
+This module provides a minimal interface to interact with Skymaps, as defined in the
+`data format for gamma ray astronomy specification
+<https://gamma-astro-data-formats.readthedocs.io/en/latest/skymaps/healpix/index.html#skymap-hdu>`_.
+
+.. autoclass:: cdshealpix.skymap.Skymap
+    :members:
+
 .. _cdshealpix: https://github.com/cds-astro/cds-healpix-python

docs/conf.py

@@ -65,6 +65,8 @@
 ]
 default_role = "py:obj"
 numpydoc_class_members_toctree = False
+numpydoc_show_class_members = False
+autosummary_generate = False
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
 bibtex_bibfiles = ["references.bib"]

docs/examples/examples.rst

@@ -42,11 +42,23 @@ Zone search
 -----------
 
 In this example, we get the ``ipix`` and ``depth`` in a zone and plot them by combining
-:external:ref:`~cdshealpix.nested.vertices`_ with :external:ref:`~matplotlib.path.Polygon`_
+`cdshealpix.nested.vertices` with `matplotlib.path.Polygon`
 
 .. plot:: examples/search_methods/zone_search.py
     :include-source:
 
+Skymaps
+=======
+
+The skymap sub-module allows to manipulate easily all-sky skymaps in the nested ordering
+and implicit schema.
+The class can be instantiated either from a fits file, with `Skymap.from_fits`, or
+directly with a numpy `numpy.array` containing the values associated to each HEALPix
+pixel.
+
+.. plot:: examples/skymaps/skymaps.py
+    :include-source:
+
 Notebook examples
 =================
 

docs/examples/skymaps/skymaps.py

@@ -0,0 +1,7 @@
+"""Read and plots a quick preview of a skymap in a FITS file."""
+
+from cdshealpix.skymap import Skymap
+
+skymap = Skymap.from_fits("skymap.fits")
+print(skymap.depth)
+skymap.quick_plot()

python/cdshealpix/skymap/skymap.py

@@ -10,13 +10,36 @@
 from pathlib import Path
 from typing import Union
 
+import matplotlib.pyplot as plt
+import numpy as np
+
 
 class Skymap:
     """A Skymap, containing values to associate to healpix cells."""
 
     def __init__(self, values):
         self.values = values
 
+    @property
+    def depth(self):
+        """The depth of the skymap.
+
+        Avoids the costly log calculation.
+
+        Returns
+        -------
+        `int`
+            The depth of the skymap.
+
+        Examples
+        --------
+        >>> from cdshealpix.skymap import Skymap
+        >>> map = Skymap.from_array([0]*12)
+        >>> map.depth
+        0
+        """
+        return cdshealpix.depth_skymap(self.values)
+
     @classmethod
     def from_fits(cls, path: Union[str, Path]):
         """Read a skymap in the nested schema from a FITS file.
@@ -34,17 +57,94 @@ def from_fits(cls, path: Union[str, Path]):
 
         Returns
         -------
-        `numpy.array`
-            The map in a numpy array. Its dtype is inferred from the fits header.
+        `Skymap`
+            A skymap. Its values are in a numpy array which data type in inferred from
+            the FITS header.
         """
         return cls(cdshealpix.read_skymap(str(path)))
 
+    @classmethod
+    def from_array(cls, values):
+        """Create a skymap from an array.
+
+        Parameters
+        ----------
+        values : `numpy.array`
+            An array-like object. It should be one-dimensional, and its length should be
+            the number of cells in a HEALPix order.
+            It should be in the nested ordering (not tested).
+
+        Returns
+        -------
+        `Skymap`
+            A skymap object.
+
+        Examples
+        --------
+        >>> from cdshealpix.skymap import Skymap
+        >>> import numpy as np
+        >>> skymap =Skymap.from_array(np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11], dtype=np.uint8))
+        """
+        # only makes a copy if it was not C-contiguous in the first place
+        values = np.ascontiguousarray(values)
+        if values.ndim != 1:
+            raise ValueError(
+                "Skymap values should be one-dimensional. Got an array of "
+                f"shape {values.shape}."
+            )
+        n = int(len(values) / 12)
+        # test if it is a power of two (1000 & 0111 = 0000)
+        if n & (n - 1) != 0 or n == 0:
+            raise ValueError(
+                "The length of values should be a valid number of cells in "
+                "a given HEALPix order, i.e something like 12, 48, 192... "
+                f"Got '{len(values)}'."
+            )
+
+        if values.dtype not in (
+            np.float64,
+            np.float32,
+            np.int64,
+            np.int32,
+            np.int16,
+            np.uint8,
+            float,
+            int,
+        ):
+            raise ValueError(
+                f"The accepted types are f64, f32, i64, i32, u8. Got '{values.dtype}'."
+            )
+        return cls(values)
+
     def to_fits(self, path):
         """Write a Skymap in a fits file.
 
         Parameters
         ----------
-        path : str, pathlib.Path
+        path : `str`, `pathlib.Path`
             The file's path.
         """
         cdshealpix.write_skymap(self.values, str(path))
+
+    def quick_plot(self, *, size=256, convert_to_gal=True, path=None):
+        """Preview a skymap in the Mollweide projection.
+
+        Parameters
+        ----------
+        size : `int`, optional
+            The size of the plot in the y-axis in pixels.
+            It fixes the resolution of the image. By default 256
+        convert_to_gal : `bool`, optional
+            Should the image be converted into a galactic frame? by default True
+        path : `str` or `pathlib.Path`, optional
+            If different from none, the image will not only be displayed, but also saved
+            at the given location. By default None
+        """
+        img = cdshealpix.pixels_skymap(self.values, size, convert_to_gal)
+        fig = plt.imshow(img)
+        plt.axis("off")
+        fig.axes.get_xaxis().set_visible(False)
+        fig.axes.get_yaxis().set_visible(False)
+        if path:
+            plt.savefig(path, bbox_inches="tight", pad_inches=0, transparent=True)
+        plt.show()

python/cdshealpix/tests/test_skymaps.py

@@ -2,8 +2,10 @@
 from tempfile import NamedTemporaryFile
 
 import numpy as np
+import pytest
 
 from ..skymap import Skymap
+from .. import cdshealpix
 
 path_to_test_skymap = Path(__file__).parent.resolve() / "resources" / "skymap.fits"
 
@@ -20,3 +22,24 @@ def test_read_write_read_conservation():
         skymap.to_fits(fp.name)
         skymap2 = Skymap.from_fits(fp.name)
         assert all(skymap.values == skymap2.values)
+
+
+def test_plot():
+    skymap = Skymap.from_array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11])
+    img = cdshealpix.pixels_skymap(skymap.values, 256, True)
+    assert img.shape == (256, 512, 4)
+
+
+def test_depth():
+    n = 12
+    skymap = Skymap.from_array(np.zeros(12 * 4**n))
+    assert skymap.depth == n
+
+
+def test_from_array():
+    with pytest.raises(ValueError, match="The length of values should be*"):
+        Skymap.from_array(np.zeros(3))
+    with pytest.raises(ValueError, match="The accepted types are*"):
+        Skymap.from_array(["a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l"])
+    with pytest.raises(ValueError, match="Skymap values should be one-dimensional*"):
+        Skymap.from_array([[1, 2, 3], [1, 2, 3]])

src/lib.rs

@@ -34,6 +34,10 @@ fn cdshealpix(_py: Python, m: &Bound<'_, PyModule>) -> PyResult<()> {
     .unwrap();
   m.add_function(wrap_pyfunction!(skymap_functions::write_skymap, m)?)
     .unwrap();
+  m.add_function(wrap_pyfunction!(skymap_functions::pixels_skymap, m)?)
+    .unwrap();
+  m.add_function(wrap_pyfunction!(skymap_functions::depth_skymap, m)?)
+    .unwrap();
 
   // wrapper of to_ring and from_ring
   #[pyfn(m)]