Coverage for python/lsst/drp/tasks/metadetection_shear.py: 19%
277 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-04 01:05 -0700
1# This file is part of drp_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/>.
22from __future__ import annotations
24__all__ = (
25 "MetadetectionProcessingError",
26 "MetadetectionShearConfig",
27 "MetadetectionShearTask",
28)
30from collections.abc import Collection, Mapping, Sequence
31from itertools import product
32from typing import Any, ClassVar
34import esutil as eu
35import numpy as np
36import pyarrow as pa
37from metadetect.lsst.masking import apply_apodized_bright_masks_mbexp, apply_apodized_edge_masks_mbexp
38from metadetect.lsst.metacal_exposures import STEP as SHEAR_STEP
39from metadetect.lsst.metadetect import MetadetectTask
40from metadetect.lsst.util import extract_multiband_coadd_data
42import lsst.geom as geom
43import lsst.pipe.base.connectionTypes as cT
44from lsst.afw.image import ExposureF
45from lsst.afw.table import SimpleCatalog
46from lsst.cell_coadds import MultipleCellCoadd, StitchedCoadd
47from lsst.daf.butler import DataCoordinate, DatasetRef
48from lsst.meas.algorithms import LoadReferenceObjectsConfig, ReferenceObjectLoader
49from lsst.meas.base import FullIdGenerator, SkyMapIdGeneratorConfig
50from lsst.pex.config import ConfigField, ConfigurableField, Field, FieldValidationError, ListField
51from lsst.pipe.base import (
52 AlgorithmError,
53 AnnotatedPartialOutputsError,
54 InputQuantizedConnection,
55 InvalidQuantumError,
56 NoWorkFound,
57 OutputQuantizedConnection,
58 PipelineTask,
59 PipelineTaskConfig,
60 PipelineTaskConnections,
61 QuantumContext,
62 Struct,
63)
64from lsst.pipe.base.connectionTypes import BaseInput, Output
65from lsst.pipe.tasks.coaddBase import makeSkyInfo
66from lsst.skymap import BaseSkyMap, Index2D
69class MetadetectionProcessingError(AlgorithmError):
70 """Exception raised when metadetection processing fails."""
72 @property
73 def metadata(self) -> dict:
74 return {}
77class MetadetectionShearConnections(PipelineTaskConnections, dimensions={"patch"}):
78 """Definitions of inputs and outputs for MetadetectionShearTask."""
80 input_coadds = cT.Input(
81 "deep_coadd_cell_predetection",
82 storageClass="MultipleCellCoadd",
83 doc="Per-band deep coadds.",
84 multiple=True,
85 dimensions={"patch", "band"},
86 )
88 sky_map = cT.Input(
89 doc="Cell-based skymap defining the patch structure.",
90 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
91 storageClass="SkyMap",
92 dimensions=("skymap",),
93 )
95 ref_cat = cT.PrerequisiteInput(
96 doc="Reference catalog used to mask bright objects.",
97 name="the_monster_20250219",
98 storageClass="SimpleCatalog",
99 dimensions=("skypix",),
100 deferLoad=True,
101 multiple=True,
102 )
104 metadetect_catalog = cT.Output(
105 "object_shear_patch",
106 storageClass="ArrowTable",
107 doc="Output catalog with all quantities measured inside the metacalibration loop.",
108 multiple=False,
109 dimensions={"patch"},
110 )
112 metadetect_schema = cT.InitOutput(
113 "object_shear_schema",
114 storageClass="ArrowSchema",
115 doc="Schema of the output catalog.",
116 )
118 config: MetadetectionShearConfig
120 def __init__(self, *, config=None):
121 super().__init__(config=config)
123 if not config:
124 return None
126 if not config.do_mask_bright_objects:
127 del self.ref_cat
129 def adjustQuantum(
130 self,
131 inputs: dict[str, tuple[BaseInput, Collection[DatasetRef]]],
132 outputs: dict[str, tuple[Output, Collection[DatasetRef]]],
133 label: str,
134 data_id: DataCoordinate,
135 ) -> tuple[
136 Mapping[str, tuple[BaseInput, Collection[DatasetRef]]],
137 Mapping[str, tuple[Output, Collection[DatasetRef]]],
138 ]:
139 # Docstring inherited.
140 # This is a hook for customizing what is input and output to each
141 # invocation of the task as early as possible, which we override here
142 # to make sure we have exactly the required bands, no more, no less.
143 connection, original_input_coadds = inputs["input_coadds"]
144 bands_missing = set(self.config.photometry_bands)
145 adjusted_input_coadds = []
146 for ref in original_input_coadds:
147 if ref.dataId["band"] in self.config.photometry_bands:
148 adjusted_input_coadds.append(ref)
149 bands_missing.remove(ref.dataId["band"])
150 if missing_shear_bands := bands_missing.intersection(self.config.metadetect.shear_bands):
151 raise NoWorkFound(f"Required bands {missing_shear_bands} not present for {label}@{data_id}).")
152 adjusted_inputs = {"input_coadds": (connection, adjusted_input_coadds)}
153 inputs.update(adjusted_inputs)
154 super().adjustQuantum(inputs, outputs, label, data_id)
155 return adjusted_inputs, {}
158class MetadetectionShearConfig(PipelineTaskConfig, pipelineConnections=MetadetectionShearConnections):
159 """Configuration definition for MetadetectionShearTask."""
161 metadetect = ConfigurableField(
162 target=MetadetectTask,
163 doc="Configuration for metadetection.",
164 )
166 photometry_bands = ListField[str](
167 "Bands expected to be present. Cells with one or more of these bands "
168 "missing will be skipped. Bands other than those listed here will "
169 "not be processed.",
170 default=["g", "r", "i", "z"],
171 )
173 do_mask_bright_objects = Field[bool](
174 doc="Mask bright objects in coadds?",
175 default=False,
176 )
178 ref_loader = ConfigField(
179 dtype=LoadReferenceObjectsConfig,
180 doc="Reference object loader used for bright-object masking.",
181 )
183 ref_loader_filter_name = Field[str](
184 "Filter name from ref_loader used for bright-object masking.",
185 default="monster_DES_r",
186 )
188 border = Field[int](
189 "Border to apply to single cell images, if skymap has no cell borders",
190 default=50,
191 )
193 id_generator = SkyMapIdGeneratorConfig.make_field()
195 def setDefaults(self):
196 super().setDefaults()
197 self.metadetect.shear_bands = ["r", "i", "z"]
198 self.metadetect.metacal.types = ["noshear", "1p", "1m", "2p", "2m"]
200 def validate(self):
201 super().validate()
202 if (shear_bands := self.metadetect.shear_bands) is not None and not set(shear_bands).issubset(
203 self.photometry_bands
204 ):
205 raise FieldValidationError(
206 self.__class__.metadetect,
207 self,
208 "photometry_bands must be a list of bands that is a superset of metadetect.shear_bands",
209 )
212class MetadetectionShearTask(PipelineTask):
213 """A PipelineTask that measures shear using metadetection."""
215 _DefaultName: ClassVar[str] = "metadetectionShear"
216 ConfigClass: ClassVar[type[MetadetectionShearConfig]] = MetadetectionShearConfig
218 config: MetadetectionShearConfig
220 def __init__(self, *, initInputs: dict[str, Any] | None = None, **kwargs: Any):
221 super().__init__(initInputs=initInputs, **kwargs)
222 self.metadetect_schema = self.make_metadetect_schema(self.config)
223 self.makeSubtask("metadetect")
225 @classmethod
226 def make_metadetect_schema(cls, config: MetadetectionShearConfig) -> pa.Schema:
227 """Construct a PyArrow Schema for this task's main output catalog.
229 Parameters
230 ----------
231 config : `MetadetectionShearConfig`
232 Configuration that may be used to control details of the schema.
234 Returns
235 -------
236 object_schema : `pyarrow.Schema`
237 Schema for the object catalog produced by this task. Each field's
238 metadata should include both a 'doc' entry and a 'unit' entry.
239 """
240 pa_schema = pa.schema(
241 [
242 # Fields from pipeline bookkeeping.
243 pa.field(
244 "shearObjectId",
245 pa.int64(),
246 nullable=False,
247 metadata={
248 "doc": (
249 "Unique identifier for a ShearObject, specific "
250 "to a single metacalibration counterfactual image."
251 ),
252 "unit": "",
253 },
254 ),
255 pa.field(
256 "tract",
257 pa.int64(),
258 nullable=False,
259 metadata={
260 "doc": "ID of the tract on which this measurement was made.",
261 "unit": "",
262 },
263 ),
264 pa.field(
265 "patch",
266 pa.int64(),
267 nullable=False,
268 metadata={
269 "doc": "ID of the patch within the tract on which this measurement was made.",
270 "unit": "",
271 },
272 ),
273 pa.field(
274 "cell_x",
275 pa.int32(),
276 nullable=False,
277 metadata={
278 "doc": "Column of the cell within the patch on which this measurement was made.",
279 "unit": "",
280 },
281 ),
282 pa.field(
283 "cell_y",
284 pa.int32(),
285 nullable=False,
286 metadata={
287 "doc": "Row of the cell within the patch on which this measurement was made.",
288 "unit": "",
289 },
290 ),
291 pa.field(
292 "is_cell_inner",
293 pa.bool_(),
294 nullable=False,
295 metadata={
296 "doc": "Whether this object is within the inner region of the cell.",
297 "unit": "",
298 },
299 ),
300 pa.field(
301 "is_patch_inner",
302 pa.bool_(),
303 nullable=False,
304 metadata={
305 "doc": "Whether this object is within the inner region of the patch.",
306 "unit": "",
307 },
308 ),
309 pa.field(
310 "is_tract_inner",
311 pa.bool_(),
312 nullable=False,
313 metadata={
314 "doc": "Whether this object is within the inner region of the tract.",
315 "unit": "",
316 },
317 ),
318 pa.field(
319 "is_primary",
320 pa.bool_(),
321 nullable=False,
322 metadata={
323 "doc": "Whether this object is to be selected after de-duplication",
324 "unit": "",
325 },
326 ),
327 # Fields from metadetection (generic).
328 pa.field(
329 "metaStep",
330 pa.string(),
331 nullable=False,
332 metadata={
333 "doc": (
334 "Type of artificial shear applied to image. "
335 "One of: 'ns', '1p', '1m', '2p', '2m'."
336 ),
337 "unit": "",
338 },
339 ),
340 pa.field(
341 "image_flags",
342 pa.int32(),
343 nullable=False,
344 metadata={
345 "doc": "Flags for the image on which this measurement was made.",
346 "unit": "",
347 },
348 ),
349 pa.field(
350 "x",
351 pa.float32(),
352 nullable=False,
353 metadata={
354 "doc": "Centroid (tract, x-axis) of the detected ShearObject.",
355 "unit": "",
356 },
357 ),
358 pa.field(
359 "y",
360 pa.float32(),
361 nullable=False,
362 metadata={
363 "doc": "Centroid (tract, y-axis) of the detected ShearObject.",
364 "unit": "",
365 },
366 ),
367 pa.field(
368 "ra",
369 pa.float64(),
370 nullable=False,
371 metadata={
372 "doc": "Detected Right Ascension of the ShearObject.",
373 "unit": "degrees",
374 },
375 ),
376 pa.field(
377 "dec",
378 pa.float64(),
379 nullable=False,
380 metadata={
381 "doc": "Detected Declination of the ShearObject.",
382 "unit": "degrees",
383 },
384 ),
385 # Original PSF measurements
386 pa.field(
387 "psfOriginal_flags",
388 pa.int32(),
389 nullable=False,
390 metadata={
391 "doc": "Flags for the original PSF measurement.",
392 "unit": "",
393 },
394 ),
395 pa.field(
396 "psfOriginal_e1",
397 pa.float32(),
398 nullable=False,
399 metadata={
400 "doc": "Distortion-style e1 of the original PSF from adaptive moments.",
401 "unit": "",
402 },
403 ),
404 pa.field(
405 "psfOriginal_e2",
406 pa.float32(),
407 nullable=False,
408 metadata={
409 "doc": "Distortion-style e2 of the original PSF from adaptive moments.",
410 "unit": "",
411 },
412 ),
413 pa.field(
414 "psfOriginal_T",
415 pa.float32(),
416 nullable=False,
417 metadata={
418 "doc": "Trace (<x^2> + <y^2>) measurement of the original PSF from adaptive moments.",
419 "unit": "arcseconds squared",
420 },
421 ),
422 pa.field(
423 "bmask_flags",
424 pa.int32(),
425 nullable=False,
426 metadata={
427 "doc": "`bmask` flags for the ShearObject",
428 "unit": "",
429 },
430 ),
431 pa.field(
432 "ormask_flags",
433 pa.int32(),
434 nullable=False,
435 metadata={
436 "doc": "`ored` mask flags for the ShearObject",
437 "unit": "",
438 },
439 ),
440 pa.field(
441 "mfrac",
442 pa.float32(),
443 nullable=False,
444 metadata={
445 "doc": "Gaussian-weighted masked fraction for the ShearObject.",
446 "unit": "",
447 },
448 ),
449 # Fields that come only from gauss algorithm.
450 # Reconvolved PSF measurements (gauss)
451 pa.field(
452 "gauss_psfReconvolved_flags",
453 pa.int32(),
454 nullable=False,
455 metadata={
456 "doc": "Flags for reconvolved PSF (measured with gauss algorithm).",
457 "unit": "",
458 },
459 ),
460 pa.field(
461 "gauss_psfReconvolved_g1",
462 pa.float32(),
463 nullable=False,
464 metadata={
465 "doc": "Reduced-shear g1 of the reconvolved PSF (measured with gauss algorithm).",
466 "unit": "",
467 },
468 ),
469 pa.field(
470 "gauss_psfReconvolved_g2",
471 pa.float32(),
472 nullable=False,
473 metadata={
474 "doc": "Reduced-shear g2 of the reconvolved PSF (measured with gauss algorithm).",
475 "unit": "",
476 },
477 ),
478 pa.field(
479 "gauss_psfReconvolved_T",
480 pa.float32(),
481 nullable=False,
482 metadata={
483 "doc": (
484 "Trace (<x^2> + <y^2>) of the reconvolved PSF (measured with gauss algorithm)."
485 ),
486 "unit": "arcseconds squared",
487 },
488 ),
489 # Object measurements (gauss algorithm).
490 pa.field(
491 "gauss_g1",
492 pa.float32(),
493 nullable=False,
494 metadata={
495 "doc": (
496 "Reduced-shear g1 measurement of the ShearObject "
497 "(measured with gauss algorithm)."
498 ),
499 "unit": "",
500 },
501 ),
502 pa.field(
503 "gauss_g2",
504 pa.float32(),
505 nullable=False,
506 metadata={
507 "doc": (
508 "Reduced-shear g2 measurement of the ShearObject "
509 "(measured with gauss algorithm)."
510 ),
511 "unit": "",
512 },
513 ),
514 pa.field(
515 "gauss_g1_g1_Cov",
516 pa.float32(),
517 nullable=False,
518 metadata={
519 "doc": (
520 "Auto-covariance of g1 measurement of the ShearObject "
521 "(measured with gauss algorithm)."
522 ),
523 "unit": "",
524 },
525 ),
526 pa.field(
527 "gauss_g1_g2_Cov",
528 pa.float32(),
529 nullable=False,
530 metadata={
531 "doc": (
532 "Cross-covariance of g1 and g2 measurement of the ShearObject "
533 "(measured with gauss algorithm)."
534 ),
535 "unit": "",
536 },
537 ),
538 pa.field(
539 "gauss_g2_g2_Cov",
540 pa.float32(),
541 nullable=False,
542 metadata={
543 "doc": (
544 "Auto-covariance of g2 measurement of the ShearObject "
545 "(measured with gauss algorithm)."
546 ),
547 "unit": "",
548 },
549 ),
550 ],
551 metadata={
552 "shear_step": str(SHEAR_STEP),
553 "shear_bands": "".join(sorted(config.metadetect.shear_bands)),
554 },
555 )
557 for alg_name in ("gauss", "pgauss"):
558 pa_schema = pa_schema.append(
559 pa.field(
560 f"{alg_name}_snr",
561 pa.float32(),
562 nullable=False,
563 metadata={
564 "doc": (
565 "Signal-to-noise ratio measure of the ShearObject "
566 f"(measured with {alg_name} algorithm)."
567 ),
568 "unit": "",
569 },
570 ),
571 )
572 pa_schema = pa_schema.append(
573 pa.field(
574 f"{alg_name}_T",
575 pa.float32(),
576 nullable=False,
577 metadata={
578 "doc": (
579 "Trace (<x^2> + <y^2>) measurement of the ShearObject "
580 f"(measured with {alg_name} algorithm)."
581 ),
582 "unit": "arcseconds squared",
583 },
584 ),
585 )
586 pa_schema = pa_schema.append(
587 pa.field(
588 f"{alg_name}_TErr",
589 pa.float32(),
590 nullable=False,
591 metadata={
592 "doc": (
593 "Uncertainty in the trace measurement of the ShearObject "
594 f"(measured with {alg_name} algorithm)."
595 ),
596 "unit": "arcseconds squared",
597 },
598 ),
599 )
600 pa_schema = pa_schema.append(
601 pa.field(
602 f"{alg_name}_shape_flags",
603 pa.int32(),
604 nullable=False,
605 metadata={
606 "doc": (
607 "Flags for the second order moments measurement of the ShearObject "
608 f"(measured with {alg_name} algorithm)."
609 ),
610 "unit": "",
611 },
612 ),
613 )
614 pa_schema = pa_schema.append(
615 pa.field(
616 f"{alg_name}_object_flags",
617 pa.int32(),
618 nullable=False,
619 metadata={
620 "doc": f"Flags for the ShearObject measurement (measured with {alg_name} algorithm).",
621 "unit": "",
622 },
623 ),
624 )
625 pa_schema = pa_schema.append(
626 pa.field(
627 f"{alg_name}_flags",
628 pa.int32(),
629 nullable=False,
630 metadata={
631 "doc": f"Overall flags for {alg_name} measurement algorithm.",
632 "unit": "",
633 },
634 ),
635 )
637 # Per-band quantities, typically fluxes and associated quantites.
638 for b in config.photometry_bands:
639 pa_schema = pa_schema.append(
640 pa.field(
641 f"{b}_{alg_name}Flux_flags",
642 pa.int32(),
643 nullable=False,
644 metadata={
645 "doc": f"Flags set for flux in {b} band measured with {alg_name} algorithm.",
646 "unit": "",
647 },
648 ),
649 )
650 pa_schema = pa_schema.append(
651 pa.field(
652 f"{b}_{alg_name}Flux",
653 pa.float32(),
654 nullable=b not in config.metadetect.shear_bands,
655 metadata={
656 "doc": f"Flux in {b} band (measured with {alg_name} algorithm).",
657 "unit": "",
658 },
659 ),
660 )
661 pa_schema = pa_schema.append(
662 pa.field(
663 f"{b}_{alg_name}FluxErr",
664 pa.float32(),
665 nullable=b not in config.metadetect.shear_bands,
666 metadata={
667 "doc": f"Flux uncertainty in {b} band (measured with {alg_name} algorithm).",
668 "unit": "",
669 },
670 ),
671 )
673 return pa_schema
675 def validate_skymap_config(self, skymap_config: BaseSkyMap.ConfigClass) -> None:
676 if not skymap_config.tractBuilder.name == "cells":
677 raise InvalidQuantumError("MetadetectionShearTask requires a cell-based skymap.")
679 cell_config = skymap_config.tractBuilder.active
680 if (self.config.border == 0 and cell_config.cellBorder == 0) or (
681 self.config.border > 0 and cell_config.cellBorder > 0
682 ):
683 raise InvalidQuantumError(
684 "MetadetectionShearTask requires a positive border to be set either in the skymap config "
685 "or in the task config (but not in both)."
686 )
688 if self.config.border:
689 # In case cellInnerDimensions are different in different directions
690 # (which is rare in practice), take the min of it
691 cell_inner_dimensions = min(cell_config.cellInnerDimensions)
692 if self.config.border > (
693 max_border_value := cell_config.numCellsInPatchBorder * cell_inner_dimensions
694 ):
695 raise InvalidQuantumError(
696 "The border value is too large for the skymap configuration. "
697 f"Maximum border value is {max_border_value}."
698 )
700 # Ensure that the amount by which we withdraw inwards does not
701 # create gaps at tract boundaries.
702 # tractOverlap is specified in degrees.
703 if (
704 skymap_config.tractOverlap * 3600 / skymap_config.pixelScale
705 < self.count_cells_along_edges(skymap_config) * cell_inner_dimensions + cell_config.cellBorder
706 ):
707 raise InvalidQuantumError(
708 "The tract overlap is insufficient given the borders. "
709 "This will result in missed regions between adjacent tracts."
710 )
712 @staticmethod
713 def count_cells_along_edges(skymap_config: BaseSkyMap.ConfigClass) -> int:
714 """Count the number of cells along the edges to skip processing."""
715 if skymap_config.tractBuilder["cells"].cellBorder > 0:
716 return 0
718 return skymap_config.tractBuilder["cells"].numCellsInPatchBorder
720 def runQuantum(
721 self,
722 qc: QuantumContext,
723 inputRefs: InputQuantizedConnection,
724 outputRefs: OutputQuantizedConnection,
725 ) -> None:
726 # Docstring inherited.
728 # Get the skyMap for this quantum
729 sky_map = qc.get(inputRefs.sky_map)
731 self.validate_skymap_config(sky_map.config)
733 id_generator = self.config.id_generator.apply(qc.quantum.dataId)
735 if self.config.do_mask_bright_objects:
736 ref_loader = ReferenceObjectLoader(
737 dataIds=[ref.datasetRef.dataId for ref in inputRefs.ref_cat],
738 refCats=[qc.get(ref) for ref in inputRefs.ref_cat],
739 name=self.config.connections.ref_cat,
740 config=self.config.ref_loader,
741 log=self.log,
742 )
743 ref_cat = ref_loader.loadRegion(
744 qc.quantum.dataId.region, filterName=self.config.ref_loader_filter_name
745 )
746 else:
747 ref_cat = None
749 # Read the coadds and put them in the order defined by
750 # config.photometry_bands (note that each MultipleCellCoadd object also
751 # knows its own band, if that's needed).
753 coadds_by_band = {
754 ref.dataId["band"]: qc.get(ref)
755 for ref in inputRefs.input_coadds
756 if ref.dataId["band"] in self.config.photometry_bands
757 }
759 try:
760 outputs = self.run(
761 patch_coadds=coadds_by_band,
762 id_generator=id_generator,
763 sky_map=sky_map,
764 ref_cat=ref_cat,
765 )
766 except AlgorithmError as err:
767 # We know there are no actual outputs in this case, but this is
768 # still the right exception to raise (it's just badly named).
769 raise AnnotatedPartialOutputsError.annotate(err, self, log=self.log) from err
770 qc.put(outputs, outputRefs)
772 def run(
773 self,
774 *,
775 patch_coadds: Mapping[str, MultipleCellCoadd],
776 id_generator: FullIdGenerator,
777 sky_map: BaseSkyMap,
778 ref_cat: SimpleCatalog | None,
779 ) -> Struct:
780 """Run metadetection on a patch.
782 Parameters
783 ----------
784 patch_coadds : `~collections.abc.Mapping` [ \
785 `~lsst.cell_coadds.MultipleCellCoadd` ]
786 Per-band, per-patch coadds, in the order specified by
787 `MetadetectionShearConfig.photometry_bands`.
788 id_generator : `~lsst.meas.base.FullIdGenerator`
789 Generator for object IDs and to seed the random number generator.
790 sky_map : `~lsst.skymap.BaseSkyMap`
791 Sky map to use for determining the patch boundaries.
792 ref_cat : `lsst.afw.table.SimpleCatalog`, optional
793 Reference catalog to use when masking bright stars.
795 Returns
796 -------
797 results : `lsst.pipe.base.Struct`
798 Structure with the following attributes:
800 - ``metadetect_catalog`` [ `pyarrow.Table` ]: the output object
801 catalog for the patch, with schema equal to `metadetect_schema`.
802 """
803 seed = id_generator.catalog_id
804 self.rng = np.random.RandomState(seed)
805 idstart = 0
807 sky_info = makeSkyInfo(
808 sky_map,
809 tractId=id_generator.data_id.tract.id,
810 patchId=id_generator.data_id.patch.id,
811 )
812 # Since cells have no borders, we cannot apply the metacal
813 # procedure to the edge cells in the same way as inner cells.
814 # We can do so for all that cells that are marked as borders,
815 # but it has to be at least 1.
816 if (num_edge_cells_skip := sky_map.config.tractBuilder.active.numCellsInPatchBorder) < 1:
817 raise InvalidQuantumError("No border cells found in the skymap configuration.")
819 dilate_by = self.config.border or sky_map.config.tractBuilder.active.cellBorder
821 grid = patch_coadds[self.config.metadetect.shear_bands[0]].grid
822 # Undo any padding that was applied when creating the grid.
823 # We only support metadetection on non-overlapping single cells.
824 # lsst_cells_v1 has overlapping cells, so we need to undo the padding.
825 grid = grid.from_bbox_shape(bbox=grid.bbox, shape=grid.shape, padding=0)
827 nx_cells, ny_cells = grid.shape
828 single_cell_tables: list[pa.Table] = []
829 for nx, ny in product(
830 range(num_edge_cells_skip, nx_cells - num_edge_cells_skip),
831 range(num_edge_cells_skip, ny_cells - num_edge_cells_skip),
832 ):
833 cell_id = Index2D(nx, ny)
834 bbox = grid.bbox_of(cell_id).dilatedBy(dilate_by)
835 cell_coadds = [patch_coadd.stitch(bbox) for patch_coadd in patch_coadds.values()]
836 self.log.debug("Processing cell %s %s", nx, ny)
838 try:
839 res = self.process_cell(cell_coadds, cell_id=cell_id)
840 except Exception as e:
841 self.log.error("Failed to process cell %s %s: %s", nx, ny, e)
842 continue
844 if len(res) > 0:
845 res["id"] = id_generator.arange(idstart, idstart + len(res))
846 # TODO: Avoid back and forth conversion between array and dict.
847 da = self._dictify(
848 res,
849 tract=id_generator.data_id.tract.id,
850 patch=id_generator.data_id.patch.id,
851 )
853 self._fill_is_inner(da, sky_info.patchInfo[cell_id].inner_bbox, "is_cell_inner")
854 self._fill_is_inner(da, sky_info.patchInfo.inner_bbox, "is_patch_inner")
855 # Calculating is_tract_inner follows a different logic with
856 # regions instead of boxes.
857 da["is_tract_inner"] = sky_info.tractInfo.inner_sky_region.contains(
858 da["ra"] * geom.degrees,
859 da["dec"] * geom.degrees,
860 )
861 da["is_primary"] = da["is_cell_inner"] & da["is_patch_inner"] & da["is_tract_inner"]
863 table = pa.Table.from_pydict(da, self.metadetect_schema)
865 single_cell_tables.append(table)
866 idstart += len(res)
868 if not single_cell_tables:
869 raise MetadetectionProcessingError("No objects found in any cell")
871 # TODO: DM-53796 De-duplicate objects before concatenation.
872 return Struct(
873 metadetect_catalog=pa.concat_tables(single_cell_tables),
874 )
876 def process_cell(
877 self,
878 cell_coadds: Sequence[StitchedCoadd],
879 cell_id: Index2D,
880 ) -> pa.Table:
881 """Run metadetection on a single cell.
883 Parameters
884 ----------
885 cell_coadds : `~collections.abc.Sequence` [ \
886 `~lsst.cell_coadds.StitchedCoadd` ]
887 Per-band, per-cell coadds, in the order specified by
888 `MetadetectionShearConfig.photometry_bands`.
889 cell_id : `~lsst.skymap.Index2D`
890 The cell ID for the cell being processed.
892 Returns
893 -------
894 metadetect_catalog : `pyarrow.Table`
895 Output object catalog for the cell, with schema equal to
896 `metadetect_schema`.
897 """
899 coadd_data = self._cell_to_coadd_data(cell_coadds)
900 # TODO get bright star etc. info as input
901 bright_info = []
903 apply_apodized_edge_masks_mbexp(**coadd_data)
905 if len(bright_info) > 0:
906 apply_apodized_bright_masks_mbexp(bright_info=bright_info, **coadd_data)
908 mask_frac = _get_mask_frac(
909 coadd_data["mfrac_mbexp"],
910 trim_pixels=0,
911 )
913 res = self.metadetect.run(rng=self.rng, **coadd_data)
915 comb_res = _make_comb_data(
916 cell_coadd=cell_coadds[0],
917 res=res,
918 mask_frac=mask_frac,
919 bands=[cell_coadd.band for cell_coadd in cell_coadds],
920 cell_id=cell_id,
921 )
923 return comb_res
925 @staticmethod
926 def _cell_to_coadd_data(cell_coadds: Sequence[StitchedCoadd]):
927 coadd_data_list = []
928 for cell_coadd in cell_coadds:
929 coadd_data = {}
930 coadd_data["coadd_exp"] = cell_coadd.asExposure()
931 coadd_data["coadd_noise_exp"] = cell_coadd.asExposure(noise_index=0)
932 coadd_data["coadd_mfrac_exp"] = ExposureF(coadd_data["coadd_exp"], deep=True)
933 coadd_data["coadd_mfrac_exp"].image = cell_coadd.mask_fractions
934 coadd_data_list.append(coadd_data)
936 return extract_multiband_coadd_data(coadd_data_list)
938 def _dictify(self, data, tract: int, patch: int):
939 output = {}
940 # TODO: Move this to a better location after DP2.
941 mapping = {
942 "bmask": "bmask_flags",
943 "cell_x": "cell_x",
944 "cell_y": "cell_y",
945 "col": "x",
946 "col_diff": "x_offset", # dropped.
947 "dec": "dec",
948 "gauss_flags": "gauss_flags",
949 "gauss_g_1": "gauss_g1",
950 "gauss_g_2": "gauss_g2",
951 "gauss_g_cov_11": "gauss_g1_g1_Cov",
952 "gauss_g_cov_12": "gauss_g1_g2_Cov", # same as 21.
953 "gauss_g_cov_22": "gauss_g2_g2_Cov",
954 "gauss_obj_flags": "gauss_object_flags",
955 "gauss_psf_flags": "gauss_psfReconvolved_flags",
956 "gauss_psf_g_1": "gauss_psfReconvolved_g1",
957 "gauss_psf_g_2": "gauss_psfReconvolved_g2",
958 "gauss_psf_T": "gauss_psfReconvolved_T",
959 "gauss_s2n": "gauss_snr",
960 "gauss_T": "gauss_T",
961 "gauss_T_err": "gauss_TErr",
962 "gauss_T_flags": "gauss_shape_flags",
963 "gauss_T_ratio": "gauss_T_ratio", # dropped.
964 "id": "shearObjectId",
965 "mfrac": "mfrac",
966 "ormask": "ormask_flags",
967 "pgauss_flags": "pgauss_flags",
968 "pgauss_obj_flags": "pgauss_object_flags",
969 "pgauss_s2n": "pgauss_snr",
970 "pgauss_T": "pgauss_T",
971 "pgauss_T_err": "pgauss_TErr",
972 "pgauss_T_flags": "pgauss_shape_flags",
973 "pgauss_T_ratio": "pgauss_T_ratio", # dropped.
974 "psfrec_flags": "psfOriginal_flags",
975 "psfrec_g_1": "psfOriginal_e1",
976 "psfrec_g_2": "psfOriginal_e2",
977 "psfrec_T": "psfOriginal_T",
978 "ra": "ra",
979 "row": "y",
980 "row_diff": "y_offset", # dropped.
981 "shear_type": "metaStep",
982 "stamp_flags": "image_flags",
983 }
985 for b, alg_name in product(self.config.photometry_bands, ("gauss", "pgauss")):
986 mapping[f"{alg_name}_band_flux_{b}"] = f"{b}_{alg_name}Flux"
987 mapping[f"{alg_name}_band_flux_err_{b}"] = f"{b}_{alg_name}FluxErr"
988 mapping[f"{alg_name}_band_flux_flags_{b}"] = f"{b}_{alg_name}Flux_flags"
990 for name in mapping:
991 if name in data.dtype.names:
992 output[mapping.get(name, name)] = data[name]
993 else:
994 if "flags" in name.lower():
995 output[mapping.get(name, name)] = np.ones_like(data["id"], dtype=np.int32)
996 else:
997 output[mapping.get(name, name)] = np.ones_like(data["id"], dtype=np.float32)
998 output[mapping.get(name, name)] *= np.nan
1000 output["tract"] = tract * np.ones_like(data["id"], dtype=np.int64)
1001 output["patch"] = patch * np.ones_like(data["id"], dtype=np.int32)
1003 return output
1005 @staticmethod
1006 def _fill_is_inner(
1007 di: dict,
1008 bbox: geom.Box2I,
1009 key: str,
1010 ) -> None:
1011 # Convert it to a Box2D object so the full pixels are considered, not
1012 # just the pixel centers.
1013 bbox = geom.Box2D(bbox)
1014 xmin = bbox.getMinX()
1015 ymin = bbox.getMinY()
1016 xmax = bbox.getMaxX()
1017 ymax = bbox.getMaxY()
1018 di[key] = (di["x"] >= xmin) & (di["x"] <= xmax) & (di["y"] >= ymin) & (di["y"] <= ymax)
1021def _make_comb_data(
1022 cell_coadd,
1023 res,
1024 mask_frac,
1025 bands,
1026 cell_id,
1027):
1028 idinfo = cell_coadd.identifiers
1030 copy_dt = [
1031 # we will copy out of arrays to these
1032 ("psfrec_g_1", "f4"),
1033 ("psfrec_g_2", "f4"),
1034 ("gauss_psf_g_1", "f4"),
1035 ("gauss_psf_g_2", "f4"),
1036 ("gauss_g_1", "f4"),
1037 ("gauss_g_2", "f4"),
1038 ("gauss_g_cov_11", "f4"),
1039 ("gauss_g_cov_12", "f4"),
1040 ("gauss_g_cov_21", "f4"),
1041 ("gauss_g_cov_22", "f4"),
1042 ]
1044 for b in bands:
1045 copy_dt.append(("gauss_band_flux_flags_%s" % b, "i4"))
1046 copy_dt.append(("gauss_band_flux_%s" % b, "f4"))
1047 copy_dt.append(("gauss_band_flux_err_%s" % b, "f4"))
1048 copy_dt.append(("pgauss_band_flux_flags_%s" % b, "i4"))
1049 copy_dt.append(("pgauss_band_flux_%s" % b, "f4"))
1050 copy_dt.append(("pgauss_band_flux_err_%s" % b, "f4"))
1052 add_dt = [
1053 ("id", "u8"),
1054 ("tract", "u4"),
1055 ("patch_x", "u1"),
1056 ("patch_y", "u1"),
1057 ("cell_x", "u1"),
1058 ("cell_y", "u1"),
1059 ("shear_type", "U2"),
1060 ("mask_frac", "f4"),
1061 ("primary", bool),
1062 ] + copy_dt
1064 if not hasattr(res, "keys"):
1065 res = {"noshear": res}
1067 dlist = []
1068 for stype in res.keys():
1069 data = res[stype]
1070 if data is not None:
1071 newdata = eu.numpy_util.add_fields(data, add_dt)
1072 newdata["psfrec_g_1"] = newdata["psfrec_g"][:, 0]
1073 newdata["psfrec_g_2"] = newdata["psfrec_g"][:, 1]
1075 newdata["gauss_psf_g_1"] = newdata["gauss_psf_g"][:, 0]
1076 newdata["gauss_psf_g_2"] = newdata["gauss_psf_g"][:, 1]
1077 newdata["gauss_g_1"] = newdata["gauss_g"][:, 0]
1078 newdata["gauss_g_2"] = newdata["gauss_g"][:, 1]
1080 newdata["gauss_g_cov_11"] = newdata["gauss_g_cov"][:, 0, 0]
1081 newdata["gauss_g_cov_12"] = newdata["gauss_g_cov"][:, 0, 1]
1082 newdata["gauss_g_cov_21"] = newdata["gauss_g_cov"][:, 1, 0]
1083 newdata["gauss_g_cov_22"] = newdata["gauss_g_cov"][:, 1, 1]
1085 # To-do make compatible with a single band better than this.
1086 if len(bands) > 1:
1087 for i, b in enumerate(bands):
1088 newdata["gauss_band_flux_flags_%s" % b] = newdata["gauss_band_flux_flags"][:, i]
1089 newdata["gauss_band_flux_%s" % b] = newdata["gauss_band_flux"][:, i]
1090 newdata["gauss_band_flux_err_%s" % b] = newdata["gauss_band_flux_err"][:, i]
1091 newdata["pgauss_band_flux_flags_%s" % b] = newdata["pgauss_band_flux_flags"][:, i]
1092 newdata["pgauss_band_flux_%s" % b] = newdata["pgauss_band_flux"][:, i]
1093 newdata["pgauss_band_flux_err_%s" % b] = newdata["pgauss_band_flux_err"][:, i]
1094 newdata["gauss_band_flux_flags_%s" % b] = newdata["gauss_band_flux_flags"][:, i]
1095 newdata["gauss_band_flux_%s" % b] = newdata["gauss_band_flux"][:, i]
1096 newdata["gauss_band_flux_err_%s" % b] = newdata["gauss_band_flux_err"][:, i]
1097 else:
1098 b = bands[0]
1099 newdata["gauss_band_flux_flags_%s" % b] = newdata["gauss_band_flux_flags"]
1100 newdata["gauss_band_flux_%s" % b] = newdata["gauss_band_flux"]
1101 newdata["gauss_band_flux_err_%s" % b] = newdata["gauss_band_flux_err"]
1102 newdata["pgauss_band_flux_flags_%s" % b] = newdata["pgauss_band_flux_flags"]
1103 newdata["pgauss_band_flux_%s" % b] = newdata["pgauss_band_flux"]
1104 newdata["pgauss_band_flux_err_%s" % b] = newdata["pgauss_band_flux_err"]
1105 newdata["gauss_band_flux_flags_%s" % b] = newdata["gauss_band_flux_flags"]
1106 newdata["gauss_band_flux_%s" % b] = newdata["gauss_band_flux"]
1107 newdata["gauss_band_flux_err_%s" % b] = newdata["gauss_band_flux_err"]
1109 newdata["tract"] = idinfo.tract
1110 newdata["patch_x"] = idinfo.patch.x
1111 newdata["patch_y"] = idinfo.patch.y
1112 newdata["cell_x"] = cell_id.x
1113 newdata["cell_y"] = cell_id.y
1115 if stype == "noshear":
1116 newdata["shear_type"] = "ns"
1117 else:
1118 newdata["shear_type"] = stype
1120 dlist.append(newdata)
1122 if len(dlist) > 0:
1123 output = eu.numpy_util.combine_arrlist(dlist)
1124 else:
1125 output = []
1127 return output
1130def _get_mask_frac(mfrac_mbexp, trim_pixels=0):
1131 """
1132 get the average mask frac for each band and then return the max of those
1133 """
1135 mask_fracs = []
1136 for mfrac_exp in mfrac_mbexp:
1137 mfrac = mfrac_exp.image.array
1138 dim = mfrac.shape[0]
1139 mfrac = mfrac[
1140 trim_pixels : dim - trim_pixels - 1,
1141 trim_pixels : dim - trim_pixels - 1,
1142 ]
1143 mask_fracs.append(mfrac.mean())
1145 return max(mask_fracs)