102. Visit and tract metadata

102. Visit and tract metadata¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
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 the visits and tracts expected to for DP2.

LSST data products: Preliminary lists of visit and tract metadata for DP2.

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¶

Preliminary files of visit and tract metadata for Data Preview 2 (DP2) are made available with this tutorial. The purpose of these files and this tutorial is to enable early exploration and visualization of the visits and tracts that are expected to be included in DP2. When DP2 is released, these files will be deprecated, and the visit and tract metadata will be available via queryable catalogs. A new set of tutorials will accompany the DP2 release.

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

Tracts: One tract is one square region of LSST's all-sky tesselation ("skymap"), $~1.66$ deg per side. Tracts are subdivided into 100 overlapping patches, and one deep coadd image is created per patch.

For DP2, the total sky area covered by DP2 visits is greater than the area covered by DP2 tracts. The total sky area over 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.

Caveat: These temporary, preliminary metadata files include all visits and tracts that were initially processed, and might include those that do not ultimately pass data validation or produce science-quality data products.

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
from astropy.wcs import WCS
from regions import PolygonSkyRegion
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/'

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¶

The file dp2_visits_table.ecsv contains a preliminary, temporary table of visits that will be included in DP2.

Column names, descriptions, and units.

• visitId: A unique long integer identifier, composed of the year-month-day and a sequential image number.
• band: The LSST filter. One of $ugrizy$.
• ra: Right Ascension, in degrees.
• dec: Declination, in degrees.
• skyRotation: Camera rotation in degrees east of north.
• azimuth: Telescope pointing azimuth, in degrees.
• altitude: Telescope pointing altitude, in degrees.
• zenithDistance: Telescope pointing distance from zenith, in degrees.
• airmass: Airmass of the telescope pointing.
• expTime: Exposure time, in seconds.
• expMidptMJD: Time at the middle of the exposure, in modified Julian date (MJD).
• obsStartMJD: Time at the start of the exposure, in MJD.
• mean_seeing: Initial estimate of seeing, averaged over all detectors, in arcseconds.
• mean_maglim: Initial estimate of PSF magnitude limit, averaged over all detectors, in mags.

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.

In [5]:

visits_table = ascii.read(file_path + 'dp2_visits_table_with_iq.ecsv')
print('Number of visits: ', len(visits_table))
visits_table[:5]

Number of visits:  28698

Out[5]:

Table length=5

visitId	band	ra	dec	skyRotation	azimuth	altitude	zenithDistance	airmass	expTime	expMidptMJD	obsStartMJD	mean_seeing	mean_maglim
int64	str1	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64
2025042400172	r	185.443060811076	5.98102138504802	146.113088289607	1.9614440076983797	53.8990970549479	36.1009029450521	1.23704795139999	29.999915599823	60790.11722739659	60790.11705378597	1.3717809522436122	23.90346908569336
2025042400173	r	187.349245626471	8.23604638508455	144.18450391305	4.59267468721633	51.5610973958556	38.4389026041444	1.27595524322485	30.0000872612	60790.1177404311	60790.117566819485	1.360464795099889	23.798002243041992
2025042400174	r	186.443750064719	7.16885705640284	138.496645426916	2.93858842208796	52.6885249328851	37.3114750671149	1.2566323773581	30.0009806156158	60790.11824723747	60790.11807362069	1.357174265291119	23.803224563598633
2025042400175	r	185.542702196852	6.0999042479778	132.754022698538	1.19423996978959	53.7910062641441	36.2089937358559	1.23874981794006	30.0010826587677	60790.118757722914	60790.118584105534	1.2905303859739274	23.877490997314453
2025042400176	r	187.450364693337	8.35450149584232	130.813892358106	3.8613155396456	51.471714859197	38.528285140803	1.27753357366488	30.0010211467743	60790.11927441511	60790.11910079809	1.390931531603405	23.81049156188965

2.2. Exposure times¶

The exposure times in column expTime 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 [6]:

xvals = []
clrs = []
lbls = []
for filt in filter_names:
    tx = np.where(visits_table['band'] == filt)[0]
    if len(tx) > 0:
        xvals.append(visits_table['expTime'][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.ylim([-500, 27000])
plt.xlabel('Exposure time [s]')
plt.ylabel('Number of visits')
plt.show()
del xvals, clrs, lbls

Figure 1: 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, it is unclear which other filters have visits with shorter or longer exposure times. 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 [7]:

tx1 = np.where((visits_table['expTime'] < 29))[0]
tx2 = np.where((visits_table['expTime'] >= 29) & (visits_table['expTime'] <= 31))[0]
tx3 = np.where((visits_table['expTime'] > 31))[0]
print(f"{'exposure time': <15} {' mean': <4}  {'std': <3}")
print(f"{'<29 s ': <15} {np.round(np.mean(visits_table['expTime'][tx1]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['expTime'][tx1]), 1): 3.1f}")
print(f"{' 30 s ': <15} {np.round(np.mean(visits_table['expTime'][tx2]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['expTime'][tx2]), 1): 3.1f}")
print(f"{'>31 s ': <15} {np.round(np.mean(visits_table['expTime'][tx3]), 1): 4.1f}  "
      f"{np.round(np.std(visits_table['expTime'][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['band'][tx1] == filt)[0]))
    vals2.append(len(np.where(visits_table['band'][tx2] == filt)[0]))
    vals3.append(len(np.where(visits_table['band'][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     14    0     85    49    0
 30 s           63    4058  4848  7870  5760  3944
>31 s           1901  94    8     4     0     0

2.3. MJD distribution¶

Plot the cumulative distribution of the modified Julian dates (MJDs) at the midpoint time of all exposures (use the expMidptMJD column).

Define a few "notable dates" to mark on the figure (see caption), then create a plot of the cumulative number of visits over time.

In [8]:

notable_mjds = [np.floor(np.min(visits_table['expMidptMJD'])),
                60857, 60881, 60973,
                np.floor(np.max(visits_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, 20000, s[0:10], rotation=90)
    del t, s
plt.plot(np.sort(visits_table['expMidptMJD']), 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 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 (Oct 25 2025).

Plot the binned cumulative distribution of MJDs for each filter.

In [9]:

fig = plt.figure(figsize=(6, 4))
for f, filt in enumerate(filter_names):
    tx = np.where(visits_table['band'] == filt)[0]
    plt.hist(visits_table['expMidptMJD'][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 3: The binned cumulative distribution of visit MJDs by filter (see figure legend).

2.4. Airmass distribution¶

Plot the cumulative distribution of airmasses for all visits (use the airmass column).

In [10]:

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['band'] == filt)[0]
    plt.hist(visits_table['airmass'][tx], histtype='step', bins=75,
             range=(1.0, 3.0), 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. Initial IQ estimates¶

Plot the initial estimates of seeing and magnitude limit, averaged over all detectors, per visit per filter.

In [12]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for f, filt in enumerate(filter_names):
    tx = np.where(visits_table['band'] == filt)[0]
    ax[0].hist(visits_table['mean_seeing'][tx], histtype='step', bins=75,
               cumulative=True, log=True,
               linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    ax[1].hist(visits_table['mean_maglim'][tx], histtype='step', bins=75,
               cumulative=True, log=True,
               linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    del tx
ax[1].legend(loc='upper left', ncol=2)
ax[0].set_xlabel('Mean Seeing [arcsec]')
ax[1].set_xlabel('Mean Magnitude Limit [mag]')
ax[0].set_ylabel('Number of visits per filter')
plt.subplots_adjust(wspace=0)
plt.suptitle('Binned cumulative distribution of initial visit image quality estimates, by filter')
plt.show()

Figure 6: The binned cumulative distribution of initial estimates for the mean seeing (left) and mean magnitude limit (right) per visit (averaged over all detectors), by filter.

2.6. Sky distribution¶

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

In [13]:

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

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 7: A 2D histogram illustrating the distribution of DP2 visits on the sky (for all filters, combined).

Print the approximate total sky area covered by visits.

In [15]:

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 DP2 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:  1633
Approximate area of all DP2 visits:  15550  square degrees.

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

In [16]:

show_filter = 'g'
use_colormap = 'Greens'
fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
tx = np.where(visits_table['band'] == 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 8: Similar to Figure 6, but for $g$-band visits only.

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 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 ~28,000 visits of DP2. If the matplotlib widget for interactive plots has been enabled, the plot will still display but it will not zoom smoothly.

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. Tracts¶

The file dp2_tracts_table.ecsv contains a preliminary, temporary table of tracts that will have at least one patch with a deep coadd image in DP2.

Column names, descriptions, and units.

• tractId: A unique integer identifier; the tract number.
• filters_list: A string of up to six characters containing the filters with a deep coadd image in the tract (e.g., gi, ugr, ugrizy).
• vertices_deg: A list of the sky coordinates for the tract's four corners. For north-up east-left alignment (RA increases left), the order of the corners are: south-east (lower-left), south-west (lower-right), north-west (upper-right), north-east (upper-left).
• psf_maglim: The $5\sigma$ PSF magnitude limit in the $r$ and $i$ bands, in AB magnitudes.
• psf_size: The PSF characteristic width as computed from the determinant radius, in the $r$ and $i$ bands, in pixels.

The two image quality (IQ) parameters, the magnitude limit and PSF size, are evaluated from the survey property maps at the center of the tract:

• deepCoadd_psf_maglim_consolidated_map_weighted_mean
• deepCoadd_psf_size_consolidated_map_weighted_mean

The psf_size is the $\sigma$ in pixels, and can be approximately converted to arcseconds with a pixel scale of $\sim0.2$ arcsec/pixel. To convert the psf_size to a FWHM in arcseconds, use: $FWHM = 0.2 \sigma \times 2 \sqrt{2 \ln(2)} = 2.3548 \times 0.2 \sigma$

3.1. Read the tracts table¶

Read the prepared file of tract metadata. Print the total number of rows (tracts) and an upper limit on the total sky area covered.

In [18]:

tracts_table = ascii.read(file_path + 'dp2_tracts_table_with_iq.ecsv')
print('Number of deep coadd tracts: ', len(tracts_table))
print('Maximum sky area covered by deep coadd tracts: ', int(len(tracts_table) * 1.66**2),
      ' square degrees.')

Number of deep coadd tracts:  2193
Maximum sky area covered by deep coadd tracts:  6043  square degrees.

Caveat: The maximum sky area covered by deep coadd tracts is an upper limit for two reasons. One, all tracts overlap and these overlap areas have not been accounted for. Two, a tract will be included in the file even if only one of its 100 patches has a deep coadd image.

Print the number of tracts covered by all six filters, $ugrizy$, and the sky area covered. This is also an upper limit.

In [19]:

tx = np.where(tracts_table['filters_list'] == 'ugrizy')[0]
print('Number of deep coadd tracts with ugrizy: ', len(tx))
print('Maximum sky area covered: ', int(len(tx) * 1.66**2), ' square degrees.')
del tx

Number of deep coadd tracts with ugrizy:  1244
Maximum sky area covered:  3427  square degrees.

Display the first five rows of the table.

In [20]:

tracts_table[:5]

Out[20]:

Table length=5

tractId	filters_list	vertices_deg	psf_maglim_r	psf_size_r	psf_maglim_i	psf_size_i
int64	str6	float64[4,2]	float64	float64	float64	float64
2054	gi	1.3578396241683794 .. -50.481666831124706	nan	nan	nan	nan
2055	gi	3.695501961830713 .. -50.481666831124706	nan	nan	22.74056826049059	3.4460980395425973
2056	gi	6.033164299493053 .. -50.481666831124706	nan	nan	22.548426447647294	3.485872864308098
2057	gi	8.370826637155393 .. -50.4816668311247	nan	nan	22.704788096643394	3.643065677369538
2058	gi	10.708488974817731 .. -50.481666831124706	nan	nan	22.43389261906287	3.4671865722720625

3.2. Coverage by filter¶

Print the total number (and fraction) of tracts with deep coadd images for each filter.

In [21]:

for filt in filter_names:
    tally = 0
    for flist in tracts_table['filters_list']:
        if flist.find(filt) >= 0:
            tally += 1
    print(filt, tally, np.round(tally/len(tracts_table), 2))
del tally

u 1416 0.65
g 1799 0.82
r 2010 0.92
i 2031 0.93
z 1891 0.86
y 1596 0.73

Some tracts are not fully filled by deep coadd images, and have a nan value for the PSF magnitude limit or PSF size. Print the total number (and fraction) of $r$- and $i$-band tracts that have sufficient data for the PSF magnitude limit to be estimated.

In [22]:

for filt in ['r', 'i']:
    tally = 0
    tally2 = 0
    for f, flist in enumerate(tracts_table['filters_list']):
        if flist.find(filt) >= 0:
            tally += 1
            if np.isfinite(tracts_table['psf_maglim_' + filt][f]):
                tally2 += 1
    print(filt, tally, tally2, np.round(tally2/tally, 2))
del tally

r 2010 1536 0.76
i 2031 1594 0.78

3.3. IQ parameters¶

Plot histograms of the tract image quality (IQ) parameters provided in the table.

In [23]:

fig, ax = plt.subplots(1, 2, figsize=(10, 4), sharey=True)
for f, filt in enumerate(['r', 'i']):
    tx1 = np.where(np.isfinite(tracts_table['psf_maglim_' + filt]))[0]
    tx2 = np.where(np.isfinite(tracts_table['psf_size_' + filt]))[0]
    ax[0].hist(tracts_table['psf_maglim_' + filt][tx1], histtype='step',
               bins=75, cumulative=True, log=True,
               linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    ax[1].hist(0.2 * 2.3548 * tracts_table['psf_size_' + filt][tx2], histtype='step',
               bins=75, cumulative=True, log=True,
               linestyle=filter_linestyles[filt], color=filter_colors[filt], label=filt)
    del tx1, tx2
ax[1].legend(loc='upper left', ncol=2)
ax[0].set_xlabel('PSF Magnitude Limit [mag]')
ax[1].set_xlabel('PSF FWHM [arcsec]')
ax[0].set_ylabel('Number of tracts per filter')
plt.subplots_adjust(wspace=0)
plt.suptitle('Binned cumulative distribution of tract image quality parameters')
plt.show()

Figure 9: The binned cumulative distribution of tract IQ parameters, the PSF magnitude limit (left) and the PSF FWHM (right) in the $r$- and $i$-filters.

3.4. Sky distribution¶

Use the tract vertices to draw one polygon per tract.

In [24]:

fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
for vertices in tracts_table['vertices_deg']:
    ras = vertices.transpose()[0]
    decs = vertices.transpose()[1]
    sp.draw_polygon(ras, decs, edgecolor='black', alpha=0.5, linewidth=0.5, facecolor=None)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()
del sp

Figure 10: All of the DP2 tracts for which deep coadd images will exist, in any filter.

Create a plot that is similar to the one above, but for only tracts that will have an $r$-band deep coadd image.

In [25]:

show_filter = 'r'
fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
for i, vertices in enumerate(tracts_table['vertices_deg']):
    flist = str(tracts_table['filters_list'][i])
    if flist.find(show_filter) >= 0:
        ras = vertices.transpose()[0]
        decs = vertices.transpose()[1]
        sp.draw_polygon(ras, decs, edgecolor=filter_colors[show_filter],
                        alpha=0.5, linewidth=0.5, facecolor=None)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()
del sp, show_filter

Figure 11: The DP2 tracts for which $r$-band deep coadd images will exist.

Create a plot that is similar to the one above, but only for tracts with $r$-band for which the $r$-band magnitude limit is $>23$ mag.

In [26]:

show_filter = 'r'
fig, ax = plt.subplots(figsize=(8, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
for i, vertices in enumerate(tracts_table['vertices_deg']):
    flist = str(tracts_table['filters_list'][i])
    if flist.find(show_filter) >= 0:
        if (tracts_table['psf_maglim_r'][i] > 23.0):
            ras = vertices.transpose()[0]
            decs = vertices.transpose()[1]
            sp.draw_polygon(ras, decs, edgecolor=filter_colors[show_filter],
                            alpha=0.5, linewidth=0.5, facecolor=None)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()
del sp, show_filter

Figure 12: The DP2 tracts for which $r$-band deep coadd images will exist, and have a magnitude limit fainter than 23 mag.

4. Is my object in DP2?¶

Recall the caveat from Section 1: these temporary, preliminary metadata files include all visits and tracts that were processed, and might include those that do not, ultimately, pass data validation or produce science-quality data products.

Caveat: coordinates that are near the edge of a visit or tract might seem to be included in DP2 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 DP2.

In [27]:

target1_ra, target1_dec = 300.0, -15.0
target2_ra, target2_dec = 300.0, +15.0

Convert the targets' coordinates to the astropy SkyCoord class.

In [28]:

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

4.1. Visits¶

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

In [29]:

radius = 1.75

Convert the visits' coordinates to the astropy SkyCoord class.

In [30]:

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

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

In [31]:

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

Sum up the number of visits that are within the LSSTCam radius of each target, and print the result.

In [32]:

tx1 = np.where(s1.degree < radius)[0]
tx2 = np.where(s2.degree < radius)[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:  49
  target 2:  0

As expected, target 1 has overlapping visits and target 2 does not.

In [33]:

del radius, vc, s1, s2

4.2. Tracts¶

For the first tract in the table, get the vertices (corners), convert them to SkyCoord, and create a PolygonSkyRegion for the tract.

In [34]:

row = 0
vertices = tracts_table['vertices_deg'][row]
ras = vertices.transpose()[0]
decs = vertices.transpose()[1]
vertices_sc = SkyCoord(ras, decs, unit='deg')
region_sky = PolygonSkyRegion(vertices=vertices_sc)
print(region_sky)
del row, vertices, ras, decs, vertices_sc

Region: PolygonSkyRegion
vertices: <SkyCoord (ICRS): (ra, dec) in deg
    [(  1.35783962, -52.14776467), (358.64206988, -52.14776363),
     (358.6905163 , -50.48166585), (  1.30939644, -50.48166683)]>

Define a third target coordinate that is within the first tract's vertices, and convert it to a SkyCoord.

In [35]:

target3_ra, target3_dec = 359.0, -51.0
c3 = SkyCoord(target3_ra, target3_dec, unit="deg")

Use the contains method for the region_sky to check if any of targets 1, 2, or 3 are within the first tract. This method requires that a generic World Coordinate System (WCS) be passed, so create one first.

The first two coordinates will not be, but the third was especially defined to be within the first tract. The output will be "False, False, True".

In [36]:

wcs = WCS(naxis=2)
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
print(region_sky.contains(c1, wcs))
print(region_sky.contains(c2, wcs))
print(region_sky.contains(c3, wcs))

[False]
[False]
[ True]

For each of the three coordinates, check whether any tracts overlap. Tally the number of overlapping tracts and print the number.

For targets 1, 2, and 3, the resulting tally will be 1, 0, and 1.

In [37]:

for c, coord in enumerate([c1, c2, c3]):
    tally = 0
    for row in range(len(tracts_table)):
        vertices = tracts_table['vertices_deg'][row]
        ras = vertices.transpose()[0]
        decs = vertices.transpose()[1]
        vertices_sc = SkyCoord(ras, decs, unit='deg')
        region_sky = PolygonSkyRegion(vertices=vertices_sc)
        if region_sky.contains(coord, wcs):
            tally += 1
        del vertices, ras, decs, vertices_sc, region_sky
    print('For target', c+1, 'the number of overlapping tracts =', tally)
del tally

For target 1 the number of overlapping tracts = 1
For target 2 the number of overlapping tracts = 0
For target 3 the number of overlapping tracts = 1

In [38]:

del c1, c2, c3

4.3. Sky distribution¶

Combine the visits and tracts maps into a single plot, and show the locations of the three targets.

In [39]:

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='Blues', vmin=0)
sp.draw_colorbar(label='number of visits', fontsize=10, pad=0.01, shrink=0.5)
for vertices in tracts_table['vertices_deg']:
    ras = vertices.transpose()[0]
    decs = vertices.transpose()[1]
    sp.draw_polygon(ras, decs, edgecolor='black', alpha=0.5, linewidth=0.5, facecolor=None)
    del ras, decs
sp.ax.plot(target1_ra, target1_dec, '*', ms=10, color='yellow')
sp.ax.plot(target2_ra, target2_dec, '^', ms=8, color='green')
sp.ax.plot(target3_ra, target3_dec, 's', ms=6, color='magenta')
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()
del sp, vras, vdecs

Figure 13: A 2D histogram illustrating the distribution of DP2 visits on the sky (for all filters, combined), similar to Figure 6 but plotted in blue-scale. Boxes are drawn to represent each skymap tract that will have at least one patch with a deep coadd image (similar to Figure 8). Overplotted are the coordinates of the three target coordinates (1 as a yellow star; 2 as a green triangle; 3 as a magenta square).

In [40]:

del target1_ra, target2_ra, target3_ra, target1_dec, target2_dec, target3_dec

5. Which DP1 fields are in DP2?¶

Define the DP1 field centers.

In [41]:

dp1_field_names = ['47 Tuc globular cluster',
                   'Low Ecliptic Latitude Field',
                   'Fornax Dwarf Spheroidal Galaxy',
                   'Extended Chandra Deep Field South',
                   'Euclid Deep Field South',
                   'Low Galactic Latitude Field',
                   'Seagull Nebula']
dp1_field_ras = [6.02, 37.86, 40.00, 53.13, 59.10, 95.00, 106.23]
dp1_field_decs = [-72.08, 6.98, -34.45, -28.10, -48.73, -25.00, -10.51]

5.1. Visits¶

Print the number of visits, in any filter, that overlap the central coordinates of each DP1 field.

In [42]:

radius = 1.75
vc = SkyCoord(visits_table['ra'], visits_table['dec'], unit="deg")
for i in range(len(dp1_field_names)):
    field_coord = SkyCoord(dp1_field_ras[i], dp1_field_decs[i], unit="deg")
    separations = vc.separation(field_coord)
    tx = np.where(separations.degree < radius)[0]
    print('%-35s %4i' % (dp1_field_names[i], len(tx)))
    del field_coord, separations, tx
del radius, vc

47 Tuc globular cluster                8
Low Ecliptic Latitude Field            0
Fornax Dwarf Spheroidal Galaxy         3
Extended Chandra Deep Field South    294
Euclid Deep Field South              175
Low Galactic Latitude Field           11
Seagull Nebula                         3

5.2. Tracts¶

Print whether there is a deep coadd tract that overlaps the central coordinates of each DP1 field.

In [43]:

for i in range(len(dp1_field_names)):
    field_coord = SkyCoord(dp1_field_ras[i], dp1_field_decs[i], unit="deg")
    overlap = False
    for row in range(len(tracts_table)):
        vertices = tracts_table['vertices_deg'][row]
        ras = vertices.transpose()[0]
        decs = vertices.transpose()[1]
        vertices_sc = SkyCoord(ras, decs, unit='deg')
        region_sky = PolygonSkyRegion(vertices=vertices_sc)
        if region_sky.contains(field_coord, wcs):
            overlap = True
        del vertices, ras, decs, vertices_sc, region_sky
    print('%-35s %5s' % (dp1_field_names[i], overlap))
    del field_coord, overlap

47 Tuc globular cluster             False
Low Ecliptic Latitude Field         False
Fornax Dwarf Spheroidal Galaxy      False
Extended Chandra Deep Field South    True
Euclid Deep Field South              True
Low Galactic Latitude Field         False
Seagull Nebula                      False

5.3. Sky distribution¶

Create a similar plot as in Figure 10, but with the locations of the DP1 fields marked.

In [44]:

dp1_field_symbols = ['o', 's', 'p', '*', '^', 'v', 'x']
dp1_field_symsizes = [6, 6, 6, 8, 6, 6, 8]
dp1_field_colors = ['red', 'lightseagreen', 'darkviolet',
                    'magenta', 'yellow', 'lime', 'darkorange']
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='Blues', vmin=0)
for vertices in tracts_table['vertices_deg']:
    ras = vertices.transpose()[0]
    decs = vertices.transpose()[1]
    sp.draw_polygon(ras, decs, edgecolor='black', alpha=0.5, linewidth=0.5, facecolor=None)
    del ras, decs
for i in range(len(dp1_field_names)):
    sp.ax.plot(dp1_field_ras[i], dp1_field_decs[i], dp1_field_symbols[i], ms=dp1_field_symsizes[i],
               color=dp1_field_colors[i], label=dp1_field_names[i])
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
sp.ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
plt.show()
del sp, vras, vdecs

Figure 14: Similar to Figure 10, this plot shows the distribution of DP2 visits (blue-scale), with boxes drawn for DP2 tracts, and with the locations of the DP1 fields marked.

In [ ]: