Coverage for python/lsst/pipe/tasks/finalizeCharacterization.py: 45%
562 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-26 09:39 +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"""Task to run a finalized image characterization, using additional data.
23"""
25__all__ = [
26 'FinalizeCharacterizationConnections',
27 'FinalizeCharacterizationConfig',
28 'FinalizeCharacterizationTask',
29 'FinalizeCharacterizationDetectorConnections',
30 'FinalizeCharacterizationDetectorConfig',
31 'FinalizeCharacterizationDetectorTask',
32 'ConsolidateFinalizeCharacterizationDetectorConnections',
33 'ConsolidateFinalizeCharacterizationDetectorConfig',
34 'ConsolidateFinalizeCharacterizationDetectorTask',
35]
37import astropy.table
38import astropy.units as u
39import numpy as np
40import esutil
41from smatch.matcher import Matcher
43import lsst.geom as geom
44from lsst.geom import LinearTransform
45import lsst.pex.config as pexConfig
46import lsst.pipe.base as pipeBase
47import lsst.daf.base as dafBase
48import lsst.afw.table as afwTable
49from lsst.afw.geom import Quadrupole
50import lsst.meas.algorithms as measAlg
51import lsst.meas.extensions.piff.piffPsfDeterminer # noqa: F401
52from lsst.meas.algorithms import MeasureApCorrTask
53from lsst.meas.base import SingleFrameMeasurementTask, ApplyApCorrTask
54from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
56from .reserveIsolatedStars import ReserveIsolatedStarsTask
57from lsst.obs.base.utils import TableVStack
60class FinalizeCharacterizationConnectionsBase(
61 pipeBase.PipelineTaskConnections,
62 dimensions=('instrument', 'visit',),
63 defaultTemplates={},
64):
65 isolated_star_cats = pipeBase.connectionTypes.Input(
66 doc=('Catalog of isolated stars with average positions, number of associated '
67 'sources, and indexes to the isolated_star_sources catalogs.'),
68 name='isolated_star_presource_associations',
69 storageClass='ArrowAstropy',
70 dimensions=('instrument', 'tract', 'skymap'),
71 deferLoad=True,
72 multiple=True,
73 )
74 isolated_star_sources = pipeBase.connectionTypes.Input(
75 doc=('Catalog of isolated star sources with sourceIds, and indexes to the '
76 'isolated_star_cats catalogs.'),
77 name='isolated_star_presources',
78 storageClass='ArrowAstropy',
79 dimensions=('instrument', 'tract', 'skymap'),
80 deferLoad=True,
81 multiple=True,
82 )
83 fgcm_standard_star = pipeBase.connectionTypes.Input(
84 doc=('Catalog of fgcm for color corrections, and indexes to the '
85 'isolated_star_cats catalogs.'),
86 name='fgcm_standard_star',
87 storageClass='ArrowAstropy',
88 dimensions=('instrument', 'tract', 'skymap'),
89 deferLoad=True,
90 multiple=True,
91 )
93 def __init__(self, *, config=None):
94 super().__init__(config=config)
95 if config is None:
96 return None
97 # Keep fgcm_standard_star if using chromatic PSF OR if adding FGCM photometry
98 needs_fgcm = (config.psf_determiner['piff'].useColor
99 or config.do_add_fgcm_photometry)
100 if not needs_fgcm:
101 self.inputs.remove("fgcm_standard_star")
104class FinalizeCharacterizationConnections(
105 FinalizeCharacterizationConnectionsBase,
106 dimensions=('instrument', 'visit',),
107 defaultTemplates={},
108):
109 srcs = pipeBase.connectionTypes.Input(
110 doc='Source catalogs for the visit',
111 name='src',
112 storageClass='SourceCatalog',
113 dimensions=('instrument', 'visit', 'detector'),
114 deferLoad=True,
115 multiple=True,
116 deferGraphConstraint=True,
117 )
118 calexps = pipeBase.connectionTypes.Input(
119 doc='Calexps for the visit',
120 name='calexp',
121 storageClass='ExposureF',
122 dimensions=('instrument', 'visit', 'detector'),
123 deferLoad=True,
124 multiple=True,
125 )
126 finalized_psf_ap_corr_cat = pipeBase.connectionTypes.Output(
127 doc=('Per-visit finalized psf models and aperture corrections. This '
128 'catalog uses detector id for the id and are sorted for fast '
129 'lookups of a detector.'),
130 name='finalized_psf_ap_corr_catalog',
131 storageClass='ExposureCatalog',
132 dimensions=('instrument', 'visit'),
133 )
134 finalized_src_table = pipeBase.connectionTypes.Output(
135 doc=('Per-visit catalog of measurements for psf/flag/etc.'),
136 name='finalized_src_table',
137 storageClass='ArrowAstropy',
138 dimensions=('instrument', 'visit'),
139 )
142class FinalizeCharacterizationDetectorConnections(
143 FinalizeCharacterizationConnectionsBase,
144 dimensions=('instrument', 'visit', 'detector',),
145 defaultTemplates={},
146):
147 src = pipeBase.connectionTypes.Input(
148 doc='Source catalog for the visit/detector.',
149 name='src',
150 storageClass='SourceCatalog',
151 dimensions=('instrument', 'visit', 'detector'),
152 )
153 calexp = pipeBase.connectionTypes.Input(
154 doc='Calibrated exposure for the visit/detector.',
155 name='calexp',
156 storageClass='ExposureF',
157 dimensions=('instrument', 'visit', 'detector'),
158 )
159 finalized_psf_ap_corr_detector_cat = pipeBase.connectionTypes.Output(
160 doc=('Per-visit/per-detector finalized psf models and aperture corrections. This '
161 'catalog uses detector id for the id.'),
162 name='finalized_psf_ap_corr_detector_catalog',
163 storageClass='ExposureCatalog',
164 dimensions=('instrument', 'visit', 'detector'),
165 )
166 finalized_src_detector_table = pipeBase.connectionTypes.Output(
167 doc=('Per-visit/per-detector catalog of measurements for psf/flag/etc.'),
168 name='finalized_src_detector_table',
169 storageClass='ArrowAstropy',
170 dimensions=('instrument', 'visit', 'detector'),
171 )
174class FinalizeCharacterizationConfigBase(
175 pipeBase.PipelineTaskConfig,
176 pipelineConnections=FinalizeCharacterizationConnectionsBase,
177):
178 """Configuration for FinalizeCharacterizationBaseTask."""
179 source_selector = sourceSelectorRegistry.makeField(
180 doc="How to select sources",
181 default="science"
182 )
183 id_column = pexConfig.Field(
184 doc='Name of column in isolated_star_sources with source id.',
185 dtype=str,
186 default='sourceId',
187 )
188 reserve_selection = pexConfig.ConfigurableField(
189 target=ReserveIsolatedStarsTask,
190 doc='Task to select reserved stars',
191 )
192 make_psf_candidates = pexConfig.ConfigurableField(
193 target=measAlg.MakePsfCandidatesTask,
194 doc='Task to make psf candidates from selected stars.',
195 )
196 psf_determiner = measAlg.psfDeterminerRegistry.makeField(
197 'PSF Determination algorithm',
198 default='piff'
199 )
200 measurement = pexConfig.ConfigurableField(
201 target=SingleFrameMeasurementTask,
202 doc='Measure sources for aperture corrections'
203 )
204 measure_ap_corr = pexConfig.ConfigurableField(
205 target=MeasureApCorrTask,
206 doc="Subtask to measure aperture corrections"
207 )
208 apply_ap_corr = pexConfig.ConfigurableField(
209 target=ApplyApCorrTask,
210 doc="Subtask to apply aperture corrections"
211 )
212 do_add_sky_moments = pexConfig.Field(
213 dtype=bool,
214 default=False,
215 doc="Add second moments in sky (RA/Dec) and alt/az coordinates to the output table.",
216 )
217 do_add_fgcm_photometry = pexConfig.Field(
218 dtype=bool,
219 default=False,
220 doc="Add FGCM photometry for all bands to the output table.",
221 )
222 fgcmPhotometryBands = pexConfig.ListField(
223 dtype=str,
224 default=['g', 'r', 'i', 'z', 'y'],
225 doc="List of bands to include for FGCM photometry (telescope-dependent).",
226 )
228 def setDefaults(self):
229 super().setDefaults()
231 source_selector = self.source_selector['science']
232 source_selector.setDefaults()
234 # We use the source selector only to select out flagged objects
235 # and signal-to-noise. Isolated, unresolved sources are handled
236 # by the isolated star catalog.
238 source_selector.doFlags = True
239 source_selector.doSignalToNoise = True
240 source_selector.doFluxLimit = False
241 source_selector.doUnresolved = False
242 source_selector.doIsolated = False
244 source_selector.signalToNoise.minimum = 50.0
245 source_selector.signalToNoise.maximum = 1000.0
247 source_selector.signalToNoise.fluxField = 'base_GaussianFlux_instFlux'
248 source_selector.signalToNoise.errField = 'base_GaussianFlux_instFluxErr'
250 source_selector.flags.bad = ['base_PixelFlags_flag_edge',
251 'base_PixelFlags_flag_nodata',
252 'base_PixelFlags_flag_interpolatedCenter',
253 'base_PixelFlags_flag_saturatedCenter',
254 'base_PixelFlags_flag_crCenter',
255 'base_PixelFlags_flag_bad',
256 'base_PixelFlags_flag_interpolated',
257 'base_PixelFlags_flag_saturated',
258 'slot_Centroid_flag',
259 'base_GaussianFlux_flag']
261 # Configure aperture correction to select only high s/n sources (that
262 # were used in the psf modeling) to avoid background problems when
263 # computing the aperture correction map.
264 self.measure_ap_corr.sourceSelector = 'science'
266 ap_selector = self.measure_ap_corr.sourceSelector['science']
267 # We do not need to filter flags or unresolved because we have used
268 # the filtered isolated stars as an input
269 ap_selector.doFlags = False
270 ap_selector.doUnresolved = False
272 import lsst.meas.modelfit # noqa: F401
273 import lsst.meas.extensions.photometryKron # noqa: F401
274 import lsst.meas.extensions.convolved # noqa: F401
275 import lsst.meas.extensions.gaap # noqa: F401
276 import lsst.meas.extensions.shapeHSM # noqa: F401
278 # Set up measurement defaults
279 self.measurement.plugins.names = [
280 'base_FPPosition',
281 'base_PsfFlux',
282 'base_GaussianFlux',
283 'modelfit_DoubleShapeletPsfApprox',
284 'modelfit_CModel',
285 'ext_photometryKron_KronFlux',
286 'ext_convolved_ConvolvedFlux',
287 'ext_gaap_GaapFlux',
288 'ext_shapeHSM_HsmShapeRegauss',
289 'ext_shapeHSM_HsmSourceMoments',
290 'ext_shapeHSM_HsmPsfMoments',
291 'ext_shapeHSM_HsmSourceMomentsRound',
292 'ext_shapeHSM_HigherOrderMomentsSource',
293 'ext_shapeHSM_HigherOrderMomentsPSF',
294 ]
295 self.measurement.slots.modelFlux = 'modelfit_CModel'
296 self.measurement.plugins['ext_convolved_ConvolvedFlux'].seeing.append(8.0)
297 self.measurement.plugins['ext_gaap_GaapFlux'].sigmas = [
298 0.5,
299 0.7,
300 1.0,
301 1.5,
302 2.5,
303 3.0
304 ]
305 self.measurement.plugins['ext_gaap_GaapFlux'].doPsfPhotometry = True
306 self.measurement.slots.shape = 'ext_shapeHSM_HsmSourceMoments'
307 self.measurement.slots.psfShape = 'ext_shapeHSM_HsmPsfMoments'
308 self.measurement.plugins['ext_shapeHSM_HsmShapeRegauss'].deblendNChild = ""
310 # TODO: Remove in DM-44658, streak masking to happen only in ip_diffim
311 # Keep track of which footprints contain streaks
312 self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['STREAK']
313 self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['STREAK']
315 # Turn off slot setting for measurement for centroid and shape
316 # (for which we use the input src catalog measurements)
317 self.measurement.slots.centroid = None
318 self.measurement.slots.apFlux = None
319 self.measurement.slots.calibFlux = None
321 names = self.measurement.plugins['ext_convolved_ConvolvedFlux'].getAllResultNames()
322 self.measure_ap_corr.allowFailure += names
323 names = self.measurement.plugins["ext_gaap_GaapFlux"].getAllGaapResultNames()
324 self.measure_ap_corr.allowFailure += names
327class FinalizeCharacterizationConfig(
328 FinalizeCharacterizationConfigBase,
329 pipelineConnections=FinalizeCharacterizationConnections,
330):
331 pass
334class FinalizeCharacterizationDetectorConfig(
335 FinalizeCharacterizationConfigBase,
336 pipelineConnections=FinalizeCharacterizationDetectorConnections,
337):
338 pass
341class FinalizeCharacterizationTaskBase(pipeBase.PipelineTask):
342 """Run final characterization on exposures."""
343 ConfigClass = FinalizeCharacterizationConfigBase
344 _DefaultName = 'finalize_characterization_base'
346 def __init__(self, **kwargs):
347 super().__init__(**kwargs)
349 # We defer setting the schema and subtasks until
350 # the run method to save I/O pressure on one schema file
351 # that otherwise needs to be loaded by thousands of tasks
352 # simultaneously.
353 self.schema = None
355 self.makeSubtask('reserve_selection')
356 self.makeSubtask('source_selector')
357 self.makeSubtask('make_psf_candidates')
358 self.makeSubtask('psf_determiner')
360 # Only log warning and fatal errors from the source_selector
361 self.source_selector.log.setLevel(self.source_selector.log.WARN)
362 self.isPsfDeterminerPiff = False
363 if isinstance(self.psf_determiner, lsst.meas.extensions.piff.piffPsfDeterminer.PiffPsfDeterminerTask):
364 self.isPsfDeterminerPiff = True
366 def _make_output_schema_mapper(self, input_schema):
367 """Make the schema mapper from the input schema to the output schema.
369 Parameters
370 ----------
371 input_schema : `lsst.afw.table.Schema`
372 Input schema.
374 Returns
375 -------
376 mapper : `lsst.afw.table.SchemaMapper`
377 Schema mapper
378 output_schema : `lsst.afw.table.Schema`
379 Output schema (with alias map)
380 """
381 mapper = afwTable.SchemaMapper(input_schema)
382 mapper.addMinimalSchema(afwTable.SourceTable.makeMinimalSchema())
383 mapper.addMapping(input_schema['slot_Centroid_x'].asKey())
384 mapper.addMapping(input_schema['slot_Centroid_y'].asKey())
386 # The aperture fields may be used by the psf determiner.
387 aper_fields = input_schema.extract('base_CircularApertureFlux_*')
388 for field, item in aper_fields.items():
389 mapper.addMapping(item.key)
391 # The following two may be redundant, but then the mapping is a no-op.
392 # Note that the slot_CalibFlux mapping will copy over any
393 # normalized compensated fluxes that are used for calibration.
394 apflux_fields = input_schema.extract('slot_ApFlux_*')
395 for field, item in apflux_fields.items():
396 mapper.addMapping(item.key)
398 calibflux_fields = input_schema.extract('slot_CalibFlux_*')
399 for field, item in calibflux_fields.items():
400 mapper.addMapping(item.key)
402 mapper.addMapping(
403 input_schema[self.config.source_selector.active.signalToNoise.fluxField].asKey(),
404 'calib_psf_selection_flux')
405 mapper.addMapping(
406 input_schema[self.config.source_selector.active.signalToNoise.errField].asKey(),
407 'calib_psf_selection_flux_err')
409 output_schema = mapper.getOutputSchema()
411 output_schema.addField(
412 'calib_psf_candidate',
413 type='Flag',
414 doc=('set if the source was a candidate for PSF determination, '
415 'as determined from FinalizeCharacterizationTask.'),
416 )
417 output_schema.addField(
418 'calib_psf_reserved',
419 type='Flag',
420 doc=('set if source was reserved from PSF determination by '
421 'FinalizeCharacterizationTask.'),
422 )
423 output_schema.addField(
424 'calib_psf_used',
425 type='Flag',
426 doc=('set if source was used in the PSF determination by '
427 'FinalizeCharacterizationTask.'),
428 )
429 output_schema.addField(
430 'visit',
431 type=np.int64,
432 doc='Visit number for the sources.',
433 )
434 output_schema.addField(
435 'detector',
436 type=np.int32,
437 doc='Detector number for the sources.',
438 )
439 output_schema.addField(
440 'psf_color_value',
441 type=np.float32,
442 doc="Color used in PSF fit."
443 )
444 output_schema.addField(
445 'psf_color_type',
446 type=str,
447 size=10,
448 doc="Color used in PSF fit."
449 )
450 output_schema.addField(
451 'psf_max_value',
452 type=np.float32,
453 doc="Maximum value in the star image used to train PSF.",
454 doReplace=True,
455 )
457 if self.config.do_add_sky_moments:
458 # Sky coordinate moments for source shape (arcsec^2)
459 output_schema.addField(
460 'shape_Iuu',
461 type=np.float32,
462 doc="Second moment of source shape along RA axis in sky coordinates (arcsec^2).",
463 )
464 output_schema.addField(
465 'shape_Ivv',
466 type=np.float32,
467 doc="Second moment of source shape along Dec axis in sky coordinates (arcsec^2).",
468 )
469 output_schema.addField(
470 'shape_Iuv',
471 type=np.float32,
472 doc="Cross-term of source shape second moments in sky coordinates (arcsec^2).",
473 )
475 # Sky coordinate moments for PSF shape (arcsec^2)
476 output_schema.addField(
477 'psfShape_Iuu',
478 type=np.float32,
479 doc="Second moment of PSF shape along RA axis in sky coordinates (arcsec^2).",
480 )
481 output_schema.addField(
482 'psfShape_Ivv',
483 type=np.float32,
484 doc="Second moment of PSF shape along Dec axis in sky coordinates (arcsec^2).",
485 )
486 output_schema.addField(
487 'psfShape_Iuv',
488 type=np.float32,
489 doc="Cross-term of PSF shape second moments in sky coordinates (arcsec^2).",
490 )
492 # Alt/Az coordinate moments for source shape (arcsec^2)
493 output_schema.addField(
494 'shape_Ialtalt',
495 type=np.float32,
496 doc="Second moment of source shape along altitude axis (arcsec^2).",
497 )
498 output_schema.addField(
499 'shape_Iazaz',
500 type=np.float32,
501 doc="Second moment of source shape along azimuth axis (arcsec^2).",
502 )
503 output_schema.addField(
504 'shape_Ialtaz',
505 type=np.float32,
506 doc="Cross-term of source shape second moments in alt/az coordinates (arcsec^2).",
507 )
509 # Alt/Az coordinate moments for PSF shape (arcsec^2)
510 output_schema.addField(
511 'psfShape_Ialtalt',
512 type=np.float32,
513 doc="Second moment of PSF shape along altitude axis (arcsec^2).",
514 )
515 output_schema.addField(
516 'psfShape_Iazaz',
517 type=np.float32,
518 doc="Second moment of PSF shape along azimuth axis (arcsec^2).",
519 )
520 output_schema.addField(
521 'psfShape_Ialtaz',
522 type=np.float32,
523 doc="Cross-term of PSF shape second moments in alt/az coordinates (arcsec^2).",
524 )
526 if self.config.do_add_fgcm_photometry:
527 # FGCM photometry for all configured bands
528 for band_name in self.config.fgcmPhotometryBands:
529 output_schema.addField(
530 f'fgcm_mag_{band_name}',
531 type=np.float32,
532 doc=f"FGCM standard star magnitude in {band_name} band (mag).",
533 )
535 alias_map = input_schema.getAliasMap()
536 alias_map_output = afwTable.AliasMap()
537 alias_map_output.set('slot_Centroid', alias_map.get('slot_Centroid'))
538 alias_map_output.set('slot_ApFlux', alias_map.get('slot_ApFlux'))
539 alias_map_output.set('slot_CalibFlux', alias_map.get('slot_CalibFlux'))
541 output_schema.setAliasMap(alias_map_output)
543 return mapper, output_schema
545 def _make_selection_schema_mapper(self, input_schema):
546 """Make the schema mapper from the input schema to the selection schema.
548 Parameters
549 ----------
550 input_schema : `lsst.afw.table.Schema`
551 Input schema.
553 Returns
554 -------
555 mapper : `lsst.afw.table.SchemaMapper`
556 Schema mapper
557 selection_schema : `lsst.afw.table.Schema`
558 Selection schema (with alias map)
559 """
560 mapper = afwTable.SchemaMapper(input_schema)
561 mapper.addMinimalSchema(input_schema)
563 selection_schema = mapper.getOutputSchema()
565 selection_schema.setAliasMap(input_schema.getAliasMap())
567 selection_schema.addField(
568 'psf_color_value',
569 type=np.float32,
570 doc="Color used in PSF fit."
571 )
572 selection_schema.addField(
573 'psf_color_type',
574 type=str,
575 size=10,
576 doc="Color used in PSF fit."
577 )
578 selection_schema.addField(
579 'psf_max_value',
580 type=np.float32,
581 doc="Maximum value in the star image used to train PSF.",
582 doReplace=True,
583 )
585 return mapper, selection_schema
587 def _compute_sky_moments(self, wcs, x, y, ixx, iyy, ixy):
588 """Compute second moments in sky coordinates from pixel moments.
590 Transforms the second moments tensor from pixel coordinates to sky
591 coordinates (RA/Dec) using the local WCS CD matrix at each source
592 position.
594 Parameters
595 ----------
596 wcs : `lsst.afw.geom.SkyWcs`
597 The WCS of the exposure.
598 x : `numpy.ndarray`
599 X pixel coordinates of the sources.
600 y : `numpy.ndarray`
601 Y pixel coordinates of the sources.
602 ixx : `numpy.ndarray`
603 Second moment Ixx in pixel coordinates (pixels^2).
604 iyy : `numpy.ndarray`
605 Second moment Iyy in pixel coordinates (pixels^2).
606 ixy : `numpy.ndarray`
607 Second moment Ixy in pixel coordinates (pixels^2).
609 Returns
610 -------
611 iuu : `numpy.ndarray`
612 Second moment along RA axis in sky coordinates (arcsec^2).
613 ivv : `numpy.ndarray`
614 Second moment along Dec axis in sky coordinates (arcsec^2).
615 iuv : `numpy.ndarray`
616 Cross-term of second moments in sky coordinates (arcsec^2).
617 """
618 n_sources = len(x)
619 iuu = np.full(n_sources, np.nan, dtype=np.float32)
620 ivv = np.full(n_sources, np.nan, dtype=np.float32)
621 iuv = np.full(n_sources, np.nan, dtype=np.float32)
623 # Conversion factor from degrees^2 to arcsec^2
624 # (CD matrix is in degrees/pixel per FITS standard)
625 deg2_to_arcsec2 = 3600.0 ** 2
627 for i in range(n_sources):
628 # Skip sources with invalid moments
629 if not np.isfinite(ixx[i]) or not np.isfinite(iyy[i]) or not np.isfinite(ixy[i]):
630 continue
632 center = geom.Point2D(x[i], y[i])
633 cd_matrix = wcs.getCdMatrix(center)
635 cd11 = cd_matrix[0, 0]
636 cd12 = cd_matrix[0, 1]
637 cd21 = cd_matrix[1, 0]
638 cd22 = cd_matrix[1, 1]
640 # Transform moments: M_sky = CD * M_pixel * CD^T
641 # Iuu = CD11*(Ixx*CD11 + Ixy*CD12) + CD12*(Ixy*CD11 + Iyy*CD12)
642 iuu[i] = (cd11 * (ixx[i] * cd11 + ixy[i] * cd12)
643 + cd12 * (ixy[i] * cd11 + iyy[i] * cd12))
644 # Ivv = CD21*(Ixx*CD21 + Ixy*CD22) + CD22*(Ixy*CD21 + Iyy*CD22)
645 ivv[i] = (cd21 * (ixx[i] * cd21 + ixy[i] * cd22)
646 + cd22 * (ixy[i] * cd21 + iyy[i] * cd22))
647 # Iuv = (CD11*Ixx + CD12*Ixy)*CD21 + (CD11*Ixy + CD12*Iyy)*CD22
648 iuv[i] = ((cd11 * ixx[i] + cd12 * ixy[i]) * cd21
649 + (cd11 * ixy[i] + cd12 * iyy[i]) * cd22)
651 # Convert from degrees^2 to arcsec^2
652 iuu[i] *= deg2_to_arcsec2
653 ivv[i] *= deg2_to_arcsec2
654 iuv[i] *= deg2_to_arcsec2
656 return iuu, ivv, iuv
658 def _rotate_moments_to_altaz(self, iuu, ivv, iuv, parallactic_angle):
659 """Rotate second moments from RA/Dec to Alt/Az coordinates.
661 Rotates the second moments tensor from equatorial (RA/Dec) coordinates
662 to horizontal (Alt/Az) coordinates using the parallactic angle.
664 Parameters
665 ----------
666 iuu : `numpy.ndarray`
667 Second moment along RA axis in sky coordinates (arcsec^2).
668 ivv : `numpy.ndarray`
669 Second moment along Dec axis in sky coordinates (arcsec^2).
670 iuv : `numpy.ndarray`
671 Cross-term of second moments in sky coordinates (arcsec^2).
672 parallactic_angle : `lsst.geom.Angle`
673 The parallactic angle (from visitInfo.boresightParAngle).
675 Returns
676 -------
677 ialtalt : `numpy.ndarray`
678 Second moment along altitude axis (arcsec^2).
679 iazaz : `numpy.ndarray`
680 Second moment along azimuth axis (arcsec^2).
681 ialtaz : `numpy.ndarray`
682 Cross-term of second moments in alt/az coordinates (arcsec^2).
683 """
685 n_sources = len(iuu)
686 ialtalt = np.full(n_sources, np.nan, dtype=np.float32)
687 iazaz = np.full(n_sources, np.nan, dtype=np.float32)
688 ialtaz = np.full(n_sources, np.nan, dtype=np.float32)
690 # Create the rotation transformation for the parallactic angle
691 rot_transform = LinearTransform.makeRotation(parallactic_angle)
693 for i in range(n_sources):
694 # Skip sources with invalid moments
695 if not np.isfinite(iuu[i]) or not np.isfinite(ivv[i]) or not np.isfinite(iuv[i]):
696 continue
698 # Create a Quadrupole from the sky moments and rotate it
699 shape = Quadrupole(iuu[i], ivv[i], iuv[i])
700 rot_shape = shape.transform(rot_transform)
702 ialtalt[i] = rot_shape.getIxx()
703 iazaz[i] = rot_shape.getIyy()
704 ialtaz[i] = rot_shape.getIxy()
706 return ialtalt, iazaz, ialtaz
708 def concat_isolated_star_cats(
709 self,
710 band,
711 isolated_star_cat_dict,
712 isolated_star_source_dict,
713 visit=None,
714 detector=None,
715 ):
716 """
717 Concatenate isolated star catalogs and make reserve selection.
719 Parameters
720 ----------
721 band : `str`
722 Band name. Used to select reserved stars.
723 isolated_star_cat_dict : `dict`
724 Per-tract dict of isolated star catalog handles.
725 isolated_star_source_dict : `dict`
726 Per-tract dict of isolated star source catalog handles.
727 visit : `int`, optional
728 Visit to down-select sources.
729 detector : `int`, optional
730 Detector to down-select sources. Will only be used if visit
731 is also set.
733 Returns
734 -------
735 isolated_table : `astropy.table.Table` (N,)
736 Table of isolated stars, with indexes to isolated sources.
737 Returns None if there are no usable isolated catalogs.
738 isolated_source_table : `astropy.table.Table` (M,)
739 Table of isolated sources, with indexes to isolated stars.
740 Returns None if there are no usable isolated catalogs.
741 """
742 isolated_tables = []
743 isolated_sources = []
744 merge_cat_counter = 0
745 merge_source_counter = 0
747 source_columns = [self.config.id_column, 'obj_index']
748 if visit is not None:
749 source_columns.append('visit')
750 if detector is not None:
751 source_columns.append('detector')
753 for tract in isolated_star_cat_dict:
754 astropy_cat = isolated_star_cat_dict[tract].get()
755 astropy_source = isolated_star_source_dict[tract].get(
756 parameters={'columns': source_columns}
757 )
759 # Cut the isolated star sources to this visit/detector.
760 sources_downselected = False
761 if visit is not None:
762 sources_downselected = True
763 these_sources = (astropy_source['visit'] == visit)
765 if these_sources.sum() == 0: 765 ↛ 766line 765 didn't jump to line 766 because the condition on line 765 was never true
766 self.log.info('No sources found for visit %d in tract %d.', visit, tract)
767 continue
769 if detector is not None:
770 these_sources &= (astropy_source['detector'] == detector)
771 if these_sources.sum() == 0:
772 self.log.info(
773 'No sources found for visit %d, detector %d in tract %d.',
774 visit,
775 detector,
776 tract,
777 )
778 continue
780 astropy_source = astropy_source[these_sources]
782 # Cut isolated star table to those observed in this band, and adjust indexes
783 # We must use all the stars in the table in the band to ensure consistent
784 # reserve star selection below.
785 (use_band,) = (astropy_cat[f'nsource_{band}'] > 0).nonzero()
787 if len(use_band) == 0:
788 # There are no sources in this band in this tract.
789 self.log.info("No sources found in %s band in tract %d.", band, tract)
790 continue
792 # With the following matching:
793 # table_source[b] <-> table_cat[use_band[a]]
794 a, b = esutil.numpy_util.match(use_band, np.asarray(astropy_source['obj_index']))
796 # Update indexes and cut to band-selected stars/sources
797 astropy_source['obj_index'][b] = a
798 if sources_downselected:
799 # The isolated star catalog is built such that each star
800 # in the catalog is matched to multiple sources. But due
801 # to the way it is constructed of isolated stars there
802 # is at most one unique source on any visit associated
803 # with a single isolated star in the catalog. Therefore,
804 # everything is guaranteed to already be uniquely
805 # matched here because we did visit (and detector)
806 # down-selection above.
807 astropy_cat[f'source_cat_index_{band}'][use_band][a] = b
808 else:
809 # If there was no down-selection then each star has
810 # multiple sources.
811 _, index_new = np.unique(a, return_index=True)
812 astropy_cat[f'source_cat_index_{band}'][use_band] = index_new
814 # After the following cuts, the catalogs have the following properties:
815 # - table_cat only contains isolated stars that have at least one source
816 # in ``band``.
817 # - table_source only contains ``band`` sources.
818 # - The slice table_cat["source_cat_index_{band}"]: table_cat["source_cat_index_{band}"]
819 # + table_cat["nsource_{band}]
820 # applied to table_source will give all the sources associated with the star.
821 # - For each source, table_source["obj_index"] points to the index of the associated
822 # isolated star.
823 astropy_source = astropy_source[b]
824 astropy_cat = astropy_cat[use_band]
826 # Add reserved flag column to tables
827 astropy_cat['reserved'] = False
828 astropy_source['reserved'] = False
830 # Get reserve star flags
831 astropy_cat['reserved'][:] = self.reserve_selection.run(
832 len(astropy_cat),
833 extra=f'{band}_{tract}',
834 )
835 astropy_source['reserved'][:] = astropy_cat['reserved'][astropy_source['obj_index']]
837 # Offset indexes to account for tract merging
838 astropy_cat[f'source_cat_index_{band}'] += merge_source_counter
839 astropy_source['obj_index'] += merge_cat_counter
841 isolated_tables.append(astropy_cat)
842 isolated_sources.append(astropy_source)
844 merge_cat_counter += len(astropy_cat)
845 merge_source_counter += len(astropy_source)
847 if len(isolated_tables) > 0: 847 ↛ 851line 847 didn't jump to line 851 because the condition on line 847 was always true
848 isolated_table = astropy.table.vstack(isolated_tables, metadata_conflicts='silent')
849 isolated_source_table = astropy.table.vstack(isolated_sources, metadata_conflicts='silent')
850 else:
851 isolated_table = None
852 isolated_source_table = None
854 return isolated_table, isolated_source_table
856 def add_src_colors(self, srcCat, fgcmCat, band):
858 needColor = self.isPsfDeterminerPiff and fgcmCat is not None
859 needColor &= self.config.psf_determiner['piff'].useColor
861 if needColor:
863 raSrc = (srcCat['coord_ra'] * u.radian).to(u.degree).value
864 decSrc = (srcCat['coord_dec'] * u.radian).to(u.degree).value
866 with Matcher(raSrc, decSrc) as matcher:
867 idx, idxSrcCat, idxColorCat, d = matcher.query_radius(
868 fgcmCat["ra"],
869 fgcmCat["dec"],
870 1. / 3600.0,
871 return_indices=True,
872 )
874 magStr1 = self.psf_determiner.config.color[band][0]
875 magStr2 = self.psf_determiner.config.color[band][2]
876 colors = fgcmCat[f'mag_{magStr1}'] - fgcmCat[f'mag_{magStr2}']
878 for idSrc, idColor in zip(idxSrcCat, idxColorCat):
879 srcCat[idSrc]['psf_color_value'] = colors[idColor]
880 srcCat[idSrc]['psf_color_type'] = f"{magStr1}-{magStr2}"
882 def add_fgcm_photometry(self, srcCat, fgcmCat):
883 """Add FGCM photometry for all configured bands to the source catalog.
885 Parameters
886 ----------
887 srcCat : `lsst.afw.table.SourceCatalog`
888 Source catalog to add photometry to.
889 fgcmCat : `numpy.ndarray`
890 FGCM standard star catalog with ra, dec, and mag_* columns.
891 """
892 # Initialize all FGCM magnitude columns to NaN
893 for band_name in self.config.fgcmPhotometryBands:
894 output_col = f'fgcm_mag_{band_name}'
895 srcCat[output_col] = np.nan
897 if fgcmCat is None:
898 return
900 raSrc = (srcCat['coord_ra'] * u.radian).to(u.degree).value
901 decSrc = (srcCat['coord_dec'] * u.radian).to(u.degree).value
903 with Matcher(raSrc, decSrc) as matcher:
904 idx, idxSrcCat, idxFgcmCat, d = matcher.query_radius(
905 fgcmCat["ra"],
906 fgcmCat["dec"],
907 1. / 3600.0,
908 return_indices=True,
909 )
911 # Add FGCM magnitudes for all configured bands
912 for band_name in self.config.fgcmPhotometryBands:
913 mag_col = f'mag_{band_name}'
914 output_col = f'fgcm_mag_{band_name}'
915 if mag_col in fgcmCat.dtype.names:
916 for idSrc, idFgcm in zip(idxSrcCat, idxFgcmCat):
917 srcCat[idSrc][output_col] = fgcmCat[mag_col][idFgcm]
919 def compute_psf_and_ap_corr_map(self, visit, detector, exposure, src,
920 isolated_source_table, fgcm_standard_star_cat):
921 """Compute psf model and aperture correction map for a single exposure.
923 Parameters
924 ----------
925 visit : `int`
926 Visit number (for logging).
927 detector : `int`
928 Detector number (for logging).
929 exposure : `lsst.afw.image.ExposureF`
930 src : `lsst.afw.table.SourceCatalog`
931 isolated_source_table : `np.ndarray` or `astropy.table.Table`
932 fgcm_standard_star_cat : `np.ndarray`
934 Returns
935 -------
936 psf : `lsst.meas.algorithms.ImagePsf`
937 PSF Model
938 ap_corr_map : `lsst.afw.image.ApCorrMap`
939 Aperture correction map.
940 measured_src : `lsst.afw.table.SourceCatalog`
941 Updated source catalog with measurements, flags and aperture corrections.
942 """
943 # Extract footprints from the input src catalog for noise replacement.
944 footprints = SingleFrameMeasurementTask.getFootprintsFromCatalog(src)
946 # Apply source selector (s/n, flags, etc.)
947 good_src = self.source_selector.selectSources(src)
948 if sum(good_src.selected) == 0: 948 ↛ 958line 948 didn't jump to line 958 because the condition on line 948 was always true
949 self.log.warning('No good sources remain after cuts for visit %d, detector %d',
950 visit, detector)
951 return None, None, None
953 # Cut down input src to the selected sources
954 # We use a separate schema/mapper here than for the output/measurement catalog because of
955 # clashes between fields that were previously run and those that need to be rerun with
956 # the new psf model. This may be slightly inefficient but keeps input
957 # and output values cleanly separated.
958 selection_mapper, selection_schema = self._make_selection_schema_mapper(src.schema)
960 selected_src = afwTable.SourceCatalog(selection_schema)
961 selected_src.reserve(good_src.selected.sum())
962 selected_src.extend(src[good_src.selected], mapper=selection_mapper)
964 # The calib flags have been copied from the input table,
965 # and we reset them here just to ensure they aren't propagated.
966 selected_src['calib_psf_candidate'] = np.zeros(len(selected_src), dtype=bool)
967 selected_src['calib_psf_used'] = np.zeros(len(selected_src), dtype=bool)
968 selected_src['calib_psf_reserved'] = np.zeros(len(selected_src), dtype=bool)
970 # Find the isolated sources and set flags
971 matched_src, matched_iso = esutil.numpy_util.match(
972 selected_src['id'],
973 np.asarray(isolated_source_table[self.config.id_column]),
974 )
975 if len(matched_src) == 0:
976 self.log.warning(
977 "No candidates from matched isolate stars for visit=%s, detector=%s "
978 "(this is probably the result of an earlier astrometry failure).",
979 visit, detector,
980 )
981 return None, None, None
983 matched_arr = np.zeros(len(selected_src), dtype=bool)
984 matched_arr[matched_src] = True
985 selected_src['calib_psf_candidate'] = matched_arr
987 reserved_arr = np.zeros(len(selected_src), dtype=bool)
988 reserved_arr[matched_src] = np.asarray(isolated_source_table['reserved'][matched_iso])
989 selected_src['calib_psf_reserved'] = reserved_arr
991 selected_src = selected_src[selected_src['calib_psf_candidate']].copy(deep=True)
993 # Make the measured source catalog as well, based on the selected catalog.
994 measured_src = afwTable.SourceCatalog(self.schema)
995 measured_src.reserve(len(selected_src))
996 measured_src.extend(selected_src, mapper=self.schema_mapper)
998 # We need to copy over the calib_psf flags because they were not in the mapper
999 measured_src['calib_psf_candidate'] = selected_src['calib_psf_candidate']
1000 measured_src['calib_psf_reserved'] = selected_src['calib_psf_reserved']
1001 if exposure.filter.hasBandLabel():
1002 band = exposure.filter.bandLabel
1003 else:
1004 band = None
1005 self.add_src_colors(selected_src, fgcm_standard_star_cat, band)
1006 self.add_src_colors(measured_src, fgcm_standard_star_cat, band)
1008 # Add FGCM photometry for all bands to the output catalog
1009 if self.config.do_add_fgcm_photometry:
1010 self.add_fgcm_photometry(measured_src, fgcm_standard_star_cat)
1012 # Select the psf candidates from the selection catalog
1013 try:
1014 psf_selection_result = self.make_psf_candidates.run(selected_src, exposure=exposure)
1015 _ = self.make_psf_candidates.run(measured_src, exposure=exposure)
1016 except Exception as e:
1017 self.log.exception('Failed to make PSF candidates for visit %d, detector %d: %s',
1018 visit, detector, e)
1019 return None, None, measured_src
1021 psf_cand_cat = psf_selection_result.goodStarCat
1023 # Make list of psf candidates to send to the determiner
1024 # (omitting those marked as reserved)
1025 psf_determiner_list = [cand for cand, use
1026 in zip(psf_selection_result.psfCandidates,
1027 ~psf_cand_cat['calib_psf_reserved']) if use]
1028 flag_key = psf_cand_cat.schema['calib_psf_used'].asKey()
1030 try:
1031 psf, cell_set = self.psf_determiner.determinePsf(exposure,
1032 psf_determiner_list,
1033 self.metadata,
1034 flagKey=flag_key)
1035 except Exception as e:
1036 self.log.exception('Failed to determine PSF for visit %d, detector %d: %s',
1037 visit, detector, e)
1038 return None, None, measured_src
1039 # Verify that the PSF is usable by downstream tasks
1040 sigma = psf.computeShape(psf.getAveragePosition(), psf.getAverageColor()).getDeterminantRadius()
1041 if np.isnan(sigma):
1042 self.log.warning('Failed to determine psf for visit %d, detector %d: '
1043 'Computed final PSF size is NAN.',
1044 visit, detector)
1045 return None, None, measured_src
1047 # Set the psf in the exposure for measurement/aperture corrections.
1048 exposure.setPsf(psf)
1050 # At this point, we need to transfer the psf used flag from the selection
1051 # catalog to the measurement catalog.
1052 matched_selected, matched_measured = esutil.numpy_util.match(
1053 selected_src['id'],
1054 measured_src['id']
1055 )
1056 measured_used = np.zeros(len(measured_src), dtype=bool)
1057 measured_used[matched_measured] = selected_src['calib_psf_used'][matched_selected]
1058 measured_src['calib_psf_used'] = measured_used
1060 # Next, we do the measurement on all the psf candidate, used, and reserved stars.
1061 # We use the full footprint list from the input src catalog for noise replacement.
1062 try:
1063 self.measurement.run(measCat=measured_src, exposure=exposure, footprints=footprints)
1064 except Exception as e:
1065 self.log.warning('Failed to make measurements for visit %d, detector %d: %s',
1066 visit, detector, e)
1067 return psf, None, measured_src
1069 # Compute sky and alt/az coordinate moments for source and PSF shapes.
1070 if self.config.do_add_sky_moments:
1071 wcs = exposure.getWcs()
1072 x = np.array(measured_src['slot_Centroid_x'])
1073 y = np.array(measured_src['slot_Centroid_y'])
1075 # Source shape sky moments
1076 shape_xx = np.array(measured_src['slot_Shape_xx'])
1077 shape_yy = np.array(measured_src['slot_Shape_yy'])
1078 shape_xy = np.array(measured_src['slot_Shape_xy'])
1079 shape_iuu, shape_ivv, shape_iuv = self._compute_sky_moments(
1080 wcs, x, y, shape_xx, shape_yy, shape_xy
1081 )
1082 measured_src['shape_Iuu'] = shape_iuu
1083 measured_src['shape_Ivv'] = shape_ivv
1084 measured_src['shape_Iuv'] = shape_iuv
1086 # PSF shape sky moments
1087 psf_xx = np.array(measured_src['slot_PsfShape_xx'])
1088 psf_yy = np.array(measured_src['slot_PsfShape_yy'])
1089 psf_xy = np.array(measured_src['slot_PsfShape_xy'])
1090 psf_iuu, psf_ivv, psf_iuv = self._compute_sky_moments(
1091 wcs, x, y, psf_xx, psf_yy, psf_xy
1092 )
1093 measured_src['psfShape_Iuu'] = psf_iuu
1094 measured_src['psfShape_Ivv'] = psf_ivv
1095 measured_src['psfShape_Iuv'] = psf_iuv
1097 # Rotate sky moments to alt/az coordinates using the parallactic angle.
1098 visit_info = exposure.info.getVisitInfo()
1099 parallactic_angle = visit_info.boresightParAngle
1101 # Source shape alt/az moments
1102 shape_ialtalt, shape_iazaz, shape_ialtaz = self._rotate_moments_to_altaz(
1103 shape_iuu, shape_ivv, shape_iuv, parallactic_angle
1104 )
1105 measured_src['shape_Ialtalt'] = shape_ialtalt
1106 measured_src['shape_Iazaz'] = shape_iazaz
1107 measured_src['shape_Ialtaz'] = shape_ialtaz
1109 # PSF shape alt/az moments
1110 psf_ialtalt, psf_iazaz, psf_ialtaz = self._rotate_moments_to_altaz(
1111 psf_iuu, psf_ivv, psf_iuv, parallactic_angle
1112 )
1113 measured_src['psfShape_Ialtalt'] = psf_ialtalt
1114 measured_src['psfShape_Iazaz'] = psf_iazaz
1115 measured_src['psfShape_Ialtaz'] = psf_ialtaz
1117 # And finally the ap corr map.
1118 try:
1119 ap_corr_map = self.measure_ap_corr.run(exposure=exposure,
1120 catalog=measured_src).apCorrMap
1121 except Exception as e:
1122 self.log.warning('Failed to compute aperture corrections for visit %d, detector %d: %s',
1123 visit, detector, e)
1124 return psf, None, measured_src
1126 # Need to merge the original normalization aperture correction map.
1127 ap_corr_map_input = exposure.apCorrMap
1128 for key in ap_corr_map_input:
1129 if key not in ap_corr_map:
1130 ap_corr_map[key] = ap_corr_map_input[key]
1132 self.apply_ap_corr.run(catalog=measured_src, apCorrMap=ap_corr_map)
1134 return psf, ap_corr_map, measured_src
1137class FinalizeCharacterizationTask(FinalizeCharacterizationTaskBase):
1138 """Run final characterization on full visits."""
1139 ConfigClass = FinalizeCharacterizationConfig
1140 _DefaultName = 'finalize_characterization'
1142 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1143 input_handle_dict = butlerQC.get(inputRefs)
1145 band = butlerQC.quantum.dataId['band']
1146 visit = butlerQC.quantum.dataId['visit']
1148 src_dict_temp = {handle.dataId['detector']: handle
1149 for handle in input_handle_dict['srcs']}
1150 calexp_dict_temp = {handle.dataId['detector']: handle
1151 for handle in input_handle_dict['calexps']}
1152 isolated_star_cat_dict_temp = {handle.dataId['tract']: handle
1153 for handle in input_handle_dict['isolated_star_cats']}
1154 isolated_star_source_dict_temp = {handle.dataId['tract']: handle
1155 for handle in input_handle_dict['isolated_star_sources']}
1157 src_dict = {detector: src_dict_temp[detector] for
1158 detector in sorted(src_dict_temp.keys())}
1159 calexp_dict = {detector: calexp_dict_temp[detector] for
1160 detector in sorted(calexp_dict_temp.keys())}
1161 isolated_star_cat_dict = {tract: isolated_star_cat_dict_temp[tract] for
1162 tract in sorted(isolated_star_cat_dict_temp.keys())}
1163 isolated_star_source_dict = {tract: isolated_star_source_dict_temp[tract] for
1164 tract in sorted(isolated_star_source_dict_temp.keys())}
1166 needs_fgcm = (self.config.psf_determiner['piff'].useColor
1167 or self.config.do_add_fgcm_photometry)
1168 if needs_fgcm:
1169 fgcm_standard_star_dict_temp = {handle.dataId['tract']: handle
1170 for handle in input_handle_dict['fgcm_standard_star']}
1171 fgcm_standard_star_dict = {tract: fgcm_standard_star_dict_temp[tract] for
1172 tract in sorted(fgcm_standard_star_dict_temp.keys())}
1173 else:
1174 fgcm_standard_star_dict = None
1176 struct = self.run(
1177 visit,
1178 band,
1179 isolated_star_cat_dict,
1180 isolated_star_source_dict,
1181 src_dict,
1182 calexp_dict,
1183 fgcm_standard_star_dict=fgcm_standard_star_dict,
1184 )
1186 butlerQC.put(struct.psf_ap_corr_cat, outputRefs.finalized_psf_ap_corr_cat)
1187 butlerQC.put(struct.output_table, outputRefs.finalized_src_table)
1189 def run(self, visit, band, isolated_star_cat_dict, isolated_star_source_dict,
1190 src_dict, calexp_dict, fgcm_standard_star_dict=None):
1191 """
1192 Run the FinalizeCharacterizationTask.
1194 Parameters
1195 ----------
1196 visit : `int`
1197 Visit number. Used in the output catalogs.
1198 band : `str`
1199 Band name. Used to select reserved stars.
1200 isolated_star_cat_dict : `dict`
1201 Per-tract dict of isolated star catalog handles.
1202 isolated_star_source_dict : `dict`
1203 Per-tract dict of isolated star source catalog handles.
1204 src_dict : `dict`
1205 Per-detector dict of src catalog handles.
1206 calexp_dict : `dict`
1207 Per-detector dict of calibrated exposure handles.
1208 fgcm_standard_star_dict : `dict`
1209 Per tract dict of fgcm isolated stars.
1211 Returns
1212 -------
1213 struct : `lsst.pipe.base.struct`
1214 Struct with outputs for persistence.
1216 Raises
1217 ------
1218 NoWorkFound
1219 Raised if the selector returns no good sources.
1220 """
1221 # Check if we have the same inputs for each of the
1222 # src_dict and calexp_dict.
1223 src_detectors = set(src_dict.keys())
1224 calexp_detectors = set(calexp_dict.keys())
1226 if src_detectors != calexp_detectors: 1226 ↛ 1227line 1226 didn't jump to line 1227 because the condition on line 1226 was never true
1227 detector_keys = sorted(src_detectors.intersection(calexp_detectors))
1228 self.log.warning(
1229 "Input src and calexp have mismatched detectors; "
1230 "running intersection of %d detectors.",
1231 len(detector_keys),
1232 )
1233 else:
1234 detector_keys = sorted(src_detectors)
1236 # We do not need the isolated star table in this task.
1237 # However, it is used in tests to confirm consistency of indexes.
1238 _, isolated_source_table = self.concat_isolated_star_cats(
1239 band,
1240 isolated_star_cat_dict,
1241 isolated_star_source_dict,
1242 visit=visit,
1243 )
1245 if isolated_source_table is None: 1245 ↛ 1246line 1245 didn't jump to line 1246 because the condition on line 1245 was never true
1246 raise pipeBase.NoWorkFound(f'No good isolated sources found for any detectors in visit {visit}')
1248 exposure_cat_schema = afwTable.ExposureTable.makeMinimalSchema()
1249 exposure_cat_schema.addField('visit', type='L', doc='Visit number')
1251 metadata = dafBase.PropertyList()
1252 metadata.add("COMMENT", "Catalog id is detector id, sorted.")
1253 metadata.add("COMMENT", "Only detectors with data have entries.")
1255 psf_ap_corr_cat = afwTable.ExposureCatalog(exposure_cat_schema)
1256 psf_ap_corr_cat.setMetadata(metadata)
1258 measured_src_tables = []
1259 measured_src_table = None
1261 if fgcm_standard_star_dict is not None: 1261 ↛ 1272line 1261 didn't jump to line 1272 because the condition on line 1261 was always true
1263 fgcm_standard_star_cat = []
1265 for tract in fgcm_standard_star_dict:
1266 astropy_fgcm = fgcm_standard_star_dict[tract].get()
1267 table_fgcm = np.asarray(astropy_fgcm)
1268 fgcm_standard_star_cat.append(table_fgcm)
1270 fgcm_standard_star_cat = np.concatenate(fgcm_standard_star_cat)
1271 else:
1272 fgcm_standard_star_cat = None
1274 self.log.info("Running finalizeCharacterization on %d detectors.", len(detector_keys))
1275 for detector in detector_keys:
1276 self.log.info("Starting finalizeCharacterization on detector ID %d.", detector)
1277 src = src_dict[detector].get()
1278 exposure = calexp_dict[detector].get()
1280 # Initialize schema, mappers, and subtasks if necessary.
1281 if self.schema is None: 1281 ↛ 1282line 1281 didn't jump to line 1282 because the condition on line 1281 was never true
1282 self.schema_mapper, self.schema = self._make_output_schema_mapper(src.schema)
1284 self.makeSubtask("measurement", schema=self.schema)
1285 self.makeSubtask("measure_ap_corr", schema=self.schema)
1286 self.makeSubtask("apply_ap_corr", schema=self.schema)
1288 psf, ap_corr_map, measured_src = self.compute_psf_and_ap_corr_map(
1289 visit,
1290 detector,
1291 exposure,
1292 src,
1293 isolated_source_table,
1294 fgcm_standard_star_cat,
1295 )
1297 # And now we package it together...
1298 if measured_src is not None: 1298 ↛ 1275line 1298 didn't jump to line 1275 because the condition on line 1298 was always true
1299 record = psf_ap_corr_cat.addNew()
1300 record['id'] = int(detector)
1301 record['visit'] = visit
1302 if psf is not None: 1302 ↛ 1304line 1302 didn't jump to line 1304 because the condition on line 1302 was always true
1303 record.setPsf(psf)
1304 if ap_corr_map is not None: 1304 ↛ 1307line 1304 didn't jump to line 1307 because the condition on line 1304 was always true
1305 record.setApCorrMap(ap_corr_map)
1307 measured_src['visit'][:] = visit
1308 measured_src['detector'][:] = detector
1310 measured_src_tables.append(measured_src.asAstropy())
1312 if len(measured_src_tables) > 0: 1312 ↛ 1315line 1312 didn't jump to line 1315 because the condition on line 1312 was always true
1313 measured_src_table = astropy.table.vstack(measured_src_tables, join_type='exact')
1315 if measured_src_table is None: 1315 ↛ 1316line 1315 didn't jump to line 1316 because the condition on line 1315 was never true
1316 raise pipeBase.NoWorkFound(f'No good sources found for any detectors in visit {visit}')
1318 return pipeBase.Struct(
1319 psf_ap_corr_cat=psf_ap_corr_cat,
1320 output_table=measured_src_table,
1321 )
1324class FinalizeCharacterizationDetectorTask(FinalizeCharacterizationTaskBase):
1325 """Run final characterization per detector."""
1326 ConfigClass = FinalizeCharacterizationDetectorConfig
1327 _DefaultName = 'finalize_characterization_detector'
1329 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1330 input_handle_dict = butlerQC.get(inputRefs)
1332 band = butlerQC.quantum.dataId['band']
1333 visit = butlerQC.quantum.dataId['visit']
1334 detector = butlerQC.quantum.dataId['detector']
1336 isolated_star_cat_dict_temp = {handle.dataId['tract']: handle
1337 for handle in input_handle_dict['isolated_star_cats']}
1338 isolated_star_source_dict_temp = {handle.dataId['tract']: handle
1339 for handle in input_handle_dict['isolated_star_sources']}
1341 isolated_star_cat_dict = {tract: isolated_star_cat_dict_temp[tract] for
1342 tract in sorted(isolated_star_cat_dict_temp.keys())}
1343 isolated_star_source_dict = {tract: isolated_star_source_dict_temp[tract] for
1344 tract in sorted(isolated_star_source_dict_temp.keys())}
1346 needs_fgcm = (self.config.psf_determiner['piff'].useColor
1347 or self.config.do_add_fgcm_photometry)
1348 if needs_fgcm:
1349 fgcm_standard_star_dict_temp = {handle.dataId['tract']: handle
1350 for handle in input_handle_dict['fgcm_standard_star']}
1351 fgcm_standard_star_dict = {tract: fgcm_standard_star_dict_temp[tract] for
1352 tract in sorted(fgcm_standard_star_dict_temp.keys())}
1353 else:
1354 fgcm_standard_star_dict = None
1356 struct = self.run(
1357 visit,
1358 band,
1359 detector,
1360 isolated_star_cat_dict,
1361 isolated_star_source_dict,
1362 input_handle_dict['src'],
1363 input_handle_dict['calexp'],
1364 fgcm_standard_star_dict=fgcm_standard_star_dict,
1365 )
1367 butlerQC.put(
1368 struct.psf_ap_corr_cat,
1369 outputRefs.finalized_psf_ap_corr_detector_cat,
1370 )
1371 butlerQC.put(
1372 struct.output_table,
1373 outputRefs.finalized_src_detector_table,
1374 )
1376 def run(self, visit, band, detector, isolated_star_cat_dict, isolated_star_source_dict,
1377 src, exposure, fgcm_standard_star_dict=None):
1378 """
1379 Run the FinalizeCharacterizationDetectorTask.
1381 Parameters
1382 ----------
1383 visit : `int`
1384 Visit number. Used in the output catalogs.
1385 band : `str`
1386 Band name. Used to select reserved stars.
1387 detector : `int`
1388 Detector number.
1389 isolated_star_cat_dict : `dict`
1390 Per-tract dict of isolated star catalog handles.
1391 isolated_star_source_dict : `dict`
1392 Per-tract dict of isolated star source catalog handles.
1393 src : `lsst.afw.table.SourceCatalog`
1394 Src catalog.
1395 exposure : `lsst.afw.image.Exposure`
1396 Calexp exposure.
1397 fgcm_standard_star_dict : `dict`
1398 Per tract dict of fgcm isolated stars.
1400 Returns
1401 -------
1402 struct : `lsst.pipe.base.struct`
1403 Struct with outputs for persistence.
1405 Raises
1406 ------
1407 NoWorkFound
1408 Raised if the selector returns no good sources.
1409 """
1410 # Initialize schema, mappers, and subtasks if necessary.
1411 if self.schema is None: 1411 ↛ 1412line 1411 didn't jump to line 1412 because the condition on line 1411 was never true
1412 self.schema_mapper, self.schema = self._make_output_schema_mapper(src.schema)
1414 self.makeSubtask("measurement", schema=self.schema)
1415 self.makeSubtask("measure_ap_corr", schema=self.schema)
1416 self.makeSubtask("apply_ap_corr", schema=self.schema)
1418 # We do not need the isolated star table in this task.
1419 # However, it is used in tests to confirm consistency of indexes.
1420 _, isolated_source_table = self.concat_isolated_star_cats(
1421 band,
1422 isolated_star_cat_dict,
1423 isolated_star_source_dict,
1424 visit=visit,
1425 detector=detector,
1426 )
1428 if isolated_source_table is None: 1428 ↛ 1429line 1428 didn't jump to line 1429 because the condition on line 1428 was never true
1429 raise pipeBase.NoWorkFound(f'No good isolated sources found for any detectors in visit {visit}')
1431 exposure_cat_schema = afwTable.ExposureTable.makeMinimalSchema()
1432 exposure_cat_schema.addField('visit', type='L', doc='Visit number')
1434 metadata = dafBase.PropertyList()
1435 metadata.add("COMMENT", "Catalog id is detector id, sorted.")
1436 metadata.add("COMMENT", "Only one detector with data has an entry.")
1438 psf_ap_corr_cat = afwTable.ExposureCatalog(exposure_cat_schema)
1439 psf_ap_corr_cat.setMetadata(metadata)
1441 self.log.info("Starting finalizeCharacterization on detector ID %d.", detector)
1443 if fgcm_standard_star_dict is not None: 1443 ↛ 1453line 1443 didn't jump to line 1453 because the condition on line 1443 was always true
1444 fgcm_standard_star_cat = []
1446 for tract in fgcm_standard_star_dict:
1447 astropy_fgcm = fgcm_standard_star_dict[tract].get()
1448 table_fgcm = np.asarray(astropy_fgcm)
1449 fgcm_standard_star_cat.append(table_fgcm)
1451 fgcm_standard_star_cat = np.concatenate(fgcm_standard_star_cat)
1452 else:
1453 fgcm_standard_star_cat = None
1455 psf, ap_corr_map, measured_src = self.compute_psf_and_ap_corr_map(
1456 visit,
1457 detector,
1458 exposure,
1459 src,
1460 isolated_source_table,
1461 fgcm_standard_star_cat,
1462 )
1464 # And now we package it together...
1465 measured_src_table = None
1466 if measured_src is not None: 1466 ↛ 1480line 1466 didn't jump to line 1480 because the condition on line 1466 was always true
1467 record = psf_ap_corr_cat.addNew()
1468 record['id'] = int(detector)
1469 record['visit'] = visit
1470 if psf is not None: 1470 ↛ 1472line 1470 didn't jump to line 1472 because the condition on line 1470 was always true
1471 record.setPsf(psf)
1472 if ap_corr_map is not None: 1472 ↛ 1475line 1472 didn't jump to line 1475 because the condition on line 1472 was always true
1473 record.setApCorrMap(ap_corr_map)
1475 measured_src['visit'][:] = visit
1476 measured_src['detector'][:] = detector
1478 measured_src_table = measured_src.asAstropy()
1480 if measured_src_table is None: 1480 ↛ 1481line 1480 didn't jump to line 1481 because the condition on line 1480 was never true
1481 raise pipeBase.NoWorkFound(f'No good sources found for visit {visit} / detector {detector}')
1483 return pipeBase.Struct(
1484 psf_ap_corr_cat=psf_ap_corr_cat,
1485 output_table=measured_src_table,
1486 )
1489class ConsolidateFinalizeCharacterizationDetectorConnections(
1490 pipeBase.PipelineTaskConnections,
1491 dimensions=('instrument', 'visit',),
1492):
1493 finalized_psf_ap_corr_detector_cats = pipeBase.connectionTypes.Input(
1494 doc='Per-visit/per-detector finalized psf models and aperture corrections.',
1495 name='finalized_psf_ap_corr_detector_catalog',
1496 storageClass='ExposureCatalog',
1497 dimensions=('instrument', 'visit', 'detector'),
1498 multiple=True,
1499 deferLoad=True,
1500 )
1501 finalized_src_detector_tables = pipeBase.connectionTypes.Input(
1502 doc=('Per-visit/per-detector catalog of measurements for psf/flag/etc.'),
1503 name='finalized_src_detector_table',
1504 storageClass='ArrowAstropy',
1505 dimensions=('instrument', 'visit', 'detector'),
1506 multiple=True,
1507 deferLoad=True,
1508 )
1509 finalized_psf_ap_corr_cat = pipeBase.connectionTypes.Output(
1510 doc=('Per-visit finalized psf models and aperture corrections. This '
1511 'catalog uses detector id for the id and are sorted for fast '
1512 'lookups of a detector.'),
1513 name='finalized_psf_ap_corr_catalog',
1514 storageClass='ExposureCatalog',
1515 dimensions=('instrument', 'visit'),
1516 )
1517 finalized_src_table = pipeBase.connectionTypes.Output(
1518 doc=('Per-visit catalog of measurements for psf/flag/etc.'),
1519 name='finalized_src_table',
1520 storageClass='ArrowAstropy',
1521 dimensions=('instrument', 'visit'),
1522 )
1525class ConsolidateFinalizeCharacterizationDetectorConfig(
1526 pipeBase.PipelineTaskConfig,
1527 pipelineConnections=ConsolidateFinalizeCharacterizationDetectorConnections,
1528):
1529 pass
1532class ConsolidateFinalizeCharacterizationDetectorTask(pipeBase.PipelineTask):
1533 """Consolidate per-detector finalize characterization catalogs."""
1534 ConfigClass = ConsolidateFinalizeCharacterizationDetectorConfig
1535 _DefaultName = 'consolidate_finalize_characterization_detector'
1537 def runQuantum(self, butlerQC, inputRefs, outputRefs):
1538 input_handle_dict = butlerQC.get(inputRefs)
1540 psf_ap_corr_detector_dict_temp = {
1541 handle.dataId['detector']: handle
1542 for handle in input_handle_dict['finalized_psf_ap_corr_detector_cats']
1543 }
1544 src_detector_table_dict_temp = {
1545 handle.dataId['detector']: handle
1546 for handle in input_handle_dict['finalized_src_detector_tables']
1547 }
1549 psf_ap_corr_detector_dict = {
1550 detector: psf_ap_corr_detector_dict_temp[detector]
1551 for detector in sorted(psf_ap_corr_detector_dict_temp.keys())
1552 }
1553 src_detector_table_dict = {
1554 detector: src_detector_table_dict_temp[detector]
1555 for detector in sorted(src_detector_table_dict_temp.keys())
1556 }
1558 result = self.run(
1559 psf_ap_corr_detector_dict=psf_ap_corr_detector_dict,
1560 src_detector_table_dict=src_detector_table_dict,
1561 )
1563 butlerQC.put(result.psf_ap_corr_cat, outputRefs.finalized_psf_ap_corr_cat)
1564 butlerQC.put(result.output_table, outputRefs.finalized_src_table)
1566 def run(self, *, psf_ap_corr_detector_dict, src_detector_table_dict):
1567 """
1568 Run the ConsolidateFinalizeCharacterizationDetectorTask.
1570 Parameters
1571 ----------
1572 psf_ap_corr_detector_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
1573 Dictionary of input exposure catalogs, keyed by detector id.
1574 src_detector_table_dict : `dict` [`int`, `lsst.daf.butler.DeferredDatasetHandle`]
1575 Dictionary of input source tables, keyed by detector id.
1577 Returns
1578 -------
1579 result : `lsst.pipe.base.struct`
1580 Struct with the following outputs:
1581 ``psf_ap_corr_cat``: Consolidated exposure catalog
1582 ``src_table``: Consolidated source table.
1583 """
1584 if not len(psf_ap_corr_detector_dict): 1584 ↛ 1585line 1584 didn't jump to line 1585 because the condition on line 1584 was never true
1585 raise pipeBase.NoWorkFound("No inputs found.")
1587 if not np.all( 1587 ↛ 1591line 1587 didn't jump to line 1591 because the condition on line 1587 was never true
1588 np.asarray(psf_ap_corr_detector_dict.keys())
1589 == np.asarray(src_detector_table_dict.keys())
1590 ):
1591 raise ValueError(
1592 "Input psf_ap_corr_detector_dict and src_detector_table_dict must have the same keys",
1593 )
1595 psf_ap_corr_cat = None
1596 for detector_id, handle in psf_ap_corr_detector_dict.items():
1597 if psf_ap_corr_cat is None:
1598 psf_ap_corr_cat = handle.get()
1599 else:
1600 psf_ap_corr_cat.append(handle.get().find(detector_id))
1602 # Make sure it is a contiguous catalog.
1603 psf_ap_corr_cat = psf_ap_corr_cat.copy(deep=True)
1605 src_table = TableVStack.vstack_handles(src_detector_table_dict.values())
1607 return pipeBase.Struct(
1608 psf_ap_corr_cat=psf_ap_corr_cat,
1609 output_table=src_table,
1610 )