209.1. Object table flags

209.1. Object table flags¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: r30.0.11
Last verified to run: 2026-10-28
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: To characterize the flag columns of the Object table and measure which objects a set of flag cuts removes, in fields of different depth and stellar density.

LSST data products: Object

Packages: numpy, matplotlib, lsst.rsp.RSPDiscovery, lsst.utils.plotting

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

A flag column is a boolean column that records a condition found during processing: a measurement failed, a pixel in the footprint carried a mask bit, or a special case was identified. A value of True (1 in ADQL) means the condition occurred.

The flag columns are defined in the catalog flag columns pages on dp2.lsst.io. Consult those pages for what each flag means.

Related tutorials: The contents of the Object table are described in the 200-level tutorial on the Object table. The image mask planes behind the pixelFlags_* columns are described in the 200-level tutorial on pixel mask planes. The TAP service is introduced in the 100-level "How to" tutorials.

1.1. Import packages¶

Import numpy for arrays and matplotlib for plots.

From the lsst package, import RSPDiscovery for the data services and set_rubin_plotstyle for the Rubin plot style.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from lsst.rsp import RSPDiscovery
from lsst.utils.plotting import set_rubin_plotstyle

1.2. Define parameters and functions¶

1.2.1. Instantiate the Table Access Protocol (TAP) service¶

Instantiate RSPDiscovery with the DP2 data release and get the TAP client.

In [2]:

discovery = RSPDiscovery("dp2")
service = discovery.get_tap_client()

1.2.2. Set the plot style and the band¶

Set up the default Rubin plotstyle.

Set the band used throughout the notebook.

In [3]:

set_rubin_plotstyle()

band = 'i'

1.2.3. Define functions¶

Four functions are used more than once and are defined here. Single-use code is defined where it is used.

Function: run_query

Submit an Astronomical Data Query Language (ADQL) query as an asynchronous TAP job, wait for it, and return the results as an astropy table.

Asynchronous jobs are the best practice for TAP queries. Deleting the job releases the resources it holds on the server.

In [4]:

def run_query(query):
    """
    Submit an ADQL query as an asynchronous TAP job and return the results.

    Parameters
    ----------
    query: str
        The ADQL query to execute.

    Returns
    -------
    table: astropy.table.Table
        The results of the query.
    """
    job = service.submit_job(query)
    job.run()
    job.wait(phases=['COMPLETED', 'ERROR'])
    print('Job phase is', job.phase)
    if job.phase == 'ERROR':
        job.raise_if_error()
    assert job.phase == 'COMPLETED'
    table = job.fetch_result().to_table()
    job.delete()
    return table

Function: unpack_columns

Return a dictionary of the table's columns as plain numpy arrays, with missing entries filled.

TAP returns masked columns, and a comparison involving a masked entry returns a masked value rather than True or False, which can drop rows from a selection with no error raised. Filling the columns once makes the treatment explicit: boolean columns take False, meaning the condition was not recorded, and numeric columns take NaN, which numpy.isfinite detects.

In [5]:

def unpack_columns(table):
    """
    Return the columns of a table as plain arrays, with masked entries
      filled.

    Parameters
    ----------
    table: astropy.table.Table
        The table to unpack.

    Returns
    -------
    data: dict
        Column names mapped to numpy arrays. Boolean columns are filled
        with False, string columns with a question mark, and numeric
        columns with NaN.
    """
    data = {}
    for name in table.colnames:
        column = table[name]
        if column.dtype.kind == 'b':
            data[name] = np.asarray(np.ma.filled(column, False)).astype(bool)
        elif column.dtype.kind in 'SU':
            data[name] = np.asarray(np.ma.filled(column, '?')).astype(str)
        else:
            data[name] = np.asarray(
                np.ma.filled(column.astype(float), np.nan))
    return data

Function: build_cuts

Return the set of flag cuts recommended for the Object table in the DP2 flag use examples, for one band, as a dictionary of keep masks.

A keep mask is a boolean array with one entry per row, True where the row passes that cut. The surviving sample is the logical AND of all the masks. The cuts are kept separate, rather than combined into one expression, so that each can be examined and removed independently.

In [6]:

def build_cuts(data, band):
    """
    Return the set of flag cuts recommended for the Object table in the
      DP2 flag use examples, for a single band.

    Parameters
    ----------
    data: dict
        Column names mapped to arrays, as returned by unpack_columns.
    band: str
        The band to build the cuts for, one of u, g, r, i, z, or y.

    Returns
    -------
    cuts: dict
        Cut labels mapped to boolean keep masks, which are True where an
        object passes that cut.
    """
    prefix = band + '_'
    return {
        prefix + 'inputCount > 0': data[prefix + 'inputCount'] > 0,
        prefix + 'psfFlux_flag = 0': ~data[prefix + 'psfFlux_flag'],
        prefix + 'invalidPsfFlag = 0': ~data[prefix + 'invalidPsfFlag'],
        prefix + 'pixelFlags_saturatedCenter = 0':
            ~data[prefix + 'pixelFlags_saturatedCenter'],
        prefix + 'pixelFlags_interpolatedCenter = 0':
            ~data[prefix + 'pixelFlags_interpolatedCenter'],
    }

Function: report_cut_impact

Print the fraction of a sample that each cut removes, and return the combined keep mask.

Two fractions are printed per cut. removed is the fraction that cut removes counted on its own. only this cut is the fraction it removes that every other cut in the set would have kept, so a value of zero means the cut is redundant with the rest of the set. The two differ because flag conditions overlap, which also means the removed column does not add up to the total fraction removed: an object rejected by three cuts is counted by all three.

In [7]:

def report_cut_impact(keep_masks, label='sample'):
    """
    Print the fraction of a sample that each cut removes.

    Parameters
    ----------
    keep_masks: dict
        Cut labels mapped to boolean keep masks, which are True where an
        object passes that cut.
    label: str
        A name for the sample, used in the summary line.

    Returns
    -------
    all_pass: numpy.ndarray
        Boolean array, True where an object passes every cut.
    """
    n_total = len(keep_masks[list(keep_masks)[0]])
    all_pass = np.logical_and.reduce(list(keep_masks.values()))

    print('{}: {} objects, {} kept, {:.2f}% removed'.format(
        label, n_total, np.sum(all_pass),
        100.0 * np.sum(~all_pass) / n_total))
    print('{:<38s} {:>10s} {:>15s}'.format('cut', 'removed', 'only this cut'))

    for name, mask in keep_masks.items():
        others = [m for other, m in keep_masks.items() if other != name]
        if others:
            pass_others = np.logical_and.reduce(others)
        else:
            pass_others = np.ones(n_total, dtype=bool)
        print('{:<38s} {:>10.4f} {:>15.4f}'.format(
            name, np.sum(~mask) / n_total,
            np.sum(~mask & pass_others) / n_total))

    return all_pass

2. Retrieve an unfiltered sample of objects¶

2.1. Define the sky region¶

Define a cone of radius 0.05 degrees on the Extended Chandra Deep Field South (ECDFS), one of the DP2 Deep Drilling Fields.

In [8]:

ecdfs_ra = 53.0
ecdfs_dec = -28.1
radius = 0.05

2.2. Query the Object table without any flag cuts¶

List the per-band column names without their band prefix in per_band, then prepend {band}_ to each one to build the query. Only the prefix changes when the band changes, so the list itself never has to be edited.

The query applies no flag cuts, so that the cuts can be built and switched off in python afterwards.

The list includes {band}_pixelFlags_edge, {band}_pixelFlags_sensor_edge, and {band}_pixelFlags_sensor_edgeCenter, which are not expected to be set at all, since the DP2 deep coadds are cell-based and so are not affected by the sensor edges of the input visit images. They are queried here in order to confirm that, in Section 3.5.

The list includes {band}_pixelFlags_nodata, which reports that no pixel data was available in this band for the object. An object appears in the Object table if it was detected in any band, so an object can have no data at all in a given band and still have a row.

In [9]:

per_band = ['psfMag', 'psfMagErr', 'cModelMag', 'inputCount',
            'i_flag', 'psfFlux_flag', 'cModel_flag', 'invalidPsfFlag',
            'pixelFlags_saturatedCenter', 'pixelFlags_saturated',
            'pixelFlags_interpolatedCenter', 'pixelFlags_interpolated',
            'pixelFlags_crCenter', 'pixelFlags_nodata',
            'pixelFlags_sensor_edge', 'pixelFlags_sensor_edgeCenter',
            'pixelFlags_bad', 'pixelFlags_edge', 'pixelFlags_suspect',
            'pixelFlags_suspectCenter', 'pixelFlags_offimage']

columns = ['coord_ra', 'coord_dec', 'refBand'] \
    + ['{}_{}'.format(band, name) for name in per_band]

query = "SELECT " + ", ".join(columns) + " " \
        "FROM dp2.Object " \
        "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " \
        "CIRCLE('ICRS', {}, {}, {})) = 1".format(ecdfs_ra, ecdfs_dec, radius)
print(query)

SELECT coord_ra, coord_dec, refBand, i_psfMag, i_psfMagErr, i_cModelMag, i_inputCount, i_i_flag, i_psfFlux_flag, i_cModel_flag, i_invalidPsfFlag, i_pixelFlags_saturatedCenter, i_pixelFlags_saturated, i_pixelFlags_interpolatedCenter, i_pixelFlags_interpolated, i_pixelFlags_crCenter, i_pixelFlags_nodata, i_pixelFlags_sensor_edge, i_pixelFlags_sensor_edgeCenter, i_pixelFlags_bad, i_pixelFlags_edge, i_pixelFlags_suspect, i_pixelFlags_suspectCenter, i_pixelFlags_offimage FROM dp2.Object WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), CIRCLE('ICRS', 53.0, -28.1, 0.05)) = 1

Run the query, unpack the columns, and report how many objects were returned.

In [10]:

objects = unpack_columns(run_query(query))
n_objects = len(objects['coord_ra'])
print(n_objects, 'objects retrieved, with no flag cuts applied')

Job phase is COMPLETED

2614 objects retrieved, with no flag cuts applied

2.3. Check the coadd coverage of the sample¶

Check the coverage before reading any flag.

refBand gives the band each object's forced measurements are referenced to.

In [11]:

input_count = objects['{}_inputCount'.format(band)]
psf_mag = objects['{}_psfMag'.format(band)]
nodata = objects['{}_pixelFlags_nodata'.format(band)]

print('{}_inputCount: min {:.0f}, median {:.0f}, max {:.0f}'.format(
    band, np.nanmin(input_count), np.nanmedian(input_count),
    np.nanmax(input_count)))
print('fraction with pixelFlags_nodata set: {:.4f}'.format(
    np.sum(nodata) / n_objects))
print('fraction with a finite psfMag:       {:.4f}'.format(
    np.sum(np.isfinite(psf_mag)) / n_objects))

print('refBand:', ', '.join(
    ['{} {}'.format(value, int(np.sum(objects['refBand'] == value)))
     for value in sorted(set(objects['refBand']))]))

i_inputCount: min 60, median 72, max 79
fraction with pixelFlags_nodata set: 0.0226
fraction with a finite psfMag:       0.9858
refBand: g 129, i 2157, r 54, u 73, y 101, z 100

A median input count in the tens, a nodata fraction of a few percent, and a finite magnitude for nearly every object mean the cone is well inside the coadd footprint and the flags below describe the objects.

A median input count of one or two, with a large nodata fraction, means the opposite. In that case reduce the radius or move the center before going further, because every occupancy below will be dominated by objects that have no data in this band.

2.4. Measure how often each flag is set¶

Measure the occupancy of each flag: the fraction of objects for which it is set. This separates the flags that will drive a selection from those that are set for almost nothing, and it is the most direct characterization of the flag columns in a given sample.

The general failure flags of the three measurements queried here, {band}_psfFlux_flag, {band}_cModel_flag and {band}_i_flag, appear in this list alongside the pixel flags. They are not interchangeable, because each reports on a different measurement. Require the one that matches the quantity used in the analysis.

In {band}_i_flag, the i in the middle is not a band: it is the symbol for the second moments of the object's light distribution, stored in the {band}_ixx, {band}_iyy, and {band}_ixy columns. The second moments measure how the light is spread about the object's center, and so describe its size and elongation. So {band}_i_flag reports a failure of that measurement, and it reads i_i_flag for the $i$ band and u_i_flag for the $u$ band.

In [12]:

flag_names = [name for name in columns if 'flag' in name.lower()]

print('{:<38s} {:>10s} {:>10s}'.format('flag column', 'n set', 'fraction'))
for name in flag_names:
    n_set = np.sum(objects[name])
    print('{:<38s} {:>10d} {:>10.5f}'.format(name, n_set,
                                             n_set / n_objects))

flag column                                 n set   fraction
i_i_flag                                      641    0.24522
i_psfFlux_flag                                247    0.09449
i_cModel_flag                                  91    0.03481
i_invalidPsfFlag                                0    0.00000
i_pixelFlags_saturatedCenter                   39    0.01492
i_pixelFlags_saturated                         90    0.03443
i_pixelFlags_interpolatedCenter                20    0.00765
i_pixelFlags_interpolated                     324    0.12395
i_pixelFlags_crCenter                          18    0.00689
i_pixelFlags_nodata                            59    0.02257
i_pixelFlags_sensor_edge                        0    0.00000
i_pixelFlags_sensor_edgeCenter                  0    0.00000
i_pixelFlags_bad                                0    0.00000
i_pixelFlags_edge                               0    0.00000
i_pixelFlags_suspect                            0    0.00000
i_pixelFlags_suspectCenter                      0    0.00000
i_pixelFlags_offimage                           0    0.00000

Plot the same fractions, sorted, on a logarithmic axis. The fractions span several orders of magnitude.

Flags that are set for no object are plotted at half of one object, so that they stay visible.

In [13]:

fractions = np.array([np.sum(objects[name]) / n_objects
                      for name in flag_names])
order = np.argsort(fractions)
floor = 0.5 / n_objects

fig, ax = plt.subplots(figsize=(6, 6))
ax.barh(np.arange(len(flag_names)),
        np.clip(fractions[order], floor, None))
ax.set_yticks(np.arange(len(flag_names)))
ax.set_yticklabels([flag_names[i] for i in order], fontsize=7)
ax.set_xscale('log')
ax.set_xlim(floor, 1.0)
ax.set_xlabel('fraction of objects with the flag set')
plt.tight_layout()
plt.show()

Figure 1: Horizontal bar chart of the fraction of objects for which each flag is set, one bar per flag, sorted with the most frequently set flag at the top. The horizontal axis is logarithmic and spans one object in the sample up to all objects. Bars at the left-hand limit are flags set for no object in this sample. The flags at the top of the chart dominate any selection built from these columns.

3. Characterize a set of flag cuts¶

3.1. Apply the baseline set of cuts¶

Build the set of cuts for an analysis that uses the PSF flux, using build_cuts from Section 1.2.3. These are the cuts recommended for the Object table in the DP2 flag use examples.

cut	what it protects against
{band}_inputCount > 0	the coadd has no input image at this position in this band
{band}_psfFlux_flag = 0	the PSF flux measurement failed
{band}_invalidPsfFlag = 0	there is no usable PSF model at the position
{band}_pixelFlags_saturatedCenter = 0	a saturated pixel sits in the central 3x3 box
{band}_pixelFlags_interpolatedCenter = 0	a central pixel was interpolated, not measured

Pass the set to report_cut_impact.

In [14]:

cuts = build_cuts(objects, band)
keep = report_cut_impact(cuts, label='Object, {} band'.format(band))

Object, i band: 2614 objects, 2334 kept, 10.71% removed
cut                                       removed   only this cut
i_inputCount > 0                           0.0000          0.0000
i_psfFlux_flag = 0                         0.0945          0.0853
i_invalidPsfFlag = 0                       0.0000          0.0000
i_pixelFlags_saturatedCenter = 0           0.0149          0.0054
i_pixelFlags_interpolatedCenter = 0        0.0077          0.0065

Read the only this cut column first. A cut whose value is zero removes nothing that the rest of the set would not remove anyway, so dropping it would not change the sample.

That is common among pixel flags, because the conditions are correlated. An object with a saturated central pixel usually has interpolated central pixels too, since saturated pixels are interpolated over.

Both columns depend on the field. A cut that is redundant in a Deep Drilling Field, where every position has many coadd inputs, can dominate in a region observed a few times (see Section 4).

3.2. Visualize which objects the cuts remove¶

Plot the kept and removed objects in the magnitude-uncertainty plane and on the sky.

In [15]:

mag_err = objects['{}_psfMagErr'.format(band)]
ra = objects['coord_ra']
dec = objects['coord_dec']

mag_range = (14.0, 30.0)
shown = np.isfinite(psf_mag) & np.isfinite(mag_err) \
    & (psf_mag > mag_range[0]) & (psf_mag < mag_range[1])
print('{} of {} objects fall outside the plotted magnitude range'.format(
    np.sum(~shown), n_objects))

75 of 2614 objects fall outside the plotted magnitude range

In [16]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4))

kept_shown = keep & shown
removed_shown = ~keep & shown
ax1.plot(psf_mag[kept_shown], mag_err[kept_shown], 'o', ms=3, mew=0,
         alpha=0.6, color='grey',
         label='kept ({})'.format(np.sum(kept_shown)))
ax1.plot(psf_mag[removed_shown], mag_err[removed_shown], 'x', ms=4, mew=1.0,
         alpha=0.8, label='removed ({})'.format(np.sum(removed_shown)))
ax1.set_xlabel('{}_psfMag (mag)'.format(band))
ax1.set_ylabel('{}_psfMagErr (mag)'.format(band))
ax1.set_yscale('log')
ax1.set_xlim(*mag_range)
ax1.legend(loc='upper left', fontsize=8)

ax2.plot(ra[keep], dec[keep], 'o', ms=2.5, mew=0, alpha=0.5,
         color='grey')
ax2.plot(ra[~keep], dec[~keep], 'x', ms=4, mew=1.0, alpha=0.8)
ax2.set_xlabel('Right Ascension (deg)')
ax2.set_ylabel('Declination (deg)')
ax2.invert_xaxis()

plt.tight_layout()
plt.show()

Figure 2: Two scatter plots. At left, PSF magnitude on the horizontal axis against its uncertainty on a logarithmic vertical axis, for the objects kept by the baseline set, drawn as small round points, and the objects it removes, drawn as crosses. The axes are clipped, because a few objects with near-zero or negative flux have magnitudes far off scale. At right, the same two populations in Right Ascension and Declination. Whether the removed objects concentrate at bright magnitudes, at faint magnitudes, or at particular positions depends on which cut dominates the set.

3.3. Compare the magnitude distribution of the objects each cut removes¶

Break the removals down by cut.

In [17]:

bins = np.arange(15.0, 28.5, 0.25)
finite = np.isfinite(psf_mag)
styles = ['-', '--', '-.', ':', (0, (3, 1, 1, 1))]

fig, ax = plt.subplots(figsize=(6, 4))
ax.hist(psf_mag[finite], bins=bins, color='lightgrey', label='all objects')
for index, (name, cut) in enumerate(cuts.items()):
    removed = psf_mag[~cut & finite]
    if len(removed) > 0:
        ax.hist(removed, bins=bins, histtype='step', lw=1.5,
                ls=styles[index % len(styles)],
                label='removed by {} ({})'.format(name, len(removed)))
ax.set_xlabel('{}_psfMag (mag)'.format(band))
ax.set_ylabel('number of objects')
ax.set_yscale('log')
ax.legend(fontsize=6)
plt.tight_layout()
plt.show()

Figure 3: Histogram of PSF magnitude on the horizontal axis against number of objects on a logarithmic vertical axis. The filled grey histogram is the full sample. Each outlined histogram, in a different line style, is the objects removed by one cut.

Note that the objects removed by the saturation cut are not necessarily the brightest objects in the sample. A faint object next to a bright one can have a saturated pixel from its neighbor among its central pixels. Additionally, for a saturated object the saturated pixel values are not propagated into the coadd, so the core is missing and the measured flux underestimates the real flux, placing the object at a fainter recorded magnitude.

3.4. Add stricter cuts¶

To switch a cut off, drop it from the set; to add one, add it. Define the stricter, science-case-dependent cuts from the flag use examples page and measure what they add.

cut	what it adds
{band}_pixelFlags_crCenter = 0	rejects a cosmic ray at the center
{band}_pixelFlags_interpolated = 0	rejects any interpolated pixel in the footprint
{band}_i_flag = 0	rejects a failed second-moments measurement

The second is the strict counterpart of the interpolatedCenter cut already in the baseline set: it rejects an object if any pixel anywhere in the footprint was interpolated. For an extended object with a large footprint that is far more aggressive.

The third shows what a flag for an unused measurement costs when added anyway.

In [18]:

stricter_cuts = {
    '{}_pixelFlags_crCenter = 0'.format(band):
        ~objects['{}_pixelFlags_crCenter'.format(band)],
    '{}_pixelFlags_interpolated = 0'.format(band):
        ~objects['{}_pixelFlags_interpolated'.format(band)],
    '{}_i_flag = 0'.format(band): ~objects['{}_i_flag'.format(band)],
}

keep_strict = report_cut_impact(
    {**cuts, **stricter_cuts},
    label='Object, {} band, baseline + stricter'.format(band))

Object, i band, baseline + stricter: 2614 objects, 1736 kept, 33.59% removed
cut                                       removed   only this cut
i_inputCount > 0                           0.0000          0.0000
i_psfFlux_flag = 0                         0.0945          0.0004
i_invalidPsfFlag = 0                       0.0000          0.0000
i_pixelFlags_saturatedCenter = 0           0.0149          0.0038
i_pixelFlags_interpolatedCenter = 0        0.0077          0.0000
i_pixelFlags_crCenter = 0                  0.0069          0.0000
i_pixelFlags_interpolated = 0              0.1239          0.0796
i_i_flag = 0                               0.2452          0.1389

Compare the magnitude distribution of the surviving sample under both selections, and with no cuts.

In [19]:

variants = {
    'no cuts': np.ones(n_objects, dtype=bool),
    'baseline': keep,
    'baseline + stricter': keep_strict,
}

styles = ['-', '--', '-.']

fig, ax = plt.subplots(figsize=(6, 4))
for index, (name, mask) in enumerate(variants.items()):
    ax.hist(psf_mag[mask & finite], bins=bins, histtype='step', lw=1.8,
            ls=styles[index % len(styles)],
            label='{} ({})'.format(name, np.sum(mask)))
ax.set_xlabel('{}_psfMag (mag)'.format(band))
ax.set_ylabel('number of objects')
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

Figure 4: Histogram of PSF magnitude on the horizontal axis against number of objects on the vertical axis, for three selections drawn in different line styles: no flag cuts, the baseline set, and the baseline set plus the stricter cuts. The vertical gap between two curves at a given magnitude is the completeness cost of the extra cuts there. A selection that removes a magnitude-dependent fraction changes the shape of the distribution, not only its normalization.

4. The same cuts in other fields¶

Repeat the measurement in three additional DP2 fields.

field	RA, Dec (deg)	galactic latitude (deg)	$i$-band visits	character
DDF ECDFS	53.0, -28.1	-54.6	96	deep extragalactic Deep Drilling Field
New Horizons	289.4, -20.2	-14.6	116	low galactic latitude, denser in stars
Rubin SV 280 -48	280.1, -48.0	-18.1	31	fewer coadd inputs

Coordinates and visit counts are from the observations page on dp2.lsst.io, which lists every DP2 field.

Each pair of fields isolates one variable. ECDFS and New Horizons have almost the same number of coadd inputs but differ in stellar density, so comparing them isolates crowding. New Horizons and Rubin SV 280 -48 sit at similar galactic latitude but differ by a factor of a few in coadd inputs, so comparing them isolates depth.

Query each field with the same column list and no flag cuts, and report the coverage of each alongside its object count, as in Section 2.3, so that a difference in flag occupancy can be attributed to crowding or to depth rather than to missing data.

In [20]:

fields = {
    'DDF ECDFS': (53.0, -28.1),
    'New Horizons': (289.4, -20.2),
    'Rubin SV 280 -48': (280.1, -48.0),
}

field_data = {}
for name, (field_ra, field_dec) in fields.items():
    field_query = "SELECT " + ", ".join(columns) + " " \
        "FROM dp2.Object " \
        "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " \
        "CIRCLE('ICRS', {}, {}, {})) = 1".format(
            field_ra, field_dec, radius)
    table = run_query(field_query)
    if len(table) == 0:
        print(name, 'returned no objects and is skipped')
        continue
    field_data[name] = unpack_columns(table)
    print('{:<20s} {:>7d} objects, median inputCount {:>4.0f}, '
          'nodata {:.4f}'.format(
              name, len(table),
              np.nanmedian(field_data[name]['{}_inputCount'.format(band)]),
              np.mean(field_data[name]['{}_pixelFlags_nodata'.format(band)])))

Job phase is COMPLETED

DDF ECDFS               2614 objects, median inputCount   72, nodata 0.0226

Job phase is COMPLETED

New Horizons            3808 objects, median inputCount   73, nodata 0.0481

Job phase is COMPLETED

Rubin SV 280 -48        1896 objects, median inputCount   26, nodata 0.0992

Compare the occupancy of the baseline flags across the fields, with two footprint-wide pixel flags added.

{band}_pixelFlags_saturated and {band}_pixelFlags_interpolated are good proxies for how many bright objects a field holds, since saturated pixels are interpolated over.

In [21]:

compare = ['inputCount = 0', 'psfFlux_flag', 'invalidPsfFlag', 'i_flag',
           'pixelFlags_saturatedCenter', 'pixelFlags_interpolatedCenter',
           'pixelFlags_saturated', 'pixelFlags_interpolated']

occupancy = {}
for name, data in field_data.items():
    values = {'inputCount = 0':
              data['{}_inputCount'.format(band)] <= 0}
    for flag in compare[1:]:
        values[flag] = data['{}_{}'.format(band, flag)]
    occupancy[name] = {key: np.mean(value) for key, value in values.items()}

header = '{:<38s}'.format('condition set (fraction of objects)')
for name in occupancy:
    header = header + '{:>18s}'.format(name)
print(header)
for flag in compare:
    row = '{:<38s}'.format(flag)
    for name in occupancy:
        row = row + '{:>18.5f}'.format(occupancy[name][flag])
    print(row)

condition set (fraction of objects)            DDF ECDFS      New Horizons  Rubin SV 280 -48
inputCount = 0                                   0.00000           0.00000           0.00000
psfFlux_flag                                     0.09449           0.04044           0.07015
invalidPsfFlag                                   0.00000           0.00000           0.00000
i_flag                                           0.24522           0.07589           0.09546
pixelFlags_saturatedCenter                       0.01492           0.02705           0.00897
pixelFlags_interpolatedCenter                    0.00765           0.00578           0.00897
pixelFlags_saturated                             0.03443           0.10163           0.04378
pixelFlags_interpolated                          0.12395           0.20851           0.15559

Plot the same numbers as grouped bars, with a different hatch per field so the fields are distinguishable without color.

In [22]:

hatches = ['', '///', '...', 'xxx']
positions = np.arange(len(compare))
width = 0.8 / max(len(occupancy), 1)
floor = 1.0e-3

fig, ax = plt.subplots(figsize=(7, 5))
for offset, (name, values) in enumerate(occupancy.items()):
    heights = np.array([values[flag] for flag in compare])
    ax.barh(positions + offset * width, np.clip(heights, floor, None),
            height=width, hatch=hatches[offset % len(hatches)], label=name)
ax.set_yticks(positions + 0.4 - width / 2)
ax.set_yticklabels(compare, fontsize=8)
ax.set_xscale('log')
ax.set_xlim(floor, 1.0)
ax.set_xlabel('fraction of objects with the condition set')
ax.invert_yaxis()
ax.legend(fontsize=8)
plt.tight_layout()
plt.show()

Figure 5: Grouped horizontal bar chart of the fraction of objects for which each condition is set, with one group of bars per condition and one bar per field, each field drawn with a different hatch pattern. The horizontal axis is logarithmic. Bars at the left-hand limit correspond to fractions below one in a thousand, which for samples of this size is a few objects at most, including none at all. The saturation and interpolation flags track how many bright objects a field holds, and are expected to be set more often in the field at lower galactic latitude. The measurement-failure flags respond to crowding and depth together, and inputCount = 0 responds to coadd coverage alone.

Apply the baseline set of cuts to each field and compare how much of each sample survives.

In [23]:

for name, data in field_data.items():
    print()
    report_cut_impact(build_cuts(data, band), label=name)

DDF ECDFS: 2614 objects, 2334 kept, 10.71% removed
cut                                       removed   only this cut
i_inputCount > 0                           0.0000          0.0000
i_psfFlux_flag = 0                         0.0945          0.0853
i_invalidPsfFlag = 0                       0.0000          0.0000
i_pixelFlags_saturatedCenter = 0           0.0149          0.0054
i_pixelFlags_interpolatedCenter = 0        0.0077          0.0065

New Horizons: 3808 objects, 3568 kept, 6.30% removed
cut                                       removed   only this cut
i_inputCount > 0                           0.0000          0.0000
i_psfFlux_flag = 0                         0.0404          0.0331
i_invalidPsfFlag = 0                       0.0000          0.0000
i_pixelFlags_saturatedCenter = 0           0.0270          0.0186
i_pixelFlags_interpolatedCenter = 0        0.0058          0.0026

Rubin SV 280 -48: 1896 objects, 1744 kept, 8.02% removed
cut                                       removed   only this cut
i_inputCount > 0                           0.0000          0.0000
i_psfFlux_flag = 0                         0.0701          0.0649
i_invalidPsfFlag = 0                       0.0000          0.0000
i_pixelFlags_saturatedCenter = 0           0.0090          0.0053
i_pixelFlags_interpolatedCenter = 0        0.0090          0.0026

5. Exercises for the learner¶

1. Change band in Section 1.2.2 from i to u and re-run the notebook. Every query, cut, and axis label follows. The $u$ band is the shallowest of the six, so expect a smaller median inputCount in Section 2.3, the measurement-failure flags to be set for a larger fraction of objects, and a smaller surviving sample in Section 3.1.

2. Replace psfFlux_flag with cModel_flag in build_cuts in Section 1.2.3, as suits an analysis using the CModel fluxes, and re-run Section 3.1. Compare the removed fraction of the two flags, using {band}_cModelMag in place of {band}_psfMag in the figures.

3. Add another field from the observations page on dp2.lsst.io to the fields dictionary in Section 4 and re-run that section. The M49 field, at RA, Dec of 186.3 deg, 6.9 deg, has more $i$-band visits than any field used here, so expect inputCount = 0 to be set for a smaller fraction of objects there than in Rubin SV 280 -48.

In [ ]: