diff --git a/pyproject.toml b/pyproject.toml
index 008439c..67f7d49 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -22,7 +22,7 @@ dependencies = [
     "camb",
     "clmm",
     "colorama",
-    "cs_util>=0.1.5",
+    "cs_util>=0.1.9",
     "emcee",
     "getdist @ git+https://github.com/benabed/getdist.git@upper_triangle_whisker",
     "h5py",
diff --git a/src/sp_validation/plots.py b/src/sp_validation/plots.py
index 3f7f460..b16910e 100644
--- a/src/sp_validation/plots.py
+++ b/src/sp_validation/plots.py
@@ -19,6 +19,7 @@
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from cs_util import plots
+from cs_util.plots import FootprintPlotter
 from lenspack.geometry.projections.gnom import radec2xy
 
 from sp_validation.plot_style import *
@@ -399,278 +400,6 @@ def plot_binned(
     plots.savefig(f"{key}_binned.png")
 
 
-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.