Coverage for python/lsst/analysis/ap/ppdb.py: 95%
254 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-22 15:23 +0000
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-22 15:23 +0000
1# This file is part of analysis_ap.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""Load PPDB catalogs (DiaObjects, DiaSources, DiaForcedSources) through the
23Rubin Science Platform TAP service.
25These tools deliberately use the scientist-facing TAP interface (ADQL over
26``pyvo``) rather than a direct database connection, so that exercising them
27also exercises the interface scientists will use once the PPDB goes live.
29The prototype PPDB currently lives at the ``data-int`` RSP environment
30(``https://data-int.lsst.cloud/api/ppdbtap``). Tooling running elsewhere
31(e.g. the USDF RSP) reaches it cross-environment with a ``data-int`` RSP
32token carrying the ``read:tap`` scope, supplied via the ``RSP_TOKEN``
33environment variable.
35The PPDB ``DiaObject`` table holds only the latest version of each object,
36so it has one row per ``diaObjectId``. ``DiaSource`` and ``DiaForcedSource``
37are append-only.
39The production PPDB will allow cone searches only on ``DiaObject``. Sources
40and forced sources for a region must therefore be loaded object-first:
41cone-search DiaObjects, then load the sources of those objects. The
42``load_sources_for_region`` / ``load_forced_sources_for_region`` methods do
43this for you and are the encouraged entry points. ``load_sources`` /
44``load_forced_sources`` require a ``diaObjectId`` (scalar or iterable) and
45take no spatial argument. Direct, object-unrestricted access to the source
46tables is available only through the explicit, prototype-only
47``load_sources_by_cone`` / ``load_forced_sources_by_cone`` (spatial) and
48``load_sources_by_time_window`` / ``load_forced_sources_by_time_window``
49(temporal) methods, which warn when used.
51All loaders return `astropy.table.Table` objects. Their columns are returned
52in SDM-schema order, since the TAP service returns ``SELECT *`` columns
53alphabetically; passing an explicit ``columns`` list overrides this with the
54order you request.
55"""
57__all__ = [
58 "PpdbTap",
59 "DiaObjectLightCurve",
60 "region_from_exposure",
61 "DEFAULT_PPDB_TAP_URL",
62 "DEFAULT_PADDING_ARCSEC",
63 "DEFAULT_ROW_LIMIT",
64]
66import dataclasses
67import functools
68import logging
69import numbers
70import os
71import time
72from importlib.resources import files
74import astropy.table
75import astropy.units
76import pyvo
77import requests
78from astropy.coordinates import SkyCoord
80from lsst.pipe.tasks.schemaUtils import readSdmSchemaFile
82from .imageUtils import _pixel_to_sky, as_masked_image
84_log = logging.getLogger(__name__)
86# TAP endpoint of the prototype PPDB (the ``data-int`` RSP environment).
87DEFAULT_PPDB_TAP_URL = "https://data-int.lsst.cloud/api/ppdbtap"
89# Default padding (arcseconds) loaded around an exposure footprint.
90DEFAULT_PADDING_ARCSEC = 5.0
92# Default cap on rows returned by a single query, to guard against
93# accidentally pulling the whole (multi-million row) PPDB. Pass ``limit=None``
94# to disable.
95DEFAULT_ROW_LIMIT = 100000
97# Valid LSST band names; used to validate band filters before they are
98# interpolated into an ADQL string.
99_VALID_BANDS = frozenset("ugrizy")
101# The PPDB tables share the APDB SDM schema (``sdm_schemas/apdb.yaml``); its
102# column order is the canonical, scientifically-grouped order. Loaders use it
103# to reorder results, because the TAP service returns ``SELECT *`` columns
104# alphabetically (its ``tap_schema.columns`` has no ``column_index``). Both
105# helpers are memoized, so the schema file is parsed at most once.
108@functools.cache
109def _sdm_schema():
110 """Parse and return the APDB SDM schema (shared by the PPDB tables)."""
111 path = os.fspath(files("lsst.sdm.schemas").joinpath("apdb.yaml"))
112 return readSdmSchemaFile(path)
115@functools.cache
116def _sdm_column_order(table_name):
117 """Return the SDM-schema column order for a PPDB table, or ``()``.
119 Parameters
120 ----------
121 table_name : `str`
122 Unqualified table name, e.g. ``"DiaSource"``.
124 Returns
125 -------
126 columns : `tuple` [`str`]
127 Column names in schema order, or an empty tuple if the table is not
128 in the SDM schema.
129 """
130 table_def = _sdm_schema().get(table_name)
131 if table_def is None: 131 ↛ 132line 131 didn't jump to line 132 because the condition on line 131 was never true
132 return ()
133 return tuple(column.name for column in table_def.columns)
136def region_from_exposure(exposure, padding=DEFAULT_PADDING_ARCSEC):
137 """Compute a cone (ra, dec, radius) covering an exposure's sky footprint.
139 The cone is the smallest circle centered on the exposure's bounding-box
140 center that contains all four corners, expanded by ``padding``
141 arcseconds. A circumscribing circle (rather than the exact polygon) keeps
142 the resulting ADQL robust: ``CIRCLE`` is universally supported by TAP
143 services, whereas ``POLYGON`` support and winding conventions vary.
145 Parameters
146 ----------
147 exposure : `lsst.images.MaskedImage` or `lsst.afw.image.Exposure`
148 Image whose sky projection and bounding box define the region, such
149 as an `lsst.images.VisitImage`. A legacy exposure is converted with
150 `as_masked_image`.
151 padding : `float`, optional
152 Margin to add around the footprint, in arcseconds.
154 Returns
155 -------
156 ra : `float`
157 Cone center right ascension, in degrees (ICRS).
158 dec : `float`
159 Cone center declination, in degrees (ICRS).
160 radius : `float`
161 Cone radius, in degrees.
163 Raises
164 ------
165 ValueError
166 If the image has no WCS, or ``padding`` is negative.
167 """
168 if padding < 0:
169 raise ValueError(f"padding must be non-negative, got {padding}")
170 image = as_masked_image(exposure)
171 if image.sky_projection is None:
172 raise ValueError("Image has no WCS; cannot compute a sky region.")
173 bbox = image.bbox
174 x_min, x_max = bbox.x.start, bbox.x.stop - 1
175 y_min, y_max = bbox.y.start, bbox.y.stop - 1
176 # The center first, then the centers of the four corner pixels.
177 ra, dec = _pixel_to_sky(image.sky_projection,
178 [0.5*(x_min + x_max), x_min, x_max, x_max, x_min],
179 [0.5*(y_min + y_max), y_min, y_min, y_max, y_max])
180 points = SkyCoord(ra, dec, unit="deg")
181 radius = points[0].separation(points[1:]).max() + padding*astropy.units.arcsec
182 return float(ra[0]), float(dec[0]), float(radius.deg)
185@dataclasses.dataclass
186class DiaObjectLightCurve:
187 """The full PPDB record for a single diaObject.
189 Attributes
190 ----------
191 diaObjectId : `int`
192 Identifier of the diaObject.
193 dia_object : `astropy.table.Row`
194 The DiaObject row.
195 dia_sources : `astropy.table.Table`
196 All DiaSources associated with the object, time-ordered.
197 dia_forced_sources : `astropy.table.Table`
198 All DiaForcedSources for the object, time-ordered.
199 """
201 diaObjectId: int
202 dia_object: "astropy.table.Row"
203 dia_sources: astropy.table.Table
204 dia_forced_sources: astropy.table.Table
206 @property
207 def n_sources(self):
208 """Number of associated DiaSources (`int`)."""
209 return len(self.dia_sources)
211 @property
212 def n_forced_sources(self):
213 """Number of associated DiaForcedSources (`int`)."""
214 return len(self.dia_forced_sources)
216 def __repr__(self):
217 return (f"DiaObjectLightCurve(diaObjectId={self.diaObjectId}, "
218 f"n_sources={self.n_sources}, "
219 f"n_forced_sources={self.n_forced_sources})")
222class PpdbTap:
223 """Load PPDB catalogs through the RSP TAP service as astropy Tables.
225 Parameters
226 ----------
227 url : `str`, optional
228 TAP endpoint to query. Defaults to the prototype PPDB at ``data-int``.
229 token : `str`, optional
230 RSP token with the ``read:tap`` scope. If not given (and ``service``
231 is also not given), the ``RSP_TOKEN`` environment variable is used.
232 service : `pyvo.dal.TAPService`, optional
233 Pre-built TAP service. When supplied, ``url`` and ``token`` are
234 ignored. Intended for use with unit tests.
235 log : `logging.Logger`, optional
236 Logger to use; defaults to the module logger.
238 Notes
239 -----
240 The token is read only from the environment or the constructor argument
241 and is never logged.
242 """
244 def __init__(self, url=DEFAULT_PPDB_TAP_URL, *, token=None,
245 service=None, log=None):
246 self.url = url
247 self.log = log if log is not None else _log
248 if service is not None:
249 self._service = service
250 return
251 if token is None:
252 token = os.environ.get("RSP_TOKEN")
253 if not token:
254 raise RuntimeError(
255 "No RSP token available for the PPDB TAP service. Generate "
256 "a token with the 'read:tap' scope on data-int "
257 "(https://data-int.lsst.cloud) and set it in the RSP_TOKEN "
258 "environment variable, or pass token=...")
259 session = requests.Session()
260 session.headers["Authorization"] = f"Bearer {token}"
261 self._service = pyvo.dal.TAPService(url, session=session)
263 # ------------------------------------------------------------------
264 # Query execution
265 # ------------------------------------------------------------------
266 def run_query(self, adql, *, maxrec=None, order_table=None):
267 """Run an ADQL query and return the result as an astropy Table.
269 Parameters
270 ----------
271 adql : `str`
272 The ADQL query to execute.
273 maxrec : `int`, optional
274 Server-side cap on the number of records returned.
275 order_table : `str`, optional
276 If given (e.g. ``"DiaSource"``), reorder the result columns into
277 that table's SDM-schema order. The loaders pass this so that
278 ``SELECT *`` results, which the TAP service returns alphabetically,
279 come back in the schema's logical column order. Columns not in the
280 schema are appended in their original order.
282 Returns
283 -------
284 table : `astropy.table.Table`
285 The query result.
286 """
287 self.log.debug("Executing ADQL query: %s", adql)
288 start = time.perf_counter()
289 results = self._service.run_async(adql, maxrec=maxrec)
290 table = results.to_table()
291 if order_table is not None:
292 table = self._order_columns(table, order_table)
293 elapsed = time.perf_counter() - start
294 self.log.info("Query returned %d rows in %.3f s", len(table), elapsed)
295 return table
297 @staticmethod
298 def _order_columns(table, table_name):
299 """Reorder a result table's columns into SDM-schema order.
301 Every result column is expected to be defined in the SDM schema for
302 ``table_name`` (the PPDB does not add columns outside it), so columns
303 are returned in schema order. Returns the table unchanged if the
304 schema is unknown or the order already matches.
306 Raises
307 ------
308 RuntimeError
309 If the result contains a column not present in the SDM schema,
310 which signals an unexpected schema mismatch.
311 """
312 schema_order = _sdm_column_order(table_name)
313 if not schema_order: 313 ↛ 314line 313 didn't jump to line 314 because the condition on line 313 was never true
314 return table
315 schema_set = set(schema_order)
316 unknown = [c for c in table.colnames if c not in schema_set]
317 if unknown:
318 raise RuntimeError(
319 f"ppdb.{table_name} returned columns not in the SDM schema: "
320 f"{unknown}")
321 new_order = [c for c in schema_order if c in table.colnames]
322 if new_order == table.colnames:
323 return table
324 return table[new_order]
326 # ------------------------------------------------------------------
327 # ADQL construction helpers
328 # ------------------------------------------------------------------
329 @staticmethod
330 def _select_clause(columns, limit):
331 """Build the ``SELECT [TOP n] cols`` prefix of a query."""
332 cols = "*" if not columns else ", ".join(columns)
333 top = f"TOP {int(limit)} " if limit is not None else ""
334 return f"SELECT {top}{cols}"
336 @staticmethod
337 def _cone_clause(ra, dec, radius):
338 """Build an ADQL cone-search predicate (all angles in degrees)."""
339 return (f"CONTAINS(POINT('ICRS', ra, dec), "
340 f"CIRCLE('ICRS', {float(ra)}, {float(dec)}, "
341 f"{float(radius)})) = 1")
343 def _resolve_cone(self, ra, dec, radius, exposure, padding):
344 """Resolve spatial arguments into a (ra, dec, radius) cone or None.
346 Either an ``exposure`` or an explicit ``ra``/``dec``/``radius`` triple
347 may be given, but not both. Returns None when no spatial constraint
348 was requested. Explicit ``radius`` is in degrees; ``padding`` (used
349 only with ``exposure``) is in arcseconds.
350 """
351 if exposure is not None:
352 if any(v is not None for v in (ra, dec, radius)):
353 raise ValueError("Specify either exposure= or ra/dec/radius, "
354 "not both.")
355 return region_from_exposure(exposure, padding=padding)
356 if ra is None and dec is None and radius is None:
357 return None
358 if None in (ra, dec, radius):
359 raise ValueError("ra, dec, and radius must all be provided for a "
360 "cone search.")
361 return float(ra), float(dec), float(radius)
363 @staticmethod
364 def _band_clause(bands):
365 """Build an ADQL ``band IN (...)`` predicate, validating each band."""
366 quoted = []
367 for band in bands:
368 if band not in _VALID_BANDS:
369 raise ValueError(f"Unknown band {band!r}; expected one of "
370 f"{sorted(_VALID_BANDS)}.")
371 quoted.append(f"'{band}'")
372 return f"band IN ({', '.join(quoted)})"
374 def _finish(self, table, limit, table_name):
375 """Warn on probable truncation and return the table unchanged."""
376 if limit is not None and len(table) == limit:
377 self.log.warning(
378 "%s query returned exactly the row limit (%d); results are "
379 "likely truncated. Narrow the query or raise limit.",
380 table_name, limit)
381 return table
383 # ------------------------------------------------------------------
384 # DiaObject loaders
385 # ------------------------------------------------------------------
386 def load_objects(self, *, ra=None, dec=None, radius=None, exposure=None,
387 padding=DEFAULT_PADDING_ARCSEC, columns=None,
388 limit=DEFAULT_ROW_LIMIT):
389 """Load DiaObjects, optionally within a spatial region.
391 Parameters
392 ----------
393 ra, dec, radius : `float`, optional
394 Cone-search center and radius, in degrees. All three must be
395 given together; mutually exclusive with ``exposure``.
396 exposure : `lsst.images.MaskedImage` or `lsst.afw.image.Exposure`, \
397 optional
398 Load objects covering this exposure's footprint (see
399 `region_from_exposure`). Mutually exclusive with
400 ``ra``/``dec``/``radius``.
401 padding : `float`, optional
402 Margin around the exposure footprint, in arcseconds (only used
403 with ``exposure``).
404 columns : `list` [`str`], optional
405 Columns to select; defaults to all columns.
406 limit : `int`, optional
407 Maximum number of rows to return; None means no limit.
409 Returns
410 -------
411 objects : `astropy.table.Table`
412 The matching DiaObjects.
413 """
414 cone = self._resolve_cone(ra, dec, radius, exposure, padding)
415 where = f" WHERE {self._cone_clause(*cone)}" if cone is not None else ""
416 adql = (f"{self._select_clause(columns, limit)} FROM ppdb.DiaObject"
417 f"{where} ORDER BY diaObjectId")
418 order_table = "DiaObject" if columns is None else None
419 table = self._finish(self.run_query(adql, order_table=order_table),
420 limit, "DiaObject")
421 if len(table) and "diaObjectId" in table.colnames:
422 n_unique = len(set(table["diaObjectId"].tolist()))
423 if n_unique != len(table):
424 self.log.warning(
425 "DiaObject query returned %d rows but only %d unique "
426 "diaObjectIds; the PPDB may contain duplicate current "
427 "versions.",
428 len(table), n_unique)
429 return table
431 def load_object(self, diaObjectId, *, columns=None):
432 """Load a single DiaObject.
434 Parameters
435 ----------
436 diaObjectId : `int`
437 Identifier of the object to load.
438 columns : `list` [`str`], optional
439 Columns to select; defaults to all columns.
441 Returns
442 -------
443 object : `astropy.table.Row`
444 The current DiaObject row.
446 Raises
447 ------
448 RuntimeError
449 If the object is not in the PPDB.
450 """
451 adql = (f"{self._select_clause(columns, None)} FROM ppdb.DiaObject "
452 f"WHERE diaObjectId = {int(diaObjectId)}")
453 order_table = "DiaObject" if columns is None else None
454 table = self.run_query(adql, order_table=order_table)
455 if len(table) == 0:
456 raise RuntimeError(
457 f"diaObjectId={diaObjectId} not found in ppdb.DiaObject")
458 if len(table) > 1: 458 ↛ 459line 458 didn't jump to line 459 because the condition on line 458 was never true
459 self.log.warning(
460 "diaObjectId=%s has %d rows; returning the first.",
461 diaObjectId, len(table))
462 return table[0]
464 # ------------------------------------------------------------------
465 # DiaSource / DiaForcedSource loaders
466 #
467 # The production PPDB will permit cone searches only on DiaObject. The
468 # encouraged way to retrieve sources for a region is therefore object-
469 # first: cone-search DiaObjects (`load_*_for_region`), then load their
470 # sources by diaObjectId. `load_sources`/`load_forced_sources` require a
471 # diaObjectId and take no spatial argument. Direct, object-unrestricted
472 # source-table access is confined to the explicit, prototype-only
473 # `load_*_by_cone` and `load_*_by_time_window` methods.
474 # ------------------------------------------------------------------
475 def _source_filter_clauses(self, bands, mjd_begin, mjd_end):
476 """Build the time/band filter clauses shared by the source loaders."""
477 clauses = []
478 if mjd_begin is not None:
479 clauses.append(f"midpointMjdTai >= {float(mjd_begin)}")
480 if mjd_end is not None:
481 clauses.append(f"midpointMjdTai < {float(mjd_end)}")
482 if bands:
483 clauses.append(self._band_clause(bands))
484 return clauses
486 def _run_source_query(self, table_name, clauses, columns, limit):
487 """Run a single time-ordered SELECT against a source table."""
488 where = (" WHERE " + " AND ".join(clauses)) if clauses else ""
489 adql = (f"{self._select_clause(columns, limit)} FROM ppdb.{table_name}"
490 f"{where} ORDER BY midpointMjdTai")
491 order_table = table_name if columns is None else None
492 return self._finish(self.run_query(adql, order_table=order_table),
493 limit, table_name)
495 def _load_source_like(self, table_name, diaObjectId, extra_clauses,
496 columns, limit, id_chunk_size):
497 """Dispatch a source load by one id or many ids.
499 ``diaObjectId`` is required: source tables are loaded by object. An
500 iterable (including the empty list a region search may produce) takes
501 the chunked ``IN`` path; a scalar takes a single ``=`` query.
502 """
503 if diaObjectId is None:
504 raise ValueError(
505 "diaObjectId is required; DiaSources and DiaForcedSources are "
506 "loaded by object. Use load_*_for_region for a spatial "
507 "region, or the prototype-only load_*_by_cone / "
508 "load_*_by_time_window methods for direct source-table access.")
509 if isinstance(diaObjectId, numbers.Integral):
510 clauses = [f"diaObjectId = {int(diaObjectId)}"] + list(extra_clauses)
511 return self._run_source_query(table_name, clauses, columns, limit)
512 return self._load_sources_for_ids(table_name, diaObjectId,
513 extra_clauses, columns, limit,
514 id_chunk_size)
516 def _load_sources_for_ids(self, table_name, ids, extra_clauses, columns,
517 limit, id_chunk_size):
518 """Load sources for many diaObjectIds via chunked ``IN`` queries.
520 A cone search can return more diaObjectIds than fit in a single
521 ``diaObjectId IN (...)`` clause, so the ids are split into chunks
522 of ``id_chunk_size`` and the per-chunk results are concatenated and
523 re-sorted by time.
524 ``limit`` caps the total rows returned across all chunks.
525 """
526 if id_chunk_size < 1: 526 ↛ 527line 526 didn't jump to line 527 because the condition on line 526 was never true
527 raise ValueError("id_chunk_size must be >= 1")
528 ids = list(dict.fromkeys(int(i) for i in ids))
529 if not ids or (limit is not None and limit <= 0):
530 # Empty result, but issue a query so columns/dtypes are populated.
531 return self._run_source_query(table_name, ["1 = 0"], columns, None)
533 order_table = table_name if columns is None else None
534 pages = []
535 total = 0
536 truncated = False
537 for start in range(0, len(ids), id_chunk_size):
538 if limit is not None and total >= limit: 538 ↛ 539line 538 didn't jump to line 539 because the condition on line 538 was never true
539 truncated = True
540 break
541 chunk = ids[start:start + id_chunk_size]
542 remaining = None if limit is None else limit - total
543 id_list = ", ".join(str(i) for i in chunk)
544 clauses = [f"diaObjectId IN ({id_list})"] + list(extra_clauses)
545 where = " WHERE " + " AND ".join(clauses)
546 adql = (f"{self._select_clause(columns, remaining)} "
547 f"FROM ppdb.{table_name}{where} ORDER BY midpointMjdTai")
548 page = self.run_query(adql, order_table=order_table)
549 total += len(page)
550 pages.append(page)
552 table = (pages[0] if len(pages) == 1
553 else astropy.table.vstack(pages, metadata_conflicts="silent"))
554 if len(table) and "midpointMjdTai" in table.colnames: 554 ↛ 556line 554 didn't jump to line 556 because the condition on line 554 was always true
555 table.sort("midpointMjdTai")
556 if limit is not None and (truncated or total == limit): 556 ↛ 557line 556 didn't jump to line 557 because the condition on line 556 was never true
557 self.log.warning(
558 "%s query reached the row limit (%d) while loading sources for "
559 "%d objects; results are likely truncated. Raise limit or "
560 "narrow the query.", table_name, limit, len(ids))
561 return table
563 def _direct_source_query(self, table_name, region_method, kind, clauses,
564 columns, limit):
565 """Warn, then run a prototype-only direct source-table query.
567 ``load_sources``/``load_forced_sources`` require a diaObjectId; these
568 direct (object-unrestricted) source-table searches will not exist in
569 the production PPDB, so callers reach them only through the explicit
570 ``load_*_by_cone`` / ``load_*_by_time_window`` methods, which route
571 here and emit this warning.
572 """
573 self.log.warning(
574 "Direct %s %s search is a prototype-only capability and will not "
575 "be available in the production PPDB; prefer %s().",
576 table_name, kind, region_method)
577 return self._run_source_query(table_name, clauses, columns, limit)
579 def _load_by_cone(self, table_name, region_method, *, ra, dec, radius,
580 exposure, padding, bands, mjd_begin, mjd_end, columns,
581 limit):
582 """Direct cone search on a source table (prototype-only)."""
583 cone = self._resolve_cone(ra, dec, radius, exposure, padding)
584 if cone is None:
585 raise ValueError("A cone search requires ra/dec/radius or an "
586 "exposure.")
587 clauses = [self._cone_clause(*cone)]
588 clauses += self._source_filter_clauses(bands, mjd_begin, mjd_end)
589 return self._direct_source_query(table_name, region_method, "cone",
590 clauses, columns, limit)
592 def _load_by_time_window(self, table_name, region_method, *, mjd_begin,
593 mjd_end, bands, columns, limit):
594 """Direct time-window scan of a source table (prototype-only)."""
595 if mjd_begin is None and mjd_end is None:
596 raise ValueError("A time-window search requires mjd_begin and/or "
597 "mjd_end.")
598 clauses = self._source_filter_clauses(bands, mjd_begin, mjd_end)
599 return self._direct_source_query(table_name, region_method,
600 "time-window", clauses, columns,
601 limit)
603 def load_sources(self, *, diaObjectId, mjd_begin=None, mjd_end=None,
604 bands=None, columns=None, limit=DEFAULT_ROW_LIMIT,
605 id_chunk_size=1000):
606 """Load DiaSources for a given object (or objects), by time and band.
608 DiaSources are loaded by object: ``diaObjectId`` is required. This
609 loader has no spatial argument, because the production PPDB does not
610 permit searches on the source tables. To load sources for a region,
611 use `load_sources_for_region` (object-first); for direct, object-
612 unrestricted prototype-only access, use `load_sources_by_cone` or
613 `load_sources_by_time_window`.
615 Parameters
616 ----------
617 diaObjectId : `int` or iterable of `int`
618 Restrict to sources of this object, or of any of these objects
619 (issued as chunked ``diaObjectId IN (...)`` queries). Required.
620 mjd_begin, mjd_end : `float`, optional
621 Restrict to ``mjd_begin <= midpointMjdTai < mjd_end`` (TAI).
622 bands : `list` [`str`], optional
623 Restrict to these bands (subset of ``ugrizy``).
624 columns : `list` [`str`], optional
625 Columns to select; defaults to all columns.
626 limit : `int`, optional
627 Maximum number of rows to return; None means no limit.
628 id_chunk_size : `int`, optional
629 Number of diaObjectIds per ``IN`` query when an iterable is given.
631 Returns
632 -------
633 sources : `astropy.table.Table`
634 The matching DiaSources, ordered by ``midpointMjdTai``.
636 Raises
637 ------
638 ValueError
639 If ``diaObjectId`` is None.
640 """
641 extra = self._source_filter_clauses(bands, mjd_begin, mjd_end)
642 return self._load_source_like("DiaSource", diaObjectId, extra,
643 columns, limit, id_chunk_size)
645 def load_forced_sources(self, *, diaObjectId, mjd_begin=None,
646 mjd_end=None, bands=None, columns=None,
647 limit=DEFAULT_ROW_LIMIT, id_chunk_size=1000):
648 """Load DiaForcedSources, with the same filters as `load_sources`.
650 Like `load_sources`, ``diaObjectId`` is required and there is no
651 spatial argument; use `load_forced_sources_for_region`, or the
652 prototype-only `load_forced_sources_by_cone` /
653 `load_forced_sources_by_time_window`.
655 Returns
656 -------
657 forced_sources : `astropy.table.Table`
658 The matching DiaForcedSources, ordered by ``midpointMjdTai``.
660 Raises
661 ------
662 ValueError
663 If ``diaObjectId`` is None.
664 """
665 extra = self._source_filter_clauses(bands, mjd_begin, mjd_end)
666 return self._load_source_like("DiaForcedSource", diaObjectId, extra,
667 columns, limit, id_chunk_size)
669 def load_sources_for_region(self, *, ra=None, dec=None, radius=None,
670 exposure=None, padding=DEFAULT_PADDING_ARCSEC,
671 mjd_begin=None, mjd_end=None, bands=None,
672 columns=None, limit=DEFAULT_ROW_LIMIT,
673 object_limit=DEFAULT_ROW_LIMIT,
674 id_chunk_size=1000):
675 """Load DiaSources for a region, the production-compatible way.
677 Performs the two-step workflow the production PPDB requires: cone-
678 search DiaObjects (the only table that allows spatial search), then
679 load the DiaSources of those objects.
681 Parameters
682 ----------
683 ra, dec, radius : `float`, optional
684 Cone center and radius, in degrees; mutually exclusive with
685 ``exposure``. The radius is applied to *object* positions.
686 exposure : `lsst.images.MaskedImage` or `lsst.afw.image.Exposure`, \
687 optional
688 Load sources for objects covering this exposure's footprint.
689 padding : `float`, optional
690 Margin around the exposure footprint, in arcseconds.
691 mjd_begin, mjd_end, bands, columns : optional
692 Passed through to the source query (see `load_sources`).
693 limit : `int`, optional
694 Cap on returned sources; None means no limit.
695 object_limit : `int`, optional
696 Cap on the DiaObjects found by the cone search.
697 id_chunk_size : `int`, optional
698 diaObjectIds per ``IN`` query when fetching sources.
700 Returns
701 -------
702 sources : `astropy.table.Table`
703 DiaSources of the objects in the region, ordered by time.
705 Notes
706 -----
707 Because this is object-first, DiaSources not associated with a
708 returned DiaObject (e.g. unassociated or ssObject-only detections)
709 are not included. To find those, use `load_sources_by_cone`.
710 """
711 ids = self._region_object_ids(
712 ra, dec, radius, exposure, padding, object_limit)
713 return self.load_sources(diaObjectId=ids, mjd_begin=mjd_begin,
714 mjd_end=mjd_end, bands=bands, columns=columns,
715 limit=limit, id_chunk_size=id_chunk_size)
717 def load_forced_sources_for_region(self, *, ra=None, dec=None,
718 radius=None, exposure=None,
719 padding=DEFAULT_PADDING_ARCSEC,
720 mjd_begin=None, mjd_end=None, bands=None,
721 columns=None, limit=DEFAULT_ROW_LIMIT,
722 object_limit=DEFAULT_ROW_LIMIT,
723 id_chunk_size=1000):
724 """Load DiaForcedSources for a region (object-first).
726 Same two-step workflow and caveats as `load_sources_for_region`.
728 Returns
729 -------
730 forced_sources : `astropy.table.Table`
731 DiaForcedSources of the objects in the region, ordered by time.
732 """
733 ids = self._region_object_ids(
734 ra, dec, radius, exposure, padding, object_limit)
735 return self.load_forced_sources(
736 diaObjectId=ids, mjd_begin=mjd_begin, mjd_end=mjd_end, bands=bands,
737 columns=columns, limit=limit, id_chunk_size=id_chunk_size)
739 def _region_object_ids(self, ra, dec, radius, exposure, padding,
740 object_limit):
741 """Cone-search DiaObjects and return their diaObjectIds."""
742 objects = self.load_objects(ra=ra, dec=dec, radius=radius,
743 exposure=exposure, padding=padding,
744 columns=["diaObjectId"],
745 limit=object_limit)
746 if not len(objects):
747 return []
748 return objects["diaObjectId"].tolist()
750 def load_sources_by_cone(self, *, ra=None, dec=None, radius=None,
751 exposure=None, padding=DEFAULT_PADDING_ARCSEC,
752 mjd_begin=None, mjd_end=None, bands=None,
753 columns=None, limit=DEFAULT_ROW_LIMIT):
754 """Directly cone-search the DiaSource table (prototype-only).
756 .. warning::
758 The production PPDB will not allow spatial searches on the source
759 tables. This works only against the prototype and is intended for
760 debugging/validation (e.g. finding DiaSources with no associated
761 DiaObject). For the production-compatible region workflow use
762 `load_sources_for_region`.
764 Parameters mirror `load_sources_for_region` minus ``object_limit``.
766 Returns
767 -------
768 sources : `astropy.table.Table`
769 DiaSources within the cone, ordered by ``midpointMjdTai``.
770 """
771 return self._load_by_cone(
772 "DiaSource", "load_sources_for_region", ra=ra, dec=dec,
773 radius=radius, exposure=exposure, padding=padding, bands=bands,
774 mjd_begin=mjd_begin, mjd_end=mjd_end, columns=columns, limit=limit)
776 def load_forced_sources_by_cone(self, *, ra=None, dec=None, radius=None,
777 exposure=None,
778 padding=DEFAULT_PADDING_ARCSEC,
779 mjd_begin=None, mjd_end=None, bands=None,
780 columns=None, limit=DEFAULT_ROW_LIMIT):
781 """Directly cone-search the DiaForcedSource table (prototype-only).
783 .. warning::
785 Not supported by the production PPDB; prototype/debugging only.
786 Prefer `load_forced_sources_for_region`.
788 Returns
789 -------
790 forced_sources : `astropy.table.Table`
791 DiaForcedSources within the cone, ordered by ``midpointMjdTai``.
792 """
793 return self._load_by_cone(
794 "DiaForcedSource", "load_forced_sources_for_region", ra=ra,
795 dec=dec, radius=radius, exposure=exposure, padding=padding,
796 bands=bands, mjd_begin=mjd_begin, mjd_end=mjd_end, columns=columns,
797 limit=limit)
799 def load_sources_by_time_window(self, *, mjd_begin=None, mjd_end=None,
800 bands=None, columns=None,
801 limit=DEFAULT_ROW_LIMIT):
802 """Load DiaSources in a time window, any object (prototype-only).
804 .. warning::
806 The production PPDB will not allow object-unrestricted searches on
807 the source tables. This works only against the prototype and is
808 intended for debugging (e.g. inspecting everything processed in an
809 MJD range). For the production-compatible workflow use
810 `load_sources_for_region`.
812 Parameters
813 ----------
814 mjd_begin, mjd_end : `float`, optional
815 Time window on ``midpointMjdTai`` (TAI), as
816 ``mjd_begin <= midpointMjdTai < mjd_end``. At least one bound is
817 required (an unbounded source-table scan is not permitted).
818 bands : `list` [`str`], optional
819 Restrict to these bands (subset of ``ugrizy``).
820 columns : `list` [`str`], optional
821 Columns to select; defaults to all columns.
822 limit : `int`, optional
823 Maximum number of rows to return; None means no limit.
825 Returns
826 -------
827 sources : `astropy.table.Table`
828 DiaSources in the window, ordered by ``midpointMjdTai``.
830 Raises
831 ------
832 ValueError
833 If neither ``mjd_begin`` nor ``mjd_end`` is given.
834 """
835 return self._load_by_time_window(
836 "DiaSource", "load_sources_for_region", mjd_begin=mjd_begin,
837 mjd_end=mjd_end, bands=bands, columns=columns, limit=limit)
839 def load_forced_sources_by_time_window(self, *, mjd_begin=None,
840 mjd_end=None, bands=None,
841 columns=None,
842 limit=DEFAULT_ROW_LIMIT):
843 """Load DiaForcedSources in a time window, any object (prototype-only).
845 .. warning::
847 Not supported by the production PPDB; prototype/debugging only.
848 Prefer `load_forced_sources_for_region`. At least one of
849 ``mjd_begin`` / ``mjd_end`` is required.
851 Returns
852 -------
853 forced_sources : `astropy.table.Table`
854 DiaForcedSources in the window, ordered by ``midpointMjdTai``.
856 Raises
857 ------
858 ValueError
859 If neither ``mjd_begin`` nor ``mjd_end`` is given.
860 """
861 return self._load_by_time_window(
862 "DiaForcedSource", "load_forced_sources_for_region",
863 mjd_begin=mjd_begin, mjd_end=mjd_end, bands=bands, columns=columns,
864 limit=limit)
866 def load_sources_for_object(self, diaObjectId, *, columns=None,
867 limit=DEFAULT_ROW_LIMIT):
868 """Load all DiaSources associated with one diaObject.
870 Returns
871 -------
872 sources : `astropy.table.Table`
873 The object's DiaSources, ordered by ``midpointMjdTai``.
874 """
875 return self.load_sources(diaObjectId=diaObjectId, columns=columns,
876 limit=limit)
878 def load_forced_sources_for_object(self, diaObjectId, *, columns=None,
879 limit=DEFAULT_ROW_LIMIT):
880 """Load all DiaForcedSources associated with one diaObject.
882 Returns
883 -------
884 forced_sources : `astropy.table.Table`
885 The object's DiaForcedSources, ordered by ``midpointMjdTai``.
886 """
887 return self.load_forced_sources(diaObjectId=diaObjectId,
888 columns=columns, limit=limit)
890 def load_source(self, diaSourceId, *, columns=None):
891 """Load a single DiaSource by id.
893 Raises
894 ------
895 RuntimeError
896 If the source does not exist.
897 """
898 return self._load_one("DiaSource", "diaSourceId", diaSourceId, columns)
900 def load_forced_source(self, diaForcedSourceId, *, columns=None):
901 """Load a single DiaForcedSource by id.
903 Raises
904 ------
905 RuntimeError
906 If the forced source does not exist.
907 """
908 return self._load_one("DiaForcedSource", "diaForcedSourceId",
909 diaForcedSourceId, columns)
911 def _load_one(self, table_name, id_column, id_value, columns):
912 """Load exactly one row from a table by id, raising if absent."""
913 adql = (f"{self._select_clause(columns, None)} FROM ppdb.{table_name} "
914 f"WHERE {id_column} = {int(id_value)}")
915 order_table = table_name if columns is None else None
916 table = self.run_query(adql, order_table=order_table)
917 if len(table) == 0:
918 raise RuntimeError(
919 f"{id_column}={id_value} not found in ppdb.{table_name}")
920 return table[0]
922 # ------------------------------------------------------------------
923 # Light-curve assembly
924 # ------------------------------------------------------------------
925 def load_light_curve(self, diaObjectId):
926 """Assemble the full PPDB record for one diaObject.
928 Loads the DiaObject plus all of its DiaSources and
929 DiaForcedSources, time-ordered.
931 Parameters
932 ----------
933 diaObjectId : `int`
934 Identifier of the object.
936 Returns
937 -------
938 light_curve : `DiaObjectLightCurve`
939 The object and its associated source time series.
941 Raises
942 ------
943 RuntimeError
944 If the object is not in the PPDB.
945 """
946 dia_object = self.load_object(diaObjectId)
947 dia_sources = self.load_sources_for_object(diaObjectId)
948 dia_forced_sources = self.load_forced_sources_for_object(diaObjectId)
949 return DiaObjectLightCurve(
950 diaObjectId=int(diaObjectId),
951 dia_object=dia_object,
952 dia_sources=dia_sources,
953 dia_forced_sources=dia_forced_sources)