203.2. The skyMap

203.2. The skyMap¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: large
LSST Science Pipelines version: r30.0.11
Last verified to run: 2026-09-15
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Understand how to access and use the skymap.

LSST data products: SkyMap.

Packages: skyproj, lsst.daf.butler, lsst.geom, lsst.sphgeom

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 demonstrates how to load and work with the skyMap. A given skyMap contains the definition of the tesselation of the sky into tracts, patches, and cells that are the fundamental blocks used for coaddition and object measurement. The skyMap object contains methods that can be used to find tracts, patches, and cells and their boundaries given an input coordinate.

Related tutorials: Information about tract/patch geometry can also be retrieved using the CoaddPatches table demonstrated in a 201-series tutorial. Detailed explorations of Object tables and DeepCoadd images are also in the 201 level. See also the 301-level overview notebooks for more details about the DP2 survey geometry.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery). Import the skyproj package for plotting astronomical skymaps (skyproj.readthedocs.io).

From the lsst package, import the RSPDiscovery class for accessing the RSP Table Access Protocol (TAP) service. Import the butler, geom and sphgeom packages, and plotting utilities from the LSST Science Pipelines (pipelines.lsst.io).

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import skyproj

from lsst.rsp import RSPDiscovery
from lsst.daf.butler import Butler
import lsst.geom as geom
import lsst.sphgeom as sphgeom
from lsst.utils.plotting import set_rubin_plotstyle

1.2. Define parameters and functions¶

Define parameters to use colorblind-friendly colors with matplotlib.

In [2]:

set_rubin_plotstyle()

Instantiate the butler, and point it to the DP2 repo and collection.

In [3]:

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

Instantiate RSPDiscovery with the DP2 release, and create an instance of the TAP service.

In [4]:

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

Define the central coordinates of a field to use throughout. For later use, create an lsst.geom.SpherePoint object at that coordinate.

In [5]:

ra = 290.0
dec = -20.8

coord = geom.SpherePoint(ra*geom.degrees, dec*geom.degrees)

2. Retrieve and examine the skyMap¶

The skymap defines the all-sky tesselation of tracts and patches for the deep coadd images.

Use the Butler's get_dataset_type function to show that the dimension for the skyMap dataset type is just skymap, and that the returned format for the skymap is SkyMap.

Show the schema for the skyMap.

In [6]:

dataset_type = butler.get_dataset_type('skyMap')
print(f"dataset_type: {dataset_type}\n")
for dimension in dataset_type.dimensions.data_coordinate_keys:
    print('dimension = ', dimension)
    print(butler.dimensions[dimension].schema)
    print(' ')

dataset_type: DatasetType('skyMap', {skymap}, SkyMap)

dimension =  skymap
skymap:
  name: string
  hash: hash
    A hash of the skymap's parameters.
  tract_max: int
    Maximum ID for tracts in this skymap, exclusive.
  patch_nx_max: int
    Number of patches in the x direction in each tract.
  patch_ny_max: int
    Number of patches in the y direction in each tract.

The skyMap DatasetType requires specifying the skymap that was used to create it. See what keys are available for the skymap dimension using query_dimension_records.

In [7]:

butler.query_dimension_records('skymap')

Out[7]:

[skymap.RecordClass(name='lsst_cells_v2', hash=b'\xa5\xf1\x1eV\x15\x05\xdck$\xec\x8fN\xf5G\xbcY.\x82\x99\xda', tract_max=18938, patch_nx_max=10, patch_ny_max=10)]

There is only the lsst_cells_v2 skymap for DP2.

Retrieve the skyMap from the butler, specifying lsst_cells_v2 in the skymap argument.

In [8]:

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

Uncomment all lines in the following cell to see all methods associated with the returned skyMap.

In [9]:

# for tmp in dir(skymap):
#     if tmp[0] != '_':
#         print(tmp)

2.1. Tracts and patches via coordinates¶

2.1.1. Nearest tract/patch to coordinate¶

Find the tract and patch whose center is nearest to the coordinate specified in Section 1. This method expects the input to be a list of SpherePoint objects (note that the list can contain multiple entries, in which case it will return the closest tract/patch for each item in the list).

In [10]:

closest_tract_patch = skymap.findClosestTractPatchList([coord])
closest_tract_patch

Out[10]:

[(TractInfo(id=6324, ctrCoord=[0.30976315626767525, -0.8866455017580095, -0.34337521930298853]),
  (PatchInfo(index=Index2D(x=0, y=0), innerBBox=(minimum=(0, 0), maximum=(2999, 2999)), outerBBox=(minimum=(-150, -150), maximum=(3149, 3149)), cellInnerDimensions=(150, 150), cellBorder=0, numCellsPerPatchInner=22),))]

This method returns the xy indices of the patch. Extract the PatchInfo object from the list that is returned, and use the getSequentialIndex to get the patch number.

In [11]:

patchinfo = closest_tract_patch[0][1][0]
print(f"patch ID: {patchinfo.getSequentialIndex()}")

patch ID: 0

Find the tract/patch overlapping a coordinate with a different method. This time, use skyMap.findTract, then findPatch from the results of the tract query; this will return only the tract and patch IDs.

In [12]:

my_tract = skymap.findTract(coord)
my_patch = my_tract.findPatch(coord)
print(my_tract.tract_id, my_patch.getSequentialIndex())

6324 0

2.1.2. All tracts overlapping coordinate¶

Find all tracts overlapping the coordinate specified in Section 1. This method expects the input to be a SpherePoint object.

In [13]:

alltracts = skymap.findAllTracts(coord)
alltracts

Out[13]:

[TractInfo(id=6096, ctrCoord=[0.3214247684691223, -0.8726543995830115, -0.36764169663272]),
 TractInfo(id=6324, ctrCoord=[0.30976315626767525, -0.8866455017580095, -0.34337521930298853]),
 TractInfo(id=6325, ctrCoord=[0.3339708415166301, -0.8778137249925514, -0.34337521930298853])]

Note that the result from findClosestTractPatchList is one of the results from findAllTracts, but that "alltracts" contains multiple entries. This is because tracts and patches overlap at their edges.

2.1.3. Tracts containing a list of coordinates¶

Create arrays of RA, Dec coordinates and search for tracts containing those coordinates using findTractIdArray. Specify degrees=True to tell the method that the input coordinates are in degrees (if not specified, they are assumed to be in radians).

In [14]:

ra_array = np.array([ra, ra+5.0, ra+12.0])
dec_array = np.array([dec, dec+2.0, dec+4.0])

tract_ids = skymap.findTractIdArray(ra_array, dec_array, degrees=True)

print(tract_ids)

[6324 6558 6795]

The findTractIdArray returns a single result for each input coordinate. If the coordinate lies in a region of overlap between tracts, the id of the nearest tract will be returned.

2.1.4. Tracts/patches within a region¶

Create a rectangular region based on four corner coordinates. Convert these to a list of SpherePoint objects, then find the tract/patch info for all patches overlapping the region.

In [15]:

width = 0.5*geom.degrees

coords = []
coeffs = [(0.5, 0.5), (0.5, -0.5), (-0.5, -0.5), (-0.5, 0.5)]
for coeff_ra, coeff_dec in coeffs:
    coords.append(geom.SpherePoint(ra*geom.degrees+(width*coeff_ra),
                                   dec*geom.degrees+(width*coeff_dec)))

print(coords)

tractpatch_in_region = skymap.findTractPatchList(coords)

[SpherePoint(290.25*degrees, -20.55*degrees), SpherePoint(290.25*degrees, -21.05*degrees), SpherePoint(289.75000000000006*degrees, -21.05*degrees), SpherePoint(289.75000000000006*degrees, -20.55*degrees)]

In [16]:

for i in range(len(tractpatch_in_region)):
    npatches = len(tractpatch_in_region[i][1])
    print(f"{npatches} patches in tract {tractpatch_in_region[i][0].getId()}.")

6 patches in tract 6325.
8 patches in tract 6096.
9 patches in tract 6324.

The region encompasses a total of 23 patches distributed across 3 tracts.

2.1.5. RA, Dec range of tract¶

Print the min/max RA, Dec of tract 6325.

In [17]:

skymap.getRaDecRange(6325)

Out[17]:

(Angle(290.04366812227079, degrees),
 Angle(291.61572052401743, degrees),
 Angle(-20.826446280991735, degrees),
 Angle(-19.338842975206607, degrees))

Print the vertices (corners) of the tract. These will appear as SpherePoint objects that contain RA, Dec pairs.

In [18]:

my_tract = skymap.generateTract(6325)
my_tract.vertex_list

Out[18]:

[SpherePoint(291.72161721625366*degrees, -20.91357632950237*degrees),
 SpherePoint(289.93771197612944*degrees, -20.913576020859665*degrees),
 SpherePoint(289.9471465352892*degrees, -19.247226759672763*degrees),
 SpherePoint(291.71218328574196*degrees, -19.247227041683743*degrees)]

3. Plotting tract/patch boundaries¶

3.1. Plot all DP2 tracts¶

Query the CoaddPatches table from the TAP service to retrieve information about all DP2 tracts and patches. Then extract a list of unique tracts.

In [19]:

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

In [20]:

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

Number of tracts:  2191

Plot the boundaries of all DP2 tracts on a map of the sky using the skyproj package. Generate the tract coordinates using the skyMap.generateTract method, then extract the vertices and convert them to RA, Dec coordinates.

In [21]:

fig, ax = plt.subplots(figsize=(10, 6))
sp = skyproj.McBrydeSkyproj(ax=ax)
for tract in tracts:
    info = skymap.generateTract(tract)
    vert_list = info.getVertexList()
    ras, decs = zip(*[[vert.getRa().asDegrees(), vert.getDec().asDegrees()] for vert in vert_list])
    sp.draw_polygon(ras, decs, edgecolor='darkgrey', alpha=1, linewidth=0.5, facecolor=None)
sp.ax.set_xlabel("Right Ascension", fontsize=14)
sp.ax.set_ylabel("Declination", fontsize=14)
plt.show()

Figure 1: All-sky spatial plot (in equatorial RA, Dec coordinates) showing the 2191 tracts with data in DP2 as small gray squares.

3.2. Plot patches in a small sky area¶

Query the butler for all r-band deep_coadd images within a 1.5 degree radius of a given position. Extract the tract and patch IDs for these from the returned DatasetRefs.

In [22]:

radius = 1.5
region = sphgeom.Region.from_ivoa_pos(f"CIRCLE {ra} {dec} {radius}")
query = "patch.region OVERLAPS :region AND band=:band"
bind_params = {'region': region, 'band': 'r'}
coadd_refs = butler.query_datasets("deep_coadd",
                                   where=query,
                                   bind=bind_params)

tractpatch_list = [(ref.dataId['tract'], ref.dataId['patch']) for ref in coadd_refs]

Use generateTract, getPatchInfo, and the WCS to extract the corner coordinates of each patch in RA, Dec. Assign each unique tract its own color, and label the patches.

In [23]:

tracts_selected = np.unique([tract for tract, patch in tractpatch_list])
tract_colors = {}
for i in range(len(tracts_selected)):
    tract_colors[tracts_selected[i]] = f"C{i}"

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
for tract, patch in tractpatch_list:
    tmp_tract = skymap.generateTract(tract)
    tmp_patch = tmp_tract.getPatchInfo(patch)
    ras = []
    decs = []
    for corn in tmp_patch.outer_bbox.getCorners():
        ras.append(tmp_patch.wcs.pixelToSky(corn.x, corn.y).getRa().asDegrees())
        decs.append(tmp_patch.wcs.pixelToSky(corn.x, corn.y).getDec().asDegrees())
    ras.append(ras[0])
    decs.append(decs[0])

    color = tract_colors[tract]
    plt.plot(ras, decs, color=color, linewidth=0.5, alpha=1, label='__none__')
    plt.text(np.max(ras)-0.05, np.min(decs)+0.05, f"{patch}", fontsize='xx-small', color='Gray')

for tract in tract_colors.keys():
    plt.plot(0.0, 0.0, color=tract_colors[tract], linewidth=1, alpha=1, label=f"{tract}")
ax.set_xlabel("Right Ascension", fontsize=12)
ax.set_ylabel("Declination", fontsize=12)
plt.legend(ncols=5)
ax.set_xlim([ra-1.3*radius, ra+1.3*radius])
ax.set_ylim([dec-1.3*radius, dec+1.3*radius])
ax.invert_xaxis()
plt.tight_layout()
plt.show()

Figure 2: Spatial plot (in equatorial RA, Dec coordinates) showing the boundaries of all patches within the 1.5-degree search radius. Each of the 11 unique tracts is assigned a different color, and patch numbers are printed within each patch.

Note that one could also plot the patch boundaries using the s_region from the CoaddPatches table.

4. Cells¶

DP2 coadds are "cell-based coadds" (see the 200-level tutorial on "Deep coadd images"). Within each patch, the cell geometry can be accessed in similar ways to the patch info within a tract.

Generate tract/patch info, then explore the cell information.

In [24]:

my_tract = skymap.generateTract(6325)
my_patch = my_tract.getPatchInfo(37)

Show that the patch contains a grid of 22x22 cells.

In [25]:

my_patch.num_cells

Out[25]:

Index2D(x=22, y=22)

Show that the cells are 150 pixels on a side.

In [26]:

my_patch.cell_inner_dimensions

Out[26]:

Extent2I(150, 150)

Get the information about the cell with index 37.

In [27]:

cell37info = my_patch.getCellInfo(37)
cell37info

Out[27]:

CellInfo(index=Index2D(x=15, y=1), innerBBox=(minimum=(23100, 9000), maximum=(23249, 9149)), outerBBox=(minimum=(23100, 9000), maximum=(23249, 9149)))

Cells in DP2 do not overlap, so the innerBBox and outerBBox are the same.

Retrieve the cell "sequential index" (37) using the x, y values from the CellInfo in two different ways.

In [28]:

my_patch.getSequentialCellIndexFromPair(cell37info.index)

Out[28]:

37

In [29]:

my_patch.getSequentialCellIndexFromPair((15, 1))

Out[29]:

37

Extract the coordinates of the cell boundary (i.e., the corners of the inner bounding box) in pixel coordinates.

In [30]:

xy_corners_cell37 = cell37info.inner_bbox.getCorners()
print(xy_corners_cell37)

[Point2I(23100, 9000), Point2I(23249, 9000), Point2I(23249, 9149), Point2I(23100, 9149)]

Use the wcs associated with the cell to convert these corners to RA, Dec coordinates.

In [31]:

radec_corners_cell37 = [cell37info.wcs.pixelToSky(corn.x, corn.y) for corn in xy_corners_cell37]
radec_corners_cell37

Out[31]:

[SpherePoint(290.3494931331977*degrees, -20.415260784222344*degrees),
 SpherePoint(290.3406613181411*degrees, -20.415236363027145*degrees),
 SpherePoint(290.34068720154846*degrees, -20.406959158780108*degrees),
 SpherePoint(290.34951854920126*degrees, -20.406983569224508*degrees)]

In [ ]: