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