diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8368a56..6428619 100644
--- a/.pre-commit-config.yaml
+++ b/.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
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 69a12f7..367a5db 100644
--- a/CHANGELOG.md
+++ b/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
 
diff --git a/docs/api.rst b/docs/api.rst
index b899e5f..12f01ed 100644
--- a/docs/api.rst
+++ b/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
diff --git a/docs/conf.py b/docs/conf.py
index 1959d34..305afd9 100644
--- a/docs/conf.py
+++ b/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"]
diff --git a/docs/examples/examples.rst b/docs/examples/examples.rst
index c26a06f..b45e8f3 100644
--- a/docs/examples/examples.rst
+++ b/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
 =================
 
diff --git a/docs/examples/skymaps/skymap.fits b/docs/examples/skymaps/skymap.fits
new file mode 100644
index 0000000..db19a0f
Binary files /dev/null and b/docs/examples/skymaps/skymap.fits differ
diff --git a/docs/examples/skymaps/skymaps.py b/docs/examples/skymaps/skymaps.py
new file mode 100644
index 0000000..42e532d
--- /dev/null
+++ b/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()
diff --git a/docs/examples/skymaps/ztf.fits b/docs/examples/skymaps/ztf.fits
deleted file mode 100644
index f60ab06..0000000
Binary files a/docs/examples/skymaps/ztf.fits and /dev/null differ
diff --git a/python/cdshealpix/skymap/skymap.py b/python/cdshealpix/skymap/skymap.py
index 0710515..bf16fe7 100644
--- a/python/cdshealpix/skymap/skymap.py
+++ b/python/cdshealpix/skymap/skymap.py
@@ -10,6 +10,9 @@
 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."""
@@ -17,6 +20,26 @@ class Skymap:
     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()
diff --git a/python/cdshealpix/tests/test_skymaps.py b/python/cdshealpix/tests/test_skymaps.py
index 9d8114f..b24844d 100644
--- a/python/cdshealpix/tests/test_skymaps.py
+++ b/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]])
diff --git a/src/lib.rs b/src/lib.rs
index 5c42042..8d27664 100644
--- a/src/lib.rs
+++ b/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)]
diff --git a/src/skymap_functions.rs b/src/skymap_functions.rs
index 85e5cf2..560c759 100644
--- a/src/skymap_functions.rs
+++ b/src/skymap_functions.rs
@@ -16,14 +16,10 @@ use healpix::{
   depth_from_n_hash_unsafe,
   nested::map::{
     fits::write::write_implicit_skymap_fits,
-    img::{to_skymap_img_default, Val},
+    img::{to_skymap_img_default, PosConversion, Val},
     skymap::{ImplicitSkyMapArrayRef, SkyMap, SkyMapEnum},
-    HHash,
   },
 };
-use mapproj::pseudocyl::mol::Mol;
-
-use crate::cdshealpix;
 
 #[pyfunction]
 #[pyo3(pass_module)]
@@ -121,6 +117,7 @@ pub fn pixels_skymap<'py>(
   module: &Bound<'py, PyModule>,
   values: SupportedArray<'py>,
   image_size: u16,
+  convert_to_gal: bool,
 ) -> PyResult<Bound<'py, PyArray3<u8>>> {
   let n_hash: u64 = match &values {
     SupportedArray::F64(values) => values.as_array().shape()[0] as u64,
@@ -135,42 +132,109 @@ pub fn pixels_skymap<'py>(
   match values {
     SupportedArray::F64(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
       skymap_ref_to_img(
-        &ImplicitSkyMapArrayRef::<'_,u64,f64>::new(depth, v),
+        &ImplicitSkyMapArrayRef::<'_, u64, f64>::new(depth, v),
+        image_size,
+        convert_to_gal,
+        module.py(),
+      )
+    }),
+    SupportedArray::I64(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
+      skymap_ref_to_img(
+        &ImplicitSkyMapArrayRef::<'_, u64, i64>::new(depth, v),
+        image_size,
+        convert_to_gal,
+        module.py(),
+      )
+    }),
+    SupportedArray::F32(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
+      skymap_ref_to_img(
+        &ImplicitSkyMapArrayRef::<'_, u64, f32>::new(depth, v),
+        image_size,
+        convert_to_gal,
+        module.py(),
+      )
+    }),
+    SupportedArray::I32(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
+      skymap_ref_to_img(
+        &ImplicitSkyMapArrayRef::<'_, u64, i32>::new(depth, v),
+        image_size,
+        convert_to_gal,
+        module.py(),
+      )
+    }),
+    SupportedArray::I16(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
+      skymap_ref_to_img(
+        &ImplicitSkyMapArrayRef::<'_, u64, i16>::new(depth, v),
+        image_size,
+        convert_to_gal,
+        module.py(),
+      )
+    }),
+    SupportedArray::U8(values) => values.as_slice().map_err(|e| e.into()).and_then(|v| {
+      skymap_ref_to_img(
+        &ImplicitSkyMapArrayRef::<'_, u64, u8>::new(depth, v),
         image_size,
+        convert_to_gal,
         module.py(),
       )
     }),
-    /*SupportedArray::I64(values) => ImplicitSkyMapArrayRef::new(depth, values.as_array()),
-    SupportedArray::F32(values) => ImplicitSkyMapArrayRef::new(depth, values.as_slice()),
-    SupportedArray::I32(values) => ImplicitSkyMapArrayRef::new(depth, values.as_slice()),
-    SupportedArray::I16(values) => ImplicitSkyMapArrayRef::new(depth, values.as_slice()),
-    SupportedArray::U8(values) => ImplicitSkyMapArrayRef::new(depth, values.as_slice()),*/
-    _ => todo!(),
   }
 }
 
 fn skymap_ref_to_img<'py, 'a, S>(
   skymap: &'a S,
   image_size: u16,
+  convert_to_gal: bool,
   py: Python<'py>,
 ) -> PyResult<Bound<'py, PyArray3<u8>>>
 where
   S: SkyMap<'a> + 'a,
   S::ValueType: Val,
 {
-  let vec = to_skymap_img_default(
-    skymap,
-    (image_size << 1, image_size),
-    None,
-    None,
-    None,
-    None,
-    None,
-  )
-  .map_err(|e| PyValueError::new_err(e.to_string()))?;
-  PyArray1::from_slice_bound(py, vec.as_slice()).reshape(Ix3(
-    image_size as usize,
-    (image_size << 1) as usize,
-    4_usize,
-  ))
+  if convert_to_gal {
+    let vec = to_skymap_img_default(
+      skymap,
+      (image_size << 1, image_size),
+      None,
+      None,
+      Some(PosConversion::EqMap2GalImg),
+      None,
+      None,
+    )
+    .map_err(|e| PyValueError::new_err(e.to_string()))?;
+    PyArray1::from_slice_bound(py, vec.as_slice()).reshape(Ix3(
+      image_size as usize,
+      (image_size << 1) as usize,
+      4_usize,
+    ))
+  } else {
+    let vec = to_skymap_img_default(
+      skymap,
+      (image_size << 1, image_size),
+      None,
+      None,
+      None,
+      None,
+      None,
+    )
+    .map_err(|e| PyValueError::new_err(e.to_string()))?;
+    PyArray1::from_slice_bound(py, vec.as_slice()).reshape(Ix3(
+      image_size as usize,
+      (image_size << 1) as usize,
+      4_usize,
+    ))
+  }
+}
+
+#[pyfunction]
+pub fn depth_skymap(values: SupportedArray) -> u8 {
+  let n_hash: u64 = match &values {
+    SupportedArray::F64(values) => values.as_array().shape()[0] as u64,
+    SupportedArray::I64(values) => values.as_array().shape()[0] as u64,
+    SupportedArray::F32(values) => values.as_array().shape()[0] as u64,
+    SupportedArray::I32(values) => values.as_array().shape()[0] as u64,
+    SupportedArray::I16(values) => values.as_array().shape()[0] as u64,
+    SupportedArray::U8(values) => values.as_array().shape()[0] as u64,
+  };
+  depth_from_n_hash_unsafe(n_hash)
 }