Coverage for python/lsst/pipe/tasks/diff_matched_tract_catalog.py: 64%
190 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-25 09:10 +0000
1# This file is part of pipe_tasks.
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__all__ = [
23 'DiffMatchedTractCatalogConfig', 'DiffMatchedTractCatalogTask', 'MatchedCatalogFluxesConfig',
24]
26import lsst.afw.geom as afwGeom
27from lsst.meas.astrom.matcher_probabilistic import ConvertCatalogCoordinatesConfig
28from lsst.meas.astrom.match_probabilistic_task import radec_to_xy
29import lsst.pex.config as pexConfig
30import lsst.pipe.base as pipeBase
31import lsst.pipe.base.connectionTypes as cT
32from lsst.skymap import BaseSkyMap
33from lsst.daf.butler import DatasetProvenance
35import astropy.table
36import astropy.units as u
37import numpy as np
38from smatch.matcher import sphdist
39from typing import Sequence
42def is_sequence_set(x: Sequence):
43 return len(x) == len(set(x))
46DiffMatchedTractCatalogBaseTemplates = {
47 "name_input_cat_ref": "truth_summary",
48 "name_input_cat_target": "objectTable_tract",
49 "name_skymap": BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
50}
53class DiffMatchedTractCatalogConnections(
54 pipeBase.PipelineTaskConnections,
55 dimensions=("tract", "skymap"),
56 defaultTemplates=DiffMatchedTractCatalogBaseTemplates,
57):
58 cat_ref = cT.Input(
59 doc="Reference object catalog to match from",
60 name="{name_input_cat_ref}",
61 storageClass="ArrowAstropy",
62 dimensions=("tract", "skymap"),
63 deferLoad=True,
64 )
65 cat_target = cT.Input(
66 doc="Target object catalog to match",
67 name="{name_input_cat_target}",
68 storageClass="ArrowAstropy",
69 dimensions=("tract", "skymap"),
70 deferLoad=True,
71 )
72 skymap = cT.Input(
73 doc="Input definition of geometry/bbox and projection/wcs for coadded exposures",
74 name="{name_skymap}",
75 storageClass="SkyMap",
76 dimensions=("skymap",),
77 )
78 cat_match_ref = cT.Input(
79 doc="Reference match catalog with indices of target matches",
80 name="match_ref_{name_input_cat_ref}_{name_input_cat_target}",
81 storageClass="ArrowAstropy",
82 dimensions=("tract", "skymap"),
83 deferLoad=True,
84 )
85 cat_match_target = cT.Input(
86 doc="Target match catalog with indices of references matches",
87 name="match_target_{name_input_cat_ref}_{name_input_cat_target}",
88 storageClass="ArrowAstropy",
89 dimensions=("tract", "skymap"),
90 deferLoad=True,
91 )
92 columns_match_target = cT.Input(
93 doc="Target match catalog columns",
94 name="match_target_{name_input_cat_ref}_{name_input_cat_target}.columns",
95 storageClass="ArrowColumnList",
96 dimensions=("tract", "skymap"),
97 )
98 cat_matched = cT.Output(
99 doc="Catalog with reference and target columns for joined sources",
100 name="matched_{name_input_cat_ref}_{name_input_cat_target}",
101 storageClass="ArrowAstropy",
102 dimensions=("tract", "skymap"),
103 )
105 def __init__(self, *, config=None):
106 if config.refcat_sharding_type != "tract":
107 if config.refcat_sharding_type == "none":
108 old = self.cat_ref
109 self.cat_ref = cT.Input(
110 doc=old.doc,
111 name=old.name,
112 storageClass=old.storageClass,
113 dimensions=(),
114 deferLoad=old.deferLoad,
115 )
116 else:
117 raise NotImplementedError(f"{config.refcat_sharding_type=} not implemented")
118 if config.target_sharding_type != "tract":
119 if config.target_sharding_type == "none":
120 old = self.cat_target
121 self.cat_target = cT.Input(
122 doc=old.doc,
123 name=old.name,
124 storageClass=old.storageClass,
125 dimensions=(),
126 deferLoad=old.deferLoad,
127 )
128 else:
129 raise NotImplementedError(f"{config.target_sharding_type=} not implemented")
132class MatchedCatalogFluxesConfig(pexConfig.Config):
133 column_ref_flux = pexConfig.Field(
134 dtype=str,
135 doc='Reference catalog flux column name',
136 )
137 columns_target_flux = pexConfig.ListField(
138 dtype=str,
139 listCheck=is_sequence_set,
140 doc="List of target catalog flux column names",
141 )
142 columns_target_flux_err = pexConfig.ListField(
143 dtype=str,
144 listCheck=is_sequence_set,
145 doc="List of target catalog flux error column names",
146 )
148 # this should be an orderedset
149 @property
150 def columns_in_ref(self) -> list[str]:
151 return [self.column_ref_flux]
153 # this should also be an orderedset
154 @property
155 def columns_in_target(self) -> list[str]:
156 columns = [col for col in self.columns_target_flux]
157 columns.extend(col for col in self.columns_target_flux_err if col not in columns)
158 return columns
161class DiffMatchedTractCatalogBaseConfig(pexConfig.Config):
162 column_match_candidate_ref = pexConfig.Field[str](
163 default='match_candidate',
164 doc='The column name for the boolean field identifying reference objects'
165 ' that were used for matching',
166 optional=True,
167 )
168 column_match_candidate_target = pexConfig.Field[str](
169 default='match_candidate',
170 doc='The column name for the boolean field identifying target objects'
171 ' that were used for matching',
172 optional=True,
173 )
174 column_matched_prefix_ref = pexConfig.Field[str](
175 default='refcat_',
176 doc='The prefix for matched columns copied from the reference catalog',
177 )
178 column_matched_prefix_target = pexConfig.Field[str](
179 default='',
180 doc='The prefix for matched columns copied from the target catalog',
181 )
182 include_unmatched = pexConfig.Field[bool](
183 default=False,
184 doc='Whether to include unmatched rows in the matched table',
185 )
186 filter_on_match_candidate = pexConfig.Field[bool](
187 default=False,
188 doc='Whether to use provided column_match_candidate_[ref/target] to'
189 ' exclude rows from the output table. If False, any provided'
190 ' columns will be copied instead.'
191 )
192 prefix_best_coord = pexConfig.Field[str](
193 default=None,
194 doc="A string prefix for ra/dec coordinate columns generated from the reference coordinate if "
195 "available, and target otherwise. Ignored if None or include_unmatched is False.",
196 optional=True,
197 )
199 @property
200 def columns_in_ref(self) -> list[str]:
201 columns_all = [self.coord_format.column_ref_coord1, self.coord_format.column_ref_coord2]
202 for column_lists in (
203 (
204 self.columns_ref_copy,
205 ),
206 (x.columns_in_ref for x in self.columns_flux.values()),
207 ):
208 for column_list in column_lists:
209 columns_all.extend(column_list)
211 return list({column: None for column in columns_all}.keys())
213 @property
214 def columns_in_target(self) -> list[str]:
215 columns_all = [self.coord_format.column_target_coord1, self.coord_format.column_target_coord2]
216 if self.coord_format.coords_ref_to_convert is not None: 216 ↛ 217line 216 didn't jump to line 217 because the condition on line 216 was never true
217 columns_all.extend(col for col in self.coord_format.coords_ref_to_convert.values()
218 if col not in columns_all)
219 for column_lists in (
220 (
221 self.columns_target_coord_err,
222 self.columns_target_select_false,
223 self.columns_target_select_true,
224 self.columns_target_copy,
225 ),
226 (x.columns_in_target for x in self.columns_flux.values()),
227 ):
228 for column_list in column_lists:
229 columns_all.extend(col for col in column_list if col not in columns_all)
230 return columns_all
232 columns_flux = pexConfig.ConfigDictField(
233 doc="Configs for flux columns for each band",
234 keytype=str,
235 itemtype=MatchedCatalogFluxesConfig,
236 default={},
237 )
238 columns_ref_mag_to_nJy = pexConfig.DictField[str, str](
239 doc='Reference table AB mag columns to convert to nJy flux columns with new names',
240 default={},
241 )
242 columns_ref_copy = pexConfig.ListField[str](
243 doc='Reference table columns to copy into cat_matched',
244 default=[],
245 listCheck=is_sequence_set,
246 )
247 columns_target_coord_err = pexConfig.ListField[str](
248 doc='Target table coordinate columns with standard errors (sigma)',
249 default=["coord_raErr", "coord_decErr"],
250 listCheck=lambda x: (len(x) == 2) and (x[0] != x[1]),
251 )
252 columns_target_copy = pexConfig.ListField[str](
253 doc='Target table columns to copy into cat_matched',
254 default=('patch',),
255 listCheck=is_sequence_set,
256 )
257 columns_target_mag_to_nJy = pexConfig.DictField[str, str](
258 doc='Target table AB mag columns to convert to nJy flux columns with new names',
259 default={},
260 )
261 columns_target_select_true = pexConfig.ListField[str](
262 doc='Target table columns to require to be True for selecting sources',
263 default=[],
264 listCheck=is_sequence_set,
265 )
266 columns_target_select_false = pexConfig.ListField[str](
267 doc='Target table columns to require to be False for selecting sources',
268 default=[],
269 listCheck=is_sequence_set,
270 )
271 coord_format = pexConfig.ConfigField[ConvertCatalogCoordinatesConfig](
272 doc="Configuration for coordinate conversion",
273 )
274 refcat_sharding_type = pexConfig.ChoiceField[str](
275 doc="The type of sharding (spatial splitting) for the reference catalog",
276 allowed={"tract": "Tract-based shards", "none": "No sharding at all"},
277 default="tract",
278 )
279 target_sharding_type = pexConfig.ChoiceField[str](
280 doc="The type of sharding (spatial splitting) for the target catalog",
281 allowed={"tract": "Tract-based shards", "none": "No sharding at all"},
282 default="tract",
283 )
285 def validate(self):
286 super().validate()
288 errors = []
290 for columns_mag, columns_in, name_columns_copy in (
291 (self.columns_ref_mag_to_nJy, self.columns_in_ref, "columns_ref_copy"),
292 (self.columns_target_mag_to_nJy, self.columns_in_target, "columns_target_copy"),
293 ):
294 columns_copy = getattr(self, name_columns_copy)
295 for column_old, column_new in columns_mag.items(): 295 ↛ 296line 295 didn't jump to line 296 because the loop on line 295 never started
296 if column_old not in columns_in:
297 errors.append(
298 f"{column_old=} key in self.columns_mag_to_nJy not found in {columns_in=}; did you"
299 f" forget to add it to self.{name_columns_copy}={columns_copy}?"
300 )
301 if column_new in columns_copy:
302 errors.append(
303 f"{column_new=} value found in self.{name_columns_copy}={columns_copy}"
304 f" this will cause a collision. Please choose a different name."
305 )
306 if errors: 306 ↛ 307line 306 didn't jump to line 307 because the condition on line 306 was never true
307 raise ValueError("\n".join(errors))
310class DiffMatchedTractCatalogConfig(
311 DiffMatchedTractCatalogBaseConfig,
312 pipeBase.PipelineTaskConfig,
313 pipelineConnections=DiffMatchedTractCatalogConnections,
314):
315 refcat_sharding_type = pexConfig.ChoiceField[str](
316 doc="The type of sharding (spatial splitting) for the reference catalog",
317 allowed={"tract": "Tract-based shards", "none": "No sharding at all"},
318 default="tract",
319 )
320 target_sharding_type = pexConfig.ChoiceField[str](
321 doc="The type of sharding (spatial splitting) for the target catalog",
322 allowed={"tract": "Tract-based shards", "none": "No sharding at all"},
323 default="tract",
324 )
327class DiffMatchedTractCatalogTaskBase(pipeBase.Task):
328 """Load subsets of matched catalogs and output a merged catalog of matched sources.
329 """
330 ConfigClass = DiffMatchedTractCatalogBaseConfig
332 def run(
333 self,
334 catalog_ref: astropy.table.Table,
335 catalog_target: astropy.table.Table,
336 catalog_match_ref: astropy.table.Table,
337 catalog_match_target: astropy.table.Table,
338 wcs: afwGeom.SkyWcs = None,
339 ) -> pipeBase.Struct:
340 """Load matched reference and target (measured) catalogs, measure summary statistics, and output
341 a combined matched catalog with columns from both inputs.
343 Parameters
344 ----------
345 catalog_ref : `astropy.table.Table`
346 A reference catalog to diff objects/sources from.
347 catalog_target : `astropy.table.Table`
348 A target catalog to diff reference objects/sources to.
349 catalog_match_ref : `astropy.table.Table`
350 A catalog with match indices of target sources and selection flags
351 for each reference source.
352 catalog_match_target : `astropy.table.Table`
353 A catalog with selection flags for each target source.
354 wcs : `lsst.afw.image.SkyWcs`
355 A coordinate system to convert catalog positions to sky coordinates,
356 if necessary.
358 Returns
359 -------
360 retStruct : `lsst.pipe.base.Struct`
361 A struct with output_ref and output_target attribute containing the
362 output matched catalogs.
363 """
364 # Would be nice if this could refer directly to ConfigClass
365 config: DiffMatchedTractCatalogConfig = self.config
367 # Strip any provenance from tables before merging to prevent
368 # warnings from conflicts being issued by astropy.utils.merge during
369 # vstack or hstack calls.
370 DatasetProvenance.strip_provenance_from_flat_dict(catalog_ref.meta)
371 DatasetProvenance.strip_provenance_from_flat_dict(catalog_target.meta)
372 DatasetProvenance.strip_provenance_from_flat_dict(catalog_match_ref.meta)
373 DatasetProvenance.strip_provenance_from_flat_dict(catalog_match_target.meta)
375 # It would be nice to make this a Selector but those are
376 # only available in analysis_tools for now
377 select_ref, select_target = (
378 (catalog[column] if column else np.ones(len(catalog), dtype=bool))
379 for catalog, column in (
380 (catalog_match_ref, self.config.column_match_candidate_ref),
381 (catalog_match_target, self.config.column_match_candidate_target),
382 )
383 )
384 # Add additional selection criteria for target sources beyond those for matching
385 # (not recommended, but can be done anyway)
386 for column in config.columns_target_select_true:
387 select_target &= catalog_target[column]
388 for column in config.columns_target_select_false:
389 select_target &= ~catalog_target[column]
391 ref, target = config.coord_format.format_catalogs(
392 catalog_ref=catalog_ref, catalog_target=catalog_target,
393 select_ref=None, select_target=select_target, wcs=wcs, radec_to_xy_func=radec_to_xy,
394 )
395 cat_ref = ref.catalog
396 cat_target = target.catalog
397 n_target = len(cat_target)
399 if not config.filter_on_match_candidate: 399 ↛ 407line 399 didn't jump to line 407 because the condition on line 399 was always true
400 for cat_add, cat_match, column in (
401 (cat_ref, catalog_match_ref, config.column_match_candidate_ref),
402 (cat_target, catalog_match_target, config.column_match_candidate_target),
403 ):
404 if column is not None: 404 ↛ 400line 404 didn't jump to line 400 because the condition on line 404 was always true
405 cat_add[column] = cat_match[column]
407 match_row = catalog_match_ref['match_row']
408 matched_ref = match_row >= 0
409 matched_row = match_row[matched_ref]
410 matched_target = np.zeros(n_target, dtype=bool)
411 matched_target[matched_row] = True
413 # Add/compute distance columns
414 coord1_target_err, coord2_target_err = config.columns_target_coord_err
415 column_dist, column_dist_err = 'match_distance', 'match_distanceErr'
416 dist = np.full(n_target, np.nan)
418 target_match_c1, target_match_c2 = (
419 np.array(coord[matched_row]) for coord in (target.coord1, target.coord2)
420 )
421 target_ref_c1, target_ref_c2 = (np.array(coord[matched_ref]) for coord in (ref.coord1, ref.coord2))
423 dist_err = np.full(n_target, np.nan)
424 dist[matched_row] = sphdist(
425 target_match_c1, target_match_c2, target_ref_c1, target_ref_c2
426 ) if config.coord_format.coords_spherical else np.hypot(
427 target_match_c1 - target_ref_c1, target_match_c2 - target_ref_c2,
428 )
429 cat_target_matched = cat_target[matched_row]
430 # This will convert a masked array to an array filled with nans
431 # wherever there are bad values (otherwise sphdist can raise)
432 c1_err, c2_err = (
433 np.ma.getdata(cat_target_matched[c_err]) for c_err in (coord1_target_err, coord2_target_err)
434 )
435 # Should probably explicitly add cosine terms if ref has errors too
436 dist_err[matched_row] = sphdist(
437 target_match_c1, target_match_c2, target_match_c1 + c1_err, target_match_c2 + c2_err
438 ) if config.coord_format.coords_spherical else np.hypot(c1_err, c2_err)
439 cat_target[column_dist], cat_target[column_dist_err] = dist, dist_err
441 # Create a matched table, preserving the target catalog's named index (if it has one)
442 cat_left = cat_target[matched_row]
443 cat_right = cat_ref[matched_ref]
444 if config.column_matched_prefix_target: 444 ↛ 445line 444 didn't jump to line 445 because the condition on line 444 was never true
445 cat_left.rename_columns(
446 list(cat_left.columns),
447 new_names=[f'{config.column_matched_prefix_target}{col}' for col in cat_left.columns],
448 )
449 if config.column_matched_prefix_ref: 449 ↛ 454line 449 didn't jump to line 454 because the condition on line 449 was always true
450 cat_right.rename_columns(
451 list(cat_right.columns),
452 new_names=[f'{config.column_matched_prefix_ref}{col}' for col in cat_right.columns],
453 )
454 cat_matched = astropy.table.hstack((cat_left, cat_right))
456 if config.include_unmatched: 456 ↛ 460line 456 didn't jump to line 460 because the condition on line 456 was never true
457 # Create an unmatched table with the same schema as the matched one
458 # ... but only for objects with no matches (for completeness/purity)
459 # and that were selected for matching (or inclusion via config)
460 cat_right = astropy.table.Table(
461 cat_ref[~matched_ref & select_ref]
462 )
463 cat_right.rename_columns(
464 cat_right.colnames,
465 [f"{config.column_matched_prefix_ref}{col}" for col in cat_right.colnames],
466 )
467 match_row_target = catalog_match_target['match_row']
468 cat_left = cat_target[~(match_row_target >= 0) & select_target]
469 cat_left.rename_columns(
470 cat_left.colnames,
471 [f"{config.column_matched_prefix_target}{col}" for col in cat_left.colnames],
472 )
473 # This may be slower than pandas but will, for example, create
474 # masked columns for booleans, which pandas does not support.
475 # See https://github.com/pandas-dev/pandas/issues/46662
476 cat_unmatched = astropy.table.vstack([cat_left, cat_right])
478 for columns_convert_base, prefix in (
479 (config.columns_ref_mag_to_nJy, config.column_matched_prefix_ref),
480 (config.columns_target_mag_to_nJy, ""),
481 ):
482 if columns_convert_base: 482 ↛ 483line 482 didn't jump to line 483 because the condition on line 482 was never true
483 columns_convert = {
484 f"{prefix}{k}": f"{prefix}{v}" for k, v in columns_convert_base.items()
485 } if prefix else columns_convert_base
486 to_convert = [cat_matched]
487 if config.include_unmatched:
488 to_convert.append(cat_unmatched)
489 for cat_convert in to_convert:
490 cat_convert.rename_columns(
491 tuple(columns_convert.keys()),
492 tuple(columns_convert.values()),
493 )
494 for column_flux in columns_convert.values():
495 cat_convert[column_flux] = u.ABmag.to(u.nJy, cat_convert[column_flux])
497 if config.include_unmatched: 497 ↛ 499line 497 didn't jump to line 499 because the condition on line 497 was never true
498 # This is probably less efficient than just doing an outer join originally; worth checking
499 cat_matched = astropy.table.vstack([cat_matched, cat_unmatched])
500 if (prefix_coord := config.prefix_best_coord) is not None:
501 columns_coord_best = (
502 f"{prefix_coord}{col_coord}" for col_coord in (
503 ("ra", "dec") if config.coord_format.coords_spherical else ("coord1", "coord2")
504 )
505 )
506 for column_coord_best, column_coord_ref, column_coord_target in zip(
507 columns_coord_best,
508 (config.coord_format.column_ref_coord1, config.coord_format.column_ref_coord2),
509 (config.coord_format.column_target_coord1, config.coord_format.column_target_coord2),
510 ):
511 column_full_ref = f'{config.column_matched_prefix_ref}{column_coord_ref}'
512 column_full_target = f'{config.column_matched_prefix_target}{column_coord_target}'
513 values = cat_matched[column_full_ref]
514 unit = values.unit
515 values_bad = np.ma.masked_invalid(values).mask
516 # Cast to an unmasked array - there will be no bad values
517 values = np.array(values)
518 values[values_bad] = cat_matched[column_full_target][values_bad]
519 cat_matched[column_coord_best] = values
520 cat_matched[column_coord_best].unit = unit
521 cat_matched[column_coord_best].description = (
522 f"Best {column_coord_best} value from {column_full_ref} if available"
523 f" else {column_full_target}"
524 )
526 retStruct = pipeBase.Struct(cat_matched=cat_matched)
527 return retStruct
530class DiffMatchedTractCatalogTask(DiffMatchedTractCatalogTaskBase, pipeBase.PipelineTask):
531 """PipelineTask version of DiffMatchedTractCatalogTaskBase."""
532 ConfigClass = DiffMatchedTractCatalogConfig
533 _DefaultName = "DiffMatchedTractCatalog"
535 def runQuantum(self, butlerQC, inputRefs, outputRefs):
536 inputs = butlerQC.get(inputRefs)
537 skymap = inputs.pop("skymap")
539 columns_match_ref = ['match_row']
540 if (column := self.config.column_match_candidate_ref) is not None:
541 columns_match_ref.append(column)
543 columns_match_target = ['match_row']
544 if (column := self.config.column_match_candidate_target) is not None and (
545 column in inputs['columns_match_target']
546 ):
547 columns_match_target.append(column)
549 outputs = self.run(
550 catalog_ref=inputs['cat_ref'].get(parameters={'columns': self.config.columns_in_ref}),
551 catalog_target=inputs['cat_target'].get(parameters={'columns': self.config.columns_in_target}),
552 catalog_match_ref=inputs['cat_match_ref'].get(parameters={'columns': columns_match_ref}),
553 catalog_match_target=inputs['cat_match_target'].get(parameters={'columns': columns_match_target}),
554 wcs=skymap[butlerQC.quantum.dataId["tract"]].wcs,
555 )
556 butlerQC.put(outputs, outputRefs)