Skykatana: Astronomical Masking Pipeline

• Skykatana is an open-source astronomical masking pipeline that builds high-resolution boolean masks using HEALPix, HEALSparse, and MOC representations for wide-field surveys.
• It integrates analytic geometry, catalog-driven masking, and on-demand distributed queries to generate science-ready masks at arcsecond resolution with low memory usage.
• The framework supports scalable pipelines for bright-star and artifact masking, random-point generation, and QA diagnostics, demonstrated with Rubin LSST and HSC datasets.

Skykatana is an open-source pipeline for constructing and combining robust spatial masks for modern wide-field astronomical surveys, with particular emphasis on the Vera Rubin Observatory and Rubin LSST use cases. It builds and combines boolean HEALPix/HEALSparse maps into science-ready masks at arcsecond resolution, represents footprints and exclusion regions with multi-order coverage maps (MOCs), and is engineered for low-memory operation on billion-pixel maps. The framework supports analytic geometry, survey property maps, and catalog-driven masking, integrates with HATS/LSDB for on-demand distributed queries, and provides random-point generation, fractional-area diagnostics, and interactive visualization for quality assurance (Lopez et al., 16 Dec 2025).

1. Scientific motivation and scope

Modern wide-field surveys require reliable spatial masks that excise areas compromised by bright-star halos and bleed trails, detector artifacts, poor observing conditions, and survey-specific geometry. Accurate angular selection functions are necessary because analyses such as galaxy clustering, cosmic shear, photometry of faint objects, and transient-rate estimation are sensitive to unmasked contamination, depth variations, and quality-driven spatial systematics (Lopez et al., 16 Dec 2025).

Skykatana is designed as a survey-agnostic masking framework. Its stated capabilities include building and combining boolean HEALPix/HEALSparse maps into science-ready masks at arcsecond resolution; representing footprints and exclusion regions using MOCs; integrating with LSDB/HATS for on-demand, distributed construction of masks from external catalogs such as Gaia in the Rubin Science Platform; implementing memory-efficient sparse, bit-packed storage and streaming I/O; and supporting logical operations, random-point generation, and fractional-area diagnostics (Lopez et al., 16 Dec 2025).

The framework is oriented toward large survey workflows rather than single-purpose mask products. Inputs may be survey property maps, catalogs such as Gaia or WISE, analytic geometry, or user-defined regions. Intermediate artifacts include MOCs, sparse boolean stages, fractional-area maps, and QA plots. Outputs include combined boolean masks at a chosen order, random catalogs within masks, compact FITS stage files, and pipeline JSON metadata. A plausible implication is that Skykatana is intended not merely as a file format converter but as a reproducible masking pipeline with persistent provenance.

2. Core representations and Boolean semantics

Skykatana is built on three complementary spherical representations: dense HEALPix indexing, sparse HEALSparse storage, and hierarchical MOC region descriptions (Lopez et al., 16 Dec 2025).

Representation	Definition	Role in Skykatana
HEALPix	Equal-area tessellation of the sphere at resolution $N_{\mathrm{side}}$ with NESTED indexing	Base pixelization and resolution model
HEALSparse	Sparse hierarchical representation built on HEALPix in which only pixels present in the mask are stored	High-resolution boolean masking with low memory
MOC	Hierarchical sets of HEALPix cells $\{(o_i,p_i)\}$ whose union approximates arbitrary sky regions at multiple orders	Visualization, chunking, and on-demand catalog queries

For HEALPix, the pixel count and pixel area are

$$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$

HEALPix is well suited to dense, full-sky scalar maps, but dense boolean arrays become memory-heavy at high order (Lopez et al., 16 Dec 2025).

Skykatana uses two HEALSparse storage modes: standard boolean arrays for sparse pixels, and bit-packed encoding with up to eight flags per byte. The latter yields an approximately $8\times$ reduction in memory and I/O while preserving exact geometry and supporting fast streaming operations. Storage is organized as one FITS table per stage, with one row per coverage pixel that has at least one valid child; each row stores bit-packed child flags in little-endian order. A stage with approximately $1.5\times 10^9$ pixels can be written or read in less than one minute and occupies approximately $600$ MB on disk (Lopez et al., 16 Dec 2025).

Boolean operations are defined across resolutions through hierarchical reprojection to a target sparse order using NESTED indexing with exact containment semantics. Skykatana evaluates union $\cup$, intersection $\cap$, and subtraction \ in streaming mode over blocks of sparse rows. When down-sampling, a parent is valid under union if any child is valid; under intersection only if all children are valid; and under subtraction when valid children of the minuend are removed wherever subtrahend children are valid. This resolution-aware semantics is central to constructing masks from heterogeneous ingredients without first materializing dense full-sky arrays.

3. Pipeline architecture and computational methods

Skykatana structures masking workflows as pipelines composed of stages, each stage being a boolean HEALSparse map with metadata (Lopez et al., 16 Dec 2025). Stages can be created from analytic geometry such as circles, ellipses, polygons, and zones; from property maps stored as HEALSparse float maps; from catalog footprints in HATS/LSDB-partitioned sources; or as derived stages from prior ones. The software interface explicitly exposes SkyMaskPipe, build_foot_mask, build_circ_mask, build_box_mask, combine, random_points, and plot_moc.

Geometry-based builders pixelize circles, ellipses, polygons, or zones directly into sparse masks using mocpy for initial MOC generation and then flatten to the requested sparse order. Property-mask builders threshold survey property maps, including coadd depth, seeing, or background variance, to create quality masks. Catalog-driven builders pixelize footprints or artifact lists by streaming sources partitioned by HATS/LSDB, where each partition is a HEALPix cell, typically in Parquet, enabling out-of-core scalability via Dask (Lopez et al., 16 Dec 2025).

The implementation emphasizes streaming and chunking. MOC breadth-first chunking keeps per-chunk area or source counts nearly constant, border padding prevents missed overlaps, and pixels are streamed into bit-packed buffers. Memory scales with active pixel count rather than total sky pixels, and logical combination has complexity

$O(N_{\mathrm{act}})$

with a small constant due to block streaming, where

$$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$

For star masking, pixelization cost scales with the sum of disc areas in pixels,

$\{(o_i,p_i)\}$0

The corresponding runtime model is

$\{(o_i,p_i)\}$1

with $\{(o_i,p_i)\}$2 reflecting catalog query plus disc rasterization cost and $\{(o_i,p_i)\}$3 reflecting streaming logical operations per active row (Lopez et al., 16 Dec 2025).

Random-point generation is defined explicitly. Uniform sampling on the sphere draws $\{(o_i,p_i)\}$4 and sets

$\{(o_i,p_i)\}$5

followed by

$\{(o_i,p_i)\}$6

For masks, Skykatana first precomputes valid cells and their areas at the sampling order, draws a pixel with probability proportional to valid area, samples uniformly within the chosen pixel, and accepts only if the high-resolution mask is True at that coordinate. If sampling from a coarser MOC, it performs accept-reject against the high-resolution sparse mask to guarantee exact containment near edges. This produces uniform randoms over the valid region at the requested resolution (Lopez et al., 16 Dec 2025).

4. Bright-star, halo, and artifact masking

A major Skykatana use case is bright-star masking calibrated from Rubin DP1 source densities around Gaia stars (Lopez et al., 16 Dec 2025). Two definitions are distinguished. For halo included, the exclusion radius is the radius where the radial density profile reaches $\{(o_i,p_i)\}$7 of field density, interpreted as the outer overdensity edge. For halo removed with blendedness_flag == False, the radius is where the profile reaches $\{(o_i,p_i)\}$8, interpreted as the onset of completeness recovery.

For halo-included masks, the radius is modeled as three exponential segments with continuity at $\{(o_i,p_i)\}$9 and $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$0:

$$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$1

For halo-removed masks, the radius is modeled as a single exponential per band:

$$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$2

The best-fit parameters measured in ECDFS are:

• Halo included, three branches per band; boundaries at 13.6 and 15.0
  • u: $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$3, $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$4, $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$5
  • g: $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$6, $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$7, $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$8
  • r: $$N_{\mathrm{pix}} = 12\,N_{\mathrm{side}}^2,\quad A_{\mathrm{pix}} = \frac{4\pi}{N_{\mathrm{pix}}}\ \mathrm{sr},\quad A_{\deg^2} = A_{\mathrm{pix}}\left(\frac{180}{\pi}\right)^2.$$9, $8\times$0, $8\times$1
  • i: $8\times$2, $8\times$3, $8\times$4
  • z: $8\times$5, $8\times$6, $8\times$7 with the note fit unstable/degenerate
  • y: $8\times$8, $8\times$9, $1.5\times 10^9$0
• Halo removed, single branch; y-band unreliable
  • u: $1.5\times 10^9$1
  • g: $1.5\times 10^9$2
  • r: $1.5\times 10^9$3
  • i: $1.5\times 10^9$4
  • z: $1.5\times 10^9$5

The no-halo $1.5\times 10^9$6-band fit also depends on the adopted completeness threshold. For $1.5\times 10^9$7, $1.5\times 10^9$8, and $1.5\times 10^9$9 completeness, the parameters are respectively $600$0, $600$1, and $600$2. Dependence on the limiting magnitude for background density is also tabulated for the no-halo $600$3 band with $600$4:

$600$5

No explicit Gaia-to-LSST-band magnitude conversions are required because the fits use Gaia $600$6 directly (Lopez et al., 16 Dec 2025).

The paper also notes clear limitations. ECDFS spans only approximately $600$7 and contains few very bright stars, so bright-end radii should be used cautiously and recalibrated as Rubin accumulates more data, especially in $600$8 band. A common misconception would be that star halos are treated as exact instrumental morphology; in fact, circular discs approximate halos, whereas spikes and ghosts require survey-specific geometry. In the HSC example, bleed trails are represented as rectangular boxes aligned with detector columns, and instrument-specific rotation is handled through per-CCD definitions in the artifact catalog (Lopez et al., 16 Dec 2025).

5. On-demand masking in Rubin and end-to-end survey workflows

Skykatana includes an on-demand star-mask builder for the Rubin Science Platform. A user-specified stage, such as the WFD footprint, is converted to a MOC up to order $600$9, producing cells

$\cup$0

This search MOC is partitioned into contiguous patches by breadth-first traversal, expanding coarse cells only when needed in order to target near-constant area or source counts per chunk. Each chunk is enlarged by a border padding $\cup$1 exceeding the brightest-star radius $\cup$2 so that overlapping masking discs are not missed if their centers lie just outside nominal chunk boundaries. Optionally, Milky Way avoidance excludes $\cup$3 with a typical range of $\cup$4–$\cup$5 deg together with an elliptical bulge, producing

$\cup$6

Gaia is then queried through LSDB/HATS with selected columns and magnitude limits, radii are assigned via $\cup$7, and discs are pixelized at order $\cup$8 and stream-merged into a single bit-packed stage (Lopez et al., 16 Dec 2025).

Two complete applications are described.

The Subaru HSC–WISE composite mask proceeds as follows: build the HSC footprint from a catalog of approximately $\cup$9 million galaxies using build_foot_mask at order_sparse=13; choose order_out=15 for final combination; build a patch-quality mask by filtering patches with PSF depth thresholds in $\cap$0, $\cap$1, and $\cap$2 such as gmag_psf_depth\>26.2, rmag_psf_depth\>25.9, and imag_psf_depth\>25.7; add bright-star halos as circles and CCD bleed trails as rectangles from HSC artifact catalogs using build_circ_mask and build_box_mask; add a WISE artifact footprint by pixelizing the AllWISE Reject table with build_foot_mask; and combine by intersecting the HSC footprint with the patch-quality mask and subtracting star circles, bleed boxes, and optionally WISE artifacts. Recommended resolution choices are order_sparse $\cap$3–$\cap$4 for moderate-density artifacts and order_out=15, corresponding to approximately $\cap$5 pixels (Lopez et al., 16 Dec 2025).

The Rubin star masks on demand use a WFD $\cap$6-band footprint as the search stage, a Gaia DR3 HATS endpoint, a magnitude range $\cap$7, and a radius function for $\cap$8 completeness, for example in $\cap$9 or $O(N_{\mathrm{act}})$0 band. In the example execution, a Dask cluster in the Rubin Science Platform with 3 workers, threads_per_worker=1, and memory_limit≈6 GiB each runs build_star_mask_online(starq, order_sparse=15, order_cov=5). The output stage, named "starmask", contains approximately $O(N_{\mathrm{act}})$1 pixels at order $O(N_{\mathrm{act}})$2 ($O(N_{\mathrm{act}})$3), is built from approximately $O(N_{\mathrm{act}})$4 Gaia stars, masks approximately $O(N_{\mathrm{act}})$5 in the illustrated configuration, and completes in approximately $O(N_{\mathrm{act}})$6 hours on standard Rubin Science Platform instances. QA products include MOC visualization, zoomed views of masking circles across the footprint, and fractional-area maps at coarser order such as $O(N_{\mathrm{act}})$7, which can then be thresholded to build masks that tolerate up to $O(N_{\mathrm{act}})$8 loss (Lopez et al., 16 Dec 2025).

6. Performance, comparative context, limitations, and released products

Skykatana’s performance claims rest on sparse, bit-packed design, streaming logical operations, and distributed chunked querying. Only active rows are stored; up to eight flags are packed per byte; and I/O throughput is sufficient for a single approximately $O(N_{\mathrm{act}})$9-pixel stage to be written or read in less than one minute at approximately $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$0 MB on disk. Parallelization is delegated to Dask, which distributes chunked queries and pixelization across workers, while breadth-first chunking yields contiguous areas and coherent boundaries and border padding prevents under-masking at chunk edges. The reported hardware and software environment includes the Rubin Science Platform on Google Cloud, Python $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$1, NumPy, Pandas, Astropy, HealSparse, MOCpy, and LSDB/HATS Gaia endpoints (Lopez et al., 16 Dec 2025).

The comparative positioning is explicit. Relative to MANGLE, exact polygonal masks support analytic area and precise boundaries but are complex and scale poorly for numerous shapes; Skykatana’s HEALPix/HEALSparse approach is described as simpler, faster, and memory-efficient for billion-pixel boolean masks and large catalogs. Relative to MOC libraries such as MOCpy, MOCs excel at set operations and visualization, whereas Skykatana extends the workflow to high-resolution boolean masking, streaming logical combination, random sampling, and HEALSparse storage while still exposing MOC views. Relative to survey-specific mask builders, Skykatana generalizes to any survey, supports on-demand catalog querying, configurable radius-magnitude relations, Milky Way avoidance, and fractional-area QA (Lopez et al., 16 Dec 2025).

Limitations are also explicit. Circular discs approximate halos, while spikes and ghosts require survey-specific geometry. Weighted or continuous maps, such as completeness fields, are identified as a natural extension. Bright-end radii require recalibration as Rubin accumulates wide-area data, especially in $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$2 band. These caveats delimit the current method: the framework is scalable and modular, but some astrophysical and instrumental masking problems remain geometry-specific or calibration-limited.

The software is publicly available at https://github.com/samotracio/skykatana under the MIT license, with documentation and notebooks included. Rubin LSST bright-star masks are distributed at https://osf.io/r5vw6. Released products include order-$$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$3 masks for WFD, NES, SCP, LMC–SMC, and SV areas, split into Gaia $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$4 ranges $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$5–$$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$6 and $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$7–$$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$8, totaling approximately $$N_{\mathrm{act}} \sim f_{\mathrm{sky}}\,12\,N_{\mathrm{side}}^2.$$9 billion pixels (Lopez et al., 16 Dec 2025).