22 """Read preprocessed bright stars and stack them to build an extended
26 from dataclasses
import dataclass
27 from typing
import List
29 from lsst.afw import image
as afwImage
41 """Single extended PSF over a focal plane region.
43 The focal plane region is defined through a list
48 extended_psf_image : `lsst.afw.image.MaskedImageF`
49 Image of the extended PSF model.
50 detector_list : `list` [`int`]
51 List of detector IDs that define the focal plane region over which this
52 extended PSF model has been built (and can be used).
54 extended_psf_image: afwImage.MaskedImageF
55 detector_list: List[int]
59 """Extended PSF model.
61 Each instance may contain a default extended PSF, a set of extended PSFs
62 that correspond to different focal plane regions, or both. At this time,
63 focal plane regions are always defined as a subset of detectors.
67 default_extended_psf : `lsst.afw.image.MaskedImageF`
68 Extended PSF model to be used as default (or only) extended PSF model.
76 """Add a new focal plane region, along wit hits extended PSF, to the
81 extended_psf_image : `lsst.afw.image.MaskedImageF`
82 Extended PSF model for the region.
84 Name of the focal plane region. Will be converted to all-uppercase.
85 detector_list : `list` [`int`]
86 List of IDs for the detectors that define the focal plane region.
88 region_name = region_name.upper()
90 raise ValueError(f
"Region name {region_name} is already used by this ExtendedPsf instance.")
92 extended_psf_image=extended_psf_image, detector_list=detector_list)
93 for det
in detector_list:
97 """Return the appropriate extended PSF.
99 If the instance contains no extended PSF defined over focal plane
100 regions, the default extended PSF will be returned regardless of
101 whether a detector ID was passed as argument.
105 detector : `int`, optional
106 Detector ID. If focal plane region PSFs are defined, is used to
107 determine which model to return.
111 extendedPsfImage : `lsst.afw.image.MaskedImageF`
112 The extended PSF model. If this instance contains extended PSFs
113 defined over focal plane regions, the extended PSF model for the
114 region that contains ``detector`` is returned. If not, the default
115 extended PSF is returned.
119 raise ValueError(
"No default extended PSF available; please provide detector number.")
126 """Returns the number of extended PSF models present in the instance.
128 Note that if the instance contains both a default model and a set of
129 focal plane region models, the length of the instance will be the
130 number of regional models, plus one (the default). This is true even
131 in the case where the default model is one of the focal plane
132 region-specific models.
140 """Returns the extended PSF for a focal plane region.
142 The region can be identified either by name, or through a detector ID.
146 region_name : `str` or `None`, optional
147 Name of the region for which the extended PSF should be retrieved.
148 Ignored if ``detector`` is provided. Must be provided if
149 ``detector`` is None.
150 detector : `int` or `None`, optional
151 If provided, returns the extended PSF for the focal plane region
152 that includes this detector.
157 Raised if neither ``detector`` nor ``regionName`` is provided.
160 if region_name
is None:
161 raise ValueError(
"One of either a regionName or a detector number must be provided.")
166 """Write this object to a file.
171 Name of file to write.
177 metadata[
"HAS_REGIONS"] =
True
180 metadata[region] = e_psf_region.detector_list
182 metadata[
"HAS_REGIONS"] =
False
183 fits_primary = afwFits.Fits(filename,
"w")
184 fits_primary.createEmpty()
185 fits_primary.writeMetadata(metadata)
186 fits_primary.closeFile()
190 default_hdu_metadata.update({
"REGION":
"DEFAULT",
"EXTNAME":
"IMAGE"})
191 self.
default_extended_psfdefault_extended_psf.image.writeFits(filename, metadata=default_hdu_metadata, mode=
"a")
192 default_hdu_metadata.update({
"REGION":
"DEFAULT",
"EXTNAME":
"MASK"})
193 self.
default_extended_psfdefault_extended_psf.mask.writeFits(filename, metadata=default_hdu_metadata, mode=
"a")
197 metadata.update({
"REGION": region,
"EXTNAME":
"IMAGE"})
198 e_psf_region.extended_psf_image.image.writeFits(filename, metadata=metadata, mode=
"a")
199 metadata.update({
"REGION": region,
"EXTNAME":
"MASK"})
200 e_psf_region.extended_psf_image.mask.writeFits(filename, metadata=metadata, mode=
"a")
203 """Alias for ``write_fits``; exists for compatibility with the Butler.
209 """Build an instance of this class from a file.
214 Name of the file to read.
217 global_metadata = afwFits.readMetadata(filename, hdu=0)
218 has_default = global_metadata.getBool(
"HAS_DEFAULT")
219 if global_metadata.getBool(
"HAS_REGIONS"):
220 focal_plane_region_names = global_metadata.getArray(
"REGION_NAMES")
222 focal_plane_region_names = []
223 f = afwFits.Fits(filename,
"r")
224 n_extensions = f.countHdus()
225 extended_psf_parts = {}
226 for j
in range(1, n_extensions):
227 md = afwFits.readMetadata(filename, hdu=j)
228 if has_default
and md[
"REGION"] ==
"DEFAULT":
229 if md[
"EXTNAME"] ==
"IMAGE":
230 default_image = afwImage.ImageF(filename, hdu=j)
231 elif md[
"EXTNAME"] ==
"MASK":
232 default_mask = afwImage.MaskX(filename, hdu=j)
234 if md[
"EXTNAME"] ==
"IMAGE":
235 extended_psf_part = afwImage.ImageF(filename, hdu=j)
236 elif md[
"EXTNAME"] ==
"MASK":
237 extended_psf_part = afwImage.MaskX(filename, hdu=j)
238 extended_psf_parts.setdefault(md[
"REGION"], {})[md[
"EXTNAME"].lower()] = extended_psf_part
241 extended_psf = cls(afwImage.MaskedImageF(default_image, default_mask))
245 if len(extended_psf_parts) != len(focal_plane_region_names):
246 raise ValueError(f
"Number of per-region extended PSFs read ({len(extended_psf_parts)}) does not "
247 "match with the number of regions recorded in the metadata "
248 f
"({len(focal_plane_region_names)}).")
250 for r_name
in focal_plane_region_names:
251 extended_psf_image = afwImage.MaskedImageF(**extended_psf_parts[r_name])
252 detector_list = global_metadata.getArray(r_name)
253 extended_psf.add_regional_extended_psf(extended_psf_image, r_name, detector_list)
259 """Alias for ``readFits``; exists for compatibility with the Butler.
265 """Configuration parameters for StackBrightStarsTask.
267 subregion_size = pexConfig.ListField(
269 doc=
"Size, in pixels, of the subregions over which the stacking will be "
270 "iteratively performed.",
273 stacking_statistic = pexConfig.ChoiceField(
275 doc=
"Type of statistic to use for stacking.",
280 "MEANCLIP":
"clipped mean",
283 num_sigma_clip = pexConfig.Field(
285 doc=
"Sigma for outlier rejection; ignored if stacking_statistic != 'MEANCLIP'.",
288 num_iter = pexConfig.Field(
290 doc=
"Number of iterations of outlier rejection; ignored if stackingStatistic != 'MEANCLIP'.",
293 bad_mask_planes = pexConfig.ListField(
295 doc=
"Mask planes that, if set, lead to associated pixels not being included in the stacking of the "
296 "bright star stamps.",
297 default=(
'BAD',
'CR',
'CROSSTALK',
'EDGE',
'NO_DATA',
'SAT',
'SUSPECT',
'UNMASKEDNAN')
299 do_mag_cut = pexConfig.Field(
301 doc=
"Apply magnitude cut before stacking?",
304 mag_limit = pexConfig.Field(
306 doc=
"Magnitude limit, in Gaia G; all stars brighter than this value will be stacked",
312 """Stack bright stars together to build an extended PSF model.
314 ConfigClass = StackBrightStarsConfig
315 _DefaultName =
"stack_bright_stars"
317 def __init__(self, initInputs=None, *args, **kwargs):
318 pipeBase.CmdLineTask.__init__(self, *args, **kwargs)
320 def _set_up_stacking(self, example_stamp):
321 """Configure stacking statistic and control from config fields.
324 stats_control.setNumSigmaClip(self.config.num_sigma_clip)
325 stats_control.setNumIter(self.config.num_iter)
326 if bad_masks := self.config.bad_mask_planes:
327 and_mask = example_stamp.mask.getPlaneBitMask(bad_masks[0])
328 for bm
in bad_masks[1:]:
329 and_mask = and_mask | example_stamp.mask.getPlaneBitMask(bm)
330 stats_control.setAndMask(and_mask)
332 return stats_control, stats_flags
334 def run(self, bss_ref_list, region_name=None):
335 """Read input bright star stamps and stack them together.
337 The stacking is done iteratively over smaller areas of the final model
338 image to allow for a great number of bright star stamps to be used.
342 bss_ref_list : `list` of
343 `lsst.daf.butler._deferredDatasetHandle.DeferredDatasetHandle`
344 List of available bright star stamps data references.
345 region_name : `str`, optional
346 Name of the focal plane region, if applicable. Only used for
347 logging purposes, when running over multiple such regions
348 (typically from `MeasureExtendedPsfTask`)
351 region_message = f
' for region "{region_name}".'
354 self.log.
info(
'Building extended PSF from stamps extracted from %d detector images%s',
355 len(bss_ref_list), region_message)
357 example_bss = bss_ref_list[0].get(datasetType=
"brightStarStamps", immediate=
True)
358 example_stamp = example_bss[0].stamp_im
360 ext_psf = afwImage.MaskedImageF(example_stamp.getBBox())
362 subregion_size =
Extent2I(*self.config.subregion_size)
363 sub_bboxes = AssembleCoaddTask._subBBoxIter(ext_psf.getBBox(), subregion_size)
365 n_subregions = int(ext_psf.getDimensions()[0]/subregion_size[0] + 1)*int(
366 ext_psf.getDimensions()[1]/subregion_size[1] + 1)
367 self.log.
info(
"Stacking will performed iteratively over approximately %d "
368 "smaller areas of the final model image.", n_subregions)
370 stats_control, stats_flags = self.
_set_up_stacking_set_up_stacking(example_stamp)
372 for jbbox, bbox
in enumerate(sub_bboxes):
374 for bss_ref
in bss_ref_list:
375 read_stars = bss_ref.get(datasetType=
"brightStarStamps", parameters={
'bbox': bbox})
376 if self.config.do_mag_cut:
377 read_stars = read_stars.selectByMag(magMax=self.config.mag_limit)
379 all_stars.extend(read_stars)
381 all_stars = read_stars
384 ext_psf.assign(coadd_sub_bbox, bbox)
389 dimensions=(
"band",
"instrument")):
390 input_brightStarStamps = pipeBase.connectionTypes.Input(
391 doc=
"Input list of bright star collections to be stacked.",
392 name=
"brightStarStamps",
393 storageClass=
"BrightStarStamps",
394 dimensions=(
"visit",
"detector"),
398 extended_psf = pipeBase.connectionTypes.Output(
399 doc=
"Extended PSF model built by stacking bright stars.",
401 storageClass=
"ExtendedPsf",
402 dimensions=(
"band",),
407 pipelineConnections=MeasureExtendedPsfConnections):
408 """Configuration parameters for MeasureExtendedPsfTask.
410 stack_bright_stars = pexConfig.ConfigurableField(
411 target=StackBrightStarsTask,
412 doc=
"Stack selected bright stars",
414 detectors_focal_plane_regions = pexConfig.DictField(
417 doc=
"Mapping from detector IDs to focal plane region names. If empty, a constant "
418 "extended PSF model is built from all selected bright stars.",
424 """Build and save extended PSF model.
426 The model is built by stacking bright star stamps, extracted and
428 `lsst.pipe.tasks.processBrightStars.ProcessBrightStarsTask`.
429 If a mapping from detector IDs to focal plane regions is provided,
430 a different extended PSF model will be built for each focal plane
431 region. If not, a single, constant extended PSF model is built using
434 ConfigClass = MeasureExtendedPsfConfig
435 _DefaultName =
"measureExtendedPsf"
437 def __init__(self, initInputs=None, *args, **kwargs):
438 pipeBase.CmdLineTask.__init__(self, *args, **kwargs)
439 self.makeSubtask(
"stack_bright_stars")
441 set(self.config.detectors_focal_plane_regions.values())}
442 for det, region
in self.config.detectors_focal_plane_regions.items():
451 """Split available sets of bright star stamps according to focal plane
457 `lsst.daf.butler._deferredDatasetHandle.DeferredDatasetHandle`
458 List of available bright star stamps data references.
461 for dataset_handle
in ref_list:
462 det_id = dataset_handle.ref.dataId[
"detector"]
466 region_name = self.config.detectors_focal_plane_regions[det_id]
468 self.log.
warning(
'Bright stars were available for detector %d, but it was missing '
469 'from the "detectors_focal_plane_regions" config field, so they will not '
470 'be used to build any of the extended PSF models', det_id)
473 region_ref_list[region_name].
append(dataset_handle)
474 return region_ref_list
477 input_data = butlerQC.get(inputRefs)
478 bss_ref_list = input_data[
'input_brightStarStamps']
480 if not self.config.detectors_focal_plane_regions:
481 self.log.
info(
"No detector groups were provided to MeasureExtendedPsfTask; computing a single, "
482 "constant extended PSF model over all available observations.")
483 output_e_psf =
ExtendedPsf(self.stack_bright_stars.
run(bss_ref_list))
487 for region_name, ref_list
in region_ref_list.items():
490 self.log.
warning(
'No valid brightStarStamps reference found for region "%s"; '
491 'skipping it.', region_name)
493 ext_psf = self.stack_bright_stars.
run(ref_list, region_name)
494 output_e_psf.add_regional_extended_psf(ext_psf, region_name,
496 output = pipeBase.Struct(extended_psf=output_e_psf)
497 butlerQC.put(output, outputRefs)
std::vector< SchemaItem< Flag > > * items
Pass parameters to a Statistics object.
Class for storing ordered metadata with comments.
def writeFits(self, filename)
def readFits(cls, filename)
detectors_focal_plane_regions
def __call__(self, detector=None)
def get_regional_extended_psf(self, region_name=None, detector=None)
def add_regional_extended_psf(self, extended_psf_image, region_name, detector_list)
def read_fits(cls, filename)
def write_fits(self, filename)
def __init__(self, default_extended_psf=None)
def select_detector_refs(self, ref_list)
def __init__(self, initInputs=None, *args, **kwargs)
def runQuantum(self, butlerQC, inputRefs, outputRefs)
def _set_up_stacking(self, example_stamp)
def run(self, bss_ref_list, region_name=None)
def __init__(self, initInputs=None, *args, **kwargs)
daf::base::PropertyList * list
daf::base::PropertySet * set
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)
std::shared_ptr< lsst::afw::image::Image< PixelT > > statisticsStack(std::vector< std::shared_ptr< lsst::afw::image::Image< PixelT >>> &images, Property flags, StatisticsControl const &sctrl=StatisticsControl(), std::vector< lsst::afw::image::VariancePixel > const &wvector=std::vector< lsst::afw::image::VariancePixel >(0))
A function to compute some statistics of a stack of Images.
Extent< int, 2 > Extent2I
def run(self, coaddExposures, bbox, wcs)