203.1. Survey property maps

203.1. Survey property maps¶

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

Learning objective: To understand how to access and visualize the survey property maps.

LSST data products: Survey property maps.

Packages: skyproj, lsst.daf.butler

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 will teach the user how to load and visualize survey property maps generated by the Rubin Science Pipelines.

Survey property maps are all-sky metrics such as total depth or mean image quality, aggregated from all the images.

The maps are available in the "healsparse" format. The healsparse package, which is not used in this notebook, provides utilities for reading and manipulating sparse healpix maps. More information can be found in the documentation "HealSparse: A sparse implementation of HEALPix".

The map names, descriptions (and units if not dimensionless) are:

• deepCoadd_exposure_time_consolidated_map_sum: Total exposure time accumulated per sky position (second)
• deepCoadd_epoch_consolidated_map_min, ...max, ...mean: Earliest, latest, and mean observation epochs (MJD)
• deepCoadd_psf_size_consolidated_map_weighted_mean: Weighted mean of PSF characteristic width as computed from the determinant radius (pixel)
• deepCoadd_psf_e1_consolidated_map_weighted_mean: Weighted mean of PSF ellipticity component e1
• deepCoadd_psf_e2_consolidated_map_weighted_mean: Weighted mean of PSF ellipticity component e2
• deepCoadd_psf_maglim_consolidated_map_weighted_mean: Weighted mean of PSF flux 5σ magnitude limit (magAB)
• deepCoadd_sky_background_consolidated_map_weighted_mean: Weighted mean of background light level from the sky (nJy)
• deepCoadd_sky_noise_consolidated_map_weighted_mean: Weighted mean of standard deviation of the sky level (nJy)
• deepCoadd_dcr_dra_consolidated_map_weighted_mean: Weighted mean of DCR-induced astrometric shift in right ascension direction, expressed as a proportionality factor
• deepCoadd_dcr_ddec_consolidated_map_weighted_mean: Weighted mean of DCR-induced astrometric shift in declination direction, expressed as a proportionality factor
• deepCoadd_dcr_e1_consolidated_map_weighted_mean: Weighted mean of DCR-induced change in PSF ellipticity (e1), expressed as a proportionality factor
• deepCoadd_dcr_e2_consolidated_map_weighted_mean: Weighted mean of DCR-induced change in PSF ellipticity (e2), expressed as a proportionality factor

Related tutorials: The 100-level tutorials on how to use the Butler and display images.

1.1. Import packages¶

Import general python packages from numpy, matplotlib, and astropy, and the LSST Science Pipelines package for the Butler. Also import the skyproj package, which provides utilities for visualizing both sparse and dense HEALPix maps, as described in the documentation "SkyProj: Sky Projections with matplotlib and PROJ".

In [1]:

import hpgeom as hpg
import numpy as np
import matplotlib.pyplot as plt
from astropy.visualization import ZScaleInterval, LinearStretch, ImageNormalize
from astropy.table import Table
import skyproj
from lsst.daf.butler import Butler

1.2. Define parameters and functions¶

Instantiate the Butler.

In [2]:

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

Define the central coordinates of a field to use throughout.

In [3]:

ra_cen = 300.0
dec_cen = -15.0

In [4]:

hpg.nside_to_resolution(32768, units='arcseconds')

Out[4]:

np.float64(6.441537022157636)

2. Data access¶

Healsparse maps are stored in two layers, a "coverage" map that is a low-resolution layer of the entire sky, and a higher-resolution map containing the data. The "resolution" of a HEALPix tesselation of the sky is set by the nside parameter (where nside is a power of 2). Thus two layers of the healsparse maps have their resolutions defined by a nside_coverage for the coverage layer and nside_sparse for the data layer.

2.1. About maps and nside¶

By default, the DP2 survey property maps have nside_coverage=32 and nside_sparse=32768. To gain some context, use the hpgeom package (already imported) to calculate the resolution in arcseconds, arcminutes, and degrees of different nside values from 5 to 15.

In [5]:

tab_list = []
exponent_range = np.arange(5, 16, 1, dtype=int)

for exp in exponent_range:
    tmpdict = {}
    nside = 2**exp
    res_arcsec = hpg.nside_to_resolution(nside, units='arcseconds')
    res_arcmin = hpg.nside_to_resolution(nside, units='arcminutes')
    res_deg = hpg.nside_to_resolution(nside, units='degrees')
    tmpdict['exponent'] = exp
    tmpdict['nside'] = nside
    tmpdict['res_arcsec'] = res_arcsec
    tmpdict['res_arcmin'] = res_arcmin
    tmpdict['res_deg'] = res_deg
    tab_list.append(tmpdict)

nside_table = Table(rows=tab_list)
nside_table['res_arcsec'].info.format = '11.2f'
nside_table['res_arcmin'].info.format = '9.3f'
nside_table['res_deg'].info.format = '11.5f'

nside_table

Out[5]:

Table length=11

exponent	nside	res_arcsec	res_arcmin	res_deg
int64	int64	float64	float64	float64
5	32	6596.13	109.936	1.83226
6	64	3298.07	54.968	0.91613
7	128	1649.03	27.484	0.45806
8	256	824.52	13.742	0.22903
9	512	412.26	6.871	0.11452
10	1024	206.13	3.435	0.05726
11	2048	103.06	1.718	0.02863
12	4096	51.53	0.859	0.01431
13	8192	25.77	0.429	0.00716
14	16384	12.88	0.215	0.00358
15	32768	6.44	0.107	0.00179

Warning: The high-resolution "sparse" data layer is stored at ~6 arcsecond resolution. Obviously, this will lead to extremely large file sizes if reading the file at its full resolution. It is thus recommended to only read the full resolution (nside=32768) if that is needed for your science case. For all-sky maps where such fine-grained information is not needed, it is recommended to use the degrade_nside parameter to retrieve maps at lower resolution.

Survey property maps are available via the Butler.

2.2. List all maps¶

Use the Butler to print a list of the available maps.

In [6]:

for dtype in sorted(butler.registry.queryDatasetTypes(expression="*consolidated_map*")):
    print(dtype.name)

deepCoadd_dcr_ddec_consolidated_map_weighted_mean
deepCoadd_dcr_dra_consolidated_map_weighted_mean
deepCoadd_dcr_e1_consolidated_map_weighted_mean
deepCoadd_dcr_e2_consolidated_map_weighted_mean
deepCoadd_epoch_consolidated_map_max
deepCoadd_epoch_consolidated_map_mean
deepCoadd_epoch_consolidated_map_min
deepCoadd_exposure_time_consolidated_map_sum
deepCoadd_psf_e1_consolidated_map_weighted_mean
deepCoadd_psf_e2_consolidated_map_weighted_mean
deepCoadd_psf_maglim_consolidated_map_weighted_mean
deepCoadd_psf_size_consolidated_map_weighted_mean
deepCoadd_sky_background_consolidated_map_weighted_mean
deepCoadd_sky_noise_consolidated_map_weighted_mean

Each of these products represents a healsparse map containing the value of an individual survey property, as described in Section 1.

Show the dataset type for one of the maps.

In [7]:

butler.get_dataset_type("deepCoadd_dcr_ddec_consolidated_map_weighted_mean")

Out[7]:

DatasetType('deepCoadd_dcr_ddec_consolidated_map_weighted_mean', {band, skymap}, HealSparseMap)

Show that the required dimensions for retrieving this map are the skymap and the band.

In [8]:

butler.get_dataset_type("deepCoadd_dcr_ddec_consolidated_map_weighted_mean").dimensions.required

Out[8]:

{band, skymap}

2.2. Retrieve a map¶

Retrieve a map from the Butler by specifying the map name and band: select the PSF magnitude limits in r-band.

Because many of the pixels tiling the sky are empty, first query the .coverage pixels (a lower-resolution map), then retrieve only the pixels that have coverage from the full map.

In [9]:

hspmap_cov = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean.coverage',
                        band='r', skymap='lsst_cells_v2')

Extract an array of only the pixels with coverage in this map.

In [10]:

cov_pixels, = np.where(hspmap_cov.coverage_mask)

Make an RA, Dec spatial plot showing the positions of these pixels, to confirm that it extracted only the DP2 regions. Use the hpgeom utilities to extract the RA, Dec values of the pixels.

In [11]:

fig = plt.figure(figsize=(9, 5))

ra, dec = hpg.pixel_to_angle(hspmap_cov.nside_coverage, cov_pixels)
plt.plot(ra, dec, 'k.', alpha=0.5)
plt.gca().invert_xaxis()
plt.xlabel('RA (deg)')
plt.ylabel('Dec (deg)')
plt.show()

Warning: Survey property maps are very large files, and can easily cause your session to run out of memory. Unless you absolutely need the full resolution, you should read the maps at degraded resolution.

Read the full set of pixels containing data at a degraded resolution of ~50 arcseconds (see Sec. 2.1) by passing the 'degrade_nside':4096 parameter when reading from the Butler.

In [12]:

hspmap = butler.get('deepCoadd_psf_maglim_consolidated_map_weighted_mean',
                    band='r', skymap='lsst_cells_v2',
                    parameters={'pixels':list(cov_pixels), 'degrade_nside':4096})

3. Manipulate the map¶

To conserve memory, HealSparse uses a dual-map approach, where a low-resolution full-sky “coverage map” is combined with a high resolution map containing the pixel data where it is available.

Show that the map has been provided as a HealSparseMap object with the requested resolution.

In [13]:

hspmap

Out[13]:

HealSparseMap: nside_coverage = 32, nside_sparse = 4096, float64, 17919510 valid pixels

3.1. Get map value at position¶

Each pixel of the healsparse map corresponds to a small region of the sky. The value of the map corresponds to the value of the survey property at that location.

To access the survey property value at a specific location or set of locations, query for the map value using the get_values_pos functionality.

Print the map value at the center of the selected field. Recall that the map is providing the mean PSF magnitude limit in $r$-band at this position.

In [14]:

print(hspmap.get_values_pos(ra_cen, dec_cen))

24.760509146978322

Extract map values at an array of locations over a small area near the center of the selected field.

In [15]:

span_dec = 0.25
span_ra = span_dec / np.cos(np.deg2rad(dec_cen))
ra = np.linspace(ra_cen-span_ra, ra_cen+span_ra, 3)
dec = np.linspace(dec_cen-span_dec, dec_cen+span_dec, 3)
print('RA: ', ra)
print('Dec: ', dec)
for d in dec:
    print(hspmap.get_values_pos(ra, d))
del ra, dec

RA:  [299.74118095 300.         300.25881905]
Dec:  [-15.25 -15.   -14.75]
[24.83969852 24.75119865 24.75570704]
[24.6544644  24.76050915 24.73502447]
[24.75652985 24.60930114 24.75569286]

A sentinel (placeholder) value is returned if the position is outside the region where the map has data.

Request the map value at a northern declination, and see the sentinel value.

In [16]:

print(hspmap.get_values_pos(180, +30))

-1.6375e+30

3.2. Get map statistics in region¶

Extract map values over a 10-degree region centered on the selected position.

In [17]:

span_dec = 10.0
span_ra = span_dec / np.cos(np.deg2rad(dec_cen))
ra = np.linspace(ra_cen-span_ra, ra_cen+span_ra, 250)
dec = np.linspace(dec_cen-span_dec, dec_cen+span_dec, 250)
x, y = np.meshgrid(ra, dec)
values = hspmap.get_values_pos(x, y)

Change map values with the sentinel value from the placeholder to "Not a Number" (NaN), so they don't contribute to the statistical measures.

In [18]:

row, col = np.where(values < 0.0)
values[row, col] = 'NaN'

Print a few statistics on the map values.

In [19]:

print(f"min, max limiting mags: {np.nanmin(values):.2F}, {np.nanmax(values):.2F}")
print(f"mean, median, std of limiting mags: {np.nanmean(values):.2F}, {np.nanmedian(values):.2F}, {np.nanstd(values):.2F}")

min, max limiting mags: 22.54, 26.44
mean, median, std of limiting mags: 24.66, 24.68, 0.35

4. Visualize the map¶

Plot the map values in the ECDFS region with the pcolormesh functionality and include a colorbar. Invert the x-axis so that the plot is east-left.

In [20]:

fig = plt.figure(figsize=(6, 5))
plt.pcolormesh(x, y, values)
plt.xlabel("Right Ascension (deg)")
plt.ylabel("Declination (deg)")
plt.colorbar(label="PSF Mag Limit (r-band)")
plt.gca().invert_xaxis()
plt.show()

Figure 1: The survey property map for PSF magnitude limit in the r-band for a field centered at (RA, Dec) = (300, -15), 20 degrees on a side.

4.1. Scale and color¶

Set the vmin and vmax values to produce a map that clearly delineates the region where the r-band image depth is greater than 25 magnitudes.

In [21]:

fig = plt.figure(figsize=(6, 5))
plt.pcolormesh(x, y, values, vmin=24.9, vmax=25.0)
plt.xlabel("Right Ascension (deg)")
plt.ylabel("Declination (deg)")
plt.colorbar(label="PSF Mag Limit (r-band)")
plt.gca().invert_xaxis()
plt.show()

Figure 2: Same data as Figure 1, but with the minimum and maximum values of the color bar set at 24.9 and 25.0 magnitudes.

Create a black-and-white plot using the "Greys" colormap (cmap).

In [22]:

fig = plt.figure(figsize=(6, 5))
plt.pcolormesh(x, y, values, cmap="Greys")
plt.xlabel("Right Ascension (deg)")
plt.ylabel("Declination (deg)")
plt.colorbar(label="PSF Mag Limit (r-band)")
plt.gca().invert_xaxis()
plt.show()

Figure 3: The same data as in Figure 1, but displayed with a greyscale colormap.

4.2. Sky projections¶

The skyproj package provides more advanced plotting capabilities.

Display the all-sky map with a grid.

In [23]:

fig, ax = plt.subplots(figsize=(9, 6))
sp = skyproj.Skyproj(ax=ax)
sp.draw_hspmap(hspmap)

cbar_ax = fig.add_axes([0.98, 0.3, 0.02, 0.4])
sp.draw_colorbar(cax=cbar_ax, label='PSF maglim (r-band)', orientation="vertical", location="right")

plt.show()
del fig, ax, sp

Figure 4: The all-sky r-band PSF depth map for Data Preview 2 shows the ~3000 square degrees of sky coverage, including deep drilling fields.

In Figure 4, notice how the deep drilling fields at lower declinations are "stretched" out with this linear grid. This grid is not a good representation of the curved sky.

Display the all-sky map with the McBrydeSkyproj. This creates a visualization using the McBryde-Thomas Flat Polar Quartic projection, centered at 0 deg RA.

In [24]:

fig, ax = plt.subplots(figsize=(12, 3))
sp = skyproj.McBrydeSkyproj(ax=ax, lon_0=0.0)
sp.draw_hspmap(hspmap)
cbar_ax = fig.add_axes([0.3, 0.98, 0.4, 0.04])
sp.draw_colorbar(cax=cbar_ax, label='PSF maglim (r-band)', orientation="horizontal", location="top")

plt.show()
del fig, ax, sp

Figure 5: The same data as Figure 4, but with a McBryde sky projection that is more appropriate for the curved plane of the sky.

5. Co-plot maps and deep coadd¶

Build a four-panel plot that shows one deep coadd image, side-by-side with three survey property maps covering the same region.

Find a deep_coadd image that overlaps the central RA and Dec of the field defined in Section 1.2.

In [25]:

query = f"band.name = 'r' AND patch.region OVERLAPS POINT({ra_cen}, {dec_cen})"
print(query)
refs = butler.query_datasets("deep_coadd",
                             where=query,
                             order_by=["patch.tract"])
print(len(refs))

band.name = 'r' AND patch.region OVERLAPS POINT(300.0, -15.0)

2

Retrieve the first image from the results.

In [26]:

deep_coadd = butler.get(refs[0])

Get the WCS and bounding box, and define the extent and corners of the image.

In [27]:

wcs = deep_coadd.sky_projection
bbox = deep_coadd.bbox

In [28]:

extent = (bbox.min.x, bbox.max.x,
          bbox.min.y, bbox.max.y)
corners = [wcs.pixel_to_sky(x=bbox.min.x, y=bbox.min.y),
           wcs.pixel_to_sky(x=bbox.min.x, y=bbox.min.y),
           wcs.pixel_to_sky(x=bbox.max.x, y=bbox.max.y),
           wcs.pixel_to_sky(x=bbox.max.x, y=bbox.min.y)]

Retrieve two more survey property maps: the weighted mean of the PSF size, and the weighted mean sky background, both in r-band.

Warning: Do not try to retrieve both of these at full resolution. That will cause your session to crash because it will run out of memory. Instead, retrieve a reduced-resolution version of the map with nside=4096.

Note that the cell reading two maps is likely to take a long time to run (perhaps a few minutes).

In [29]:

hspmap_psf = butler.get('deepCoadd_psf_size_consolidated_map_weighted_mean',
                        band='r', skymap='lsst_cells_v2',
                        parameters={'pixels':list(cov_pixels), 'degrade_nside':4096})
hspmap_sky = butler.get('deepCoadd_sky_background_consolidated_map_weighted_mean',
                        band='r', skymap='lsst_cells_v2',
                        parameters={'pixels':list(cov_pixels), 'degrade_nside':4096})

Make a new grid that matches the boundaries of the deep_coadd image.

In [30]:

lower = wcs.pixel_to_sky(x=bbox.start.x, y=bbox.start.y)
upper = wcs.pixel_to_sky(x=bbox.stop.x, y=bbox.stop.y)
ra = np.linspace(lower.ra.deg, upper.ra.deg, 150)
dec = np.linspace(lower.dec.deg, upper.dec.deg, 150)
x, y = np.meshgrid(ra, dec)

It's challenging to overlay a projected RA, Dec grid in a subplot. Instead, prepare to substitute pixels for sky coordinates as tick labels for the image.

In [31]:

xtick_locs = np.linspace(extent[0], extent[1], 5)
ytick_locs = np.linspace(extent[2], extent[3], 5)
xtick_lbls = []
ytick_lbls = []
for xt, yt in zip(xtick_locs, ytick_locs):
    temp = wcs.pixel_to_sky(x=xt, y=yt)
    xtick_lbls.append(str(np.round(temp.ra.deg, 2)))
    ytick_lbls.append(str(np.round(temp.dec.deg, 2)))

Create the four-panel figure.

In [32]:

fig, ax = plt.subplots(2, 2, figsize=(9, 6))

temp = deep_coadd.image.array.flatten()
norm = ImageNormalize(temp, interval=ZScaleInterval(),
                      stretch=LinearStretch())
im = ax[0, 0].imshow(deep_coadd.image.array, cmap='gray', norm=norm,
                     extent=extent, origin='lower')
fig.colorbar(im, ax=ax[0, 0], label="Pixel Value (r-band)")
ax[0, 0].set_xticks(xtick_locs, xtick_lbls)
ax[0, 0].set_yticks(ytick_locs, ytick_lbls)
ax[0, 0].axis('tight')

values = hspmap.get_values_pos(x, y)
pcm = ax[0, 1].pcolormesh(x, y, values, cmap='viridis')
fig.colorbar(pcm, ax=ax[0, 1], label="PSF Mag Limit (r-band)")
ax[0, 1].axis('tight')
ax[0, 1].invert_xaxis()
del pcm, values

values = hspmap_psf.get_values_pos(x, y)
pcm = ax[1, 0].pcolormesh(x, y, values, cmap='plasma')
fig.colorbar(pcm, ax=ax[1, 0], label="PSF Size (r-band)")
ax[1, 0].axis('tight')
ax[1, 0].invert_xaxis()
del pcm, values

values = hspmap_sky.get_values_pos(x, y)
pcm = ax[1, 1].pcolormesh(x, y, values, cmap='cividis')
fig.colorbar(pcm, ax=ax[1, 1], label="Sky Background (r-band)")
ax[1, 1].axis('tight')
ax[1, 1].invert_xaxis()
del pcm, values

plt.tight_layout()

del fig, ax, temp, norm

Figure 6: Four different views on the field centered at (RA, Dec) = (300, -15): the $r$-band deep coadd image pixel data (grayscale, upper left); the $r$-band PSF magnitude limit (viridis, upper right); the $r$-band PSF size (plasma, lower left); and the $r$-band sky background (cividis, lower right).

In [ ]: