301.1. DP2 overview#

301_1_DP2_overview

301.1. DP2 overview¶

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-24
Repository: github.com/lsst/tutorial-notebooks

Learning objective: Learn about the Data Preview 2 (DP2) dataset.

LSST data products: Visit and CoaddPatches tables, the skymap, and survey property maps

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

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¶

This notebook provides an overview of the observations that were obtained during commissioning and science validation with the LSST Camera (lsstcam.lsst.io), processed with v30 of the Rubin Science Pipelines (pipelines.lsst.io), and released as Data Preview 2 (DP2; dp2.lsst.io).

The code in this notebook was used to generated the DP2 observations overview webpage.

Key concepts for the DP2 dataset

Visit: One visit is one observation (one image; one exposure) with the LSSTCam, centered on a sky coordinate and obtained with a single filter and sky rotation.

Skymap, tract, and patch: One tract is one square region of LSST's all-sky tesselation, the "skymap". Tracts are $~1.66$ deg per side, identified by a unique four-digit number, and are subdivided into 100 overlapping patches. One deep coadd image is created per patch, provided sufficient input visits exist. The tract and patch numbers together with the filter uniquely identify a deep coadd image.

For DP2, the area covered by visits is greater than the area covered by deep coadd images. The sky area for which DP2 deep coadd images would be created was defined in late 2025 to include only the Wide-Fast-Deep, Small Field Survey, and Deep Drilling Field sky regions. See the Science Validation Survey Overview - Summary 20250930 webpage for more information about these regions. Additional visits outside of the deep coadd regions are included in DP2 as processed visit images only. Sources in this additional sky area appear only in the Source table (detections in the individual visit images), and not in the Object table (detections in the deep coadd images).

Related tutorials: Refer to the 100-level tutorials for how to use the RSP's butler and TAP services, and the 200-level tutorials for details on the visit table, deep coadd images, and survey property maps.

1.1. Import packages¶

Import numpy (numpy.org) for data array manipulation, matplotlib (matplotlib.org) for plotting. Import astropy (astropy.org) and skyproj (skyproj.readthedocs.io) for astronomy-specific functionality, and the healpy (healpy.readthedocs.io) and hpgeom (hpgeom.readthedocs.io) packages for dealing with HEALPix. Also import python's "garbage collector" gc in order to clear memory while running the notebook.

From the lsst package, import modules for data access: RSPDiscovery and the Butler. 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
import matplotlib.colors as mcolors
from astropy.table import Table, join
from astropy.time import Time
import skyproj
import healpy as hp
import hpgeom as hpg
import gc
from lsst.rsp import RSPDiscovery
from lsst.daf.butler import Butler
from lsst.utils.plotting import (
    get_multiband_plot_colors,
    get_multiband_plot_symbols,
    get_multiband_plot_linestyles
)

1.2. Define parameters and functions¶

Instantiate the Butler for the DP2 dataset.

In [2]:

butler = Butler('dp2', collections='dp2')

Instantiate the TAP service.

In [3]:

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

Define colors, symbols, and linestyles to represent the six LSST filters, $ugrizy$.

In [4]:

filter_colors = get_multiband_plot_colors()
filter_names = filter_colors.keys()
filter_symbols = get_multiband_plot_symbols()
filter_linestyles = get_multiband_plot_linestyles()

Uncomment and execute the following code cell to display plots interactively. Revert to non-interactive plots by executing %matplotlib inline, or by re-commenting the line, restarting the kernel and clearing all outputs, and re-executing the notebook.

In [5]:

# %matplotlib widget

Define dictionaries with the Data Preview 2 names and centers for Deep Drilling Fields (DDFs) and small survey fields (SFS). The coordinates are Right Ascension and Declination as International Celestial Reference System (ICRS) coordinates in decimal degrees.

A radius in degrees is defined with the coordinates: it is approximate, and encompasses only the boresight (field-of-view center) coordinates for all individual visits of that field. It does not encompass every detector, or every skymap patch, or every detected source or object associated with that field.

In [6]:

regions = {
    "DDF_ELAIS_S1": [9.5, -44.0, 1.0],
    "DDF_XMM_LSS": [35.6, -4.8, 1.0],
    "DDF_ECDFS": [53.0, -28.1, 1.0],
    "DDF_EDFS_a": [59.2, -49.2, 1.5],
    "DDF_EDFS_b": [63.2, -47.8, 1.5],
    "DDF_COSMOS": [150.1, 2.1, 1.0],
    "Abell_2764": [5.5, -49.0, 1.9],
    "DESI_SV3_R1": [180.5, -0.3, 0.8],
    "M49": [186.3, 6.9, 2.2],
    "Prawn": [253.5, -41.0, 2.7],
    "Trifid-Lagoon": [271.7, -23.9, 1.9],
    "New_Horizons": [289.4, -20.2, 1.5],
    "Rubin_SV_212_-7": [211.7, -7.0, 0.6],
    "Rubin_SV_216_-17": [216.1, -16.7, 0.4],
    "Rubin_SV_225_-40": [225.0, -39.5, 1.9],
    "Rubin_SV_280_-48": [280.1, -48.0, 2.0],
    "Rubin_SV_300_-41": [300.3, -41.0, 1.8],
    "Rubin_SV_320_-15": [320.2, -15.1, 4.5],
}
region_names = list(regions.keys())

Define a function that will remove a figure and all it's data, to help clear memory.

In [7]:

def remove_figure(fig):
    """
    Remove a figure to reduce memory footprint.

    Parameters
    ----------
    fig: matplotlib.figure.Figure
        Figure to be removed.

    Returns
    -------
    None
    """

    for ax in fig.get_axes():
        for im in ax.get_images():
            im.remove()
    fig.clf()
    plt.close(fig)
    gc.collect()

Define a function to truncate a colormap. This can be used to cut off, e.g., the darkest or lightest colors, to help with data visibility.

In [8]:

def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=256):
    """
    Customize a matplotlib colormap.

    Parameters
    ----------
    cmap: the matplotlib colormap
    minval: desired lower bound of colormap
    maxval: desired upper bound of colormap
    n: number of intervals (colormap resolution)

    Returns
    -------
    new_cmap: the truncated colormap
    """

    new_cmap = mcolors.LinearSegmentedColormap.from_list(
        'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
        cmap(np.linspace(minval, maxval, n)))
    return new_cmap

2. Visit metadata¶

Use the TAP service to return the visit id (visit), coordinates (ra and dec), filter (band), modified Julian date (MJD; expMidptMJD), and airmass of all visits included in DP2, from the Visit table.

In [9]:

query = """SELECT visit, ra, dec, band, expMidptMJD, airmass FROM dp2.Visit"""
job = tap_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_table = job.fetch_result().to_table()
job.delete()
del query

Job phase is COMPLETED

2.1. Sky map¶

Visualize the density of visits on an all-sky map.

Use the healpy package to print the area of one 19-sided HEALPix, to show that it is similar to the LSSTCam field of view (FOV) of 9.6 square degrees.

In [10]:

print('One 19-sided HEALPix is ', np.round(hp.nside2pixarea(19, degrees=True), 2),
      ' square degrees.')
print('It takes ', int(41253.0 / hp.nside2pixarea(19, degrees=True)),
      ' 19-sided HEALPix to cover the full sky.')

One 19-sided HEALPix is  9.52  square degrees.
It takes  4332  19-sided HEALPix to cover the full sky.

Use the skyproj package with a McBryde skyprojection to visualize the distribution of DP2 visits on the sky in HEALPix that are approximately the same size as the LSSTCam FOV (nside=19).

In [11]:

orig_cmap = plt.get_cmap('Blues')
new_cmap = truncate_colormap(orig_cmap, 0.10, 0.75)

fig, ax = plt.subplots(figsize=(12, 8))
sp = skyproj.McBrydeSkyproj(ax=ax)
vras = np.asarray(visit_table['ra'], dtype='float')
vdecs = np.asarray(visit_table['dec'], dtype='float')
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap=new_cmap)
sp.draw_colorbar(label='Number of visits (any filter)', shrink=0.5, pad=0.01)
for r, name in enumerate(region_names):
    color = 'black'
    symbol = 'o'
    if name[0:3] == 'DDF':
        color = 'darkviolet'
        symbol = 's'
    sp.ax.text(regions[name][0], regions[name][1]+4, str(r+1), color=color, fontweight='bold')
    sp.ax.plot(regions[name][0], regions[name][1], symbol, ms=10, color='None', mec=color,
               label=str(r+1)+': '+name)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.1), ncol=5, handletextpad=0)
plt.tight_layout()
plt.show()
del sp, vras, vdecs

No description has been provided for this image

Figure 1: A density map of the number of visits in blue, with Deep Drilling Fields marked with violet squares and small field survey areas marked with black circles. Field numbers are labeled in the legend above. Note that fields "Abell_2764" (7) and "DESI_SV3_R1" (8) have very few visits.

2.2. Filter distributions¶

Sum up the number of visits per filter.

Instantiate an astropy table to hold the values: field_filt_dist.

In [12]:

columns = ['Field', 'u', 'g', 'r', 'i', 'z', 'y', 'Total']
fields = np.asarray(region_names + ['All'], dtype='str')
temp = np.zeros(len(fields), dtype='int')
field_filt_dist = Table([fields, temp, temp, temp, temp, temp, temp, temp], names=columns)
del columns, fields, temp

Tally the number of visits per field, per filter, and over all filters and all visits.

In [13]:

for i, name in enumerate(region_names):
    ra, dec, rad = regions[name]
    offsets = np.sqrt((visit_table['ra']-ra)**2 + (visit_table['dec']-dec)**2)
    tx = np.where(offsets < rad)[0]
    for filt in filter_names:
        field_filt_dist[filt][i] = len(np.where(visit_table['band'][tx] == filt)[0])
    field_filt_dist['Total'][i] = len(tx)
    del ra, dec, rad, offsets, tx

i = np.where(field_filt_dist['Field'] == "All")[0]
for filt in filter_names:
    field_filt_dist[filt][i] = len(np.where(visit_table['band'] == filt)[0])
field_filt_dist['Total'][i] = len(visit_table)

Display the table.

In [14]:

field_filt_dist

Out[14]:

Table length=19

Field	u	g	r	i	z	y	Total
str16	int64	int64	int64	int64	int64	int64	int64
DDF_ELAIS_S1	57	153	128	232	141	59	770
DDF_XMM_LSS	42	15	8	47	10	10	132
DDF_ECDFS	10	61	50	96	61	11	289
DDF_EDFS_a	7	31	33	59	33	7	170
DDF_EDFS_b	5	37	33	57	33	7	172
DDF_COSMOS	105	87	173	162	110	77	714
Abell_2764	0	6	0	2	0	0	8
DESI_SV3_R1	0	0	0	3	0	0	3
M49	234	280	378	213	0	0	1105
Prawn	193	159	139	91	30	0	612
Trifid-Lagoon	231	201	129	119	9	6	695
New_Horizons	37	53	77	116	80	33	396
Rubin_SV_212_-7	0	139	187	123	0	0	449
Rubin_SV_216_-17	0	27	28	59	0	0	114
Rubin_SV_225_-40	307	524	391	381	239	104	1946
Rubin_SV_280_-48	30	29	26	30	29	0	144
Rubin_SV_300_-41	0	8	0	0	30	0	38
Rubin_SV_320_-15	12	47	104	259	211	77	710
All	1964	4166	4856	7959	5809	3944	28698

2.3. Epochs (nights)¶

Create a table with the number of visits per field, the number of unique nights (epochs) the field was observed, and the average number of visits per night.

In [15]:

columns = ['Field', 'visits', 'nights', 'mean(visits/night)']
fields = np.asarray(region_names + ['All'], dtype='str')
temp = np.zeros(len(fields), dtype='int')
field_epochs = Table([fields, temp, temp, temp], names=columns)
del columns, fields, temp

for i, name in enumerate(region_names):
    ra, dec, rad = regions[name]
    offsets = np.sqrt((visit_table['ra']-ra)**2 + (visit_table['dec']-dec)**2)
    tx = np.where(offsets < rad)[0]
    field_epochs['visits'][i] = len(tx)
    unique_nights = np.unique(np.floor(visit_table['expMidptMJD'][tx]))
    field_epochs['nights'][i] = len(unique_nights)
    field_epochs['mean(visits/night)'][i] = int(np.round(len(tx)/len(unique_nights), 0))
    del ra, dec, rad, offsets, tx, unique_nights

i = np.where(field_epochs['Field'] == "All")[0]
field_epochs['visits'][i] = len(visit_table)
unique_nights = np.unique(np.floor(visit_table['expMidptMJD']))
field_epochs['nights'][i] = len(unique_nights)
field_epochs['mean(visits/night)'][i] = int(np.round(len(visit_table)/len(unique_nights), 0))
del i, unique_nights

In [16]:

field_epochs

Out[16]:

Table length=19

Field	visits	nights	mean(visits/night)
str16	int64	int64	int64
DDF_ELAIS_S1	770	62	12
DDF_XMM_LSS	132	23	6
DDF_ECDFS	289	37	8
DDF_EDFS_a	170	51	3
DDF_EDFS_b	172	50	3
DDF_COSMOS	714	21	34
Abell_2764	8	2	4
DESI_SV3_R1	3	1	3
M49	1105	7	158
Prawn	612	7	87
Trifid-Lagoon	695	23	30
New_Horizons	396	24	16
Rubin_SV_212_-7	449	5	90
Rubin_SV_216_-17	114	2	57
Rubin_SV_225_-40	1946	28	70
Rubin_SV_280_-48	144	1	144
Rubin_SV_300_-41	38	2	19
Rubin_SV_320_-15	710	34	21
All	28698	131	219

2.4. Accumulation over time¶

Plot the cumulative number of visits over time, for all regions and filters.

In [17]:

notable_mjds = [np.floor(np.min(visit_table['expMidptMJD'])),
                60857, 60881, 60973,
                np.floor(np.max(visit_table['expMidptMJD']))]
fig = plt.figure(figsize=(6, 4))
for i, mjd in enumerate(notable_mjds):
    plt.axvline(mjd, ls='dotted', color='grey')
    t = Time(mjd, format='mjd', scale='utc')
    s = str(t.to_datetime())
    plt.text(mjd+2, 22000, s[0:10], rotation=90)
    del t, s
plt.plot(np.sort(visit_table['expMidptMJD']), np.arange(len(visit_table)),
         ls='solid', lw=1, color='black')
plt.xlabel('MJD')
plt.ylabel('Total number of visits')
plt.title('Cumulative distribution of visit MJDs')
plt.show()
del notable_mjds

Figure 2: The cumulative distribution of visit MJDs (modified Julian dates) for all visits, for all filters combined. Notable dates are marked with vertical dotted lines: the start (Apr 25 2025) and end (Jan 07 2026) of the data collection window; a phase of rapid data acquisition (July 01-25 2025); and the end of the engineering shutdown when observations recommenced (Oct 25 2025).

2.5. Airmass distribution¶

Plot the cumulative distribution of airmass, for all regions and filters.

In [18]:

fig = plt.figure(figsize=(6, 4))
xvals = np.sort(visit_table['airmass'])
yvals = np.arange(len(visit_table))
a_50 = xvals[np.where(yvals < len(visit_table)/2.0)[0][-1]]
plt.plot(xvals, yvals, ls='solid', lw=1, color='black')
plt.axvline(a_50, ls='dashed', color='lightgrey')
plt.text(a_50 + 0.02, 0, 'P50 = ' + str(np.round(a_50, 2)))
plt.xlabel('Airmass')
plt.ylabel('Total number of visits')
plt.title('Cumulative distribution of visit airmass (full)')
plt.show()
del xvals, yvals, a_50

Figure 3: The cumulative distribution of visit airmass for all visits, for all filters combined. The 50% percentile for airmass is 1.2, and is marked with a dashed grey line.

2.6. Image quality¶

The seeing and magnitude limit magLim are estimated per-detector, and are available in the VisitDetector table.

Use the TAP service to return the seeing and magnitude limits, averaged over all detectors per visit, from the VisitDetector table.

In [19]:

query = """SELECT visitId, AVG(seeing) AS mean_seeing, AVG(magLim) AS mean_maglim
           FROM dp2.VisitDetector GROUP BY visitId"""
job = tap_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_detector_table = job.fetch_result().to_table()
job.delete()
del query

Job phase is COMPLETED

Combine the mean_seeing and mean_maglim columns with the visit_table using an astropy table join, to create the joined_table. To save memory, delete the original tables and run the garbage collector.

In [20]:

joined_table = join(visit_table, visit_detector_table, keys_left="visit", keys_right="visitId")
del visit_table, visit_detector_table
gc.collect()

Out[20]:

Plot the distributions of seeing and magnitude limit.

Customize the line transparency and width for the image quality plots.

In [21]:

filter_alphas = {'u': 1, 'g': 1, 'r': 0.8, 'i': 0.6, 'z': 1, 'y': 1}
filter_lw = {'u': 1.1, 'g': 1, 'r': 1, 'i': 0.8, 'z': 1, 'y': 1.1}

Create the seeing distribution plots.

In [22]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
for filt in filter_names:
    fx = np.where(joined_table['band'] == filt)[0]
    xvals = np.sort(joined_table['mean_seeing'][fx])
    yvals = (np.arange(len(xvals), dtype='float') + 1) / len(xvals)
    ax[0].hist(joined_table['mean_seeing'][fx], 30, histtype='step', cumulative=True,
               alpha=filter_alphas[filt], lw=filter_lw[filt],
               ls=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    ax[1].plot(xvals, yvals, ls=filter_linestyles[filt], lw=filter_lw[filt],
               alpha=filter_alphas[filt], color=filter_colors[filt], label=filt)
    del fx, xvals, yvals
ax[0].set_xlim([0.6, 3.0])
ax[1].set_xlim([0.6, 3.0])
ax[0].set_xlabel('Mean seeing (PSF FWHM, arcsec)')
ax[1].set_xlabel('Mean seeing (PSF FWHM, arcsec)')
ax[0].set_ylabel('Number of visits')
ax[1].set_ylabel('Fraction of visits')
ax[1].legend(loc='lower right')
plt.show()

Figure 4: The cumulative number (left) and fraction (right) of visits as a function of mean seeing (PSF FWHM in arcseconds, averaged over all detectors), for a given filter. These plots include all regions.

Create the magnitude limits plots.

In [23]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4))
for filt in filter_names:
    fx = np.where(joined_table['band'] == filt)[0]
    xvals = np.sort(joined_table['mean_maglim'][fx])
    yvals = (np.arange(len(xvals), dtype='float') + 1) / len(xvals)
    ax[0].hist(joined_table['mean_maglim'][fx], 30, histtype='step', cumulative=True,
               alpha=filter_alphas[filt], lw=filter_lw[filt],
               ls=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    ax[1].plot(xvals, yvals, ls=filter_linestyles[filt], lw=filter_lw[filt],
               alpha=filter_alphas[filt], color=filter_colors[filt], label=filt)
    del fx, xvals, yvals
ax[0].set_xlim([20, 26])
ax[1].set_xlim([20, 26])
ax[0].set_xlabel('Mean magnitude limit')
ax[1].set_xlabel('Mean magnitude limit')
ax[0].set_ylabel('Number of visits')
ax[1].set_ylabel('Fraction of visits')
ax[1].legend(loc='lower right')
plt.show()

Figure 5: The cumulative number (left) and fraction (right) of visits as a function of mean magnitude limit (averaged over all detectors), for a given filter. These plots include all regions.

Create tables of the mean seeing and magnitude limit, averaged over all visits in a given filter, per region.

Because seeing and magnitude limit are filter-dependent, but an average over all filters is not calculated.

In [24]:

columns = ['Field', 'u', 'g', 'r', 'i', 'z', 'y']
fields = np.asarray(region_names + ['All'], dtype='str')
temp = np.zeros(len(fields), dtype='float')
field_seeing = Table([fields, temp, temp, temp, temp, temp, temp], names=columns)
field_maglim = Table([fields, temp, temp, temp, temp, temp, temp], names=columns)
del columns, fields, temp

for i, name in enumerate(region_names):
    ra, dec, rad = regions[name]
    offsets = np.sqrt((joined_table['ra']-ra)**2 + (joined_table['dec']-dec)**2)
    tx = np.where(offsets < rad)[0]
    for filt in filter_names:
        fx = np.where(joined_table['band'][tx] == filt)[0]
        if len(fx) > 0:
            field_seeing[filt][i] = np.round(np.mean(joined_table['mean_seeing'][tx[fx]]), 2)
            field_maglim[filt][i] = np.round(np.mean(joined_table['mean_maglim'][tx[fx]]), 2)
        else:
            field_seeing[filt][i] = np.nan
            field_maglim[filt][i] = np.nan
        del fx
    del ra, dec, rad, offsets, tx

i = np.where(field_seeing['Field'] == "All")[0]
for filt in filter_names:
    fx = np.where(joined_table['band'] == filt)[0]
    field_seeing[filt][i] = np.round(np.mean(joined_table['mean_seeing'][fx]), 2)
    field_maglim[filt][i] = np.round(np.mean(joined_table['mean_maglim'][fx]), 2)
    del fx

In [25]:

field_seeing

Out[25]:

Table length=19

Field	u	g	r	i	z	y
str16	float64	float64	float64	float64	float64	float64
DDF_ELAIS_S1	1.63	1.51	1.51	1.32	1.46	1.33
DDF_XMM_LSS	1.52	1.24	1.08	1.01	0.99	0.94
DDF_ECDFS	1.68	1.59	1.35	1.15	1.34	1.32
DDF_EDFS_a	1.56	1.39	1.39	1.12	1.29	0.94
DDF_EDFS_b	1.34	1.41	1.39	1.16	1.3	0.91
DDF_COSMOS	1.36	1.16	1.25	1.14	1.17	1.22
Abell_2764	nan	1.01	nan	0.85	nan	nan
DESI_SV3_R1	nan	nan	nan	1.47	nan	nan
M49	1.49	1.37	1.27	1.22	nan	nan
Prawn	1.49	1.38	1.28	1.24	1.42	nan
Trifid-Lagoon	1.22	1.15	1.15	1.03	1.1	1.55
New_Horizons	1.26	1.36	1.12	1.04	1.01	1.12
Rubin_SV_212_-7	nan	1.42	1.3	1.12	nan	nan
Rubin_SV_216_-17	nan	1.35	1.25	1.32	nan	nan
Rubin_SV_225_-40	1.5	1.42	1.27	1.24	1.31	1.25
Rubin_SV_280_-48	1.69	1.54	1.51	1.34	1.53	nan
Rubin_SV_300_-41	nan	1.86	nan	nan	1.65	nan
Rubin_SV_320_-15	1.26	1.68	1.47	1.36	1.37	1.15
All	1.36	1.3	1.27	1.18	1.21	1.19

In [26]:

field_maglim

Out[26]:

Table length=19

Field	u	g	r	i	z	y
str16	float64	float64	float64	float64	float64	float64
DDF_ELAIS_S1	23.23	23.64	23.27	23.2	22.47	21.48
DDF_XMM_LSS	23.21	23.7	23.93	23.62	23.03	21.92
DDF_ECDFS	23.01	23.52	23.17	23.42	22.57	21.68
DDF_EDFS_a	23.1	23.62	23.23	23.47	22.51	21.94
DDF_EDFS_b	23.26	23.68	23.25	23.42	22.48	21.96
DDF_COSMOS	23.39	24.39	23.54	22.98	22.58	21.57
Abell_2764	nan	24.87	nan	23.41	nan	nan
DESI_SV3_R1	nan	nan	nan	23.5	nan	nan
M49	23.29	24.3	23.85	23.44	nan	nan
Prawn	23.1	24.07	23.74	23.38	22.71	nan
Trifid-Lagoon	23.42	24.28	23.86	23.65	22.63	21.05
New_Horizons	23.17	24.41	24.11	23.63	22.94	21.88
Rubin_SV_212_-7	nan	24.23	23.83	23.71	nan	nan
Rubin_SV_216_-17	nan	24.34	24.03	23.62	nan	nan
Rubin_SV_225_-40	23.3	23.97	23.52	23.33	22.76	21.63
Rubin_SV_280_-48	22.59	23.45	23.28	23.18	22.27	nan
Rubin_SV_300_-41	nan	23.25	nan	nan	22.31	nan
Rubin_SV_320_-15	23.55	24.08	23.29	23.05	22.39	21.67
All	23.37	24.12	23.58	23.32	22.65	21.68

3. Deep coadd metadata¶

Metadata for the deep coadd images comes from the skymap and the survey property maps.

3.1. Skymap (tracts and patches)¶

The skymap of DP2 tracts and patches was defined before all the observations were obtained.

Just because a tract and patch exists does not mean it is populated with a deep coadd image. This will be demonstrated in Section 3.2.

Retrieve the DP2 patch and tract numbers (lsst_patch, lsst_tract), the center coordinates RA, Dec (s_ra, s_dec), and the polygon region defining the patch (s_region) from the CoaddPatches table.

In [27]:

query = """SELECT lsst_patch, lsst_tract, s_dec, s_ra, s_region FROM dp2.CoaddPatches"""
job = tap_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'
patches_table = job.fetch_result().to_table()
job.delete()
del query

Job phase is COMPLETED

Option to display the table.

In [28]:

# patches_table

Get the unique tract number for DP2, and print their total number.

In [29]:

tracts = np.unique(patches_table['lsst_tract'])
print('Number of tracts: ', len(tracts))

Number of tracts:  2191

Retrieve the DP2 skymap, named lsst_cells_v2, from the Butler.

In [30]:

skymap = butler.get("skyMap", skymap="lsst_cells_v2")

3.1.1. All-sky map of tracts¶

Draw one box per tract on an all-sky map using the same projection as in Figure 1.

In [31]:

fig, ax = plt.subplots(figsize=(10, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
for tract in tracts:
    info = skymap.generateTract(tract)
    vertex_list = info.getVertexList()
    ras = []
    decs = []
    for vertice in vertex_list:
        ras.append(vertice.getRa().asDegrees())
        decs.append(vertice.getDec().asDegrees())
    sp.draw_polygon(ras, decs, edgecolor='darkgrey', alpha=1, linewidth=0.5, facecolor=None)
    del info, vertex_list, ras, decs
for r, name in enumerate(region_names):
    color = 'black'
    symbol = 'o'
    if name[0:3] == 'DDF':
        color = 'darkviolet'
        symbol = 's'
    sp.ax.text(regions[name][0], regions[name][1]+5, str(r+1), color=color, fontweight="bold")
    sp.ax.plot(regions[name][0], regions[name][1], symbol, ms=10, color='None', mec=color,
               label=str(r+1)+': '+name)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.05), ncol=5, handletextpad=0)
plt.show()
del sp

Figure 6: All pre-defined DP2 tracts. Note that fields DDF XMM LSS (2) and DESI SV3 R1 (8) do not overlap with any tracts and will not have any deep coadd images.

3.1.2. Single-field map of patches¶

Draw one box per patch for the DDF ECDFS field (2 in Figure 6).

Recall from Section 1.2. that the third component of regions[<name>] is the radius within which the boresight center coordinates for all visits of the field are included. To capture all tracts that contain patches for which these visits might have contributed to a deep coadd image, add to that radius the tract's diagonal size $\sqrt{2\times1.66^2}$ degrees.

To better visualize tract and patch overlap, define a color dictionary for the three central tracts using the first colors from the seaborn-v0_8-notebook color map (reference).

In [32]:

name = "DDF_ECDFS"
ra, dec, rad = regions[name]
radius = rad + np.sqrt(2*1.66**2)
delta_ra = (ra - patches_table['s_ra'])*(np.cos(np.deg2rad(dec)))
delta_dec = dec - patches_table['s_dec']
offset = np.sqrt(delta_ra**2 + delta_dec**2)
tx = np.where(offset < radius)[0]
tract_color_dict = {4848: '#1f77b4', 4849: '#ff7f0e', 5063: '#2ca02c'}

fig, ax = plt.subplots(1, 2, figsize=(10, 5))
for x in tx:
    s_region = patches_table['s_region'][x]
    temp = s_region.split(' ')
    ras = np.asarray([temp[2], temp[4], temp[6], temp[8], temp[2]], dtype='float')
    decs = np.asarray([temp[3], temp[5], temp[7], temp[9], temp[3]], dtype='float')
    for a in [0, 1]:
        color = 'grey'
        if patches_table['lsst_tract'][x] in tract_color_dict:
            color = tract_color_dict[patches_table['lsst_tract'][x]]
        ax[a].plot(ras, decs, color=color, linewidth=0.5, alpha=1)
        ax[a].set_xlabel("Right Ascension", fontsize=12)
        ax[a].set_ylabel("Declination", fontsize=12)
    del s_region, temp, ras, decs
ax[1].set_xlim([ra-0.5, ra+0.5])
ax[1].set_ylim([dec-0.5, dec+0.5])
plt.tight_layout()
plt.show()
del name, ra, dec, rad, radius, tx
del delta_ra, delta_dec, offset

Figure 7: At left, all pre-defined patches for the EDCFS. At right, the zoom-in better shows that both tracts and individual patches within a tract overlap at their edges. The center three tracts are assigned colors just to help with visualizing the overlap.

Option to print the unique tract numbers for a given field.

In [33]:

# name = "DDF_ECDFS"
# ra, dec, rad = regions[name]
# radius = rad + np.sqrt(2*1.66**2)
# tx = np.where(np.sqrt(((ra - patches_table['s_ra'])*(np.cos(np.deg2rad(dec))))**2 +
#                       (dec - patches_table['s_dec'])**2) < radius)[0]
# tract_list = np.unique(patches_table['lsst_tract'][tx])
# print(name, tract_list)
# del name, ra, dec, rad, radius, tx, tract_list

3.2. Image quality (survey property maps)¶

Use the survey property maps to explore the PSF (point-spread function) limiting magnitude and size (FWHM; full-width half-max) for the deep coadd images.

The two maps used in this demonstration are:

PSF magnitude limit, or "depth": deepCoadd_psf_maglim_consolidated_map_weighted_mean
PSF size (FWHM) in pixels (0.2"/pix): deepCoadd_psf_size_consolidated_map_weighted_mean

Maps are available per filter.

3.2.1. All-sky map of depth (i-band)¶

Load the all-sky PSF magnitude limit map for the i-band. Degrade the resolution of the data returned, from nside_sparse = 32768 (the default) to 1024, so as to take up less space in memory.

In [34]:

map_name = "deepCoadd_psf_maglim_consolidated_map_weighted_mean"
map_band = "i"
hspmap = butler.get(map_name, band=map_band, skymap="lsst_cells_v2",
                    parameters={'degrade_nside': 1024})

Print the resolution of the coverage map (always low-resolution) and of the retrieved survey property map (will match the degrade_nside parameter passed in the code cell above).

In [35]:

default_nside_coverage = hspmap.nside_coverage
print('Coverage map resolution: ', hspmap.nside_coverage)
print('Retrieved data resolution:', hspmap.nside_sparse)

Coverage map resolution:  32
Retrieved data resolution: 1024

Visualize the all-sky PSF magnitude limit map.

In [36]:

fig, ax = plt.subplots(figsize=(14, 8))
sp = skyproj.McBrydeSkyproj(ax=ax, lon_0=0.0)
sp.draw_hspmap(hspmap)
for r, name in enumerate(region_names):
    color = 'black'
    symbol = 'o'
    if name[0:3] == 'DDF':
        color = 'darkviolet'
        symbol = 's'
    sp.ax.text(regions[name][0], regions[name][1]+4, str(r+1), color=color, fontweight="bold")
    sp.ax.plot(regions[name][0], regions[name][1], symbol, ms=10, color='None', mec=color,
               label=str(r+1)+' '+name)
temp = map_name.split("_")
label = "PSF " + temp[2] + " (" + map_band + "-band)"
sp.draw_colorbar(label=label, shrink=0.5, pad=0.01)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.legend(loc="lower center", bbox_to_anchor=(0.5, 1.1), ncol=5, handletextpad=0)
plt.show()
del sp, temp, label

Figure 8: The all-sky map of the PSF magnitude limit in the i-band. Note that fields Abell 2764 (7) and Rubin SV 300 -41 (17), which do overlap with pre-defined tracts (Fig 6), have shallow and non-existent deep coadds, respectively.

Clear space in memory.

In [37]:

remove_figure(fig)
del hspmap
gc.collect()

Out[37]:

3.2.2. Single-field map of PSF size (g-band)¶

Instead of magnitude limit, load the PSF size map for the $g$-band, and only for HEALPix pixels that overlap DDF ECDFS.

Note: the PSF size survey property map stores $\sigma$ in pixels, not FWHM in arcseconds. To convert, multiply by the pixel scale, 0.2 arcsec/pixel, and by 2.35 to convert from $\sigma$ to FWHM.

In [38]:

map_name = "deepCoadd_psf_size_consolidated_map_weighted_mean"
map_band = "g"
name = "DDF_ECDFS"
coords = regions[name]
radius = coords[2] + np.sqrt(2*1.66**2)
delta_ra = (coords[0] - patches_table['s_ra'])*(np.cos(np.deg2rad(coords[1])))
delta_dec = coords[1] - patches_table['s_dec']
offset = np.sqrt(delta_ra**2 + delta_dec**2)
tx = np.where(offset < radius)[0]
pixels = np.unique(np.array(hpg.angle_to_pixel(default_nside_coverage,
                                               patches_table['s_ra'][tx],
                                               patches_table['s_dec'][tx])))
hspmap = butler.get(map_name, band=map_band, skymap='lsst_cells_v2',
                    parameters={'pixels': list(pixels), 'degrade_nside': 1024})

Visualize the PSF FWHM for the DDF ECDFS field.

In [39]:

fig, ax = plt.subplots(figsize=(6, 6))
sp = skyproj.Skyproj(ax=ax)
sp.draw_hspmap(hspmap)
temp = map_name.split("_")
label = "PSF " + temp[2] + " (" + map_band + "-band)"
sp.draw_colorbar(label=label, shrink=0.6, pad=0.05)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.grid(False)
del temp, label
for x in tx:
    s_region = patches_table['s_region'][x]
    temp = s_region.split(' ')
    ras = np.asarray([temp[2], temp[4], temp[6], temp[8], temp[2]], dtype='float')
    decs = np.asarray([temp[3], temp[5], temp[7], temp[9], temp[3]], dtype='float')
    sp.ax.plot(ras, decs, color='black', alpha=0.5, linewidth=0.5)
    del s_region, temp, ras, decs
plt.show()
del sp

Figure 9: The map of the PSF size, $\sigma$, in pixels (0.2"/pix) for the $g$-band, for HEALPix of the survey property map that overlap the DDF ECDFS. Patch boxes are drawn in black. The FWHM in arcseconds is $\sigma$ $\times$0.2"/pixel $\times$2.35 ($=0.47\sigma$).

In [40]:

remove_figure(fig)
del name, coords, radius, tx, pixels
del hspmap
gc.collect()

Out[40]:

3.2.3. Image quality at field centers¶

Option to create tables of the image quality values, PSF magnitude limit and size, at the centers of the individual fields (regions).

These are commented-out and are optional because they take >10 minutes to execute. They are included because they were used to generated the DP2 observations overview webpage.

In [41]:

columns = ['Field', 'u', 'g', 'r', 'i', 'z', 'y']
fields = np.asarray(region_names, dtype='str')
temp = np.zeros(len(fields), dtype='float')
field_coadd_maglim = Table([fields, temp, temp, temp, temp, temp, temp], names=columns)
field_coadd_size = Table([fields, temp, temp, temp, temp, temp, temp], names=columns)
del columns, fields, temp

In [42]:

map_name = "deepCoadd_psf_maglim_consolidated_map_weighted_mean"
print("Reading psf maglim map for filter: ")
for filt in filter_names:
    print(filt)
    hspmap = butler.get(map_name, band=filt, skymap='lsst_cells_v2',
                        parameters={'degrade_nside':512})
    for i, name in enumerate(region_names):
        ra, dec, rad = regions[name]
        field_coadd_maglim[filt][i] = np.round(hspmap.get_values_pos(ra, dec), 2)
        if field_coadd_maglim[filt][i] < 0:
            field_coadd_maglim[filt][i] = np.nan
        del ra, dec, rad
    del hspmap

Reading psf maglim map for filter:
u
g
r
i
z
y

In [43]:

field_coadd_maglim

Out[43]:

Table length=18

Field	u	g	r	i	z	y
str16	float64	float64	float64	float64	float64	float64
DDF_ELAIS_S1	25.18	26.57	25.88	26.07	25.07	23.64
DDF_XMM_LSS	nan	nan	nan	nan	nan	nan
DDF_ECDFS	24.35	25.84	25.47	25.87	24.76	22.23
DDF_EDFS_a	nan	25.63	25.53	25.65	24.53	23.08
DDF_EDFS_b	23.69	25.75	25.47	25.71	24.58	23.27
DDF_COSMOS	25.77	26.73	26.35	25.94	25.0	23.69
Abell_2764	nan	25.88	nan	22.58	nan	nan
DESI_SV3_R1	nan	nan	nan	nan	nan	nan
M49	25.62	27.03	26.72	26.05	nan	nan
Prawn	25.66	26.52	26.06	25.41	24.12	nan
Trifid-Lagoon	26.22	26.95	26.13	25.29	23.41	22.18
New_Horizons	25.26	26.37	26.37	25.99	25.25	23.65
Rubin_SV_212_-7	nan	26.81	26.56	26.09	nan	nan
Rubin_SV_216_-17	nan	26.15	25.93	25.83	nan	nan
Rubin_SV_225_-40	26.26	27.26	26.75	26.36	25.52	24.06
Rubin_SV_280_-48	23.92	25.09	24.89	24.88	23.7	nan
Rubin_SV_300_-41	nan	24.41	nan	nan	23.85	nan
Rubin_SV_320_-15	24.09	24.82	24.42	24.99	24.21	22.42

In [44]:

map_name = "deepCoadd_psf_size_consolidated_map_weighted_mean"
print("Reading psf size map for filter: ")
for filt in filter_names:
    print(filt)
    hspmap = butler.get(map_name, band=filt, skymap='lsst_cells_v2',
                        parameters={'degrade_nside':1024})
    for i, name in enumerate(region_names):
        ra, dec, rad = regions[name]
        field_coadd_size[filt][i] = np.round(hspmap.get_values_pos(ra, dec), 2)
        if field_coadd_size[filt][i] < 0:
            field_coadd_size[filt][i] = np.nan
        del ra, dec, rad
    del hspmap

Reading psf size map for filter:
u
g
r
i
z
y

In [45]:

field_coadd_size

Out[45]:

Table length=18

Field	u	g	r	i	z	y
str16	float64	float64	float64	float64	float64	float64
DDF_ELAIS_S1	2.69	2.55	2.55	2.53	2.57	2.6
DDF_XMM_LSS	nan	nan	nan	nan	nan	nan
DDF_ECDFS	2.53	2.5	2.2	2.25	2.14	2.89
DDF_EDFS_a	nan	2.47	2.31	2.14	2.36	2.04
DDF_EDFS_b	2.93	2.6	2.46	2.15	2.25	2.0
DDF_COSMOS	2.76	2.31	2.49	2.18	2.37	2.38
Abell_2764	nan	2.13	nan	3.11	nan	nan
DESI_SV3_R1	nan	nan	nan	nan	nan	nan
M49	2.94	2.68	2.54	2.47	nan	nan
Prawn	2.86	2.82	2.54	2.62	2.84	nan
Trifid-Lagoon	2.5	2.36	2.46	2.16	2.14	2.28
New_Horizons	2.98	2.77	2.3	2.27	2.2	2.23
Rubin_SV_212_-7	nan	2.68	2.62	2.37	nan	nan
Rubin_SV_216_-17	nan	2.28	2.13	2.49	nan	nan
Rubin_SV_225_-40	2.82	2.87	2.66	2.59	2.6	2.6
Rubin_SV_280_-48	3.27	3.18	3.13	2.81	3.1	nan
Rubin_SV_300_-41	nan	2.39	nan	nan	3.19	nan
Rubin_SV_320_-15	2.29	3.0	2.53	2.69	2.63	2.45