301.0. DP1 overview

301.0. DP1 overview¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 1
Container Size: large
LSST Science Pipelines version: v29.1.1
Last verified to run: 2025-06-20
Repository: github.com/lsst/tutorial-notebooks

Learning objective: To gain an overview of the target fields of Data Preview 1 (DP1) data taken with LSSTComCam.

LSST data products:

Packages: lsst.daf.butler, lsst.rsp, lsst.geom

Credit: Originally developed by Andrés A. Plazas Malagón and the Rubin Community Science team. Sections 3, 4, and 5 are based on notebooks developed for the DP1 paper in github.com/lsst/rtn-095. 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¶

This notebook provides an overview and exploratory analysis of the seven target fields observed during the ComCam on-sky commissioning campaign of the Vera C. Rubin Observatory.

Between October and December 2024, ComCam was used to acquire ∼16,000 exposures for systems engineering, science performance validation, and early data processing trials (SITCOMTN-149). A subset of these exposures formed the foundation of Data Preview 1 (DP1), offering the science community an opportunity to gain experience with real LSST-like data and tools. Each field was selected to support specific validation goals and LSST science use cases, covering a broad range of sky locations, object densities, and overlap with external datasets.

Related tutorials: Other tutorials in the 301 series about ComCam fields during commissioning.

1.1. Field descriptions¶

47 Tuc
A crowded stellar field, 47 Tucanae was targeted for early testing of image quality and Active Optics System (AOS) alignment performance. The high stellar density of this field makes it ideal for testing crowded-field photometry pipelines.

Low Ecliptic Latitude Field (Rubin SV 38 7)
This low ecliptic latitude field was selected to optimize Solar System science. It was observed using a 2x2 dithered grid pattern to facilitate tracking of moving objects across a larger region. This field supports time-domain science and tests related to Solar System Object association.

Fornax Dwarf Spheroidal (dSph) Galaxy
Centered on the Fornax Dwarf Spheroidal Galaxy, this field is ideal for studying resolved stellar populations in a Local Group dwarf galaxy.

Extended Chandra Deep Field South (ECDFS)
The Extended Chandra Deep Field South (ECDFS) received the densest and most temporally consistent coverage among all fields, with hundreds of visits across all six bands. It was used as the main field for photometric calibration validation using FGCM (Forward Global Calibration Method), and includes the CalSpec star C26202.

Euclid Deep Field South (EDFS)
This field overlaps with the Euclid Deep Field South, and was selected to test deep photometry and weak lensing pipelines. Although its band coverage is more limited than ECDFS, it is highly relevant for extragalactic science and future cross-survey comparisons (DP1 paper draft, §2.1).

Low Galactic Latitude Field (Rubin SV 95 -25)
Located at low Galactic latitude, this field has a very high stellar density and is relevant for Galactic structure studies.

Seagull Nebula
Targeting a region near the Seagull Nebula, this field supports ISM and star formation studies. It has moderately high stellar density and was included to explore structured backgrounds and variations in image quality and calibration.

1.1.1. Stellar densities¶

Stellar densities vary significantly across the fields and were a major consideration during field selection. Fields such as 47 Tuc and Rubin SV 95 -25 contain high-density stellar regions, testing the limits of the photometric and astrometric pipelines under crowded conditions. In contrast, ECDFS and EDFS are extragalactic fields with sparse stellar populations, supporting deep coadditions and weak lensing science.

1.1.2. Temporal cadence¶

Temporal cadence also varies by field and was driven by science goals and operational feasibility. ECDFS had the most uniform and nightly sampling, supporting difference imaging and photometric calibration studies. Rubin SV 38 7 was revisited with a structured dither pattern across multiple nights to enable Solar System science. Other fields, such as Seagull and Fornax, received sparser but still valuable coverage across filters and nights.

1.1.3. Summary table¶

See also the figure in Section 3.3.

Field	Suggested Science Use	Bands	Stellar Density	Relative Cadence	External Data
47 Tuc	Crowded field photometry	griy	High	Sparse	GAIA
Rubin SV 38 7	Solar System objects	griz	Medium	Low	—
Fornax dSph	Resolved dwarf galaxies	gri	High	Very Sparse	—
ECDFS	Extragalactic, time domain	ugrizy	Low	Highest	HST, DECam, The Monster
EDFS	Extragalactic, weak lensing	ugrizy	Low	High	Euclid
Rubin SV 95 -25	Galactic plane science	ugrizy	Medium	Medium	—
Seagull	Star formation, ISM	ugrz	Medium	Sparse	—

1.2. 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 Table Access Protocol (TAP) service, 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), and astroplan for working with astronomical targets.

Use skyproj, a specialized sky projection package, to create consistent sky maps
(skyproj).

Use the plotting utilities in lsst.utils.plotting to access LSST-standard color schemes, symbols, and line styles for multi-band visualizations.

In [1]:

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

from lsst.rsp import get_tap_service
from lsst.daf.butler import Butler
from lsst.utils.plotting import (
    get_multiband_plot_colors,
    get_multiband_plot_symbols,
    get_multiband_plot_linestyles
)

1.3. Define parameters and functions¶

Instantiate the Butler with the appropiate DP1 repository, collection, and skymap.

In [2]:

butler = Butler('dp1', collections='LSSTComCam/DP1')
assert butler is not None

Instantiate the TAP service.

In [3]:

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

Define colors, symbols, and linestyles to represent the six LSST filters, $ugrizy$, as defined in RTN-045.

In [4]:

filter_names = ['u', 'g', 'r', 'i', 'z', 'y']
filter_colors = get_multiband_plot_colors()
filter_symbols = get_multiband_plot_symbols()
filter_linestyles = get_multiband_plot_linestyles()

Define a dictionary with each field's name and central coordinates, Right Ascension and Declination, as defined in SITCOMTN-149. The following are International Celestial Reference System (ICRS) coordinates, shared in units of decimal degrees.

In [5]:

field_centers = {
    "47 Tuc": [6.02, -72.08],
    "SV 38 7": [37.86, 6.98],
    "Fornax": [40.00, -34.45],
    "ECDFS": [53.13, -28.10],
    "EDFS": [59.10, -48.73],
    "SV 95 -25": [95.00, -25.00],
    "Seagull": [106.23, -10.51],
}
field_names = list(field_centers.keys())

2. Sky map¶

Use skyproj to visualize the fields on the celestial spehere.

In [6]:

fig, ax = plt.subplots(figsize=(8, 4))
sp = skyproj.McBrydeSkyproj(ax=ax, celestial=True, galactic=False, gridlines=True)
sp.ax.set_xlabel("Right Ascension", fontsize=14, fontweight="bold", labelpad=12)
sp.ax.set_ylabel("Declination", fontsize=15, fontweight="bold", labelpad=12)

for name, (ra, dec) in field_centers.items():
    sp.ax.plot(ra, dec, marker="o", color="red", markersize=11.5)
    x, y = sp.proj(ra, dec)
    sp.ax.annotate(name, (x, y), fontsize=10, fontweight="bold",
                   textcoords="offset points", xytext=(3.4, 3.4))

ecl_lon = np.linspace(0, 360, 1000)
ecl_lat = np.zeros_like(ecl_lon)
ecl_eq = SkyCoord(ecl_lon * u.deg, ecl_lat * u.deg,
                  frame="geocentrictrueecliptic").transform_to("icrs")
ecl_ra = ecl_eq.ra.deg
ecl_dec = ecl_eq.dec.deg

sp.ax.plot(ecl_ra, ecl_dec, color="gray", linestyle="--", lw=1.3, zorder=5)
sp.ax.text(*sp.proj(260, 50), "Ecliptic", fontsize=12, color="gray",
           ha="left", va="bottom")

l_vals = np.linspace(0, 360, 1000)
for b in [0, 10, -10]:
    gal = SkyCoord(l_vals * u.deg, b * u.deg, frame="galactic")
    eq = gal.transform_to("icrs")
    sp.ax.plot(eq.ra.deg, eq.dec.deg,
               color="black" if b == 0 else "gray",
               linestyle="-" if b == 0 else "--",
               lw=1.0 if b == 0 else 0.8, alpha=0.6)

sp.ax.gridlines.set_edgecolor("black")
sp.ax.gridlines.set_linestyle("--")
sp.ax.gridlines.set_alpha(0.3)
sp.ax.gridlines.zorder = 0
plt.tight_layout()
plt.show()

Figure 1: An all-sky view with the seven DP1 fields marked with red circles and labeled by name.

3. Visits¶

Use the TAP service to return the visit id, coordinates, band, and MJD of all visits.

In [7]:

query = """SELECT visit, ra, dec, band, expMidptMJD
        FROM dp1.Visit
        ORDER BY visit ASC"""
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'
visits = job.fetch_result().to_table()
del query, job

Job phase is COMPLETED

3.1. Filter distributions¶

Print the number of visits per filter, per field, and the total per field.

In [8]:

print('Field       visits per filter       total')
print('             u   g   r   i   z   y  ')
for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    offsets = np.sqrt((visits['ra']-ra)**2
                      + (visits['dec']-dec)**2)
    tx = np.where(offsets < 3.0)[0]
    tmp = np.zeros(6, dtype='int')
    for f, filt in enumerate(filter_names):
        tmp[f] = len(np.where(visits['band'][tx] == filt)[0])
    print('%-10s %3i %3i %3i %3i %3i %3i   %3i' %
          (name, tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], len(tx)))
    del name, ra, dec, offsets, tx, tmp

Field       visits per filter       total
             u   g   r   i   z   y
47 Tuc       0  10  32  19   0   5    66
SV 38 7      0  44  40  55  20   0   159
Fornax       0   5  25  12   0   0    42
ECDFS       43 230 237 162 153  30   855
EDFS        20  61  87  42  42  20   272
SV 95 -25   33  82  84  23  60  10   292
Seagull     10  37  43   0  10   0   100

3.2. Epochs (nights)¶

Print the number of nights the field was observed, and the average number of visits per night.

In [9]:

print('Field     # nights    mean visits/night')
for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    offsets = np.sqrt((visits['ra']-ra)**2
                      + (visits['dec']-dec)**2)
    tx = np.where(offsets < 3.0)[0]
    days = np.floor(visits['expMidptMJD'][tx])
    values, counts = np.unique(days, return_counts=True)
    tmp1 = len(values)
    tmp2 = np.mean(counts)
    print('%-10s  %3i       %2.1f' % (name, tmp1, tmp2))
    del name, ra, dec, offsets, tx
    del days, values, counts, tmp1, tmp2

Field     # nights    mean visits/night
47 Tuc        4       16.5
SV 38 7       5       31.8
Fornax        2       21.0
ECDFS        21       40.7
EDFS          9       30.2
SV 95 -25    10       29.2
Seagull       4       25.0

3.3. Temporal cadence¶

In [10]:

plt.figure(figsize=(8, 6))
spacing = 0.6

for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    offsets = np.sqrt((visits['ra']-ra)**2
                      + (visits['dec']-dec)**2)
    tx = np.where(offsets < 3.0)[0]
    for f, filt in enumerate(filter_names):
        fx = np.where(visits['band'][tx] == filt)[0]
        if len(fx) > 0:
            yval = (i * spacing) - (f * 0.12) + 0.3
            for j in range(len(fx)):
                plt.plot(visits['expMidptMJD'][tx[fx[j]]], yval,
                         filter_symbols[filt], alpha=0.5,
                         color=filter_colors[filt],
                         label=filt if (j == 0) & (i == 3) else None)

yticks = [i * spacing for i in range(len(field_names))]
plt.yticks(yticks, [t for t in field_names], fontsize=12)
plt.xlabel("Midpoint of Exposure (MJD)", fontsize=12)
plt.grid(axis='x', color='0.8', linewidth=1)
plt.legend(loc='upper left', ncol=3)
plt.ylim(-spacing / 1.5, (len(field_names) - 1) * spacing + spacing / 1.5)
plt.tight_layout()

4. Image quality¶

Use the TAP service to return the seeing and magnitude limits that were evaluated for all visit_images (all detectors) from the CcdVisit table.

In [11]:

query = """SELECT visitId, ra, dec, band, seeing, magLim
        FROM dp1.CcdVisit
        ORDER BY visitId ASC"""
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'
ccdvisits = job.fetch_result().to_table()
del query, job

Job phase is COMPLETED

4.1. Seeing¶

Calculate the average seeing over all detectors by filter, and for all filters.

In [12]:

print('Field      seeing per filter               mean')
print('              u    g    r    i    z    y  ')
for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    offsets = np.sqrt((ccdvisits['ra']-ra)**2
                      + (ccdvisits['dec']-dec)**2)
    tx = np.where(offsets < 3.0)[0]
    tmp = []
    for f, filt in enumerate(filter_names):
        fx = np.where(ccdvisits['band'][tx] == filt)[0]
        if len(fx) > 0:
            tmp.append(str(np.round(np.mean(ccdvisits['seeing'][tx[fx]]), 2)))
        else:
            tmp.append('  - ')
    tmp.append(str(np.round(np.mean(ccdvisits['seeing'][tx]), 2)))
    print('%-10s %-4s %-4s %-4s %-4s %-4s %-4s   %-4s' %
          (name, tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], tmp[6]))
    del name, ra, dec, offsets, tx

Field      seeing per filter               mean
              u    g    r    i    z    y
47 Tuc       -  1.47 1.26 1.11   -  1.34   1.25
SV 38 7      -  1.16 1.13 1.11 1.23   -    1.14
Fornax       -  1.16 0.83 0.87   -    -    0.88
ECDFS      1.43 1.17 1.09 1.01 1.0  1.07   1.1
EDFS       1.88 1.2  1.2  1.05 1.17 0.98   1.21
SV 95 -25  1.48 1.24 1.13 0.99 1.23 0.83   1.2
Seagull    1.53 1.26 1.12   -  1.21   -    1.22

4.2. Magnitude limits¶

Calculate the average magnitude limit over all detectors by filter, and for all filters.

In [13]:

print('Field      mag lim per filter                   mean')
print('              u     g     r     i     z     y  ')
for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    offsets = np.sqrt((ccdvisits['ra']-ra)**2
                      + (ccdvisits['dec']-dec)**2)
    tx = np.where(offsets < 3.0)[0]
    tmp = []
    for f, filt in enumerate(filter_names):
        fx = np.where(ccdvisits['band'][tx] == filt)[0]
        if len(fx) > 0:
            tmp.append(str(np.round(np.mean(ccdvisits['magLim'][tx[fx]]), 2)))
        else:
            tmp.append('  - ')
    tmp.append(str(np.round(np.mean(ccdvisits['magLim'][tx]), 2)))
    print('%-10s %-5s %-5s %-5s %-5s %-5s %-5s   %-5s' %
          (name, tmp[0], tmp[1], tmp[2], tmp[3], tmp[4], tmp[5], tmp[6]))
    del name, ra, dec, offsets, tx

Field      mag lim per filter                   mean
              u     g     r     i     z     y
47 Tuc       -   24.16 24.28 23.99   -   21.7    23.98
SV 38 7      -   24.45 24.19 23.74 22.9    -     23.94
Fornax       -   24.78 24.74 24.15   -     -     24.58
ECDFS      23.5  24.49 24.15 23.9  23.14 21.94   23.9
EDFS       23.3  24.57 24.15 23.87 23.18 22.15   23.84
SV 95 -25  23.54 24.63 24.36 23.95 23.12 22.17   23.98
Seagull    23.38 24.02 23.76   -   23.18   -     23.76

5. Coadded depth¶

Use the survey property maps to print the coadded depth in each band in the center of each field.

Load the PSF magnitude limits in each filter.

In [14]:

u_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='u')
g_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='g')
r_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='r')
i_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='i')
z_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='z')
y_maglim = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                      band='y')

Print the value per filter at the center of each field. Where there is no data, print "nan". Notice that all the depths for 47 Tuc are evaluated to be "nan", and the reason why is explored in Section 5.1.

In [15]:

print('Field      coadded mag lim per filter        ')
print('              u     g     r     i     z     y')
for i in range(len(field_names)):
    name = field_names[i]
    ra, dec = field_centers[name]
    lims = np.asarray([u_maglim.get_values_pos(ra, dec),
                       g_maglim.get_values_pos(ra, dec),
                       r_maglim.get_values_pos(ra, dec),
                       i_maglim.get_values_pos(ra, dec),
                       z_maglim.get_values_pos(ra, dec),
                       y_maglim.get_values_pos(ra, dec)], dtype='float')
    tx = np.where(lims < 0.0)[0]
    lims[tx] = float('NaN')
    print('%-10s %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f' %
          (name, lims[0], lims[1], lims[2], lims[3], lims[4], lims[5]))

Field      coadded mag lim per filter
              u     g     r     i     z     y
47 Tuc       nan   nan   nan   nan   nan   nan
SV 38 7      nan 25.29 25.06 24.72 23.49   nan
Fornax       nan   nan 25.08 25.31   nan   nan
ECDFS      25.18 27.36 27.09 26.50 25.61 23.57
EDFS       23.89 26.53 26.47 25.73 24.98 23.57
SV 95 -25  25.12 26.74 26.53 25.44 25.14 23.28
Seagull    24.32 25.83 25.57   nan 24.09   nan

Warning! The coadded depth in the center of 47 Tuc is not evaluated (not a number, "nan"). Proceed to the next section to understand why.

5.1. Depth maps¶

Plot the r-band depth maps for every field. Notice that 47 Tuc has a "hole" in the center, where there field was too crowded for the images to be processed for Data Preview 1. Software improvements and configuration tuning are underway to be able to handle such crowded fields.

In [16]:

for i in range(len(field_names)):
    name = field_names[i]
    ra_cen, dec_cen = field_centers[name]
    span_dec = 0.75
    if i == 1:
        span_dec = 1.00
    span_ra = span_dec / np.cos(np.deg2rad(dec_cen))
    ra = np.linspace(ra_cen-span_ra, ra_cen+span_ra, 250)
    dec = np.linspace(dec_cen-span_dec, dec_cen+span_dec, 250)
    x, y = np.meshgrid(ra, dec)
    values = r_maglim.get_values_pos(x, y)
    row, col = np.where(values < 0.0)
    values[row, col] = 'NaN'

    fig = plt.figure(figsize=(4, 3))
    plt.pcolormesh(x, y, values, vmin=24.0, vmax=27)
    plt.xlabel("Right Ascension (deg)")
    plt.ylabel("Declination (deg)")
    plt.title(name)
    plt.colorbar(label="PSF Mag Limit (r-band)")
    plt.gca().invert_xaxis()
    plt.show()

Figure 2: Seven panels showing the r-band depth survey property map for each field.

In [ ]: