204.2. The Monster reference catalog

204.2. The Monster reference catalog¶

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

Learning objective: To learn about "The Monster" reference catalog.

LSST data products: the_monster, Visit table, Source table.

Packages: numpy, matplotlib, astropy, lsst.geom, lsst.daf.butler, lsst.meas.algorithms, lsst.rsp.

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¶

To support the photometric and astrometric calibration of data from the Vera C. Rubin Observatory, an all-sky reference catalog called The Monster has been developed. This catalog provides synthetic fluxes in the LSST ugrizy bandpasses, constructed by combining and transforming data from a prioritized list of external reference catalogs.

The Monster is designed to meet the demanding calibration requirements of the LSST, ensuring high source density ($>$10 reference stars per NSIDE=256 HEALPix pixel) and full-sky coverage for all bands. It includes astrometry based on Gaia DR3, and photometry transformed into LSST-native bandpasses, minimizing the need for additional color transformations during image processing.

The catalog is organized using Hierarchical Triangular Mesh (HTM) level-7 trixels, a spatial indexing system that divides the celestial sphere into approximately 1.5 million small triangular regions, each covering about 0.05 square degrees—roughly the footprint of a single LSSTCam CCD. This hierarchical scheme enables efficient data access: rather than querying the entire sky, users can retrieve only the relevant subset of the catalog corresponding to their region of interest.

This notebook accesses The Monster from the Data Preview 2 (DP2) dp2 Butler repository at data-int.lsst.cloud, which contains data from the LSST Science Camera (LSSTCam). The workflow starts from a target sky position: the Butler's visit dimension records identify which visits overlap that position, and then two catalog tables are queried via the Table Access Protocol (TAP) service — the Visit table, to obtain the exposure epoch used to correct The Monster positions for proper motion, and the Source table, to obtain the observed stars for a cross match. Neither the Object table nor the visit images are needed.

References:

• DMTN-277: The Monster: A reference catalog with synthetic ugrizy-band fluxes for the Vera C. Rubin observatory
• The DP2 schema browser: sdm-schemas.lsst.io/lsstcam.html

Related tutorials: The DP2 tutorial on astrometric calibration uses The Monster as the astrometric reference, and also demonstrates the new proper-motion and parallax catalogs.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org).

From astropy, import modules for celestial coordinate transformations, units, and time handling (astropy.org).

From the lsst package, import the Butler and geometry utilities, the ReferenceObjectLoader (which retrieves The Monster while accounting for proper motions), and the RSPDiscovery class for accessing the RSP TAP service. (pipelines.lsst.io).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import astropy.units as u
from astropy.time import Time

import lsst.geom
from lsst.daf.butler import Butler
from lsst.rsp import RSPDiscovery
from lsst.meas.algorithms import ReferenceObjectLoader
from lsst.rsp import get_tap_service

1.2. Define parameters and functions¶

Instantiate the Butler with the DP2 repository and collection.

In [2]:

butler = Butler("dp2", collections=["dp2"])

Instantiate RSPDiscovery with the DP2 release, create an instance of the TAP service, and assert that it exists

In [3]:

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

2. Access¶

Define a target sky position (in degrees) and a band. Use a location near the center of the Extended Chandra Deep Field South (ECDFS).

In [4]:

ra_cen, dec_cen = 53.076, -28.110
my_band = 'i'

Find the visits whose footprint overlaps the target position, using the Butler's visit dimension records. This queries only visit metadata in the registry (no images or catalogs are loaded) and is the recommended way to identify visits by sky position. Take one overlapping visit in the chosen band.

In [5]:

visit_records = butler.query_dimension_records(
    "visit",
    where="band.name = :band AND visit.region OVERLAPS POINT(:ra, :dec)",
    bind={"band": my_band, "ra": ra_cen, "dec": dec_cen})
visit_ids = sorted(rec.id for rec in visit_records)
print(f"Number of {my_band}-band visits overlapping the target: {len(visit_ids)}")
visit_id = visit_ids[20]
print("Using visit:", visit_id)

Number of i-band visits overlapping the target: 100
Using visit: 2025072100533

Query the Visit table with the TAP service to get this visit's exposure midpoint. The expMidpt column is the exposure midpoint as a timestamp (TAI), and expMidptMJD is the same quantity as a Modified Julian Date, which is convenient for building an astropy Time. See the DP2 schema.

In [6]:

query = (
    "SELECT visit, expMidpt, expMidptMJD "
    "FROM dp2.Visit "
    f"WHERE visit = {visit_id}"
)
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'
visit_info = job.fetch_result().to_table()

Job phase is COMPLETED

Build the epoch from the exposure midpoint.

In [7]:

epoch = Time(float(visit_info['expMidptMJD'][0]), format='mjd', scale='tai')
print("Exposure midpoint:", visit_info['expMidpt'][0], "| epoch (MJD):", epoch.mjd)

Exposure midpoint: 2025-07-22T07:23:21.681 | epoch (MJD): 60878.30788983057

Get the tract corresponding to the field center, using the skymap.

In [8]:

my_spherePoint = lsst.geom.SpherePoint(ra_cen*lsst.geom.degrees,
                                       dec_cen*lsst.geom.degrees)
skymap = butler.get('skyMap', skymap='lsst_cells_v2')
tract = skymap.findTract(my_spherePoint)
my_tract = tract.tract_id
print('my_tract: ', my_tract)

my_tract:  5063

Find the name of The Monster catalog in the registry. The catalog is versioned by date, so the exact dataset type name may differ between data releases.

In [9]:

monster_dataset_types = sorted(
    dt.name for dt in butler.registry.queryDatasetTypes("*monster*")
    if dt.name.startswith("the_monster"))
monster_dataset_types

Out[9]:

['the_monster_20250219']

Use the most recent version of the catalog.

In [10]:

monster_name = monster_dataset_types[-1]
print(monster_name)

the_monster_20250219

Get the reference datasets for the tract above.

In [11]:

monsterRefs = butler.registry.queryDatasets(monster_name,
                                            tract=my_tract,
                                            skymap='lsst_cells_v2').expanded()

Get the data IDs containing all the htm7 IDs and the reference catalogs for each htm7 in this tract.

In [12]:

data_ids = [ref.dataId for ref in monsterRefs]
refCats = [butler.getDeferred(ref) for ref in monsterRefs]

Define a ReferenceObjectLoader object for The Monster catalog. This enables proper motions to be taken into account relative to Gaia's reference epoch, J2016.0.

In [13]:

monsterLoader = ReferenceObjectLoader(dataIds=data_ids, refCats=refCats)

Warning: It is also possible to retrieve The Monster catalog via butler.get, but it is recommended to use ReferenceObjectLoader instead to take into account proper motions.

Define a radius around the central point, and get The Monster, corrected to the visit's exposure epoch.

In [14]:

radius = 0.5
monster = monsterLoader.loadSkyCircle(my_spherePoint,
                                      radius*lsst.geom.degrees,
                                      'phot_g_mean',
                                      epoch=epoch).refCat.asAstropy()

lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader INFO: Loading reference objects from None in region bounded by [52.50913041, 53.64286959], [-28.61000172, -27.60999828] RA Dec

lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader INFO: Loaded 2812 reference objects

lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader INFO: Correcting reference catalog for proper motion to <Time object: scale='tai' format='mjd' value=60878.30788983057>

Note: the argument phot_g_mean is needed for the call to succeed, but is not relevant in this case.

Filter out NaNs.

In [15]:

monster = monster[np.isfinite(monster['coord_ra'])]

2.1. Schema¶

Use the catalog above to check the schema.

As noted in DMTN-277, for each photometric system and bandpass included in The Monster, there are three associated columns:

• monster_{system}_{band}_flux: estimated flux in the given band of the specified system.
• monster_{system}_{band}_fluxErr: uncertainty on the flux estimate.
• monster_{system}_{band}_source_flag: flag denoting which catalog was the source of the flux measurement.

In addition to the photometric data, The Monster includes astrometric and photometric measurements from Gaia DR3 for all entries: positions and their uncertainties, proper motions, parallaxes, G, BP, and RP fluxes, and astrometric flags and covariances.

In [16]:

monster.colnames

Out[16]:

['id',
 'coord_ra',
 'coord_dec',
 'phot_g_mean_flux',
 'phot_bp_mean_flux',
 'phot_rp_mean_flux',
 'phot_g_mean_fluxErr',
 'phot_bp_mean_fluxErr',
 'phot_rp_mean_fluxErr',
 'coord_raErr',
 'coord_decErr',
 'epoch',
 'pm_ra',
 'pm_dec',
 'pm_raErr',
 'pm_decErr',
 'pm_flag',
 'parallax',
 'parallaxErr',
 'parallax_flag',
 'coord_ra_coord_dec_Cov',
 'coord_ra_pm_ra_Cov',
 'coord_ra_pm_dec_Cov',
 'coord_ra_parallax_Cov',
 'coord_dec_pm_ra_Cov',
 'coord_dec_pm_dec_Cov',
 'coord_dec_parallax_Cov',
 'pm_ra_pm_dec_Cov',
 'pm_ra_parallax_Cov',
 'pm_dec_parallax_Cov',
 'astrometric_excess_noise',
 'monster_ComCam_g_flux',
 'monster_ComCam_g_fluxErr',
 'monster_ComCam_g_source_flag',
 'monster_ComCam_i_flux',
 'monster_ComCam_i_fluxErr',
 'monster_ComCam_i_source_flag',
 'monster_ComCam_r_flux',
 'monster_ComCam_r_fluxErr',
 'monster_ComCam_r_source_flag',
 'monster_ComCam_y_flux',
 'monster_ComCam_y_fluxErr',
 'monster_ComCam_y_source_flag',
 'monster_ComCam_z_flux',
 'monster_ComCam_z_fluxErr',
 'monster_ComCam_z_source_flag',
 'monster_DES_g_flux',
 'monster_DES_g_fluxErr',
 'monster_DES_g_source_flag',
 'monster_DES_i_flux',
 'monster_DES_i_fluxErr',
 'monster_DES_i_source_flag',
 'monster_DES_r_flux',
 'monster_DES_r_fluxErr',
 'monster_DES_r_source_flag',
 'monster_DES_y_flux',
 'monster_DES_y_fluxErr',
 'monster_DES_y_source_flag',
 'monster_DES_z_flux',
 'monster_DES_z_fluxErr',
 'monster_DES_z_source_flag',
 'monster_LATISS_g_flux',
 'monster_LATISS_g_fluxErr',
 'monster_LATISS_g_source_flag',
 'monster_LATISS_i_flux',
 'monster_LATISS_i_fluxErr',
 'monster_LATISS_i_source_flag',
 'monster_LATISS_r_flux',
 'monster_LATISS_r_fluxErr',
 'monster_LATISS_r_source_flag',
 'monster_LATISS_y_flux',
 'monster_LATISS_y_fluxErr',
 'monster_LATISS_y_source_flag',
 'monster_LATISS_z_flux',
 'monster_LATISS_z_fluxErr',
 'monster_LATISS_z_source_flag',
 'monster_ComCam_u_flux',
 'monster_ComCam_u_fluxErr',
 'monster_ComCam_u_source_flag',
 'monster_SDSS_u_flux',
 'monster_SDSS_u_fluxErr',
 'monster_SDSS_u_source_flag']

3. Cross match objects with The Monster¶

Extract The Monster stars' coordinates as arrays coord_ra and coord_dec.

In [17]:

coord_ra = monster['coord_ra']
coord_dec = monster['coord_dec']

Check one of the columns to verify that the coordinates are in radians.

In [18]:

coord_ra

Out[18]:

<Column name='coord_ra' dtype='float64' unit='rad' description='position in ra/dec' length=2121>

0.9172279796470385
0.9172857206573967
0.917241358157227
0.9173895268502462
0.9176169625797136
0.9172450725212975
0.9178700348322678
0.91777750572389
0.917806818245459
0.917429726764242
0.9176133822472147
0.9173647278506093
0.9174534335428495
0.9175710789646457
0.9176253412974905
0.9175708928458335
0.9180833755090329
0.9167052150014452
0.9170157621921692
0.9167874300739216
0.9173239225455009
0.91730707827123
0.9178298071981559
0.9170945497381964
0.9169522608335602
0.916899718387485
0.9173048952632232
...
0.9329637610140844
0.933579361919568
0.9319276522625346
0.9316923559555184
0.9315894398093628
0.93162899680871
0.9317768201814454
0.9326515102349993
0.9326789691436769
0.9329732779801092
0.9322004145250246
0.9321306914150354
0.9327271292123137
0.9345179627332703
0.9340013892819302
0.9341197574321265
0.9343717431462082
0.9342950542962819
0.9335529848286463
0.9340510305366722
0.9308985531768568
0.9309765564942306
0.9317373774440136
0.9314178592066111
0.9314242987554479
0.9305509596957661
0.9318240809893591

Retrieve the observed stars used for Point Spread Function (PSF) modeling in this visit, from the Source table via the TAP service. Select only the stars flagged with calib_psf_used, within the field radius.

In [19]:

query = (
    "SELECT ra, dec "
    "FROM dp2.Source "
    f"WHERE visit = {visit_id} "
    "AND calib_psf_used = 1 "
    "AND CONTAINS(POINT('ICRS', ra, dec), "
    f"CIRCLE('ICRS', {ra_cen}, {dec_cen}, {radius})) = 1"
)
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'
source_table = job.fetch_result().to_table()

Job phase is COMPLETED

In [20]:

source_ra = np.asarray(source_table['ra'])
source_dec = np.asarray(source_table['dec'])
print(len(source_ra))

1146

Crossmatch Gaia stars from The Monster with the stars used for PSF modeling. First, build SkyCoord objects (note that The Monster coordinates are in radians, and the Source coordinates are in degrees).

In [21]:

source_coords = SkyCoord(ra=source_ra * u.deg, dec=source_dec * u.deg)
monster_coords = SkyCoord(ra=coord_ra.to(u.deg), dec=coord_dec.to(u.deg))

Match each source to the closest Monster star.

In [22]:

idx, d2d, _ = source_coords.match_to_catalog_sky(monster_coords)

Filter matches within 1 arcsec.

In [23]:

max_sep = 1.0
matched = d2d.arcsecond < max_sep
print(f"Matched {matched.sum()} out of {len(source_coords)} sources within {max_sep} arcsec.")

Matched 964 out of 1146 sources within 1.0 arcsec.

Visualize the matches.

In [24]:

plt.figure(figsize=(6, 6))
plt.scatter(source_ra, source_dec, s=12, color='blue', alpha=1.0,
            label='Source catalog, used for PSF')
plt.scatter(coord_ra.to(u.deg), coord_dec.to(u.deg), s=22, color='orange',
            marker='+', alpha=1.0, label='Monster catalog')
plt.plot(source_ra[matched], source_dec[matched], 'o', ms=10, mec='black', mew=0.5,
         color='None', label='matched objects')
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.xlim([ra_cen + radius, ra_cen - radius])
plt.ylim([dec_cen - radius, dec_cen + radius])
plt.title('Source matching: Source catalog vs Monster catalog')
plt.legend(loc='upper right', framealpha=1, handletextpad=0)
plt.grid(True)
plt.tight_layout()
plt.show()

Figure 1: RA vs. Dec for the stars used to estimate the PSF in this visit (from the Source table; blue circles) and stars from The Monster reference catalog (orange crosses). Matches within 1 arcsec are circled.

3.1. Astrometric residuals¶

Calculate and visualize the residuals between the matched Source and Monster positions.

Compute residuals.

In [25]:

residuals_ra = source_coords[matched].ra - monster_coords[idx][matched].ra
residuals_dec = source_coords[matched].dec - monster_coords[idx][matched].dec

Extract sky positions.

In [26]:

ra = source_coords[matched].ra.deg
dec = source_coords[matched].dec.deg

In [27]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6),
                               sharex=True, sharey=True)

sc1 = ax1.scatter(ra, dec, c=residuals_ra.to('mas').value, s=20)
fig.colorbar(sc1, ax=ax1, label='RA Residual (mas)')
ax1.set_title("RA Residuals")
ax1.set_xlabel("RA (deg)")
ax1.set_ylabel("Dec (deg)")
ax1.invert_xaxis()

sc2 = ax2.scatter(ra, dec, c=residuals_dec.to('mas').value, s=20)
fig.colorbar(sc2, ax=ax2, label='Dec Residual (mas)')
ax2.set_title("Dec Residuals")
ax2.set_xlabel("RA (deg)")
ax2.invert_xaxis()

plt.tight_layout()
plt.show()

Figure 2: Astrometric residuals in RA (left) and Dec (right), between the matched Source positions and The Monster reference positions.

In [ ]: