Add FootprintPlotter class and hsp_map_logical_or to plots module

- Move FootprintPlotter class from sp_validation to cs_util.plots
- Add hsp_map_logical_or function for combining HealSparse maps
- Add dependencies: healpy, healsparse, skyproj
- Bump version from 0.1.7 to 0.1.8

This allows shapepipe to use cs_util.plots.FootprintPlotter instead of
requiring the separate sp_validation package.

Generated with [Claude Code](https://claude.com/claude-code)

Co-Authored-By: Claude Sonnet 4.5 <noreply@anthropic.com>

Author: martinkilbinger

cs_util/plots.py

@@ -1,4 +1,4 @@
-"""PLOTS. 
+"""PLOTS.
 
 :Name: plots.py
 
@@ -8,10 +8,17 @@
 
 """
 
+from collections import Counter
+
+import healpy as hp
+import healsparse as hsp
 import matplotlib
 import matplotlib.pylab as plt
 import matplotlib.ticker as ticker
 import numpy as np
+import skyproj
+from astropy import units as u
+from astropy.coordinates import SkyCoord
 
 
 def figure(figsize=(30, 30)):
@@ -426,3 +433,316 @@ def log_ticks(x):
                 ticks.append(val)
 
     return np.array(ticks)
+
+
+class FootprintPlotter:
+    """Class to create footprint plots.
+
+    Parameters
+    -----------
+    nside_coverage: int, optional
+        basic resolution of map; default is 32
+    nside_map:
+        fine resolution for plotting; default is 2048
+
+    """
+
+    # Dictionary storing region parameters
+    _regions = {
+        "NGC": {"ra_0": 180, "extend": [120, 270, 20, 70], "vmax": 60},
+        "SGC": {"ra_0": 15, "extend": [-20, 45, 20, 45], "vmax": 60},
+        "fullsky": {"ra_0": 150, "extend": [0, 360, -90, 90], "vmax": 60},
+    }
+
+    def __init__(self, nside_coverage=32, nside_map=2048):
+
+        self._nside_coverage = nside_coverage
+        self._nside_map = nside_map
+
+    def create_hsp_map(self, ra, dec):
+        """Create Hsp Map.
+
+        Create healsparse map.
+
+        Parameters
+        ----------
+        ra : numpy.ndarray
+            right ascension values
+        dec : numpy.ndarray
+            declination values
+
+        Returns
+        -------
+        hsp.HealSparseMap
+            map
+
+        """
+        # Create empty map
+        hsp_map = hsp.HealSparseMap.make_empty(
+            self._nside_coverage, self._nside_map, dtype=np.float32, sentinel=np.nan
+        )
+
+        # Get pixel list corresponding to coordinates
+        hpix = hp.ang2pix(self._nside_map, ra, dec, nest=True, lonlat=True)
+
+        # Get count of objects per pixel
+        pixel_counts = Counter(hpix)
+
+        # List of unique pixels
+        unique_hpix = np.array(list(pixel_counts.keys()))
+
+        # Number of objects
+        values = np.array(list(pixel_counts.values()), dtype=np.float32)
+
+        # Create maps with numbers per pixel
+        hsp_map[unique_hpix] = values
+
+        return hsp_map
+
+    def plot_area(
+        self,
+        hsp_map,
+        ra_0=0,
+        extend=[120, 270, 29, 70],
+        vmax=60,
+        projection=None,
+        outpath=None,
+        title=None,
+        colorbar=True,
+        colorbar_label="Coverage depth",
+    ):
+        """Plot Area.
+
+        Plot catalogue in an area on the sky.
+
+        Parameters
+        ----------
+        hsp_map : hsp_HealSparseMap
+            input map
+        ra_0 : float, optional
+            anchor point in R.A.; default is 0
+        extend : list, optional
+            sky region, extend=[ra_low, ra_high, dec_low, dec_high];
+            default is [120, 270, 29, 70]
+        vmax : float, optional
+            maximum pixel value to plot with color; default is 60
+        projection : skyproj.McBrydeSkyproj
+            if ``None`` (default), a new plot is created
+        outpath : str, optional
+            output path, default is ``None``
+        title : str, optional
+            print title if not ``None`` (default)
+        colorbar : bool, optional
+            add colorbar; default is ``True``
+        colorbar_label : str, optional
+            colorbar label; default is "Coverage depth"
+
+        Returns
+        --------
+        skyproj.McBrydeSkyproj
+            projection instance
+        plt.axes.Axes
+            axes instance
+
+        Raises
+        ------
+        ValueError
+            if no object found in region
+
+        """
+        if not projection:
+
+            # Create new figure and axes
+            fig, ax = plt.subplots(figsize=(10, 10))
+
+            # Create new projection
+            projection = skyproj.McBrydeSkyproj(
+                ax=ax, lon_0=ra_0, extent=extend, autorescale=True, vmax=vmax
+            )
+        else:
+            ax = None
+
+        im = None
+        try:
+            im, lon_raster, lat_raster, values_raster = projection.draw_hspmap(
+                hsp_map, lon_range=extend[0:2], lat_range=extend[2:]
+            )
+        except ValueError:
+            msg = "No object found in region to draw"
+            print(f"{msg}, continuing...")
+            # raise ValueError(msg)
+
+        projection.draw_milky_way(width=25, linewidth=1.5, color="black", linestyle="-")
+
+        # Add colorbar if requested and image was drawn
+        if colorbar and im is not None:
+            plt.colorbar(im, ax=ax if ax else projection.ax, label=colorbar_label)
+
+        if title:
+            plt.title(title, pad=5)
+
+        if outpath:
+            plt.savefig(outpath)
+
+        return projection, ax
+
+    def plot_region(self, hsp_map, region, projection=None, outpath=None, title=None, colorbar=True, colorbar_label="Coverage depth"):
+        """Plot Region.
+
+        Plot catalogue in a predefined region on the sky.
+
+        Parameters
+        ----------
+        hsp_map : hsp_HealSparseMap
+            input map
+        region : dict
+            region dictionary with keys 'ra_0', 'extend', 'vmax'
+        projection : skyproj.McBrydeSkyproj, optional
+            if ``None`` (default), a new plot is created
+        outpath : str, optional
+            output path, default is ``None``
+        title : str, optional
+            print title if not ``None`` (default)
+        colorbar : bool, optional
+            add colorbar; default is ``True``
+        colorbar_label : str, optional
+            colorbar label; default is "Coverage depth"
+
+        Returns
+        --------
+        skyproj.McBrydeSkyproj
+            projection instance
+        plt.axes.Axes
+            axes instance
+
+        """
+        return self.plot_area(
+            hsp_map,
+            region["ra_0"],
+            region["extend"],
+            region["vmax"],
+            projection=projection,
+            outpath=outpath,
+            title=title,
+            colorbar=colorbar,
+            colorbar_label=colorbar_label,
+        )
+
+    def plot_all_regions(self, hsp_map, outbase=None):
+
+        for region in self._regions:
+            if outbase:
+                outpath = f"{outbase}_{region}.png"
+            else:
+                outpath = None
+            self.plot_region(hsp_map, self._regions[region], outpath=outpath)
+
+    @classmethod
+    def hp_pixel_centers(cls, nside, nest=False):
+
+        # Get number of pixels for given nside
+        npix = hp.nside2npix(nside)
+
+        # Get pixel indices
+        pix_indices = np.arange(npix)
+
+        # Get coordinates of pixel centers
+        ra, dec = hp.pix2ang(nside, pix_indices, nest=nest, lonlat=True)
+
+        return ra, dec, npix
+
+    @classmethod
+    def plot_footprint_as_hp(cls, hsp_map, nside, outpath=None, title=None):
+
+        ra, dec, npix = cls.hp_pixel_centers(nside)
+
+        # Create an empty HEALPix map
+        m = np.full(npix, np.nan)
+
+        fig, ax = plt.subplots(figsize=(10, 10))
+
+        # Plot the HEALPix grid
+        hp.mollview(m, title=title, coord="C", notext=True, rot=(180, 0, 0))
+
+        # Define the Galactic Plane: l = [0, 360], b = 0°
+        for l0, ls in zip((-5, 0, 5), (":", "-", ":")):
+            l_values = np.linspace(0, 360, 500)  # 500 points along the plane
+            b_values = np.zeros_like(l_values)  # Galactic latitude is 0 (the plane)
+
+            # Convert (l, b) to (λ, β) - Ecliptic coordinates
+            coords = SkyCoord(
+                l=l_values * u.degree, b=b_values * u.degree, frame="galactic"
+            )
+            ecl_coords = coords.transform_to(
+                "barycentrictrueecliptic"
+            )  # Ecliptic frame
+
+            # Extract Ecliptic longitude (λ) and latitude (β)
+            lambda_ecl = ecl_coords.lon.deg  # Ecliptic longitude
+            beta_ecl = ecl_coords.lat.deg  # Ecliptic latitude
+
+            # Convert to HEALPix projection coordinates (colatitude, longitude)
+            theta = np.radians(90 - beta_ecl)  # HEALPix uses colatitude
+            phi = np.radians(lambda_ecl)  # HEALPix uses longitude
+
+            # Create a healpy Mollweide projection in Ecliptic coordinates
+            hp.projplot(
+                theta, phi, linestyle=ls, color="black", linewidth=1
+            )  # Plot the outline
+
+        # Apply mask
+        mask_values = hsp_map.get_values_pos(ra, dec, valid_mask=True, lonlat=True)
+
+        ok = np.where(mask_values == False)[0]
+        # nok = np.where(mask_values == False)[0]
+
+        hp.projscatter(ra[ok], dec[ok], lonlat=True, color="green", s=1, marker=".")
+        # hp.projscatter(ra[nok], dec[nok], lonlat=True, color="red", s=1, marker=".")
+
+        plt.tight_layout()
+
+        if outpath:
+            plt.savefig(outpath)
+
+        plt.show()
+
+
+def hsp_map_logical_or(maps, verbose=False):
+    """
+    Hsp Map Logical Or.
+
+    Logical AND of HealSparseMaps.
+
+    """
+    if verbose:
+        print("Combine all maps...")
+
+    # Ensure consistency in coverage and data type
+    nside_coverage = maps[0].nside_coverage
+    nside_sparse = maps[0].nside_sparse
+    dtype = maps[0].dtype
+
+    for m in maps:
+        # MKDEBUG TODO: Change nside if possible
+        if m.nside_coverage != nside_coverage:
+            raise ValueError(
+                f"Coverage nside={m.nside_coverage} does not match {nside_coverage}"
+            )
+        if m.dtype != dtype:
+            raise ValueError(f"Data type {m.dtype} does not match {dtype}")
+
+    # Create an empty HealSparse map
+    map_comb = hsp.HealSparseMap.make_empty(nside_coverage, nside_sparse, dtype=dtype)
+    for idx, m in enumerate(maps):
+        map_comb |= m
+
+        if verbose:
+            valid_pixels = map_comb.valid_pixels
+            n_tot = np.sum(valid_pixels)
+            n_true = np.count_nonzero(valid_pixels)
+            n_false = n_tot - n_true
+            print(
+                f"after map {idx}: frac_true={n_true / n_tot:g}, frac_false={n_false / n_tot:g}"
+            )
+
+    return map_comb

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "cs_util"
-version = "0.1.7"
+version = "0.1.8"
 description = "Utility library for CosmoStat"
 authors = [
 	{ name = "Martin Kilbinger", email = "martin.kilbinger@cea.fr" },
@@ -20,6 +20,9 @@ dependencies = [
 	"keyring>=25.2.0",
 	"pyccl>=3.0.2",
 	"camb>=1.5.9",
+	"healpy>=1.16.0",
+	"healsparse>=1.8.0",
+	"skyproj>=1.0.0",
 ]
 
 [project.optional-dependencies]