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 1
Container Size: large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2025-09-02
Repository: github.com/lsst/tutorial-notebooks

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

LSST data products: the_monster_20250219.

Packages: numpy, matplotlib, scipy, astropy, lsst.geom, lsst.daf.butler, lsst.rsp, lsst.rsp.get_tap_service, lsst.rsp.retrieve_query.

Credit: Originally developed by Andrés A. Plazas Malagón and 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 early data from the Vera C. Rubin Observatory, a new 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 structure supports scalable and fast photometric and astrometric calibrations across the sky during Rubin operations.

References:

• DMTN-277: The Monster: A reference catalog with synthetic ugrizy-band fluxes for the Vera C. Rubin observatory
• See notebooks folder in the GitHub repository associated with DMTN-277, github.com/lsst-dm/dmtn-277, for further examples.

Related tutorials: See other DP1 200-level calibration tutorials.

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; matplotlib gallery).

From the lsst package, import modules for accessing the butler and geometry utilities from the LSST Science Pipelines (pipelines.lsst.io).

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

From lsst.meas.algorithms get ReferenceObjectLoader to retrieve the Monster and account for proper motions.

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
import lsst.daf.butler as dafButler
from lsst.meas.algorithms import ReferenceObjectLoader
from lsst.rsp import get_tap_service

1.2. Define parameters and functions¶

Instantiate the butler.

In [2]:

butler = dafButler.Butler("dp1", collections="LSSTComCam/DP1")
assert butler is not None

Create an instance of the Table Access Protocol (TAP) service, and assert that it exists.

In [3]:

rsp_tap = get_tap_service("tap")
assert rsp_tap is not None

2. Access¶

Define a particular location on the sky. Use the center of the Extended Chandra Deep Field South, in degrees.

In [4]:

ra_cen, dec_cen = 53.13, -28.10

Get the tract correspoding to the coordinates.

In [5]:

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

my_tract:  5063

Get the name of "The Monster" catalog in the registry.

In [6]:

butler.registry.queryDatasetTypes('*monster*')

Out[6]:

[DatasetType('the_monster_20250219', {htm7}, SimpleCatalog)]

Get the reference datasets for the tract above.

In [7]:

monsterRefs = butler.registry.queryDatasets('the_monster_20250219',
                                            tract=my_tract,
                                            skymap='lsst_cells_v1').expanded()

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

In [8]:

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 will enable proper motions to be taken into account relative to Gaia’s reference epoch, J2016.0.

In [9]:

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.

Get a visit ID to define an epoch for the Monster.

In [10]:

my_band = 'i'

In [11]:

visitimage_refs = butler.query_datasets('visit_image',
                                        where="visit.region OVERLAPS POINT(:ra, :dec) "
                                        f"and BAND='{my_band}'",
                                        bind={'ra': ra_cen, 'dec': dec_cen},
                                        order_by=['visit', 'detector'],
                                        with_dimension_records=True,
                                        limit=1)

In [12]:

day_obs = visitimage_refs[0].dataId['day_obs']

In [13]:

print(day_obs)

20241108

Define an epoch and a radius around the central point.

In [14]:

epoch = Time({"year": 2024, "month": 11, "day": 8}, format='ymdhms')
radius = 1

Get the Monster.

In [15]:

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 [51.99635699, 54.26364301], [-29.10000172, -27.09999828] RA Dec

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

lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader INFO: Correcting reference catalog for proper motion to <Time object: scale='utc' format='ymdhms' value=(2024, 11, 8, 0, 0, 0.0)>

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

Filter out NaNs.

In [16]:

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.

The systems and their corresponding bands included in The Monster are:

• DES: g, r, i, z, y
• SDSS: u
• LATISS: g, r, i, z, y
• ComCam: u, g, r, i, z, y

In addition to the photometric data, The Monster includes astrometric and photometric measurements from Gaia DR3 for all entries. These include:

• Positions (RA, Dec) and their uncertainties
• Proper motions and their uncertainties
• Parallax and associated errors
• G, BP, and RP bands fluxes and errors
• Astrometric flags and covariances.

In [17]:

monster.colnames

Out[17]:

['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 Gaia¶

Extract the Monster catalog stars' coordinates as arrays coord_ra and coord_dec.

In [18]:

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

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

In [19]:

coord_ra

Out[19]:

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

0.9129291372703638
0.9125873066972817
0.9130009849032074
0.9132343447765262
0.9130617730487425
0.9128618848547745
0.9130670634959845
0.9130523239529361
0.9123769347596139
0.9127178279136074
0.9128095281275611
0.9128855635869073
0.9133358279061989
0.9134832603666909
0.9136056123275252
0.913702280287212
0.9138683431286071
0.913908303389697
0.913968547923707
0.914117430217921
0.9141147753975971
...
0.944179976668437
0.9442100238303095
0.9444706683536135
0.9446435659185383
0.9446534416231338
0.9439860255121967
0.9439073478023187
0.9440784148000932
0.9443543528516638
0.9452037055391516
0.9451000137772472
0.9445070062338412
0.9444262856659509
0.9445188527883461
0.9448679310243236
0.943891011252815
0.9440695436928136
0.9438101259215892
0.9436970021325495
0.9443351724866541
0.9442073156295125

In [20]:

len(coord_ra), len(coord_dec)

Out[20]:

(7463, 7463)

Get thesource catalog for the visit above. Retrieve only the sky coordinates and the flag to identify objects used in Point Spread Function modeling.

In [21]:

visit_id = visitimage_refs[0].dataId['visit']

In [22]:

query = (
    "SELECT coord_ra, coord_dec "
    "FROM dp1.Source "
    f"WHERE visit = {visit_id} "
    "AND calib_psf_used = 1 "
    "ORDER BY coord_ra "
)

In [23]:

rsp_tap = get_tap_service("tap")
job = rsp_tap.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'
psf_used_table = job.fetch_result()

Job phase is COMPLETED

Crossmatch Gaia stars from The Monster with stars used for PSF modeling. First, convert RA and DEC columns to SkyCoord objects.

In [24]:

source_coords = SkyCoord(ra=psf_used_table['coord_ra'] * u.deg,
                         dec=psf_used_table['coord_dec'] * u.deg)

monster_coords = SkyCoord(ra=coord_ra.to(u.deg),
                          dec=coord_dec.to(u.deg))

Match each object to the closest Monster source.

In [25]:

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

Filter matches within 1 arcsec.

In [26]:

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 561 out of 813 sources within 1.0 arcsec.

Visualize the matches.

In [27]:

plt.figure(figsize=(6, 6))
plt.scatter(psf_used_table['coord_ra'], psf_used_table['coord_dec'],
            s=12, color='blue', alpha=1.0,
            label='Object 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(psf_used_table['coord_ra'][matched],
         psf_used_table['coord_dec'][matched],
         'o', ms=10, mec='black', mew=0.5, color='None',
         label='matched objects')
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.ylim([-28, -27.65])
plt.xlim([53, 53.5])
plt.title('Source Matching: Object Catalog vs Monster Catalog')
plt.legend(loc='upper right', framealpha=1, handletextpad=0)
plt.gca().invert_xaxis()
plt.grid(True)
plt.tight_layout()
plt.show()

Figure 1: RA vs. Dec for objects that were used to estimate the PSF (i.e., bright, isolated stars; blue circles) and objects from the Monster reference catalog (orange crosses). Matches within 1" are circled.

3.1. Astrometric residuals¶

Calculate and visualize residuals.

Compute residuals.

In [28]:

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 [29]:

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

In [30]:

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)
cb1 = 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)
cb2 = 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) across this particular field.