One catalogue-masking path for PSF stars and galaxies (manual-mask bit 8 via sky maps)

config/calibration/mask_v2.0.yaml

@@ -32,6 +32,15 @@ params:
   sky_regions: False
   verbose: True
 
+# Sky maps for mask columns a catalogue lacks: a cut on such a column reads
+# the map (True = masked) at each object's position. Galaxy catalogues carry
+# MASK_8_Manual; ShapePipe's PSF-star catalogue does not, so stars selected
+# with this config (a cat_config psf block's mask_config/mask_cuts) take the
+# manual mask from the DR6 bit-8 map. Values are paths or arc:/vos: URLs,
+# fetched once into ~/.cache/sp_validation/masks (sp_validation.masks.map_path).
+maps:
+  MASK_8_Manual: arc:home/mhudson/masks/mask_ugriz_nside131072_n8.hsp
+
 # Masks
 ## Cuts on catalogue columns, in the v2 grammar (sp_validation.grammar)
 dat:

scripts/masking.py

@@ -7,7 +7,7 @@
 import yaml
 
 from sp_validation.io import Catalogue
-from sp_validation.masks import apply_condition, catalogue_cuts
+from sp_validation.masks import catalogue_cuts, catalogue_selection
 
 # -------------------------
 # Spatially-structured cuts: these define the survey footprint.
@@ -51,6 +51,8 @@ def apply_masks(data, mask_config, footprint_only=False):
         Dictionary parsed from the YAML mask configuration file.
         Expected structure:
             - mask_config["dat"]     : list of cuts applied to `data`
+            - mask_config["maps"]    : optional sky maps for cut columns
+              `data` lacks (see sp_validation.masks.MapColumns)
             - mask_config["metacal"] : derived-quantity parameters
               (e.g. relative size limits)
 
@@ -69,17 +71,14 @@ def apply_masks(data, mask_config, footprint_only=False):
         False indicates the object is masked (removed).
     """
 
-    # Initialize mask
-    mask = np.ones(len(data), dtype=bool)
-
-    for cut in catalogue_cuts(mask_config):
-        col = cut["col_name"]
-        if footprint_only and col not in SPATIAL_CUTS:
-            continue
-        kind = cut["kind"]
-        value = cut["value"]
-
-        mask &= apply_condition(data[col], kind, value, cut.get("ignore_bits", 0))
+    cuts = None
+    if footprint_only:
+        cuts = [
+            cut["col_name"]
+            for cut in catalogue_cuts(mask_config)
+            if cut["col_name"] in SPATIAL_CUTS
+        ]
+    mask = catalogue_selection(mask_config, data, cuts=cuts)
 
     # --- metacal relative size (skip for footprint-only) ---
     if not footprint_only:

src/sp_validation/cosmo_val/core.py

@@ -406,6 +406,13 @@ def resolve_paths_for_version(ver):
                     section["path"] = (
                         str(path) if path.is_absolute() else str(subdir / path)
                     )
+            # A psf block's mask config is a config file, not data: relative
+            # to the catalogue config, not to subdir.
+            psf = cc[ver].get("psf")
+            if isinstance(psf, dict) and psf.get("mask_config"):
+                psf["mask_config"] = str(
+                    self.catalog_config_path.parent / psf["mask_config"]
+                )
 
         resolve_paths_for_version("nz")
         processed = {"nz"}

src/sp_validation/cosmo_val/psf_systematics.py

@@ -12,6 +12,7 @@
 import matplotlib.ticker as mticker
 import numpy as np
 import treecorr
+import yaml
 from cs_util import plots as cs_plots
 from shear_psf_leakage import leakage
 from shear_psf_leakage import plots as psfleak_plots
@@ -20,6 +21,7 @@
 
 from .. import blinding, sacc_io
 from ..io import open_entry
+from ..masks import catalogue_selection
 from ..rho_tau import (
     get_rho_tau_w_cov,
     get_samples,
@@ -358,16 +360,31 @@ def _get_galaxy_mask(self, ver, tomo_bin_id):
         return gal_mask
 
     def _get_star_mask(self, ver):
-        cat_star = open_entry(self.cc[ver]["psf"])
-        PSF_flag = self.cc[ver]["psf"].get("PSF_flag")
-        star_flag = self.cc[ver]["psf"].get("star_flag")
-        if PSF_flag is not None:
-            if star_flag is not None:
-                star_mask = (cat_star[PSF_flag] == 0) & (cat_star[star_flag] == 0)
-            else:
-                star_mask = cat_star[PSF_flag] == 0
-        else:
-            star_mask = np.ones(len(cat_star), dtype=bool)
+        """Return the stars that enter rho/tau for version ``ver``.
+
+        A star is kept when the psf block's ``PSF_flag`` and ``star_flag``
+        columns, where configured, are zero, and it passes the block's
+        ``mask_config``: a mask config file (relative to the catalogue
+        config) whose ``mask_cuts`` (column names; default all) select the
+        stars as they select galaxies. Star catalogues carry no mask columns,
+        so those cuts read the config's ``maps`` at ``ra_col``/``dec_col``.
+        """
+        entry = self.cc[ver]["psf"]
+        cat_star = open_entry(entry)
+        star_mask = np.ones(len(cat_star), dtype=bool)
+        for key in ("PSF_flag", "star_flag"):
+            if entry.get(key) is not None:
+                star_mask &= cat_star[entry[key]] == 0
+        if entry.get("mask_config"):
+            with open(entry["mask_config"]) as f:
+                mask_config = yaml.safe_load(f)
+            star_mask &= catalogue_selection(
+                mask_config,
+                cat_star,
+                cuts=entry.get("mask_cuts"),
+                ra_col=entry["ra_col"],
+                dec_col=entry["dec_col"],
+            )
         return star_mask
 
     def set_params_rho_tau(self, ver, params, params_psf):

src/sp_validation/masks.py

@@ -7,13 +7,20 @@
 :Author: Martin Kilbinger
 """
 
+import functools
+import os
+import subprocess
+from pathlib import Path
+
 import healsparse as hsp
 import numpy as np
 from astropy.io import fits
 from scipy import stats
 
 _KIND_ALIASES = {"smaller_equal": "less_equal"}
 
+REMOTE_SCHEMES = ("arc:", "vos:")
+
 _KIND_OPS = {
     "equal": lambda a, v: a == v,
     "not_equal": lambda a, v: a != v,
@@ -64,6 +71,67 @@ def apply_condition(array, kind, value, ignore_bits=0):
     return _KIND_OPS[kind](array, value)
 
 
+def map_path(source):
+    """Return a local path to the sky map ``source`` names.
+
+    ``source`` is a path, or an ``arc:``/``vos:`` URL on CANFAR. A URL is
+    fetched once with ``vcp`` into the mask cache, ``$SP_VALIDATION_MASK_CACHE``
+    or ``~/.cache/sp_validation/masks``, mirroring its remote path, and read
+    from there afterwards.
+    """
+    source = str(source)
+    if not source.startswith(REMOTE_SCHEMES):
+        return source
+    cache = Path(
+        os.environ.get(
+            "SP_VALIDATION_MASK_CACHE", Path.home() / ".cache/sp_validation/masks"
+        )
+    )
+    local = cache / source.split(":", 1)[1].lstrip("/")
+    if not local.exists():
+        local.parent.mkdir(parents=True, exist_ok=True)
+        part = local.with_name(f"{local.name}.part{os.getpid()}")
+        subprocess.run(["vcp", source, str(part)], check=True)
+        os.replace(part, local)
+    return str(local)
+
+
+@functools.lru_cache(maxsize=4)
+def read_map(source):
+    """Return the healsparse map ``source`` names (see :func:`map_path`)."""
+    return hsp.HealSparseMap.read(map_path(source))
+
+
+class MapColumns:
+    """A catalogue whose missing columns are read from sky maps.
+
+    ``maps`` takes a column name to a boolean healsparse map (True = masked).
+    A column the catalogue has is read from the catalogue; a column it lacks
+    is the map's value at each object's ``ra_col``/``dec_col``. So one cut,
+    e.g. ``MASK_8_Manual equal False``, masks galaxies from their column and
+    stars, whose catalogue has no mask columns, from the map.
+    """
+
+    def __init__(self, cat, maps, ra_col="RA", dec_col="Dec"):
+        self._cat = cat
+        self._maps = maps
+        self._ra_col = ra_col
+        self._dec_col = dec_col
+
+    def __len__(self):
+        return len(self._cat)
+
+    def __getitem__(self, col):
+        if col not in self._maps:
+            return self._cat[col]
+        try:
+            return self._cat[col]
+        except (KeyError, ValueError):
+            ra, dec = (np.asarray(self._cat[c]) for c in (self._ra_col, self._dec_col))
+            values = read_map(self._maps[col]).get_values_pos(ra, dec, lonlat=True)
+            return np.asarray(values, dtype=bool)
+
+
 def correlation_matrix(masks, confidence_level=0.9):
     n_key = len(masks)
     print(n_key)
@@ -351,11 +419,15 @@ def catalogue_cuts(config):
     return config.get("dat") or []
 
 
-def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False):
+def get_masks_from_config(
+    config, dat, masks_to_apply=None, verbose=False, ra_col="RA", dec_col="Dec"
+):
     """Get Masks From Config.
 
     Return the masks of the config's ``dat`` cut list (``catalogue_cuts``),
-    evaluated on ``dat``.
+    evaluated on ``dat``. A cut on a column ``dat`` lacks reads it from the
+    config's ``maps`` (column name -> healsparse map path or ``arc:``/``vos:``
+    URL), at ``dat``'s ``ra_col``/``dec_col`` (see :class:`MapColumns`).
 
     Parameters
     ----------
@@ -368,6 +440,8 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False):
         masks to apply exclusively; if `None` (default), use all masks
     verbose : bool, optional
         verbose output if ``True``; default is ``False``
+    ra_col, dec_col : str, optional
+        position columns for map lookups; default ``RA``, ``Dec``
 
     Returns
     -------
@@ -377,6 +451,9 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False):
         list of indices for given mask column name (label)
 
     """
+    if config.get("maps"):
+        dat = MapColumns(dat, config["maps"], ra_col=ra_col, dec_col=dec_col)
+
     # List to store all mask objects
     masks = []
 
@@ -402,3 +479,34 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False):
         masks.append(my_mask)
 
     return masks, labels
+
+
+def catalogue_selection(config, dat, cuts=None, ra_col="RA", dec_col="Dec"):
+    """Return the objects of ``dat`` that pass a mask config's cuts.
+
+    The one place a mask config is evaluated on a catalogue, galaxies and PSF
+    stars alike: the AND of :func:`get_masks_from_config`'s masks.
+
+    Parameters
+    ----------
+    config : dict
+        mask config (``dat`` cut list, optional ``maps``)
+    dat : mapping
+        catalogue, indexable by column name
+    cuts : list of str, optional
+        column names of the cuts to apply; default all
+    ra_col, dec_col : str, optional
+        position columns for map lookups
+
+    Returns
+    -------
+    numpy.ndarray
+        boolean, True for objects kept
+    """
+    masks, _ = get_masks_from_config(
+        config, dat, masks_to_apply=cuts, ra_col=ra_col, dec_col=dec_col
+    )
+    keep = np.ones(len(dat), dtype=bool)
+    for my_mask in masks:
+        keep &= np.asarray(my_mask._mask, dtype=bool)
+    return keep