I want to…

Load a catalog#

Find catalogs

data.lsdb.io is the catalog browser: Gaia, ZTF, SDSS, 2MASS, DES, and more, each with its row count, size on disk, and a copy-paste URL. Catalogs are mirrored in several regions — grab the one nearest you, and prefer an s3:// endpoint over https:// when you have one.

Open a catalog

open_catalog reads the catalog's metadata and nothing else. You get back a Catalog that knows its columns, types, and partitioning, but holds no rows.

import lsdb

ztf = lsdb.open_catalog(
    "https://data.lsdb.io/hats/ztf_dr14/ztf_object",
    columns=["ra", "dec", "ps1_objid", "mean_mag_r"],
)
ztf

Pick columns

Every catalog defines a set of default columns, loaded when you don't ask for specific ones. That default exists to stop you accidentally pulling two hundred columns across the network. Name what you need and the rest is never read.

# What am I getting by default, and what else is there?
ztf.columns
ztf.all_columns

# Ask for exactly what you need. The rest is never read.
ztf = lsdb.open_catalog(
    "https://data.lsdb.io/hats/ztf_dr14/ztf_object",
    columns=["ra", "dec", "mean_mag_g", "mean_mag_r"],
)

To keep the defaults and add to them, put ... in the list. LSDB expands the ellipsis in place to the catalog's default columns, skipping any you have already named, so nothing arrives twice. One ellipsis per list.

# The defaults, plus the two columns that say which partition a row lives in.
ztf = lsdb.open_catalog(
    "https://data.lsdb.io/hats/ztf_dr14/ztf_object",
    columns=[..., "Norder", "Npix"],
)

If you really do want every column, columns="all" says so explicitly — which is the point. Loading all two hundred should be a decision, not an accident.

Go faster#

Cut it down first

Write your whole workflow against one small region of sky. Check the results look the way you expected. Then delete the search and run the same code on the full catalog. This is the habit that saves the most time while you are still figuring out what you want, because a mistake costs you seconds instead of an hour.

A region search narrows partitions and rows, so you read a handful of files rather than thousands.

# While developing: one degree of sky.
gaia = gaia.cone_search(ra=180.0, dec=10.0, radius_arcsec=3600)

# When it all works, delete that line and the same code runs on the full catalog.

There are other shapes when a cone doesn't fit: box_search, polygon_search, moc_search, and order_search.

Ask for fewer columns

Covered up in pick columns, but it belongs here too. It is the cheapest speed-up to apply — one keyword argument, no restructuring — and unselected columns are never read at all.

Start a Dask client

LSDB runs on Dask. Without a client you still get parallelism across your cores; with one you get a dashboard showing exactly what is running, how much memory each worker holds, and which task is the bottleneck.

from dask.distributed import Client

client = Client(n_workers=4, threads_per_worker=1, memory_limit="auto")
client  # displays a link to the dashboard

In a script, a context manager guarantees the cluster shuts down even if the work fails. In a notebook you usually want the opposite — one client kept alive across many cells, closed with client.close() when you are done.

with Client(n_workers=4, threads_per_worker=1, memory_limit="auto") as client:
    result = matched.compute()

Read the nearest copy

Catalogs are mirrored across regions and providers. s3:// endpoints generally outperform https://, and if you are working on a cluster that already holds a local copy, use that instead of pulling across the internet at all.

Examine the catalog#

Look at a few rows

head() is the cheap way to see real values: it reads from the first partitions that have data rather than computing everything.

gaia.head()          # a handful of rows
gaia.head(20)

Find out how big it is

Before you compute anything, ask what it would cost. estimate_size() prints a report based on metadata alone — no row data is read.

gaia.npartitions     # how many files this workflow will touch
gaia.estimate_size()

That is the number worth checking before you reach for compute().

Check range of values in a column

HATS stores per-column minimum, maximum, and null counts in the parquet footers, so LSDB can report them across the whole catalog without touching a single row of data. This returns immediately even on a billion-row catalog.

gaia.aggregate_column_statistics()      # min, max, null_count, row_count
gaia.per_partition_statistics()         # the same, broken down by partition

These come from the catalog's own recorded metadata, so on a catalog you have already filtered they describe what was there originally — LSDB warns you when that happens.

Use one partition

One partition is a normal, laptop-sized dataframe. It is equivalent to a single Healpix pixel. It is the fastest way to try something out on real data before committing to the whole sky.

one = gaia.partitions[0].compute()      # by position

gaia.get_healpix_pixels()[:3]           # or address a pixel directly
gaia.get_partition(2, 108).compute()    # healpix order 2, pixel 108

See where the data is

gaia.plot_pixels()      # the partitioning itself
gaia.plot_coverage()    # which sky this catalog covers
gaia.plot_points()      # the objects (cut the catalog down first)

Modify the catalog#

Select a region of sky

Region searches are how you cut a catalog down spatially. Each runs in two stages: a coarse filter that decides which partitions to open at all, then a fine filter on the rows inside them. Pass fine=False to keep whole partitions and skip the second stage.

gaia.cone_search(ra=180.0, dec=10.0, radius_arcsec=3600)
gaia.box_search(ra=(179.0, 181.0), dec=(9.0, 11.0))
gaia.polygon_search([(179.0, 9.0), (181.0, 9.0), (181.0, 11.0), (179.0, 11.0)])

For anything that isn't one simple shape, build a MOC — a multi-order coverage map — and search on that. MOCs support set operations, so several regions become one filter. Here are the Rubin DP1 deep fields as a single selection:

from functools import reduce
from hats.pixel_math import region_to_moc

fields = {
    "ELAISS1": (9.45, -44.02),   # ELAIS-S1
    "XMM_LSS": (35.57, -4.82),   # XMM Large Scale Structure
    "ECDFS": (52.98, -28.12),    # Extended Chandra Deep Field South
    "COSMOS": (150.11, 2.23),    # COSMOS
    "EDFS_a": (58.9, -49.32),    # Euclid Deep Field South a
    "EDFS_b": (63.6, -47.60),    # Euclid Deep Field South b
}

cone_mocs = [
    # A 4 degree cone around each field center.
    region_to_moc.cone_to_moc(ra=ra, dec=dec, radius_arcsec=4 * 3600, max_depth=10)
    for ra, dec in fields.values()
]
cone_moc = reduce(lambda a, b: a.union(b), cone_mocs)

deep_fields = gaia.moc_search(cone_moc)

Because the union is a single MOC, LSDB opens only the partitions touching any of the six fields: one pass, not six. The same shape works for intersections and differences, and region_to_moc also offers box_to_moc, polygon_to_moc and pixel_list_to_moc as starting points.

A MOC can also be applied when the catalog is opened, so the partitions outside it are never considered at all:

from lsdb.core.search.region_search import MOCSearch

gaia = lsdb.open_catalog(
    "https://data.lsdb.io/hats/gaia_dr3",
    search_filter=MOCSearch(cone_moc),
    columns=["ra", "dec", "phot_g_mean_mag"],
)

Filter rows

query() takes a pandas-style expression string and stays lazy. Column names that aren't valid Python identifiers go in backticks.

bright = gaia.query("phot_g_mean_mag < 19")

limit = 18.5
brighter = gaia.query(f"phot_g_mean_mag < {limit}")

Apply my own function

map_partitions hands your function one partition at a time as a pandas dataframe, in parallel, across the whole catalog. Dask works out the resulting columns by running your function on an empty frame; pass meta explicitly when that guess comes out wrong, or when running your function even once is expensive.

def add_flux(df):
    df["flux"] = 10 ** (-0.4 * df["phot_g_mean_mag"])
    return df

with_flux = gaia.map_partitions(add_flux)

When your function is naturally per-object rather than per-partition, map_rows is the friendlier option — it handles the packing and unpacking for you, including nested columns.

Work with nested data#

Light curves, spectra, and detection histories are many measurements per object. A nested column holds a whole sub-table inside each row, so one row stays one object no matter how many epochs it has. Plenty of catalogs ship that way already — the ALeRCE ZTF object catalog below arrives with its light curves nested, nothing to assemble. See Ask for what that means for the data model, and for the catalogs that hand you bare list columns instead.

Look inside a nested frame

lc = lsdb.open_catalog(
    "https://data.lsdb.io/hats/alerce",
    columns=["oid", "mean_ra", "mean_dec", "lc"],
)
lc.nested_columns          # ['lc']

rows = lc.head(2)
rows["lc"].nest.columns    # lc_mjd, lc_magpsf, lc_fid, ...
rows["lc"].nest.to_flat()  # one row per detection instead of per object

Two objects unpack to dozens of detections. Flattening is there when you want it, but the nested form is what keeps a million objects a million rows.

Get one value per object

This is the common shape: reduce each object's light curve to a number, and keep it beside the object. map_rows passes each row in, nested column included.

import numpy as np

mean_mags = lc.map_rows(
    lambda row: float(np.nanmean(row["lc.lc_magpsf"])),
    columns=["lc.lc_magpsf"],
    meta={"mean_mag": "float64"},
    output_names="mean_mag",
    append_columns=True,      # keep the original columns alongside the result
)

Two related shapes worth knowing: join_nested attaches another catalog's rows as a nested column, and crossmatch_nested does the same for crossmatch results — so a many-to-one match stays one row per object.

Crossmatch#

Match two catalogs

Crossmatching is positional and lazy. LSDB aligns the two catalogs' partitions, then matches within each overlapping pair — so nothing is loaded until you compute.

matched = ztf.crossmatch(
    gaia,
    radius_arcsec=1.0,
    n_neighbors=1,
)
matched.head()

Keep all the matches

A crossmatch keeps the single closest neighbour by default, so you get one row per matched object. Ask for several with n_neighbors and you get a row per pair instead — five neighbours, five rows for that object. crossmatch_nested puts those matches in a nested column instead, so the object stays one row, with the distances and any other extra columns inside it.

matches = ztf.crossmatch_nested(
    gaia,
    radius_arcsec=1.0,
    n_neighbors=5,
    nested_column_name="gaia_matches",
    require_right_margin=True,
)
matches.nested_columns      # ['gaia_matches']

This is usually what you want when a match is genuinely one-to-many — every counterpart within a radius, every epoch of a light curve — and it keeps the nested-data recipes available on the result.

Not get wrong answers

gaia = lsdb.open_catalog(
    "https://data.lsdb.io/hats/gaia_dr3/gaia",
    margin_cache="https://data.lsdb.io/hats/gaia_dr3/gaia_10arcs",
    columns=["ra", "dec", "source_id", "phot_g_mean_mag"],
)

matched = ztf.crossmatch(
    gaia,
    radius_arcsec=1.0,
    require_right_margin=True,   # raise instead of quietly under-matching
)

Keep unmatched rows

how="inner" (the default) keeps only objects that matched. how="left" keeps every row of the left catalog, filling in nulls where nothing matched — useful when the question is "which of my objects have no counterpart?"

everything = ztf.crossmatch(gaia, radius_arcsec=1.0, how="left")

# Name the columns yourself when both catalogs share column names.
labelled = ztf.crossmatch(gaia, radius_arcsec=1.0, suffixes=("_ztf", "_gaia"))

Use other xmatch algorithms

The default is a k-d tree nearest-neighbour match. Other algorithms live in a separate package, lsdb-crossmatch, and are passed in with algorithm=.

pip install lsdb-crossmatch

Magnitude difference extends the default k-d tree match with a brightness check, so a positional match that is wildly the wrong magnitude can be filtered out. It adds _dist_arcsec and _magnitude_difference to the result.

from lsdb_crossmatch.mag_difference_crossmatch import MagnitudeDifferenceCrossmatch

matched = ztf.crossmatch(
    gaia,
    algorithm=MagnitudeDifferenceCrossmatch(
        left_mag_col="mean_mag_r",
        right_mag_col="phot_g_mean_mag",
        radius_arcsec=1.0,
    ),
    require_right_margin=True,
)

NWAY is a Bayesian matcher for multiwavelength catalogs, where positional errors differ between surveys and the nearest source is not always the right one. It needs the extra nwaylib dependency, and takes per-catalog positional errors, sky areas, and magnitude priors rather than a single radius.

Save my results#

Pull it into pandas

Once the result is genuinely small, compute() gives you an ordinary in-memory dataframe to plot, fit, or hand to another library.

result = bright.cone_search(ra=180.0, dec=10.0, radius_arcsec=600).compute()
result.to_parquet("my_result.parquet")

Write a HATS catalog

When the result is too big for memory, or you want to pick up where you left off tomorrow, write it back out as a HATS catalog. It reopens with open_catalog like any other, and it is partitioned, so the next workflow starts fast.

bright.write_catalog(
    "./my_catalogs/bright_gaia",
    catalog_name="bright_gaia",
    overwrite=True,
)

again = lsdb.open_catalog("./my_catalogs/bright_gaia")