102. LSSTCam visits metadata (2026)

102. LSSTCam visits metadata 2026¶

For the Rubin Science Platform at data.lsst.cloud.
Container Size: Large
LSST Science Pipelines version: r29.2.0
Last verified to run: 2026-05-01
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Explore and visualize metadata for visits obtained Dec 2025 - Apr 2026.

LSST data products: A temporary, preliminary list of visits metadata.

Packages: skyproj

Credit: Originally developed by the Rubin Community Science team. Please cite the DOI above and 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¶

The purpose of this tutorial, and the associated metadata file, is to enable early exploration and visualization of some of the first science visits obtained between Dec 2025 and Apr 2026.

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.

Caveats: This tutorial uses a temporary, preliminary metadata file for exposures that were labeled with an "observation reason" of "science" by the scheduler. The processed visit images and related catalogs might become available when Prompt images and the Prompt Products Database (PPDB) are released to users. However, the data validation stage of the data release preparations might imposed a different start date, and/or rejected some visits as not appropriate for scientific analysis. Furthermore, Prompt-processed images are only planned to be persisted for 30 days, and already that window has passed for many visits (but catalogs will persist). As a final reminder, difference image analysis (DIA) is only possible in the limited regions with templates.

When the Prompt products are released, this static, temporary, preliminary file will be deprecated, and the visit metadata will be available in queryable catalogs which are continuously updated on ~24 hour timescales. A new set of tutorials will accompany the Prompt data release.

Related tutorials: See also the Commissioning tutorial 101 on the LSSTCam visits database, which uses a metadata file for science validation visits obtained up to the end of Sep 2025.

1.1. Import packages¶

Import standard python packages astropy, matplotlib, and numpy. Import the skyproj package (skyproj.readthedocs.io) for plotting all-sky projection plots, the healpy package (healpy.readthedocs.io) for dealing with HEALPix, and import from lsst.utils the standard colorblind-friendly colors for the LSST filters, $ugrizy$.

In [1]:

from astropy.io import ascii
from astropy.time import Time
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import skyproj
import healpy as hp
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_linestyles)

1.2. Define parameters¶

Define the colors and linestyles to use for each of the LSST filters.

In [2]:

filter_colors = get_multiband_plot_colors()
filter_names = list(filter_colors.keys())
filter_linestyles = get_multiband_plot_linestyles()

Define the path to the data files used in this tutorial.

In [3]:

file_path = '/rubin/cst_repos/tutorial-notebooks-data/data/'
file_name = file_path + 'lsstcam_visits_by_2026-04-03.ecsv'

Option: to create interactive plots with zoom-in capabilities, un-comment and execute the following cell. All plots created after this will be interactive. To return to static plots, comment-out this cell, restart the kernel and clear all outputs, and re-execute the cells of this notebook.

In [4]:

# %matplotlib widget

2. Visits¶

Column names, descriptions, and units.

• visitId: A unique long integer identifier, composed of the year-month-day and a sequential image number.
• ra: Telescope tracking* Right Ascension, in degrees.
• dec: Telescope tracking* Declination, in degrees.
• start_time: Time at the start of the exposure, in ISO format (International Atomic Time; TAI).
• filter: The LSST filter. One of $ugrizy$.
• exposure_time: Exposure time, in seconds.
• zenith: Telescope pointing distance from zenith, in degrees.
• observation_reason: The source of the visit in the Feature Based Scheduler (FBS).
• target_names: The name of the sky region(s) for the visit, as a comma-separated list if more than one named sky region overlaps the pointing.

*The ra and dec coordinates are for the center of the visit, and were taken from the telescope pointing. They are not the result of an astrometric solution.

2.1. Read the visits table¶

Read the prepared file of visit metadata, print the total number of rows (visits), and display the first five rows of the table.

Warning: The following cell produces a pink InvalidEcsvDatatypeWarning but the datetime64 format for the start_time column is actually read just fine.

In [5]:

visits_table = ascii.read(file_name)
print('Number of visits: ', len(visits_table))
visits_table[:5]

WARNING: InvalidEcsvDatatypeWarning: unexpected datatype 'datetime64[ns]' of column 'start_time' is not in allowed ECSV datatypes ('bool', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64', 'float16', 'float32', 'float64', 'float128', 'string'). Using anyway as a numpy dtype but beware since unexpected results are possible. [astropy.io.ascii.ecsv]

astropy WARNING: InvalidEcsvDatatypeWarning: unexpected datatype 'datetime64[ns]' of column 'start_time' is not in allowed ECSV datatypes ('bool', 'int8', 'int16', 'int32', 'int64', 'uint8', 'uint16', 'uint32', 'uint64', 'float16', 'float32', 'float64', 'float128', 'string'). Using anyway as a numpy dtype but beware since unexpected results are possible.

Number of visits:  32329

Out[5]:

Table length=5

visit_id	ra	dec	start_time	filter	exposure_time	zenith	observation_reason	target_names
int64	float64	float64	datetime64[ns]	str1	float64	float64	str45	str32
2025120900398	59.23116437875386	-48.95615395053968	2025-12-10T06:56:16.873295009	g	30.0	43.7781087040236	ddf_edfs_a	lowdust, ddf_edfs_a
2025120900399	63.174033033099036	-47.752618714422766	2025-12-10T06:56:56.800859571	g	30.0	41.1560290027598	ddf_edfs_b	lowdust, ddf_edfs_b
2025120900400	59.81726041841907	-49.063156298797004	2025-12-10T06:59:41.056697667	i	30.0	43.9631168958183	ddf_edfs_a	lowdust, ddf_edfs_a
2025120900401	62.65539682732036	-48.22511031319913	2025-12-10T07:00:21.007114649	i	30.0	42.1238056452539	ddf_edfs_b	lowdust, ddf_edfs_b
2025120900405	85.26677375387374	-58.08896620777662	2025-12-10T07:12:22.081417143	z	30.0	34.9897919699449	template_blob_z_33.0	lowdust

2.2. MJD distribution¶

Plot the cumulative distribution of the modified Julian dates (MJDs) at the midpoint time of all visits.

In [6]:

visits_table['mjd'] = Time(visits_table['start_time']).mjd

notable_mjds = [np.floor(np.min(visits_table['mjd'])),
                61095.0, 61110.0, 61125.0,
                np.ceil(np.max(visits_table['mjd']))]
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, 20000, s[0:10], rotation=90)
    del t, s
plt.plot(np.sort(visits_table['mjd']), np.arange(len(visits_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()

Figure 1: 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 and end dates included in the file, the date of first alerts (Feb 24 2026), and a period of engineering time when no alerts were released.

Plot the binned cumulative distribution of MJDs for each filter.

In [7]:

fig = plt.figure(figsize=(6, 4))
for f, filt in enumerate(filter_names):
    tx = np.where(visits_table['filter'] == filt)[0]
    plt.hist(visits_table['mjd'][tx], histtype='step', bins=20,
             range=(notable_mjds[0], notable_mjds[-1]), cumulative=True,
             linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    del tx
plt.legend(loc='upper left')
plt.xlabel('MJD')
plt.ylabel('Total number of visits per filter')
plt.title('Binned cumulative distribution of visit MJDs, by filter')
plt.show()

Figure 2: The binned cumulative distribution of visit MJDs by filter (see figure legend).

2.3. Exposure times¶

The exposure times in column exposure_time are in seconds. For filters $grizy$, the standard exposure time is 30 seconds, whereas for $u$-band it is 38 seconds.

Create a stacked bar plot of the exposure time distribution by filter.

In [8]:

xvals = []
clrs = []
lbls = []
for filt in filter_names:
    tx = np.where(visits_table['filter'] == filt)[0]
    if len(tx) > 0:
        xvals.append(visits_table['exposure_time'][tx])
        clrs.append(filter_colors[filt])
        lbls.append(filt)
fig = plt.figure(figsize=(6, 4))
plt.hist(xvals, 10, histtype='bar', stacked=True, facecolor=clrs, label=lbls)
plt.legend(loc='upper left')
plt.xlabel('Exposure time [s]')
plt.ylabel('Number of visits')
plt.show()
del xvals, clrs, lbls

Figure 3: Stacked bar chart of exposure times per visit, by filter.

Although it is clear from the plot above that most visits with exposure times >30s are in the $u$-band, and most with <30s are in the $z$-band, it is unclear if any other filters have a very small number visits in these bins. The bin width furthermore makes it unclear whether all exposure times are exactly the same, or if there is scatter within a bin.

For visits with 30s, <30s, and >30s exposure time, print the mean and standard deviation in the exposure time, and the number of visits per filter.

In [9]:

tx1 = np.where((visits_table['exposure_time'] < 29))[0]
tx2 = np.where((visits_table['exposure_time'] >= 29) & (visits_table['exposure_time'] <= 31))[0]
tx3 = np.where((visits_table['exposure_time'] > 31))[0]
print(f"{'exposure time': <15} {' mean': <4}  {'std': <3}")
print(f"{'<29 s ': <15} {np.round(np.mean(visits_table['exposure_time'][tx1]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['exposure_time'][tx1]), 1): 3.1f}")
print(f"{' 30 s ': <15} {np.round(np.mean(visits_table['exposure_time'][tx2]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['exposure_time'][tx2]), 1): 3.1f}")
print(f"{'>31 s ': <15} {np.round(np.mean(visits_table['exposure_time'][tx3]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['exposure_time'][tx3]), 1): 3.1f}")
print(' ')
vals1 = ['<29 s ']
vals2 = [' 30 s ']
vals3 = ['>31 s ']

for filt in filter_names:
    vals1.append(len(np.where(visits_table['filter'][tx1] == filt)[0]))
    vals2.append(len(np.where(visits_table['filter'][tx2] == filt)[0]))
    vals3.append(len(np.where(visits_table['filter'][tx3] == filt)[0]))
print(f"{'exposure time': <15} {'u': <4}  {'g': <4}  {'r': <4}  {'i': <4}  {'z': <4}  {'y': <4}")
print(f"{vals1[0]: <15} {vals1[1]: <4}  {vals1[2]: <4}  {vals1[3]: <4}  "
      f"{vals1[4]: <4}  {vals1[5]: <4}  {vals1[6]: <4}")
print(f"{vals2[0]: <15} {vals2[1]: <4}  {vals2[2]: <4}  {vals2[3]: <4}  "
      f"{vals2[4]: <4}  {vals2[5]: <4}  {vals2[6]: <4}")
print(f"{vals3[0]: <15} {vals3[1]: <4}  {vals3[2]: <4}  {vals3[3]: <4}  "
      f"{vals3[4]: <4}  {vals3[5]: <4}  {vals3[6]: <4}")

exposure time    mean  std
<29 s            15.0   0.0
 30 s            30.0   0.0
>31 s            38.0   0.0

exposure time   u     g     r     i     z     y
<29 s           0     0     66    36    306   0
 30 s           0     3171  4233  9805  7697  5962
>31 s           1053  0     0     0     0     0

The output above shows that visits are either 15, 30, or 38 seconds exactly (because the standard deviations are 0), and that only $u$-band visits are 38 seconds, and only some $y$-band visits are 15 seconds.

2.4. Airmass distribution¶

Plot the cumulative distribution of airmasses for all visits. Convert the zenith column to airmass, first.

In [10]:

visits_table['airmass'] = 1.0 / np.cos(np.deg2rad(visits_table['zenith']))

fig = plt.figure(figsize=(6, 4))
plt.plot(np.sort(visits_table['airmass']), np.arange(len(visits_table)),
         ls='solid', lw=1, color='black')
plt.xlabel('Airmass')
plt.ylabel('Total number of visits')
plt.title('Cumulative distribution of visit airmass (full)')
plt.show()

Figure 4: The cumulative distribution of visit airmass for all visits, for all filters combined.

Plot the binned cumulative distribution of airmasses for each filter.

In [11]:

fig = plt.figure(figsize=(6, 4))
for f, filt in enumerate(filter_names):
    tx = np.where(visits_table['filter'] == filt)[0]
    plt.hist(visits_table['airmass'][tx], histtype='step', bins=100,
             range=(1.0, np.max(visits_table['airmass'])), cumulative=True, log=True,
             linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    del tx
plt.legend(loc='best', ncol=2)
plt.xlim([0.98, 1.6])
plt.xlabel('Airmass')
plt.ylabel('Total number of visits per filter')
plt.title('Binned cumulative distribution of visit airmass, by filter (cropped)')
plt.show()

Figure 5: The binned cumulative distribution of visit airmass by filter (see figure legend). The x-axis has been cropped to show only airmass $\leq 1.6$.

2.5. Sky distribution¶

Use the healpy package to print the area of one 19-sided HEALPix. It is similar to the LSSTCam FOV of 9.6 square degrees.

In [12]:

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.

Create a two-dimensional (2D) histogram of the visits on the sky, for all filters together. Use 19-sided HEALPix to approximately match the area of one visit on the sky.

In [13]:

fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
vras = np.asarray(visits_table['ra'], dtype='float')
vdecs = np.asarray(visits_table['dec'], dtype='float')
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap='Greys', vmin=-10)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()

Figure 6: A 2D histogram illustrating the distribution of visits on the sky (for all filters, combined).

Print the approximate total sky area covered by visits.

In [14]:

hpx_bins = sp.draw_hpxbin(vras, vdecs, nside=19)
array = hpx_bins[0]
tx = np.where(array > 0)[0]
print('Number of HEALPix with non-zero visits: ', len(tx))
print('Approximate area of all visits: ', int(len(tx) * hp.nside2pixarea(19, degrees=True)),
      ' square degrees.')
del sp, vras, vdecs, hpx_bins, array, tx

Number of HEALPix with non-zero visits:  2236
Approximate area of all visits:  21293  square degrees.

Show a similar 2D histogram as above, but only for the g-band filter, and use the "Greens" colormap.

In [15]:

show_filter = 'g'
use_colormap = 'Greens'
fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
tx = np.where(visits_table['filter'] == show_filter)[0]
vras = np.asarray(visits_table['ra'][tx], dtype='float')
vdecs = np.asarray(visits_table['dec'][tx], dtype='float')
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap=use_colormap, vmin=-10)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()
del sp, tx, vras, vdecs, show_filter, use_colormap

Figure 7: Similar to Figure 6, but for $g$-band visits only.

In Figure 6, stripes of visits with constant declination stand out in the sky distribution. The cause of this non-LSST-like survey pattern is in their observation reason: "fbs_driven_aos_stability_test". FBS stands for "feature-based scheduler" and AOS for "active optics system". These visits were obtained by the FBS for the purpose of testing the AOS, but are still anticipated to be, potentially, scientifically useful.

Visualize only the FBS AOS testing visits.

In [16]:

use_colormap = 'Purples'
fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
tx = np.where(visits_table['observation_reason'] == 'fbs_driven_aos_stability_test')[0]
vras = np.asarray(visits_table['ra'][tx], dtype='float')
vdecs = np.asarray(visits_table['dec'][tx], dtype='float')
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap=use_colormap, vmin=-10)
sp.ax.set_xlabel("Right Ascension", fontsize=12)
sp.ax.set_ylabel("Declination", fontsize=12, labelpad=22)
plt.show()

Figure 8: Similar to Figure 6, but for visits with an observation reason related to AOS stability testing.

Instead of a 2D histogram, the cell below offers the option to mark each visit as a separate, semi-transparent marker the size of the LSSTCam field of view (FOV) on the sky (radius ~1.75 deg).

Warning: The mode of plotting in the following cell takes a few minutes to render, is not as informative as a histogram, and will not scale well as the number of visits increases. But it is possible, as there are only ~30,000 visits in this temporary file.

In [17]:

# fig, ax = plt.subplots(figsize=(8, 6))
# sp = skyproj.McBrydeSkyproj(ax=ax)
# vras = np.asarray(visits_table['ra'], dtype='float')
# vdecs = np.asarray(visits_table['dec'], dtype='float')
# for ra, dec in zip(vras, vdecs):
#     sp.ax.circle(ra, dec, 1.75, edgecolor=None, color='grey', alpha=0.1, fill=True)
# plt.show()
# del sp, vras, vdecs

3. Templates overlap¶

At the time this notebook was made, template images for alert production only existed in the deep drilling fields (DDFs).

Define the coordinates, names, and filter coverage for the DDF templates.

In [18]:

ddf_coords = SkyCoord([[150.11, 2.23], [52.98, -28.12], [9.45, -44.02],
                       [58.9, -49.32], [63.6, -47.6], [35.57, -4.82],
                       [187.4, +8]], frame="icrs", unit="deg")
ddf_names = ["COSMOS", "ECDFS", "ELAIS-S1", "EDFS_a", "EDFS_b",
             "XMM-LSS", "M49 (Virgo)"]
ddf_filters = ["ugrizy", "riz", "griz", "griz", "griz", "iz", "ugri"]

Define the radius of a visit. This is the radius of the LSSTCam's FOV in degrees.

In [19]:

radius = 1.75

Tally the number of visits for which the boresight RA, Dec are within one LSSTCam radius of the center of a DDF, for each of the filters for which that DDF has a template image.

In [20]:

ddf_tally = np.zeros((len(ddf_coords), 6), dtype='int')
for visit in visits_table:
    coord = SkyCoord(visit['ra'], visit['dec'], frame="icrs", unit="deg")
    seps = coord.separation(ddf_coords).deg
    tx = np.where(seps < radius)[0]
    if len(tx) == 1:
        temp_list = ddf_filters[tx[0]]
        if temp_list.find(visit['filter']) >= 0:
            fx = np.where(visit['filter'] == np.asarray(filter_names))[0]
            ddf_tally[tx[0], fx[0]] += 1
            del fx
        del temp_list
    del coord, seps, tx

Print the number of overlapping visits per filter for each DDF, with a - if the field has no template in that filter. Also print the total number of overlapping visits.

Notice: This tally is indicative of the number of visits that produced alerts, but alert production a specific step of Prompt processing and did not necessarily proceed (or succeed) for each overlapping visit.

In [21]:

line = f"{'field': <15} "
for f, filt in enumerate(filter_names):
    line += f"{filt: <4} "
line += f"{'sum ': <5}"
print(line)
for d, ddf_name in enumerate(ddf_names):
    temp_list = ddf_filters[d]
    line = f"{ddf_name: <15} "
    for f, filt in enumerate(filter_names):
        if temp_list.find(filt) >= 0:
            line += f"{ddf_tally[d, f]: <4} "
        else:
            line += f"{'-   ': <4} "
    line += f"{np.sum(ddf_tally[d, :]): <5}"
    print(line)

field           u    g    r    i    z    y    sum
COSMOS          45   279  261  502  258  32   1377
ECDFS           -    -    13   41   17   -    71
ELAIS-S1        -    15   12   27   9    -    63
EDFS_a          -    106  120  148  109  -    483
EDFS_b          -    105  120  150  110  -    485
XMM-LSS         -    -    -    28   7    -    35
M49 (Virgo)     0    206  147  198  -    -    551

Note that the observation_reason and target_names columns, instead of or in addition to the spatial overlap, also indicate whether a given visit was an observation of a DDF and thus may have produced alerts.

Recreate Figure 6, but overplot the locations of the template fields.

In [22]:

ddf_sym = ['o', 's', 'p', '^', '^', 'v', '*']
ddf_size = [10, 10, 10, 10, 10, 10, 15]
ddf_color = ['green', 'darkorange', 'purple', 'red', 'magenta', 'blue', 'brown']

fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
vras = np.asarray(visits_table['ra'], dtype='float')
vdecs = np.asarray(visits_table['dec'], dtype='float')
sp.draw_hpxbin(vras, vdecs, nside=19, alpha=1, cmap='Greys', vmin=-10)
for d, coord in enumerate(ddf_coords):
    sp.ax.plot(coord.ra.deg, coord.dec.deg, ddf_sym[d], ms=ddf_size[d],
               mec=ddf_color[d], color='None', mew=2, label=ddf_names[d])
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

Figure 9: Similar to Figure 6, but with the locations of the template fields marked with symbols as in the legend.

4. Was my object in a visit?¶

Recall the caveat from Section 1: this temporary, preliminary metadata file includes visits that were obtained in the specific time frame, and it cannot be guaranteed that the visit will be successfully processed and pass data validation for inclusion in the Prompt data release.

Caveat: furthermore, coordinates that are near the edge of a visit might seem to be included using the code below, but this is not a guarantee that there will be high-quality data exactly at that location.

Define two target coordinates. Based on the figures above, the first one will, but the second one will not, be included.

In [23]:

target1_ra, target1_dec = 90.0, -15.0
target2_ra, target2_dec = 90.0, +30.0

Convert the targets' coordinates to the astropy SkyCoord class.

In [24]:

c1 = SkyCoord(target1_ra, target1_dec, unit="deg")
c2 = SkyCoord(target2_ra, target2_dec, unit="deg")

Use the radius of a visit defined above.

In [25]:

print(radius)

1.75

Convert the visits' coordinates to the astropy SkyCoord class.

In [26]:

vc = SkyCoord(visits_table['ra'], visits_table['dec'], unit="deg")

Calculate the on-sky separations between each target and all of the visits.

In [27]:

s1 = vc.separation(c1)
s2 = vc.separation(c2)

Sum up the number of visits that are within the LSSTCam radius of each target.

Assert the expectation that target 1 has overlapping visits and target 2 does not, and print the result.

In [28]:

tx1 = np.where(s1.degree < radius)[0]
tx2 = np.where(s2.degree < radius)[0]
assert len(tx1) != 0
assert len(tx2) == 0
print('Number of visits containing')
print('  target 1: ', len(tx1))
print('  target 2: ', len(tx2))
del tx1, tx2

Number of visits containing
  target 1:  5
  target 2:  0