Skip to content

Quality Control Pipeline

This document designs the automated quality control system for OpenAstro — the gatekeeping layer that catches bad data before it contaminates science products. It covers every stage from raw FITS arrival through time-series combination, specifies the actual flag schema, and maps the Python implementation.

The existing schema in Details on componets.md already includes a quality_flags field (VARCHAR(64)[] in PostgreSQL). This document fills in what that field should actually contain.

Architecture Overview

QC runs in four sequential stages, each triggered by the previous:

FITS arrives at server
        │
        ▼
┌─────────────────────┐
│  Stage 1: INPUT     │  ← Header completeness, plate solve, star count,
│  VALIDATION         │    background, FWHM, saturation
└─────────┬───────────┘
          │ passes or warns
          ▼
┌─────────────────────┐
│  Stage 2: PHOTOMETRIC│  ← Comparison star scatter, airmass trend,
│  QUALITY            │    intra-exposure systematics, outlier detection
└─────────┬───────────┘
          │ passes or warns
          ▼
┌─────────────────────┐
│  Stage 3: TIME-SERIES│  ← Run nightly after combining multiple epochs
│  QUALITY            │    Gap detection, aliasing risk, red noise
└─────────┬───────────┘
          │
          ▼
     quality_flags[]
     quality_score INTEGER
     stored in observations table

Important design constraint: Stages 1 and 2 run server-side on upload and require the raw FITS file to be present. Stage 3 runs as a nightly batch job against the assembled light curve in the database — it does not need raw FITS. This has a direct implication for the upload API: clients should upload full FITS files, not just extracted photometry numbers, if they want full QC coverage.

Stage 1: Input Validation Checks

These checks run immediately when a FITS file arrives, before any photometry is computed. They answer: "Is this image even worth processing?"

1.1 FITS Header Completeness

The following headers are required. Missing ones get specific flags.

Header Keyword Requirement Flag if Missing/Bad
DATE-OBS or DATE-BEG ISO 8601 UTC, not local time HDR_NO_TIMESTAMP
EXPTIME or EXPOSURE > 0 seconds HDR_NO_EXPTIME
FILTER Non-empty string HDR_NO_FILTER
TELESCOP or INSTRUME Non-empty string HDR_NO_INSTRUMENT
SITELAT + SITELON or LAT-OBS + LONG-OBS Decimal degrees, plausible range HDR_NO_SITE_COORDS
GAIN > 0 e-/ADU HDR_NO_GAIN (warning only)
RDNOISE or RDNOIS > 0 e- HDR_NO_RDNOISE (warning only)
XBINNING + YBINNING Positive integer HDR_NO_BINNING (warning only)

WCS keywords (CTYPE1, CTYPE2, CRPIX1, CRPIX2, CRVAL1, CRVAL2, CD1_1 or CDELT1) are checked after plate solving. If the client uploaded pre-solved FITS and the WCS keywords are present but inconsistent (e.g., CRVAL1 outside [0, 360]), flag HDR_WCS_CORRUPT.

Timestamp sanity: After parsing DATE-OBS, check: - Not in the future (> now + 5 minutes tolerance for clock drift): HDR_TIMESTAMP_FUTURE - Not before 2010 (earliest realistic automated site): HDR_TIMESTAMP_ANCIENT - Server submission time minus DATE-OBS is less than EXPTIME + 3600s (allow 1 hour for upload lag): HDR_TIMESTAMP_STALE

Site coordinate cross-check: The FITS header reports site coords. Cross-check against the sites table for the submitting site_id. If they differ by more than 0.1 degrees, flag HDR_COORDS_MISMATCH — either the telescope moved (mobile setup without updating registration), or the header is wrong.

# Implementation: astropy.io.fits
from astropy.io import fits
from astropy.time import Time
import numpy as np

REQUIRED_HEADERS = ['EXPTIME', 'FILTER']
TIMESTAMP_KEYS = ['DATE-OBS', 'DATE-BEG', 'MJD-OBS']

def check_header_completeness(header: fits.Header, site_lat: float, site_lon: float) -> list[str]:
    flags = []

    # Timestamp
    has_timestamp = any(k in header for k in TIMESTAMP_KEYS)
    if not has_timestamp:
        flags.append('HDR_NO_TIMESTAMP')
    else:
        try:
            obs_time = Time(header.get('DATE-OBS') or header.get('DATE-BEG'), format='isot', scale='utc')
            now = Time.now()
            if obs_time > now + 300:      # 5 min future tolerance
                flags.append('HDR_TIMESTAMP_FUTURE')
            if obs_time.jd < 2455197.5:  # 2010-01-01 in JD
                flags.append('HDR_TIMESTAMP_ANCIENT')
        except Exception:
            flags.append('HDR_TIMESTAMP_UNPARSEABLE')

    # Required fields
    for key in REQUIRED_HEADERS:
        if key not in header:
            flags.append(f'HDR_NO_{key}')

    # Exposure time validity
    exptime = header.get('EXPTIME', header.get('EXPOSURE', 0))
    if exptime <= 0:
        flags.append('HDR_EXPTIME_ZERO')

    # Site coords
    lat = header.get('SITELAT', header.get('LAT-OBS'))
    lon = header.get('SITELON', header.get('LONG-OBS'))
    if lat is None or lon is None:
        flags.append('HDR_NO_SITE_COORDS')
    elif abs(lat - site_lat) > 0.1 or abs(lon - site_lon) > 0.1:
        flags.append('HDR_COORDS_MISMATCH')

    return flags

1.2 Plate Solve Success/Failure

Every FITS must be plate solved. The Plate Solving on Server.md document specifies ASTAP primary, astrometry.net fallback. The QC check here wraps the outcome:

Rejection threshold: If both ASTAP and astrometry.net fail, flag SOLVE_FAILED. This is an automatic hard reject for astrometry-dependent science (occultations, precise TTV timing). For photometry-only targets where the client uploaded pre-assigned target coordinates, a failed solve degrades to a warning SOLVE_FAILED_COORDS_ASSUMED.

Solve quality check: When a solve succeeds, verify it is pointing at the right field. Compute the angular separation between the solved image centre and the target's known coordinates from the targets table: - Separation > 5 arcmin: SOLVE_WRONG_FIELD — hard reject. The telescope was pointing at the wrong object. - Separation 1–5 arcmin: SOLVE_POINTING_OFFSET — warning. Possibly acceptable depending on target size relative to FOV. - Solve residuals (mean position error of matched stars) > 2 arcsec: SOLVE_POOR_RESIDUALS — warning.

Pixel scale sanity check: The solved pixel scale should match the registered site's optics. If it differs by more than 20%, the observer changed their optical setup without updating their registration. Flag SOLVE_PIXSCALE_MISMATCH.

from astropy.coordinates import SkyCoord, angular_separation
import astropy.units as u

def check_plate_solve(solve_result: dict, target_ra: float, target_dec: float,
                      registered_pixscale: float) -> list[str]:
    flags = []

    if not solve_result.get('success'):
        flags.append('SOLVE_FAILED')
        return flags

    solved_ra = solve_result['ra_deg']
    solved_dec = solve_result['dec_deg']
    solved_pixscale = solve_result.get('pixel_scale_arcsec')
    solve_residuals = solve_result.get('fit_residuals_arcsec', 0)

    # Angular separation from expected target
    sep = angular_separation(
        np.radians(solved_ra), np.radians(solved_dec),
        np.radians(target_ra), np.radians(target_dec)
    )
    sep_arcmin = np.degrees(sep) * 60

    if sep_arcmin > 5.0:
        flags.append('SOLVE_WRONG_FIELD')
    elif sep_arcmin > 1.0:
        flags.append('SOLVE_POINTING_OFFSET')

    if solve_residuals > 2.0:
        flags.append('SOLVE_POOR_RESIDUALS')

    if solved_pixscale and registered_pixscale:
        if abs(solved_pixscale - registered_pixscale) / registered_pixscale > 0.20:
            flags.append('SOLVE_PIXSCALE_MISMATCH')

    return flags

1.3 Star Count in Field

Count the number of detected point sources in the image above a threshold SNR (~5σ). This single number diagnoses multiple problems.

Thresholds:

Star count Interpretation Flag Action
< 10 Clouds, misfocus, wrong field, or camera problem STARCOUNT_CRITICAL Hard reject
10–20 Thin clouds, marginal transparency, or very sparse field STARCOUNT_LOW Warning — photometry is marginal
20–200 Normal range for typical amateur FOV None Pass
> 500 Galaxy contamination, crowded Milky Way field STARCOUNT_CROWDED Warning — aperture photometry may be contaminated by source blending
> 2000 Almost certainly a dense galactic plane field, globular cluster, or galaxy core STARCOUNT_EXTREME_CROWDING Warning — differential photometry will require PSF fitting, not aperture photometry

Reference: ZTF's image quality pipeline (Masci et al. 2019, PASP) rejects images with fewer than ~10 successfully matched catalog stars for photometric calibration. LCOGT's BANZAI pipeline requires a minimum of 5 sources for a valid catalog match; in practice they see failures below 15–20 sources as indicating clouds or bad focus. The lower bound of 10 is conservative and appropriate for amateur equipment.

Implementation note: Use photutils.detection.DAOStarFinder or sep (Python wrapper for SExtractor) for source detection. sep is significantly faster for large images — relevant if processing hundreds of submissions concurrently.

import sep
import numpy as np

def count_stars(data: np.ndarray, gain: float = 1.0, threshold_sigma: float = 5.0) -> tuple[int, list]:
    """
    Count point sources in the image using SEP (SExtractor-based).
    Returns (star_count, source_table).
    """
    data = data.astype(np.float64)
    bkg = sep.Background(data)
    data_sub = data - bkg

    # Extract sources
    sources = sep.extract(data_sub, threshold_sigma, err=bkg.globalrms)

    # Filter to point-like sources: A/B ratio (elongation) < 2.0, not too faint, not saturated
    point_sources = sources[
        (sources['a'] / sources['b'] < 2.0) &       # roughly round
        (sources['peak'] < 0.9 * (2**16 - 1)) &     # not saturated (assumes 16-bit)
        (sources['flux'] / sources['fluxerr'] > 5)   # SNR > 5
    ]

    return len(point_sources), point_sources

def check_star_count(star_count: int) -> list[str]:
    flags = []
    if star_count < 10:
        flags.append('STARCOUNT_CRITICAL')
    elif star_count < 20:
        flags.append('STARCOUNT_LOW')
    elif star_count > 500:
        flags.append('STARCOUNT_CROWDED')
    if star_count > 2000:
        flags.append('STARCOUNT_EXTREME_CROWDING')
    return flags

1.4 Background Sky Level

The sky background level (in ADU or electrons after gain correction) diagnoses transparency and twilight contamination.

Measurement method: Estimate the 2D background using photutils.background.Background2D with a median estimator on a grid. The median sky value across the image is the background level. The standard deviation of the background map is the sky noise.

Thresholds:

Absolute sky level (scaled to ADU/arcsec², instrument-dependent): The absolute ADU threshold is instrument-dependent. Instead, compare against the site's own historical distribution. For a new site with no history, use: - Sky > 30% of detector full well capacity in a single pixel: SKY_BRIGHT — twilight or heavy moonlight, photometric precision will be degraded. - Sky > 60% of detector full well: SKY_VERY_BRIGHT — hard reject for faint targets; marginal for bright targets.

A more portable (instrument-independent) check uses the estimated limiting magnitude: if the computed 5σ limiting magnitude is brighter than the target star by less than 2 mag, the sky is too bright for useful science.

Sky background spatial variation: Compute the coefficient of variation (CV = std/mean) of the background map tiles. - CV > 0.3 (30% variation): SKY_VARIABLE — patchy clouds passing through during the exposure. This is a critical flag because it produces correlated photometric noise that mimics astrophysical signals. - CV > 0.6: SKY_VERY_VARIABLE — hard reject. Clouds visible in the image itself.

Reference: The Huntsman Telescope (Macquarie/AAO) automated pipeline uses a background RMS threshold normalized by the expected sky noise to detect cloud contamination in wide-field photometry (Lao et al. 2021). The LCOGT BANZAI pipeline flags images where the background level exceeds a configurable percentile of the historical distribution for that site and filter.

from photutils.background import Background2D, MedianBackground
from astropy.stats import sigma_clipped_stats

def check_sky_background(data: np.ndarray, full_well_adu: float = 65535) -> tuple[list[str], dict]:
    """
    Measure sky background and flag problematic conditions.
    Returns (flags, metrics_dict).
    """
    flags = []

    bkg_estimator = MedianBackground()
    bkg = Background2D(data, box_size=64, filter_size=3, bkg_estimator=bkg_estimator)

    median_sky = float(np.median(bkg.background))
    sky_rms = float(np.std(bkg.background))
    sky_cv = sky_rms / median_sky if median_sky > 0 else 0
    sky_fraction = median_sky / full_well_adu

    if sky_fraction > 0.60:
        flags.append('SKY_VERY_BRIGHT')
    elif sky_fraction > 0.30:
        flags.append('SKY_BRIGHT')

    if sky_cv > 0.60:
        flags.append('SKY_VERY_VARIABLE')
    elif sky_cv > 0.30:
        flags.append('SKY_VARIABLE')

    return flags, {
        'sky_median_adu': median_sky,
        'sky_rms_adu': sky_rms,
        'sky_cv': sky_cv,
        'sky_fraction_fullwell': sky_fraction
    }

1.5 FWHM of PSF

The Full Width at Half Maximum of the Point Spread Function measures image quality: focus, seeing, and tracking.

Measurement method: Select the 20 brightest non-saturated, non-crowded stars in the field. Fit a 2D Gaussian to each using photutils.psf or astropy.modeling. Report the median FWHM and the elongation ratio (major axis / minor axis of the fitted Gaussian).

Convert FWHM from pixels to arcseconds: FWHM_arcsec = FWHM_pixels × pixel_scale. The pixel scale is taken from the plate solve result or the registered site parameters.

Thresholds (in arcseconds):

FWHM Interpretation Flag Action
< 1.0 arcsec Exceptional seeing (unlikely for amateur sites, possible from high-altitude sites) None Pass (possible lucky imaging)
1.0–3.0 arcsec Good to typical seeing for mid-latitude sites None Pass
3.0–5.0 arcsec Poor seeing or slight defocus PSF_POOR_FWHM Warning — photometric precision degraded
5.0–8.0 arcsec Bad seeing, significant defocus, or guiding problem PSF_BAD_FWHM Warning — photometry unreliable for crowded fields
> 8.0 arcsec Severely out of focus, trailing, or hardware problem PSF_REJECT_FWHM Hard reject

Reference: ZTF rejects images with FWHM > 4 arcsec (Bellm et al. 2019, PASP). LCOGT's BANZAI pipeline uses a per-site configurable FWHM threshold, typically set between 3–5 arcsec depending on the site's typical seeing. For sub-arcsecond photometry science (occultations), even the warning threshold should be treated as a reject. For bright-target photometry (V < 12), FWHM up to 8 arcsec can still produce useful differential photometry.

Elongation check: If the PSF major/minor axis ratio > 1.5, the star images are elongated — this indicates tracking failure, polar misalignment, or wind shake. Flag PSF_ELONGATED. If ratio > 2.5, flag PSF_SEVERE_TRAILING — hard reject for all precision photometry. Elongated stars have more sky background in their footprint and produce systematically noisy photometry even if the FWHM appears acceptable.

FWHM spatial variation: Measure FWHM at the centre vs. corners of the image. If corner FWHM is > 50% larger than centre FWHM, flag PSF_FIELD_CURVATURE. This indicates either optical field curvature or focus that is not optimised for the full field.

from photutils.psf import extract_stars, EPSFBuilder
from astropy.modeling import models, fitting
from astropy.nddata import NDData
from astropy.table import Table

def measure_psf_fwhm(data: np.ndarray, sources: np.ndarray,
                      pixel_scale_arcsec: float) -> tuple[list[str], dict]:
    """
    Measure PSF FWHM from bright non-saturated stars.
    Returns (flags, metrics_dict).
    """
    flags = []

    # Select 20 brightest non-saturated stars away from edges
    h, w = data.shape
    margin = 50  # pixels from edge
    good = sources[
        (sources['x'] > margin) & (sources['x'] < w - margin) &
        (sources['y'] > margin) & (sources['y'] < h - margin) &
        (sources['peak'] < 0.85 * 65535)  # not saturated
    ]
    good = good[np.argsort(good['flux'])[::-1]][:20]

    fwhms_pixels = []
    elongations = []

    for src in good:
        # Cut out a stamp around each star
        x, y = int(src['x']), int(src['y'])
        stamp = data[y-15:y+15, x-15:x+15]
        if stamp.shape != (30, 30):
            continue

        # Fit 2D Gaussian
        g_init = models.Gaussian2D(
            amplitude=stamp.max(),
            x_mean=15, y_mean=15,
            x_stddev=3, y_stddev=3
        )
        fit_g = fitting.LevMarLSQFitter()
        y_grid, x_grid = np.mgrid[0:30, 0:30]
        try:
            g = fit_g(g_init, x_grid, y_grid, stamp)
            sigma_x = abs(g.x_stddev.value)
            sigma_y = abs(g.y_stddev.value)
            fwhm = 2.355 * np.sqrt(sigma_x * sigma_y)  # geometric mean FWHM
            elongation = max(sigma_x, sigma_y) / min(sigma_x, sigma_y) if min(sigma_x, sigma_y) > 0 else 1
            fwhms_pixels.append(fwhm)
            elongations.append(elongation)
        except Exception:
            continue

    if not fwhms_pixels:
        flags.append('PSF_MEASUREMENT_FAILED')
        return flags, {}

    median_fwhm_px = float(np.median(fwhms_pixels))
    median_fwhm_arcsec = median_fwhm_px * pixel_scale_arcsec
    median_elongation = float(np.median(elongations))

    # FWHM flags
    if median_fwhm_arcsec > 8.0:
        flags.append('PSF_REJECT_FWHM')
    elif median_fwhm_arcsec > 5.0:
        flags.append('PSF_BAD_FWHM')
    elif median_fwhm_arcsec > 3.0:
        flags.append('PSF_POOR_FWHM')

    # Elongation flags
    if median_elongation > 2.5:
        flags.append('PSF_SEVERE_TRAILING')
    elif median_elongation > 1.5:
        flags.append('PSF_ELONGATED')

    return flags, {
        'fwhm_arcsec': median_fwhm_arcsec,
        'fwhm_pixels': median_fwhm_px,
        'elongation': median_elongation
    }

1.6 Saturation Check on Target Star

If the target star's peak pixel value exceeds the detector's saturation threshold, the measured flux is underestimated and the data is scientifically useless for photometry.

Measurement: Extract the peak pixel value within a 2×FWHM radius of the target's known position (from WCS).

Threshold: - Peak pixel > 90% of full well capacity: TARGET_SATURATED — hard reject. The ADU threshold depends on the instrument. Use SATURATE header keyword if present; otherwise default to 90% of (2^BITPIX − 1). For 16-bit cameras: 58,981 ADU. - Peak pixel > 75%: TARGET_NEAR_SATURATION — warning. Linearity may be compromised. The actual linearity cutoff varies by sensor: Sony IMX sensors are typically linear to ~85%; KAF sensors to ~70%.

Non-linearity note: Many CCD and CMOS sensors begin showing non-linearity well below the hard saturation point. For science-grade photometry, the safe working range is typically < 70% full well. Flag TARGET_NEAR_SATURATION at 70% to catch this.

Comparison stars: Run the same check on all comparison stars used in differential photometry. A saturated comparison star is as bad as a saturated target. Flag COMP_STAR_SATURATED if any comparison star peak > 90% full well.

Anti-saturation strategy: If TARGET_SATURATED is flagged, check whether the observer submitted data with multiple exposure times. If shorter exposures exist for the same epoch, prefer those. This is relevant for bright targets (V < 8) which many amateur cameras cannot image without saturating.

Stage 2: Photometric Quality Checks

These checks run after differential photometry has been extracted — either by the client-side software or by the server's own photometry pipeline (Section 8 of Stage 1 Data Pipeline Design.md).

2.1 Comparison Star Scatter (RMS vs. Expected Poisson Noise)

The key diagnostic: how much do the comparison stars scatter relative to how much they should scatter given Poisson statistics?

Expected noise model: For a star with flux F (in electrons), the expected per-measurement noise is: 

sigma_expected = sqrt(F + n_pix * (sky_per_pix + rdnoise²))
mag_noise_expected = 1.0857 / SNR    where SNR = F / sigma_expected

Comparison star RMS: For N comparison stars observed simultaneously, compute the RMS of their differential magnitudes (each comparison star measured relative to the ensemble of the others): 

RMS_comp = std(delta_mag for each comp star over this epoch)

Noise excess factor: 

noise_excess = RMS_comp / median(sigma_expected for comp stars)

Noise excess Interpretation Flag Action
< 1.5 Normal — scatter consistent with Poisson noise None Pass
1.5–2.5 Marginal — slight transparency variation, thin cirrus PHOT_EXCESS_NOISE Warning
2.5–5.0 Poor — clouds, guiding variation, focus drift PHOT_HIGH_NOISE Warning — include with caution
> 5.0 Very poor — obvious clouds or instrumental problem PHOT_REJECT_NOISE Hard reject

Minimum comparison star count: - < 3 comparison stars: PHOT_FEW_COMPS — warning. Single comparison stars can be variable; ensemble of 3+ is the minimum for robust differential photometry. - < 2 comparison stars: hard reject for science-critical targets (occultations, transit timing). For reconnaissance photometry, allow with warning.

Reference: The AAVSO photometry guide requires a minimum of 2 comparison stars (preferably 3+) with known magnitudes. The ExoClock network (a direct competitor/collaborator for TTV science) uses a noise excess threshold of ~2× expected Poisson noise to flag bad observations. The BANZAI pipeline at LCOGT computes the "reduced chi-squared" of the comparison star ensemble, rejecting images with chi² > 3 (corresponding roughly to a noise excess factor of ~1.7).

import numpy as np

def check_comparison_scatter(comp_star_mags: list[float], comp_star_errors: list[float]) -> tuple[list[str], dict]:
    """
    comp_star_mags: list of differential magnitudes for comparison stars
    comp_star_errors: list of expected photon noise errors for each comp star
    """
    flags = []
    n_comps = len(comp_star_mags)

    if n_comps < 2:
        flags.append('PHOT_TOO_FEW_COMPS')
        return flags, {}
    if n_comps < 3:
        flags.append('PHOT_FEW_COMPS')

    mags = np.array(comp_star_mags)
    errors = np.array(comp_star_errors)

    rms_actual = float(np.std(mags - np.mean(mags)))
    sigma_expected = float(np.median(errors))

    noise_excess = rms_actual / sigma_expected if sigma_expected > 0 else 0

    if noise_excess > 5.0:
        flags.append('PHOT_REJECT_NOISE')
    elif noise_excess > 2.5:
        flags.append('PHOT_HIGH_NOISE')
    elif noise_excess > 1.5:
        flags.append('PHOT_EXCESS_NOISE')

    return flags, {
        'comp_rms_mmag': rms_actual * 1000,
        'expected_noise_mmag': sigma_expected * 1000,
        'noise_excess_factor': noise_excess,
        'n_comp_stars': n_comps
    }

If the comparison star magnitudes show a systematic trend with airmass within a single observation session, the site's atmospheric extinction is not properly calibrated. This manifests as a fake slope in the target's light curve.

Measurement: For observations within a single night that span a significant airmass range (> 0.3 airmass units), fit a linear regression of (ensemble comparison star magnitude) vs. airmass. The expected slope is zero (after extinction correction). A non-zero slope indicates either: 1. No extinction correction was applied 2. The extinction coefficient used was wrong 3. Variable atmospheric extinction during the session (clouds)

Threshold: - Slope > 0.05 mag/airmass: PHOT_EXTINCTION_TREND — warning. The Bouguer extinction law gives ~0.15–0.3 mag/airmass for V band; if uncorrected, this slope will appear in any airmass-spanning dataset. Even partially applied, a residual > 0.05 mag/airmass is significant for precision photometry. - Slope > 0.15 mag/airmass: PHOT_SEVERE_EXTINCTION — hard reject for transit photometry (a 0.15 mag/airmass slope over a 2-hour transit spanning 0.5 airmass = 75 mmag systematic, completely swamping a 10 mmag transit depth).

from scipy import stats

def check_airmass_trend(airmass_values: list[float], delta_mags: list[float]) -> tuple[list[str], dict]:
    """
    Check for systematic trend of differential magnitude vs airmass.
    This uses comparison star magnitudes over multiple exposures within one session.
    """
    flags = []

    if len(airmass_values) < 5:
        return flags, {}  # Not enough points to fit a trend

    airmass = np.array(airmass_values)
    dmag = np.array(delta_mags)

    if np.ptp(airmass) < 0.3:
        return flags, {}  # Not enough airmass range to detect extinction slope

    slope, intercept, r_value, p_value, std_err = stats.linregress(airmass, dmag)

    if abs(slope) > 0.15:
        flags.append('PHOT_SEVERE_EXTINCTION')
    elif abs(slope) > 0.05:
        flags.append('PHOT_EXTINCTION_TREND')

    return flags, {
        'airmass_slope_mag_per_unit': float(slope),
        'airmass_slope_significance': float(abs(slope) / std_err) if std_err > 0 else 0,
        'airmass_range': float(np.ptp(airmass))
    }

2.3 Intra-Exposure Systematics (Rapid Variability in Comparison Stars)

When the comparison star ensemble shows rapid correlated variations within a single night's data, this indicates patchy clouds or other variable transparency effects. This is the single most dangerous source of false astrophysical signals in ground-based photometry.

Measurement: For each time bin within the session (e.g., every 3 consecutive exposures), compute the mean comparison star differential magnitude. A flat, near-zero sequence is good. Rapid variations indicate transparency changes.

Threshold: Compute the running standard deviation of the comparison star ensemble over 5-point windows: - Ensemble scatter > 5 mmag on timescales < 10 minutes: PHOT_RAPID_SYSTEMATICS — warning. Thin cirrus. Any astrophysical signal on similar timescales (fast variable, short-period occultation) is untrustworthy. - Ensemble scatter > 20 mmag on timescales < 10 minutes: PHOT_CLOUD_CONTAMINATION — hard reject.

Correlation check: Compute the Pearson correlation between the target's differential magnitude and the comparison ensemble magnitude. If |r| > 0.7, the target signal is strongly correlated with the comparison star systematics — i.e., it is not astrophysical. Flag PHOT_TARGET_CORRELATED_WITH_SYSTEMATICS — hard reject for science use.

from scipy.stats import pearsonr

def check_intranight_systematics(target_dmags: np.ndarray, ensemble_dmags: np.ndarray,
                                  timestamps: np.ndarray) -> tuple[list[str], dict]:
    """
    Check for rapid correlated systematics between target and comparison ensemble.
    """
    flags = []

    if len(target_dmags) < 5:
        return flags, {}

    # Running scatter of comparison ensemble
    window = 5
    running_std = []
    for i in range(len(ensemble_dmags) - window):
        running_std.append(np.std(ensemble_dmags[i:i+window]))

    max_running_std_mmag = float(np.max(running_std)) * 1000 if running_std else 0

    if max_running_std_mmag > 20:
        flags.append('PHOT_CLOUD_CONTAMINATION')
    elif max_running_std_mmag > 5:
        flags.append('PHOT_RAPID_SYSTEMATICS')

    # Correlation between target and ensemble
    if len(target_dmags) == len(ensemble_dmags) and len(target_dmags) >= 5:
        r, p = pearsonr(target_dmags, ensemble_dmags)
        if abs(r) > 0.7 and p < 0.05:
            flags.append('PHOT_TARGET_CORRELATED_WITH_SYSTEMATICS')

    return flags, {
        'max_running_comp_scatter_mmag': max_running_std_mmag,
        'target_ensemble_correlation': float(r) if 'r' in dir() else None
    }

2.4 Light Curve Outlier Detection (Sigma Clipping)

Individual data points in a light curve that deviate significantly from the local trend are outliers — potentially bad exposures not caught by earlier checks, cosmic rays, or satellite trails affecting the target aperture.

Method: Apply iterative sigma clipping using astropy.stats.sigma_clip. For a stable source (comparison star or check star), the data should be consistent. For a variable target, apply sigma clipping relative to a smoothed version of the light curve (using a running median).

Threshold: 3.5σ is the standard threshold used by AAVSO, ExoClock, and the majority of published light curve pipelines. At 3σ, there is a 0.3% chance of a genuine good data point being rejected — for a session of 100 exposures, you would expect to clip ~0.3 real points. At 4σ the false rejection rate drops to near zero but true outliers may slip through.

Recommended values: - 3.5σ clip: Flag individual exposures as POINT_SIGMA_OUTLIER. These are excluded from the ensemble average but retained in the database. - 5σ clip: Flag as POINT_EXTREME_OUTLIER. These are excluded even from human review unless specifically requested. - If > 10% of data points in a session are clipped: flag the entire session with SESSION_HIGH_CLIP_FRACTION — the underlying data quality is poor.

from astropy.stats import sigma_clip
import numpy as np

def clip_outliers(times: np.ndarray, mags: np.ndarray,
                  sigma: float = 3.5) -> tuple[np.ndarray, list[int]]:
    """
    Apply sigma clipping to a light curve.
    Returns (mask of good points, indices of rejected points).
    """
    clipped = sigma_clip(mags, sigma=sigma, maxiters=5, masked=True)
    good_mask = ~clipped.mask
    outlier_indices = list(np.where(clipped.mask)[0])

    clip_fraction = len(outlier_indices) / len(mags)
    session_flags = []
    if clip_fraction > 0.10:
        session_flags.append('SESSION_HIGH_CLIP_FRACTION')

    return good_mask, outlier_indices, session_flags

Stage 3: Time-Series Quality Checks

These run as a nightly batch job after assembling the combined light curve from multiple nights/sites. They do not require raw FITS — only the assembled light_curve_points table.

3.1 Gap Detection

A time-domain science pipeline needs to know when coverage is missing. Gaps affect period searches (aliasing), TTV completeness, and occultation chord coverage.

Method: Given a time-sorted array of observation timestamps, compute the time differences between successive points. Compare against the expected cadence from the targets table (cadence_hours).

Thresholds: - Gap > 2× expected cadence: GAP_DETECTED — informational. Note the gap in the light curve metadata. - Gap > 10× expected cadence: GAP_SIGNIFICANT — warning. Period searches may be unreliable across this gap. - Gap > 365 days: GAP_SEASON_BREAK — informational. Expected seasonal gap when the target is behind the Sun. The pipeline should not treat this as anomalous. - Total coverage fraction (fraction of target window with observations) < 50%: GAP_POOR_COVERAGE — warning. Science objectives may not be achievable.

Gap reporting format: Each gap is stored as a separate metadata record with start time, end time, and duration. This allows the dashboard to show observers which time windows need filling.

def detect_gaps(times_bjd: np.ndarray, cadence_hours: float) -> tuple[list[str], list[dict]]:
    """
    Detect gaps in a time series relative to expected cadence.
    Returns session-level flags and a list of gap records.
    """
    if len(times_bjd) < 2:
        return ['GAP_INSUFFICIENT_DATA'], []

    times = np.sort(times_bjd)
    diffs_hours = np.diff(times) * 24  # BJD in days → hours

    gaps = []
    session_flags = []

    for i, dt in enumerate(diffs_hours):
        if dt > 2 * cadence_hours:
            severity = 'significant' if dt > 10 * cadence_hours else 'normal'
            season_break = dt > 365 * 24  # > 1 year = seasonal gap
            gaps.append({
                'start_bjd': float(times[i]),
                'end_bjd': float(times[i+1]),
                'duration_hours': float(dt),
                'severity': 'season_break' if season_break else severity
            })
            if season_break:
                session_flags.append('GAP_SEASON_BREAK')
            elif dt > 10 * cadence_hours:
                session_flags.append('GAP_SIGNIFICANT')
            else:
                session_flags.append('GAP_DETECTED')

    # Overall coverage fraction
    total_span_hours = (times[-1] - times[0]) * 24
    total_gap_hours = sum(g['duration_hours'] - cadence_hours for g in gaps if g['severity'] != 'season_break')
    coverage_fraction = 1 - (total_gap_hours / total_span_hours) if total_span_hours > 0 else 0
    if coverage_fraction < 0.5:
        session_flags.append('GAP_POOR_COVERAGE')

    # Deduplicate flags
    session_flags = list(set(session_flags))

    return session_flags, gaps

3.2 Period Aliasing Risk Assessment

The 24-hour sampling cadence imposed by Earth's rotation creates aliases in period searches. For OpenAstro, the geographic distribution is supposed to eliminate daytime gaps — but in practice, a network with poor longitude coverage will still show aliasing.

Method: Compute the Lomb-Scargle window function (the power spectrum of the sampling pattern alone, without any signal). Peaks in the window function at a frequency f indicate that a real signal at frequency f±(1/day) will be confused with f.

Alias risk thresholds: - Peak window function power at 1/day and harmonics > 0.3: ALIAS_DAILY_HIGH — the 24-hour alias is strong. Period searches near 1 day, 0.5 day, 2 days must be treated with extreme suspicion. - Peak window function power > 0.5: ALIAS_DAILY_CRITICAL — the dataset is probably dominated by single-site observations and is essentially a 1-day aliased sample.

Longitude coverage diagnostic: Count the number of unique longitude bins (each 30° = 2 hours) represented in the dataset. For a period search to be alias-free: - Fewer than 3 unique longitude bins: ALIAS_POOR_LONGITUDE_COVERAGE — the network is not providing the geographic diversity advantage. - 6+ unique longitude bins spread over at least 4 time zones: no alias flag — coverage is adequate.

from astropy.timeseries import LombScargle
import numpy as np

def assess_aliasing_risk(times_bjd: np.ndarray, site_longitudes: list[float]) -> tuple[list[str], dict]:
    """
    Assess aliasing risk from window function and longitude coverage.
    """
    flags = []

    # Compute window function: LS of constant signal = sampling pattern
    ones = np.ones_like(times_bjd)
    ls = LombScargle(times_bjd, ones, fit_mean=False, center_data=False)
    freqs = np.linspace(0.1, 3.0, 3000)  # cycles/day
    power = ls.power(freqs)

    # Check power near 1/day alias
    alias_mask = np.abs(freqs - 1.0) < 0.05  # within 0.05 c/day of 1/day
    peak_alias_power = float(power[alias_mask].max()) if alias_mask.any() else 0

    if peak_alias_power > 0.5:
        flags.append('ALIAS_DAILY_CRITICAL')
    elif peak_alias_power > 0.3:
        flags.append('ALIAS_DAILY_HIGH')

    # Longitude coverage
    lon_bins = set(int(lon // 30) for lon in site_longitudes)  # 30° bins = 2-hour bins
    n_lon_bins = len(lon_bins)
    if n_lon_bins < 3:
        flags.append('ALIAS_POOR_LONGITUDE_COVERAGE')

    return flags, {
        'peak_daily_alias_power': peak_alias_power,
        'n_longitude_bins': n_lon_bins
    }

3.3 Red Noise Characterisation

Red noise — correlated, low-frequency noise — is the enemy of ground-based photometric precision. It arises from atmospheric scintillation, temperature-driven focus drift, and imperfect flat fielding. Unlike white (Poisson) noise, it does not average down as 1/√N.

Why it matters: If red noise is present, a light curve made of N exposures has an effective noise of sigma_red (not sigma_white / √N). This systematically underestimates measurement uncertainties and inflates the significance of spurious detections.

Measurement: Compute the variance of the binned light curve as a function of bin size. For pure white noise, variance decreases as 1/N (bin size). For red noise, it decreases more slowly. Fit the Allan deviation curve: 

sigma²(N) = sigma²_white / N + sigma²_red
where sigma_red is the floor that the noise cannot be averaged below.

Threshold: - sigma_red > 5 mmag: RED_NOISE_MODERATE — warning. The effective noise floor is 5 mmag regardless of how many exposures are taken. For targets with transit depths < 10 mmag (Earth-sized transits of Sun-like stars), this makes detection impossible. - sigma_red > 15 mmag: RED_NOISE_HIGH — hard warning. Only bright targets with large transit/variability amplitudes are tractable. - sigma_red is stored as a numeric column (red_noise_floor_mmag) in the site performance table — useful for site comparison.

def characterise_red_noise(times_bjd: np.ndarray, mags: np.ndarray,
                            mag_errors: np.ndarray) -> tuple[list[str], dict]:
    """
    Estimate the red noise floor using bin size vs. variance scaling.
    Uses comparison star light curve (should be constant) if available.
    """
    flags = []

    # Bin the light curve at increasing timescales
    bin_sizes = [1, 2, 4, 8, 16, 32]
    variances = []

    for bs in bin_sizes:
        n_bins = len(mags) // bs
        if n_bins < 3:
            break
        binned = [np.mean(mags[i*bs:(i+1)*bs]) for i in range(n_bins)]
        variances.append(np.var(binned))

    if len(variances) < 3:
        return flags, {}

    # Fit sigma²_white / N + sigma²_red
    from scipy.optimize import curve_fit
    def noise_model(N, sigma_white_sq, sigma_red_sq):
        return sigma_white_sq / N + sigma_red_sq

    try:
        popt, _ = curve_fit(noise_model, bin_sizes[:len(variances)], variances,
                            p0=[variances[0], 0.0001], bounds=(0, np.inf))
        sigma_red_mmag = float(np.sqrt(popt[1])) * 1000
    except Exception:
        sigma_red_mmag = 0

    if sigma_red_mmag > 15:
        flags.append('RED_NOISE_HIGH')
    elif sigma_red_mmag > 5:
        flags.append('RED_NOISE_MODERATE')

    return flags, {'red_noise_floor_mmag': sigma_red_mmag}

Stage 4: Quality Flag Schema

4.1 Integer Quality Score

The quality_score column in the observations table stores an integer from 0–3:

Score Meaning Science pipeline inclusion
0 Clean — no flags raised Included in all products
1 Warning — one or more warning-level flags Included but marked; excluded from precision products (TTV timing, occultation timing)
2 Degraded — one or more degraded-level flags Excluded from default products; available via API with include_degraded=true
3 Rejected — one or more hard-reject flags Excluded from all products automatically; stored for human review

The score is computed as the maximum severity of all raised flags.

4.2 Complete Flag Registry

Each flag has an assigned severity (W = warning, D = degraded, R = reject) and a stage.

Stage 1 — Input Validation:

Flag Severity Description
HDR_NO_TIMESTAMP R FITS header contains no parseable timestamp
HDR_TIMESTAMP_FUTURE R Observation timestamp is in the future
HDR_TIMESTAMP_ANCIENT R Timestamp before 2010-01-01 (likely parsing error)
HDR_TIMESTAMP_STALE W Upload received > 1h after expected end of exposure
HDR_TIMESTAMP_UNPARSEABLE R Timestamp string cannot be parsed
HDR_NO_EXPTIME R Exposure time not in header
HDR_EXPTIME_ZERO R Exposure time ≤ 0
HDR_NO_FILTER W Filter name not in header
HDR_NO_INSTRUMENT W No telescope/instrument identifier
HDR_NO_SITE_COORDS W Site latitude/longitude not in header
HDR_COORDS_MISMATCH W Header coords differ from registered site by > 0.1°
HDR_NO_GAIN W Detector gain not in header (photon noise estimate degraded)
HDR_NO_RDNOISE W Read noise not in header
HDR_WCS_CORRUPT W WCS keywords present but self-inconsistent
SOLVE_FAILED R Plate solve failed with both ASTAP and astrometry.net
SOLVE_FAILED_COORDS_ASSUMED W Plate solve failed; using client-supplied coordinates
SOLVE_WRONG_FIELD R Solved field centre > 5 arcmin from target
SOLVE_POINTING_OFFSET W Solved field centre 1–5 arcmin from target
SOLVE_POOR_RESIDUALS W Plate solve residuals > 2 arcsec per star
SOLVE_PIXSCALE_MISMATCH W Pixel scale differs from registered value by > 20%
STARCOUNT_CRITICAL R Fewer than 10 point sources detected
STARCOUNT_LOW W 10–20 point sources (thin clouds or poor transparency)
STARCOUNT_CROWDED W > 500 sources (crowded field; blending risk)
STARCOUNT_EXTREME_CROWDING W > 2000 sources (PSF photometry required)
SKY_VERY_BRIGHT R Background > 60% full well (twilight/bright moon)
SKY_BRIGHT W Background > 30% full well (bright conditions)
SKY_VERY_VARIABLE R Background spatial CV > 60% (clouds in image)
SKY_VARIABLE W Background spatial CV > 30% (variable transparency)
PSF_REJECT_FWHM R FWHM > 8 arcsec
PSF_BAD_FWHM D FWHM 5–8 arcsec
PSF_POOR_FWHM W FWHM 3–5 arcsec
PSF_MEASUREMENT_FAILED W Could not measure PSF from available stars
PSF_SEVERE_TRAILING R Elongation > 2.5 (severe tracking failure)
PSF_ELONGATED W Elongation 1.5–2.5 (tracking degraded)
PSF_FIELD_CURVATURE W Corner FWHM > 50% larger than centre FWHM
TARGET_SATURATED R Target peak pixel > 90% full well
TARGET_NEAR_SATURATION W Target peak pixel 70–90% full well
COMP_STAR_SATURATED R One or more comparison stars > 90% full well

Stage 2 — Photometric Quality:

Flag Severity Description
PHOT_TOO_FEW_COMPS R Fewer than 2 comparison stars
PHOT_FEW_COMPS W 2 comparison stars (minimum; prefer 3+)
PHOT_REJECT_NOISE R Comparison scatter > 5× expected Poisson noise
PHOT_HIGH_NOISE D Comparison scatter 2.5–5× expected noise
PHOT_EXCESS_NOISE W Comparison scatter 1.5–2.5× expected noise
PHOT_SEVERE_EXTINCTION R Airmass slope > 0.15 mag/airmass
PHOT_EXTINCTION_TREND W Airmass slope 0.05–0.15 mag/airmass
PHOT_CLOUD_CONTAMINATION R Running comparison scatter > 20 mmag on <10 min timescale
PHOT_RAPID_SYSTEMATICS W Running comparison scatter 5–20 mmag on <10 min timescale
PHOT_TARGET_CORRELATED_WITH_SYSTEMATICS R Target light curve Pearson r > 0.7 with comparison ensemble
POINT_SIGMA_OUTLIER D Individual data point > 3.5σ from local median
POINT_EXTREME_OUTLIER R Individual data point > 5σ from local median
SESSION_HIGH_CLIP_FRACTION D > 10% of session data points clipped

Stage 3 — Time-Series Quality:

Flag Severity Description
GAP_DETECTED W Gap > 2× expected cadence
GAP_SIGNIFICANT W Gap > 10× expected cadence
GAP_SEASON_BREAK W Gap > 365 days (expected seasonal gap — informational)
GAP_POOR_COVERAGE W Coverage fraction < 50% of target observation window
GAP_INSUFFICIENT_DATA W Fewer than 3 data points for period analysis
ALIAS_DAILY_HIGH W Window function peak power > 0.3 at 1/day alias
ALIAS_DAILY_CRITICAL D Window function peak power > 0.5 at 1/day alias
ALIAS_POOR_LONGITUDE_COVERAGE W Fewer than 3 longitude bins in dataset (geographic diversity insufficient)
RED_NOISE_MODERATE W Red noise floor 5–15 mmag
RED_NOISE_HIGH D Red noise floor > 15 mmag

Science-case-specific flags (applied by campaign-aware logic):

Flag Severity When Applied
TIMING_NO_GPS D Site has no GPS time sync; occultation or TTV campaign
TIMING_OUTSIDE_WINDOW R Observation > 5 min from event window (occultation/transit)
CROSS_SITE_OUTLIER D Magnitude deviates > 4σ from same-epoch observations at other sites
NO_COLOR_CORRECTION W B-V color correction not applied due to missing B observation

4.3 Score Computation Rules

FLAG_SEVERITY = {
    # R = 3 (reject), D = 2 (degraded), W = 1 (warning)
    # Stage 1
    'HDR_NO_TIMESTAMP': 3, 'HDR_TIMESTAMP_FUTURE': 3, 'HDR_TIMESTAMP_ANCIENT': 3,
    'HDR_TIMESTAMP_UNPARSEABLE': 3, 'HDR_NO_EXPTIME': 3, 'HDR_EXPTIME_ZERO': 3,
    'SOLVE_FAILED': 3, 'SOLVE_WRONG_FIELD': 3,
    'STARCOUNT_CRITICAL': 3, 'SKY_VERY_BRIGHT': 3, 'SKY_VERY_VARIABLE': 3,
    'PSF_REJECT_FWHM': 3, 'PSF_SEVERE_TRAILING': 3,
    'TARGET_SATURATED': 3, 'COMP_STAR_SATURATED': 3,
    # Stage 1 degraded
    'PSF_BAD_FWHM': 2,
    # Stage 1 warnings
    'HDR_TIMESTAMP_STALE': 1, 'HDR_NO_FILTER': 1, 'HDR_NO_INSTRUMENT': 1,
    'HDR_NO_SITE_COORDS': 1, 'HDR_COORDS_MISMATCH': 1, 'HDR_NO_GAIN': 1,
    'HDR_NO_RDNOISE': 1, 'HDR_WCS_CORRUPT': 1, 'SOLVE_FAILED_COORDS_ASSUMED': 1,
    'SOLVE_POINTING_OFFSET': 1, 'SOLVE_POOR_RESIDUALS': 1, 'SOLVE_PIXSCALE_MISMATCH': 1,
    'STARCOUNT_LOW': 1, 'STARCOUNT_CROWDED': 1, 'STARCOUNT_EXTREME_CROWDING': 1,
    'SKY_BRIGHT': 1, 'SKY_VARIABLE': 1, 'PSF_POOR_FWHM': 1, 'PSF_MEASUREMENT_FAILED': 1,
    'PSF_ELONGATED': 1, 'PSF_FIELD_CURVATURE': 1, 'TARGET_NEAR_SATURATION': 1,
    # Stage 2
    'PHOT_TOO_FEW_COMPS': 3, 'PHOT_REJECT_NOISE': 3, 'PHOT_SEVERE_EXTINCTION': 3,
    'PHOT_CLOUD_CONTAMINATION': 3, 'PHOT_TARGET_CORRELATED_WITH_SYSTEMATICS': 3,
    'POINT_EXTREME_OUTLIER': 3,
    'PHOT_HIGH_NOISE': 2, 'SESSION_HIGH_CLIP_FRACTION': 2, 'POINT_SIGMA_OUTLIER': 2,
    'PHOT_FEW_COMPS': 1, 'PHOT_EXCESS_NOISE': 1, 'PHOT_EXTINCTION_TREND': 1,
    'PHOT_RAPID_SYSTEMATICS': 1,
    # Stage 3
    'ALIAS_DAILY_CRITICAL': 2, 'RED_NOISE_HIGH': 2,
    'GAP_DETECTED': 1, 'GAP_SIGNIFICANT': 1, 'GAP_SEASON_BREAK': 1,
    'GAP_POOR_COVERAGE': 1, 'ALIAS_DAILY_HIGH': 1, 'ALIAS_POOR_LONGITUDE_COVERAGE': 1,
    'RED_NOISE_MODERATE': 1,
    # Science-specific
    'TIMING_NO_GPS': 2, 'CROSS_SITE_OUTLIER': 2,
    'TIMING_OUTSIDE_WINDOW': 3,
    'NO_COLOR_CORRECTION': 1,
}

def compute_quality_score(flags: list[str]) -> int:
    """Compute overall quality score as max severity of all raised flags."""
    if not flags:
        return 0
    return max(FLAG_SEVERITY.get(f, 1) for f in flags)

4.4 Science-Case Override Logic

Some flags that are merely warnings for general photometry are hard rejects for specific science cases. The pipeline applies a second pass after the base score is computed:

SCIENCE_OVERRIDES = {
    'occultation': {
        # Any timing uncertainty is a reject for occultations
        'TIMING_NO_GPS': 3,
        'PSF_POOR_FWHM': 2,  # Upgrade from W to D for occultations
        'PSF_BAD_FWHM': 3,   # Upgrade from D to R
    },
    'transit_timing': {
        # TTV work needs clean photometry
        'PHOT_EXTINCTION_TREND': 2,   # Upgrade W to D
        'PHOT_RAPID_SYSTEMATICS': 3,  # Upgrade W to R
        'TIMING_NO_GPS': 2,
    },
    'frb_followup': {
        # Fast transient — need fast readout, any trailing is disqualifying
        'PSF_ELONGATED': 3,
    }
}

def apply_science_overrides(flags: list[str], campaign_type: str) -> int:
    overrides = SCIENCE_OVERRIDES.get(campaign_type, {})
    severities = []
    for flag in flags:
        severity = overrides.get(flag, FLAG_SEVERITY.get(flag, 1))
        severities.append(severity)
    return max(severities) if severities else 0

Stage 5: Implementation Notes

5.1 Python Library Stack

Task Library Notes
FITS I/O and header parsing astropy.io.fits Standard; handles all FITS variants
WCS coordinate transforms astropy.wcs Required for pixel-to-sky and sky-to-pixel
Source detection sep (SExtractor Python) 10–100× faster than photutils for large images; use for production
Source detection (backup/debug) photutils.detection.DAOStarFinder Slower but pure Python, easier to install
Background estimation photutils.background.Background2D 2D background mesh; handles non-uniform sky
PSF fitting astropy.modeling.models.Gaussian2D For FWHM measurement from bright stars
Aperture photometry photutils.aperture CircularAperture + CircularAnnulus
Statistical tools astropy.stats.sigma_clip, scipy.stats sigma clipping, linear regression
Period analysis / window function astropy.timeseries.LombScargle Built-in; handles uneven sampling
Time conversion (UTC → BJD_TDB) astropy.time.Time Essential for TTV science
Coordinate matching astropy.coordinates.SkyCoord.match_to_catalog_sky Matching sources to APASS/Gaia
Plate solving (via CLI) subprocess wrapping ASTAP CLI ASTAP is not a Python library; call via subprocess
Curve fitting (red noise model) scipy.optimize.curve_fit For Allan deviation fitting

Install list for QC pipeline: 

pip install astropy photutils sep scipy numpy
# ASTAP: install binary separately from https://www.hnsky.org/astap.htm

sep vs photutils decision: Use sep for the production QC pipeline. The sep library is a C-level implementation of the SExtractor algorithm with a Python interface. For a 4096×4096 pixel image, sep.extract() runs in ~0.1s versus ~5s for photutils.DAOStarFinder. Given that QC runs on every uploaded image, this matters at scale.

5.2 Where QC Runs

Server-side on upload (Stages 1 and 2): The QC pipeline runs as a FastAPI background task triggered when a FITS file is received at POST /api/v1/observations. The full FITS file must be present. Use FastAPI's BackgroundTasks for this:

from fastapi import BackgroundTasks

@app.post("/api/v1/observations")
async def submit_observation(
    fits_file: UploadFile,
    metadata: ObservationMetadata,
    background_tasks: BackgroundTasks
):
    # Save FITS to temporary storage
    fits_path = await save_fits_upload(fits_file)

    # Create preliminary observation record with pending status
    obs_id = await db.create_observation(metadata, quality_score=-1)

    # Queue QC as background task (doesn't block the API response)
    background_tasks.add_task(run_qc_pipeline, obs_id, fits_path, metadata)

    return {"obs_id": str(obs_id), "status": "accepted", "qc_status": "pending"}

async def run_qc_pipeline(obs_id: UUID, fits_path: str, metadata: ObservationMetadata):
    flags = []
    metrics = {}

    with fits.open(fits_path) as hdul:
        header = hdul[0].header
        data = hdul[0].data.astype(np.float64)

    # Stage 1 checks
    flags += check_header_completeness(header, metadata.site_lat, metadata.site_lon)

    solve_result = await run_plate_solve(fits_path, metadata.target_ra, metadata.target_dec)
    flags += check_plate_solve(solve_result, metadata.target_ra, metadata.target_dec, metadata.registered_pixscale)

    star_count, sources = count_stars(data)
    flags += check_star_count(star_count)

    sky_flags, sky_metrics = check_sky_background(data)
    flags += sky_flags
    metrics.update(sky_metrics)

    pixel_scale = solve_result.get('pixel_scale_arcsec', metadata.registered_pixscale)
    psf_flags, psf_metrics = measure_psf_fwhm(data, sources, pixel_scale)
    flags += psf_flags
    metrics.update(psf_metrics)

    # Saturation check requires target pixel position from WCS
    if solve_result.get('success'):
        sat_flags = check_saturation(data, header, solve_result, metadata.target_ra, metadata.target_dec)
        flags += sat_flags

    quality_score = compute_quality_score(flags)

    # Update observation record
    await db.update_observation_qc(obs_id, flags=flags, score=quality_score, metrics=metrics)

Nightly batch (Stage 3): Run as a cron job at 08:00 UTC daily (after all night-time observations worldwide have been uploaded):

# cron: 0 8 * * * python run_timeseries_qc.py

def run_timeseries_qc():
    """Run Stage 3 QC on all active campaign targets."""
    active_targets = db.query(Target).filter(Target.active == True).all()

    for target in active_targets:
        campaign_type = target.primary_campaign_type

        # Fetch assembled light curve
        lc = db.fetch_light_curve(target.target_id, quality_score_max=2)
        if len(lc) < 3:
            continue

        times = np.array([p.bjd_tdb for p in lc])
        mags = np.array([p.magnitude for p in lc])
        mag_errors = np.array([p.mag_error for p in lc])
        site_lons = [p.site_longitude for p in lc]

        ts_flags = []
        ts_flags_gap, gaps = detect_gaps(times, target.cadence_hours)
        ts_flags += ts_flags_gap

        alias_flags, _ = assess_aliasing_risk(times, site_lons)
        ts_flags += alias_flags

        rn_flags, rn_metrics = characterise_red_noise(times, mags, mag_errors)
        ts_flags += rn_flags

        db.update_target_timeseries_flags(target.target_id, ts_flags, gaps, rn_metrics)

5.3 Database Schema Additions

The existing observations table in Details on componets.md has quality_flags VARCHAR(64)[]. Add these columns:

ALTER TABLE observations ADD COLUMN quality_score INTEGER DEFAULT -1;
-- -1 = QC pending, 0 = clean, 1 = warning, 2 = degraded, 3 = rejected

ALTER TABLE observations ADD COLUMN qc_metrics JSONB;
-- Stores numeric QC results: fwhm_arcsec, sky_median_adu, star_count,
-- noise_excess_factor, comp_rms_mmag, etc.

ALTER TABLE observations ADD COLUMN qc_completed_at TIMESTAMP;
-- When QC finished (NULL if pending)

Add a new table for time-series QC results (per target, not per observation): 

CREATE TABLE target_timeseries_qc (
    target_id UUID REFERENCES targets,
    computed_at TIMESTAMP DEFAULT NOW(),
    ts_flags VARCHAR(64)[],
    gaps JSONB,              -- list of {start_bjd, end_bjd, duration_hours, severity}
    red_noise_floor_mmag DOUBLE PRECISION,
    n_longitude_bins INTEGER,
    daily_alias_power DOUBLE PRECISION,
    coverage_fraction DOUBLE PRECISION,
    PRIMARY KEY (target_id, computed_at)
);

Add an index on quality_score since the light curve assembly query (quality_score = 0) will run very frequently: 

CREATE INDEX idx_obs_quality ON observations (quality_score);
CREATE INDEX idx_obs_target_quality ON observations (target_id, quality_score, observed_at);

5.4 Client-Side Pre-checks (Optional but Encouraged)

A subset of Stage 1 checks can run client-side before upload to save bandwidth and give observers immediate feedback. Recommended client-side checks:

  1. Header completeness — trivial, zero compute cost
  2. Star count — fast with sep, runs in < 1s
  3. FWHM — fast, gives immediate focus feedback to the observer
  4. Saturation — fast pixel max check

Checks that should remain server-side only: - Plate solving (ASTAP index files are large; expensive to distribute to all clients) - Cross-site consistency (requires database access) - Airmass trend (requires full session data) - Time-series checks (require the assembled database)

If the client reports its own QC metrics in the upload payload, the server can skip re-running those checks and trust the client values (with a random re-check rate of ~5% for audit purposes).

Summary: Threshold Reference Card

Check Warn threshold Reject threshold
FWHM > 3 arcsec > 8 arcsec
PSF elongation > 1.5 > 2.5
Star count (min) < 20 < 10
Star count (max) > 500 —
Background level > 30% full well > 60% full well
Background CV > 30% > 60%
Target saturation > 70% full well > 90% full well
Comparison scatter > 1.5× Poisson > 5× Poisson
Airmass slope > 0.05 mag/airmass > 0.15 mag/airmass
Comparison running scatter > 5 mmag / 10 min > 20 mmag / 10 min
Target-ensemble correlation — Pearson r > 0.7
Sigma clip threshold — 3.5σ (individual points)
Plate solve field offset > 1 arcmin > 5 arcmin
Plate solve residuals > 2 arcsec —
Red noise floor > 5 mmag > 15 mmag
Minimum comparison stars < 3 < 2

Relationship to Existing Vault Documents

  • Details on componets.md: The quality_flags field this document populates is defined there. The quality_score column is a new addition.
  • Stage 1 Data Pipeline Design.md: Section 7 of that document has an earlier, simpler QC list. This document supersedes and expands it — the flags defined here replace those in Section 7.
  • Plate Solving on Server.md: The plate solve result (success/failure, residuals, pixel scale, field centre) feeds directly into Section 1.2 of this document.
  • failure_analysis_guide.md: Section 5.1 ("bad photometry") and Section 5.1 ("timing errors") identify the failure modes this pipeline catches. The quality control here is the direct engineering response to those failure modes.
  • Gap Analysis.md: The open question "Quality control automation: how do we flag bad data?" in the Tech — Unresolved table is answered by this document.

Last updated: 2026-03-21 Status: Design document — not yet implemented Implementation priority: Phase 2 (50+ sites). Phase 1 can use a simplified subset: header check + star count + FWHM only. First implementation target: Stages 1 and 2 header/FWHM/star count checks as FastAPI background task on upload.