Coverage for metadetect / lsst / measure.py: 17%
310 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-05 01:45 -0800
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-05 01:45 -0800
1import logging
2import warnings
3import numpy as np
4import esutil as eu
5import ngmix
7import lsst.afw.table as afw_table
8from lsst.meas.algorithms import SourceDetectionTask
9from lsst.meas.algorithms import (
10 SourceDetectionConfig as OriginalSourceDetectionConfig,
11)
12from lsst.meas.deblender import SourceDeblendTask, SourceDeblendConfig
13from lsst.meas.base import (
14 SingleFrameMeasurementConfig,
15 SingleFrameMeasurementTask,
16)
17from lsst.pex.config import Config, ConfigurableField, Field
18from lsst.pipe.base import Task
19import lsst.afw.image as afw_image
20import lsst.geom as geom
21from lsst.pex.exceptions import (
22 InvalidParameterError,
23 LengthError,
24)
26from ..procflags import (
27 EDGE_HIT,
28 ZERO_WEIGHTS,
29 NO_ATTEMPT,
30)
31from ..fitting import (
32 fit_mbobs_gauss,
33 fit_mbobs_wavg,
34 get_wavg_output_struct,
35)
37from . import util
38from .util import ContextNoiseReplacer, get_jacobian
39from . import vis
40from .defaults import DEFAULT_THRESH
42warnings.filterwarnings('ignore', category=FutureWarning)
44LOG = logging.getLogger('lsst_measure')
47class SourceDetectionConfig(OriginalSourceDetectionConfig):
48 @property
49 def thresh(self):
50 return self.thresholdValue
52 @thresh.setter
53 def thresh(self, value):
54 self.thresholdValue = value
57SourceDetectionTask.ConfigClass = SourceDetectionConfig
60class DetectAndDeblendConfig(Config):
61 meas = ConfigurableField[SingleFrameMeasurementConfig](
62 doc="Measurement config",
63 target=SingleFrameMeasurementTask,
64 )
66 detect = ConfigurableField[SourceDetectionConfig](
67 doc="Detection config", target=SourceDetectionTask
68 )
70 deblend = ConfigurableField[SourceDeblendConfig](
71 doc="Deblend config", target=SourceDeblendTask
72 )
74 seed = Field[int](
75 doc="Random rng seed",
76 default=42,
77 )
79 def setDefaults(self):
80 super().setDefaults()
82 # defaults for measurement config
83 self.meas.plugins.names = [
84 "base_SdssCentroid",
85 "base_PsfFlux",
86 "base_SkyCoord",
87 ]
89 # set these slots to none because we aren't running these algorithms
90 self.meas.slots.apFlux = None
91 self.meas.slots.gaussianFlux = None
92 self.meas.slots.calibFlux = None
93 self.meas.slots.modelFlux = None
95 # goes with SdssShape above
96 self.meas.slots.shape = None
98 # fix odd issue where it things things are near the edge
99 self.meas.plugins['base_SdssCentroid'].binmax = 1
101 # defaults for detection config
102 # DM does not have config default stability. Set all of them
103 # explicitly
104 self.detect.minPixels = 1
105 self.detect.isotropicGrow = True
106 self.detect.combinedGrow = True
107 self.detect.nSigmaToGrow = 2.4
108 self.detect.returnOriginalFootprints = False
109 self.detect.includeThresholdMultiplier = 1.0
110 self.detect.thresholdPolarity = "positive"
111 self.detect.adjustBackground = 0.0
112 self.detect.doApplyFlatBackgroundRatio = False
113 # these are ignored since we are doing reEstimateBackground = False
114 # detection_config.background
115 # detection_config.tempLocalBackground
116 # detection_config.doTempLocalBackground
117 # detection_config.tempWideBackground
118 # detection_config.doTempWideBackground
119 self.detect.nPeaksMaxSimple = 1
120 self.detect.nSigmaForKernel = 7.0
121 self.detect.excludeMaskPlanes = util.get_detection_mask()
123 # the defaults changed from from stdev to pixel_std but
124 # we don't want that
126 self.detect.thresholdType = "stdev"
127 # our changes from defaults
128 self.detect.reEstimateBackground = False
130 self.detect.thresholdValue = DEFAULT_THRESH
132 # these will be ignored when finding the image standard deviation
133 self.detect.statsMask = util.get_stats_mask()
135 # deblend config
136 # these tasks must use the same schema and all be constructed before
137 # any other tasks using the same schema are run because schema is
138 # modified in place by tasks, and the constructor does a check that
139 # fails if we do this afterward
140 self.deblend.maxFootprintArea = 0
143class DetectAndDeblendTask(Task):
144 ConfigClass = DetectAndDeblendConfig
145 _DefaultName = "detect_and_deblend"
147 def __init__(self, *args, **kwargs):
148 super().__init__(*args, **kwargs)
149 schema = afw_table.SourceTable.makeMinimalSchema()
150 self.makeSubtask("meas", schema=schema)
151 self.makeSubtask("detect", schema=schema)
152 self.makeSubtask("deblend", schema=schema)
153 self.rng = np.random.RandomState(seed=self.config.seed)
155 def run(self, mbexp, show=False):
156 if len(mbexp.singles) > 1:
157 detexp = util.coadd_exposures(mbexp.singles)
158 if detexp is None:
159 return [], None
160 else:
161 detexp = mbexp.singles[0]
163 # background measurement within the detection code requires ExposureF
164 if not isinstance(detexp, afw_image.ExposureF):
165 detexp = afw_image.ExposureF(detexp, deep=True)
167 schema = self.deblend.schema # should be the same for all tasks
168 table = afw_table.SourceTable.make(schema)
169 result = self.detect.run(table, detexp)
171 if result is not None:
172 sources = result.sources
173 self.deblend.run(detexp, sources)
175 with ContextNoiseReplacer(
176 detexp,
177 sources,
178 self.rng,
179 config=self.meas.config.noiseReplacer,
180 ) as replacer:
181 for source in sources:
182 if source.get('deblend_nChild') != 0:
183 continue
185 source_id = source.getId()
187 with replacer.sourceInserted(source_id):
188 self.meas.callMeasure(source, detexp)
190 else:
191 sources = []
193 if show:
194 vis.show_exp(detexp, use_mpl=True, sources=sources)
196 return sources, detexp
199def detect_and_deblend(
200 mbexp,
201 rng=None,
202 thresh=DEFAULT_THRESH,
203 show=False,
204 config=None,
205):
206 """
207 run detection and deblending of peaks, as well as basic measurments such as
208 centroid. The SDSS deblender is run in order to split footprints.
210 We must combine detection and deblending in the same function because the
211 schema gets modified in place, which means we must construct the deblend
212 task at the same time as the detect task
214 Parameters
215 ----------
216 mbexp: lsst.afw.image.MultibandExposure
217 The exposures to process
218 rng: np.random.RandomState
219 Random number generator for noise replacer
220 thresh: float, optional
221 The detection threshold in units of the sky noise
222 show: bool, optional
223 If set to True, show images
225 Returns
226 -------
227 sources, detexp
228 The sources and the detection exposure
229 """
230 config_override = config if config is not None else {}
231 if thresh:
232 if 'detect' not in config_override:
233 config_override['detect'] = {}
234 config_override['detect']['thresholdValue'] = thresh
236 config = DetectAndDeblendConfig()
237 config.setDefaults()
239 util.override_config(config, config_override)
241 config.freeze()
242 config.validate()
243 task = DetectAndDeblendTask(config=config)
244 if rng is not None:
245 task.rng = rng
246 return task.run(mbexp, show)
249def measure(
250 mbexp,
251 detexp,
252 sources,
253 config,
254 rng,
255):
256 """
257 run measurements on the input exposure, given the input measurement task,
258 list of sources, and config.
260 Parameters
261 ----------
262 mbexp: lsst.afw.image.MultibandExposure
263 The exposure on which to detect and measure
264 detexp: lsst.afw.image.Exposure*
265 The detection exposure, used for bmask info
266 sources: list of sources
267 From a detection task
268 config: dict, optional
269 Configuration for the fitter, metacal, psf, detect, Entries
270 in this dict override defaults; see lsst_configs.py
271 rng: np.random.RandomState
272 Random number generator
274 Returns
275 -------
276 ndarray with results or None
277 """
279 if len(sources) == 0:
280 return None
282 pgauss_fitter = get_pgauss_fitter(config)
284 nband = len(mbexp.bands)
285 shear_band_names = config["shear_bands"] or mbexp.bands
286 if not all([sb in mbexp.bands for sb in shear_band_names]):
287 raise RuntimeError(
288 "Not all requested bands for shear are available. "
289 f"Bands `{shear_band_names}` were requested but the only "
290 f"bands available are `{mbexp.bands}`."
291 )
292 shear_bands = [
293 i for i, band in enumerate(mbexp.bands) if band in shear_band_names
294 ]
295 exp_bbox = mbexp.getBBox()
296 wcs = mbexp.singles[0].getWcs()
297 results = []
299 # bmasks will be different within the loop below due to the replacer
300 bmasks = get_bmasks(sources=sources, exposure=detexp)
302 for i, source in enumerate(sources):
303 if source.get('deblend_nChild') != 0:
304 continue
306 bmask = bmasks[i]
308 stamp_flags = 0
309 try:
310 mbobs = _get_stamp_mbobs(
311 mbexp=mbexp,
312 source=source,
313 stamp_size=config['stamp_size'],
314 )
315 except LengthError as err:
316 # This is raised when a bbox hits an edge
317 LOG.debug('%s', err)
318 stamp_flags = EDGE_HIT
319 except AllZeroWeightError as err:
320 # failure creating some observation due to all weights being zero
321 # across the stamp
322 LOG.info('%s', err)
323 stamp_flags = ZERO_WEIGHTS
325 if stamp_flags != 0:
327 this_gauss_res = get_wavg_output_struct(
328 nband=nband,
329 model='gauss',
330 shear_bands=shear_bands,
331 )
332 this_pgauss_res = get_wavg_output_struct(
333 nband=nband,
334 model='pgauss',
335 )
337 else:
339 # TODO do something with bmask_flags?
340 psf_fitter = ngmix.admom.AdmomFitter(rng=rng)
341 psf_guesser = ngmix.guessers.GMixPSFGuesser(
342 rng=rng,
343 ngauss=1,
344 guess_from_moms=True,
345 )
346 psf_runner = ngmix.runners.PSFRunner(
347 fitter=psf_fitter,
348 guesser=psf_guesser,
349 ntry=4,
350 )
352 this_gauss_res = fit_mbobs_gauss(
353 mbobs=mbobs,
354 bmask_flags=0,
355 psf_runner=psf_runner,
356 rng=rng,
357 shear_bands=shear_bands,
358 )
360 this_pgauss_res = fit_mbobs_wavg(
361 mbobs=mbobs,
362 fitter=pgauss_fitter,
363 bmask_flags=0,
364 )
366 this_res = _get_combined_struct(this_gauss_res, this_pgauss_res)
367 this_res['stamp_flags'] = stamp_flags
369 res = get_output(
370 wcs=wcs,
371 source=source,
372 res=this_res,
373 bmask=bmask,
374 stamp_size=config['stamp_size'],
375 exp_bbox=exp_bbox,
376 )
378 results.append(res)
380 if len(results) > 0:
381 results = eu.numpy_util.combine_arrlist(results)
382 else:
383 results = None
385 return results
388def get_pgauss_fitter(config):
389 return ngmix.prepsfmom.PGaussMom(fwhm=config['pgauss']['fwhm'])
392def _get_combined_struct(gauss_res, pgauss_res):
394 skip = ['pgauss_g', 'pgauss_g_cov', 'shear_bands']
395 keep_dt = [('stamp_flags', 'i4')]
397 for pdt in pgauss_res.dtype.descr:
399 n = pdt[0]
401 if 'psf' in n or n in skip:
402 continue
404 keep_dt.append(pdt)
406 out = eu.numpy_util.add_fields(gauss_res, keep_dt)
408 for n in pgauss_res.dtype.names:
409 if n in out.dtype.names and n != "shear_bands":
410 out[n] = pgauss_res[n]
412 return out
415def get_bmasks(sources, exposure):
416 """
417 get a list of all the bmasks for the sources
419 Parameters
420 ----------
421 sources: lsst.afw.table.SourceCatalog
422 The sources
423 exposure: lsst.afw.image.ExposureF
424 The exposure
426 Returns
427 -------
428 list of bmask values
429 """
430 bmasks = []
431 for source in sources:
432 bmask = get_bmask(source=source, exposure=exposure)
433 bmasks.append(bmask)
434 return bmasks
437def get_bmask(source, exposure):
438 """
439 get bmask based on original peak position
441 Parameters
442 ----------
443 sources: lsst.afw.table.SourceRecord
444 The sources
445 exposure: lsst.afw.image.ExposureF
446 The exposure
448 Returns
449 -------
450 bmask value
451 """
452 peak = source.getFootprint().getPeaks()[0]
453 orig_cen = peak.getI()
454 maskval = exposure.mask[orig_cen]
455 return maskval
458def extract_obs(exp, source):
459 """
460 convert an image object into an ngmix.Observation, including
461 a psf observation
463 parameters
464 ----------
465 imobj: lsst.afw.image.ExposureF
466 The exposure
467 source: lsst.afw.table.SourceRecord
468 The source record
470 returns
471 --------
472 obs: ngmix.Observation
473 The Observation unless all the weight are zero, in which
474 case AllZeroWeightError is raised
475 """
477 im = exp.image.array
479 wt = _extract_weight(exp)
480 if np.all(wt <= 0):
481 raise AllZeroWeightError('all weights <= 0')
483 bmask = exp.mask.array
484 jacob = _extract_jacobian_at_source(
485 exp=exp,
486 source=source,
487 )
489 orig_cen = source.getCentroid()
491 psf_im = extract_psf_image(exposure=exp, orig_cen=orig_cen)
493 # fake the psf pixel noise
494 psf_err = psf_im.max() * 0.0001
495 psf_wt = psf_im * 0 + 1.0 / psf_err**2
497 # use canonical center for the psf
498 psf_cen = (np.array(psf_im.shape) - 1.0) / 2.0
499 psf_jacob = jacob.copy()
500 psf_jacob.set_cen(row=psf_cen[0], col=psf_cen[1])
502 # we will have need of the bit names which we can only
503 # get from the mask object
504 # this is sort of monkey patching, but I'm not sure of
505 # a better solution
507 meta = {'orig_cen': orig_cen}
509 psf_obs = ngmix.Observation(
510 psf_im,
511 weight=psf_wt,
512 jacobian=psf_jacob,
513 )
514 obs = ngmix.Observation(
515 im,
516 weight=wt,
517 bmask=bmask,
518 jacobian=jacob,
519 psf=psf_obs,
520 meta=meta,
521 )
523 return obs
526def _get_stamp_mbobs(mbexp, source, stamp_size, clip=False):
527 """
528 Get a postage stamp MultibandExposure
530 Parameters
531 ----------
532 mbexp: lsst.afw.image.MultibandExposure
533 The exposures
534 source: lsst.afw.table.SourceRecord
535 The source for which to get the stamp
536 stamp_size: int
537 If sent, a bounding box is created with about this size rather than
538 using the footprint bounding box. Typically the returned size is
539 stamp_size + 1
540 clip: bool, optional
541 If set to True, clip the bbox to fit into the exposure.
543 If clip is False and the bbox does not fit, a
544 lsst.pex.exceptions.LengthError is raised
546 Only relevant if stamp_size is sent. Default False
548 Returns
549 -------
550 lsst.afw.image.ExposureF
551 """
553 bbox = _get_bbox(mbexp, source, stamp_size, clip=clip)
555 mbobs = ngmix.MultiBandObsList()
556 for band in mbexp.bands:
557 subexp = mbexp[band][bbox]
558 obs = extract_obs(
559 exp=subexp,
560 source=source,
561 )
563 obslist = ngmix.ObsList(meta={'band': band})
564 obslist.append(obs)
565 mbobs.append(obslist)
567 return mbobs
570def _get_bbox(mbexp, source, stamp_size, clip=False):
571 """
572 Get a bounding box at the location of the specified source.
574 Parameters
575 ----------
576 mbexp: lsst.afw.image.MultibandExposure
577 The exposures
578 source: lsst.afw.table.SourceRecord
579 The source for which to get the stamp
580 stamp_size: int
581 If sent, a bounding box is created with about this size rather than
582 using the footprint bounding box. Typically the returned size is
583 stamp_size + 1
584 clip: bool, optional
585 If set to True, clip the bbox to fit into the exposure.
587 If clip is False and the bbox does not fit, a
588 lsst.pex.exceptions.LengthError is raised
590 Only relevant if stamp_size is sent. Default False
592 Returns
593 -------
594 lsst.geom.Box2I
595 """
597 fp = source.getFootprint()
598 peak = fp.getPeaks()[0]
600 x_peak, y_peak = peak.getIx(), peak.getIy()
602 bbox = geom.Box2I(
603 geom.Point2I(x_peak, y_peak),
604 geom.Extent2I(1, 1),
605 )
606 bbox.grow(stamp_size // 2)
608 exp_bbox = mbexp.getBBox()
609 if clip:
610 bbox.clip(exp_bbox)
611 else:
612 if not exp_bbox.contains(bbox):
613 source_id = source.getId()
614 raise LengthError(
615 f'requested stamp size {stamp_size} for source '
616 f'{source_id} does not fit into the exposoure. '
617 f'Use clip=True to clip the bbox to fit'
618 )
620 return bbox
623def extract_psf_image(exposure, orig_cen):
624 """
625 get the psf associated with this image.
627 coadded psfs from DM are generally not square, but the coadd in cells code
628 makes them so. We will assert they are square and odd dimensions
630 Parameters
631 ----------
632 exposure: lsst.afw.image.ExposureF
633 The exposure data
634 orig_cen: lsst.geom.Point2D
635 The location at which to draw the image
637 Returns
638 -------
639 ndarray
640 """
641 try:
642 psfobj = exposure.getPsf()
643 psfim = psfobj.computeKernelImage(orig_cen).array
644 except InvalidParameterError:
645 raise MissingDataError("could not reconstruct PSF")
647 psfim = np.asarray(psfim, dtype='f4')
649 shape = psfim.shape
650 assert shape[0] == shape[1], 'require square psf images'
651 assert shape[0] % 2 != 0, 'require odd psf images'
653 return psfim
656def _extract_psf_image_fix(exposure, orig_cen):
657 """
658 get the psf associated with this image
660 coadded psfs are generally not square, so we will
661 trim it to be square and preserve the center to
662 be at the new canonical center
664 TODO: should we really trim the psf to be even? will this
665 cause a shift due being off-center?
666 """
667 try:
668 psfobj = exposure.getPsf()
669 psfim = psfobj.computeKernelImage(orig_cen).array
670 except InvalidParameterError:
671 raise MissingDataError("could not reconstruct PSF")
673 psfim = np.asarray(psfim, dtype='f4')
675 psfim = util.trim_odd_image(psfim)
676 return psfim
679def _extract_weight(exp):
680 """
681 TODO get the estimated sky variance rather than this hack
682 TODO should we zero out other bits?
684 extract a weight map
686 Areas with NO_DATA will get zero weight.
688 Because the variance map includes the full poisson variance, which
689 depends on the signal, we instead extract the median of the parts of
690 the image without NO_DATA set
692 parameters
693 ----------
694 exp: sub exposure object
695 """
697 # TODO implement bit checking
698 var_image = exp.variance.array
700 weight = var_image.copy()
702 weight[:, :] = 0
704 wuse = np.where(var_image > 0)
706 if wuse[0].size > 0:
707 medvar = np.median(var_image[wuse])
708 weight[:, :] = 1.0 / medvar
709 else:
710 print(' weight is all zero, found none that passed cuts')
712 return weight
715def _extract_jacobian_at_source(exp, source):
716 """
717 extract an ngmix.Jacobian from the image object
718 and object record
720 exp: lsst.afw.image.ExposureF
721 An exposure object
722 source: lsst.afw.table.SourceRecord
723 The source record created during detection
725 returns
726 --------
727 Jacobian: ngmix.Jacobian
728 The local jacobian
729 """
731 orig_cen = exp.getWcs().skyToPixel(source.getCoord())
733 if np.isnan(orig_cen.getY()):
734 LOG.info('falling back on integer location')
735 # fall back to integer pixel location
736 peak = source.getFootprint().getPeaks()[0]
737 orig_cen_i = peak.getI()
738 orig_cen = geom.Point2D(
739 x=orig_cen_i.getX(),
740 y=orig_cen_i.getY(),
741 )
743 return get_jacobian(exp, orig_cen)
746def get_output_dtype():
747 dt = [
748 ('stamp_size', 'i4'),
749 ('row0', 'i4'), # bbox row start
750 ('col0', 'i4'), # bbox col start
751 ('row', 'f4'), # row in image. Use row0 to get to global pixel coords
752 ('col', 'f4'), # col in image. Use col0 to get to global pixel coords
753 ('row_diff', 'f4'), # difference from peak location
754 ('col_diff', 'f4'), # difference from peak location
755 ('ra', 'f8'),
756 ('dec', 'f8'),
757 ('psfrec_flags', 'i4'), # psfrec is the original psf
758 ('psfrec_g', 'f8', 2),
759 ('psfrec_T', 'f8'),
760 # values from .mask of input exposures
761 ('bmask', 'i4'),
762 # values for ormask across all input exposures to coadd
763 ('ormask', 'i4'),
764 # fraction of images going into a pixel that were masked
765 ('mfrac', 'f4'),
766 ]
768 return dt
771def get_output_struct(res):
772 """
773 get the output struct
775 Parameters
776 ----------
777 res: ndarray
778 The result from running metadetect.fitting.fit_mbobs_wavg
780 Returns
781 -------
782 ndarray
783 Has the fields from res, with new fields added, see get_output_dtype
784 """
785 dt = get_output_dtype()
786 output = eu.numpy_util.add_fields(res, dt)
788 for subdt in dt:
789 name = subdt[0]
790 dtype = subdt[1]
792 if 'flags' in name:
793 output[name] = NO_ATTEMPT
794 elif name in ('bmask', 'ormask'):
795 output[name] = 0
796 elif dtype[0] == 'i':
797 output[name] = -9999
798 else:
799 output[name] = np.nan
801 return output
804def get_output(wcs, source, res, bmask, stamp_size, exp_bbox):
805 """
806 get the output structure, copying in results
808 The following fields are not set:
809 psfrec_flags, psfrec_g, psfrec_T
810 mfrac
812 Parameters
813 ----------
814 wcs: a stack wcs
815 The wcs with which to determine the ra, dec
816 res: ndarray
817 The result from running metadetect.fitting.fit_mbobs_wavg
818 bmask: int
819 The bmask value at the location of this object
820 stamp_size: int
821 The stamp size used for the measurement
822 exp_bbox: lsst.geom.Box2I
823 The bounding box used for measurement
825 Returns
826 -------
827 ndarray
828 Has the fields from res, with new fields added, see get_output_dtype
829 """
830 output = get_output_struct(res)
832 orig_cen = source.getCentroid()
834 skypos = wcs.pixelToSky(orig_cen)
836 peak = source.getFootprint().getPeaks()[0]
837 peak_loc = peak.getI()
839 if np.isnan(orig_cen.getY()):
840 orig_cen = peak.getCentroid()
841 cen_offset = geom.Point2D(np.nan, np.nan)
842 else:
843 cen_offset = geom.Point2D(
844 orig_cen.getX() - peak_loc.getX(),
845 orig_cen.getY() - peak_loc.getY(),
846 )
848 output['stamp_size'] = stamp_size
849 output['row0'] = exp_bbox.getBeginY()
850 output['col0'] = exp_bbox.getBeginX()
851 output['row'] = orig_cen.getY()
852 output['col'] = orig_cen.getX()
853 output['row_diff'] = cen_offset.getY()
854 output['col_diff'] = cen_offset.getX()
856 output['ra'] = skypos.getRa().asDegrees()
857 output['dec'] = skypos.getDec().asDegrees()
859 # remove DETECTED bit, it is just clutter since all detected
860 # objects have this bit set
861 detected = afw_image.Mask.getPlaneBitMask('DETECTED')
862 output['bmask'] = bmask & ~detected
864 return output
867class MissingDataError(Exception):
868 """
869 Some number was out of range
870 """
872 def __init__(self, value):
873 super().__init__(value)
874 self.value = value
876 def __str__(self):
877 return repr(self.value)
880class AllZeroWeightError(Exception):
881 """
882 Some number was out of range
883 """
885 def __init__(self, value):
886 super().__init__(value)
887 self.value = value
889 def __str__(self):
890 return repr(self.value)