310.1. Photo-z access and visualization#

310_1_Photoz_access

310.1. Photo-z access and visualization¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: v30.0.10
Last verified to run: 2026-07-31
Repository: github.com/lsst/tutorial-notebooks

Learning objective: How to access, cross-match, and visualize the LSDB-formatted photo-z.

LSST data products: Object

Packages: lsdb

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¶

Members of the Rubin Commissioning Science Unit for photometric redshifts have generated photometric redshift (photo-z; PZ) estimates for Data Preview 2 objects using a variety of PZ estimators. These PZ are being made available as LSDB-formatted files, and documented in "Photometric Redshifts for Data Preview 2" RTN-124.

LSDB (Large Scale Database) is a python tool for scalable analysis of large catalogs (query and cross-match). Built on top of Dask, LSDB uses the HATS (Hierarchical Adaptive Tiling Scheme) data format to efficiently perform spatial operations.

Algorithms:

bpz: Bayesian Photometric Redshifts (BPZ; Benítez 2000)
dnf: Directional Neighbourhood Fitting (DNF; De Vicente et al. 2016)
fzboost: FlexZBoost (Izbicki & Lee 2017)
gpz: Gaussian processes for photometric redshifts (GPz; Almosallam et al. 2016)
knn: k-Nearest Neighbors (kNN; The RAIL Team et al. 2025)
tpz: Trees for Photo-Z (Carrasco Kind & Brunner 2013)

Measurements: For each algorithm these tables carry four point estimates, z_best, z_mode, z_mean, and z_median plus $1\sigma$ and $2\sigma$ credible intervals (z_err68_low/high, z_err95_low/high).

Caveats: These estimates and their error bars are only very lightly validated. Treat them as provisional and use the results with caution. They are not necessarily representative of future PZ because the DP2 depths and filter distributions are not the same as planned for the LSST, and input spectroscopic training sets (and other PZ estimator inputs) will be improved and refined in the future.

A citable publication for the DP2 PZ is forthcoming.

This tutorial only demonstrates how to access and visualize the PZ estimates with LSDB. It does not provide instruction regarding the different PZ estimators or demonstrate how to generate PZ estimates.

Related tutorials: The 100-level tutorial on how to access LSDB-formatted files. The 100-level tutorials on the TAP (Table Access Protocol) service, the Butler, and displaying images with Firefly.

1.1. Import packages¶

Import the LSDB package to work with LSDB-formatted files, along with standard astronomy packages and LSST software for data access.

In [1]:

import lsdb
from upath import UPath
import matplotlib.pyplot as plt
from bokeh.io import output_notebook
import holoviews as hv
from holoviews.operation.datashader import datashade, dynspread
from lsst.rsp import RSPDiscovery

1.2. Define parameters and functions¶

Create an instance of the TAP service.

In [2]:

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

Set up for interactive plotting with holoviews. Order is important here: the output_notebook() command must be run after the hv.extension('bokeh') command or plots will not display.

In [3]:

hv.extension('bokeh')

ⓘ

In [4]:

output_notebook()

Loading BokehJS ...

Set the base path to the LSDB-formatted DP2 data on the RSP.

In [5]:

base_path = UPath("/rubin/lsdb_data/dp2")

2. Explore the PZ catalog¶

Open the read-only LSDB PZ catalog.

In [6]:

pz_cat = lsdb.open_catalog(base_path / "object_photoz")

Display the results.

In [7]:

pz_cat

Out[7]:

lsdb Catalog object_photoz:

	objectId	ra	dec	bpz_z_best	dnf_z_best	fzboost_z_best	gpz_z_best	knn_z_best	tpz_z_best
npartitions=1356
Order: 3, Pixel: 265	int64[pyarrow]	double[pyarrow]	double[pyarrow]	float[pyarrow]	float[pyarrow]	float[pyarrow]	float[pyarrow]	float[pyarrow]	float[pyarrow]
...	...	...	...	...	...	...	...	...	...
Order: 4, Pixel: 3070	...	...	...	...	...	...	...	...	...
Order: 6, Pixel: 49136	...	...	...	...	...	...	...	...	...

9 out of 51 available columns in the catalog have been loaded lazily, meaning no data has been read, only the catalog schema

This catalog has an estimated size of 33.7 GB

Lazily loaded catalogs is always the default for LSDB catalogs, and it means that only the metadata is loaded at first.

2.1. Column names¶

The columns of the LSDB PZ catalog include the objectId, ra, dec, and standardized outputs from the PZ estimators in the format of <estimator>_z_mode, z_mean, z_median, z_best, z_err68_low, z_err68_high, z_err95_low, and z_err95_high.

Option to display all column names.

In [8]:

# pz_cat.all_columns

2.2. Sky partitions¶

Partitions are how the LSDB-formatted files are stored, with each partition typically having about the same number of objects.

Show the partitions of pz_cat on the sky. This is not necessarily the same as a sky coverage map, but rather shows the file's polygonal partition boundaries.

In [9]:

fig = pz_cat.plot_pixels(plot_title="Sky Partition Map")

No description has been provided for this image

Figure 1: The sky partitions of the lazily-loaded PZ catalog.

2.3. Table data¶

Show the first 10 rows of the PZ catalog.

In [10]:

pz_cat.head(n=10)

Computing Catalog:   0%|          | 0/1 [00:00<?, ?it/s]

Out[10]:

	objectId	ra	dec	bpz_z_best	dnf_z_best	fzboost_z_best	gpz_z_best	knn_z_best	tpz_z_best
_healpix_29
1196750063051683238	767490864585903255	351.158498	-17.850872	2.037524	2.082143	3.07092	0.537438	0.298913	0.660856
1196750064058168017	767490864585917450	351.154374	-17.850893	0.374382	2.086247	1.159724	0.50978	0.29908	0.722733
1196750064998526636	767490864585917508	351.15485	-17.849085	1.911994	2.082051	3.068399	0.532672	0.298944	0.615579
1196750065055727242	767490864585917509	351.155907	-17.849272	0.72311	2.138241	1.013317	0.509443	0.299221	1.033516
1196750081847571868	767490864585901490	351.135249	-17.849504	0.374776	2.081838	0.896646	0.52057	0.299029	0.702268
1196750082018620727	767490864585901486	351.132396	-17.849868	0.618678	2.081431	1.092273	0.524669	0.299074	0.682498
1196750086503664774	767490864585917522	351.143259	-17.848738	0.221739	2.086576	0.939482	0.509705	0.29904	0.435762
1196750086690209664	767490864585917523	351.144107	-17.848468	0.276589	2.08648	1.047891	0.513002	0.299032	1.025774
1196750087070297366	767490864585903256	351.151689	-17.850359	0.694074	2.081557	3.07025	0.52755	0.299037	0.721395
1196750088054841611	767490864585903265	351.152676	-17.847991	1.494436	2.082051	3.067286	0.534329	0.298981	0.645686

10 rows × 9 columns

2.4. Plot N(pz)¶

$N(pz)$ refers to a histogram of the number of objects, $N$, as a function of photo-z, PZ.

To visualize how the PZ results vary by estimator, for each of the following fields plot the histogram N(pz) for $0<pz<3$ for each PZ estimator. These fields were selected for having deep coadd images in each of the six filters $ugrizy$.

In [11]:

field_centers = {
    "DDF_ELAIS_S1": [9.5, -44.0],
    "DDF_ECDFS": [53.0, -28.1],
    "DDF_EDFS_b": [63.2, -47.8],
    "DDF_COSMOS": [150.1, 2.1],
    "New_Horizons": [289.4, -20.2],
    "Rubin_SV_225_-40": [225.0, -39.5],
    "Rubin_SV_320_-15": [320.2, -15.1],
}

Retrieve the names of the PZ estimators and store them in array all_pz_codes.

In [12]:

all_pz_codes = []
for col in pz_cat.columns:
    if col.find('z_best') >= 0:
        all_pz_codes.append(col.split('_')[0])
print(all_pz_codes)

['bpz', 'dnf', 'fzboost', 'gpz', 'knn', 'tpz']

Create the histograms.

In [13]:

for field in field_centers.keys():
    coords = field_centers[field]
    pz_cat_cone = pz_cat.cone_search(ra=coords[0], dec=coords[1],
                                     radius_arcsec=2.0 * 3600.0)
    pz_df_cone = pz_cat_cone.compute()
    fig, ax = plt.subplots(2, 3, figsize=(8, 4))
    ze = 0
    for i in range(2):
        for j in range(3):
            pz_code = all_pz_codes[ze]
            tx = (pz_df_cone[pz_code+'_z_best'] > 0) & \
                 (pz_df_cone[pz_code+'_z_best'] < 3)
            ax[i, j].hist(pz_df_cone[pz_code+'_z_best'][tx],
                          bins=100, histtype='step', color='black')
            ax[i, j].set_title(all_pz_codes[ze])
            ax[i, j].set_xlim([0, 3])
            ze += 1
            del pz_code, tx
    plt.suptitle(field)
    plt.tight_layout()
    del coords, pz_cat_cone, fig

Computing Catalog:   0%|          | 0/28 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/22 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/26 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/22 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/46 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/36 [00:00<?, ?it/s]

Computing Catalog:   0%|          | 0/18 [00:00<?, ?it/s]

Figure 2: Histograms of the best PZ for each field, and each PZ estimator.

2.5. Compare estimates¶

For a given field, retrieve two PZ columns and make a plot comparing them.

For this example, use the DDF Elais S1, and compare the BZP results for best and mean PZ.

In [14]:

field = 'DDF_ELAIS_S1'
coords = field_centers[field]
col1 = 'bpz_z_best'
col2 = 'bpz_z_mean'
pz_cat = lsdb.open_catalog(base_path / "object_photoz",
                           columns=[col1, col2])
pz_cat_cone = pz_cat.cone_search(ra=coords[0], dec=coords[1],
                                 radius_arcsec=2.0 * 3600.0)
pz_df_cone = pz_cat_cone.compute()

Computing Catalog:   0%|          | 0/28 [00:00<?, ?it/s]

Option to display the table.

In [15]:

# pz_df_cone

Set up to make a dynamic, interactive, 2D histogram to compare the two selected columns.

In [16]:

points = hv.Points((pz_df_cone[col1], pz_df_cone[col2]))
boundsxy = (0, 0, 0, 0)
box = hv.streams.BoundsXY(source=points, bounds=boundsxy)
bounds = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box])
p = dynspread(datashade(points, cmap="Viridis"))
p = p.opts(width=600, height=300, padding=0.05, show_grid=True,
           xlim=(0, 3), ylim=(0, 6), xlabel=col1, ylabel=col2,
           tools=['box_select'])

Display the plot.

In [17]:

p * bounds

Out[17]:

Figure 3: The BPZ photo-$z$ results of the mean vs. best, for galaxies within 2 degrees of the center of DDF Elais S1. Use the zoom-in function (dashed box with magnifying glass) and see that the histogram's resolution changes dyamically depending on how many points are in the frame.

Clean up.

In [18]:

del pz_cat

3. Join PZ and Object subsets¶

Small pandas dataframes ($<1$ million rows, $<1$ GB dataframe in-memory) can be read in to LSDB format and then joined to another LSDB table, as long as they both contain a common identifier column such as objectId.

Use the TAP service to obtain Objects within 0.7 degrees of the center of the COSMOS field (search radius set to keep the total <1 million).

In [19]:

field = 'DDF_COSMOS'
coords = field_centers[field]
search_radius = 0.7
query = """SELECT objectId, coord_ra, coord_dec, g_cModelMag, r_cModelMag
        FROM dp2.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
        CIRCLE('ICRS', {}, {}, {})) = 1 """.format(coords[0], coords[1], search_radius)
print(query)
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'
obj_cat_df = job.fetch_result().to_table().to_pandas()
print(len(obj_cat_df))

SELECT objectId, coord_ra, coord_dec, g_cModelMag, r_cModelMag
        FROM dp2.Object
        WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec),
        CIRCLE('ICRS', 150.1, 2.1, 0.7)) = 1
Job phase is COMPLETED
608145

Option to show obj_cat_df.

In [22]:

# obj_cat_df

The lsdb.from_dataframe method requires that columns named "ra" and "dec" exist, so rename the coordinate columns and then convert to an lsdb data frame. Delete redundant versions of the table.

In [23]:

obj_cat_rename = obj_cat_df.rename(columns={'coord_ra': 'ra', 'coord_dec': 'dec'})
obj_cat_lsdb = lsdb.from_dataframe(obj_cat_rename, catalog_name="object")
del obj_cat_rename, obj_cat_df

Option to display the first 10 rows of obj_cat_lsdb.

In [24]:

# obj_cat_lsdb.head(n=10)

Retrieve the objectId, ra, dec, and gpz_z_median columns from the LSDB PZ catalog, and query for rows in the same field as defined above.

In [25]:

pz_cat = lsdb.open_catalog(base_path / "object_photoz",
                           columns=['objectId', 'ra', 'dec', 'gpz_z_median'])
pz_cat_cone = pz_cat.cone_search(ra=coords[0], dec=coords[1],
                                 radius_arcsec=search_radius * 3600.0)

Use the join method to do an "inner" join on the shared column objectId. An "inner" join means the resulting catalog, joined_cat, only contains objects that appear in both obj_cat_lsdb and pz_cat_cone.

In [26]:

joined_cat = pz_cat_cone.join(
    obj_cat_lsdb,
    left_on='objectId',
    right_on='objectId',
    suffix_method='all_columns',
)

Option to display the first 10 rows of joined_cat.

In [27]:

joined_cat.head(n=10)

Computing Catalog:   0%|          | 0/4 [00:00<?, ?it/s]

Out[27]:

	objectId_object_photoz	ra_object_photoz	dec_object_photoz	gpz_z_median_object_photoz	objectId_object	ra_object	dec_object	g_cModelMag_object	r_cModelMag_object
_healpix_29
1917229402985056399	786346045893280050	150.193101	1.406708	0.824808	786346045893280050	150.193101	1.406708	25.8074	25.1101
1917229403002778086	786346045893280048	150.192226	1.40668	0.457578	786346045893280048	150.192226	1.40668	22.2682	20.9513
1917229403287839521	786346045893280109	150.194043	1.40817	1.761897	786346045893280109	150.194043	1.40817	25.691401	26.009199
1917229410617786451	786346045893279949	150.208491	1.409257	0.670406	786346045893279949	150.208491	1.409257	24.8484	24.1968
1917229410772948629	786346045893279957	150.207499	1.409857	1.568267	786346045893279957	150.207499	1.409857	25.5828	25.353701
1917229410853823956	786346045893252779	150.210302	1.409074	1.481173	786346045893252779	150.210302	1.409074	26.095699	26.0147
1917229411093865982	786346045893279960	150.21219	1.409931	2.192985	786346045893279960	150.21219	1.409931	26.6115	25.4974
1917229411213198760	786346045893279963	150.2114	1.410271	0.733254	786346045893279963	150.2114	1.410271	23.507799	22.775299
1917229411322710413	786346733088037028	150.211674	1.41131	1.29831	786346733088037028	150.211674	1.41131	24.824499	24.802099
1917229411394609214	786346045893279965	150.209246	1.410199	1.036921	786346045893279965	150.209246	1.410199	25.4676	25.134399

10 rows × 9 columns

As an example, plot the GPZ median redshift vs. the $g-r$ cModel color for all objects.

In [28]:

points = hv.Points((joined_cat['gpz_z_median_object_photoz'],
                    joined_cat['g_cModelMag_object'] - joined_cat['r_cModelMag_object']))
boundsxy = (0, 0, 0, 0)
box = hv.streams.BoundsXY(source=points, bounds=boundsxy)
bounds = hv.DynamicMap(lambda bounds: hv.Bounds(bounds), streams=[box])
p = dynspread(datashade(points, cmap="Viridis"))
p = p.opts(width=600, height=300, padding=0.05, show_grid=True,
           xlim=(0, 3), ylim=(-4, 6), xlabel='gpz_z_median', ylabel='g-r cModelMag [mag]',
           tools=['box_select'])

In [29]:

p * bounds

Out[29]:

Figure 4: A plot made by joining the PZ catalog with the Object table, showing the GPZ median redshift vs. the $g-r$ cModel magnitude color.

In [ ]: