Coverage for metadetect / lsst / measure.py: 17%

310 statements  

« prev     ^ index     » next       coverage.py v7.13.4, created at 2026-03-23 22:27 -0700

1import logging 

2import warnings 

3import numpy as np 

4import esutil as eu 

5import ngmix 

6 

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) 

25 

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) 

36 

37from . import util 

38from .util import ContextNoiseReplacer, get_jacobian 

39from . import vis 

40from .defaults import DEFAULT_THRESH 

41 

42warnings.filterwarnings('ignore', category=FutureWarning) 

43 

44LOG = logging.getLogger('lsst_measure') 

45 

46 

47class SourceDetectionConfig(OriginalSourceDetectionConfig): 

48 @property 

49 def thresh(self): 

50 return self.thresholdValue 

51 

52 @thresh.setter 

53 def thresh(self, value): 

54 self.thresholdValue = value 

55 

56 

57SourceDetectionTask.ConfigClass = SourceDetectionConfig 

58 

59 

60class DetectAndDeblendConfig(Config): 

61 meas = ConfigurableField[SingleFrameMeasurementConfig]( 

62 doc="Measurement config", 

63 target=SingleFrameMeasurementTask, 

64 ) 

65 

66 detect = ConfigurableField[SourceDetectionConfig]( 

67 doc="Detection config", target=SourceDetectionTask 

68 ) 

69 

70 deblend = ConfigurableField[SourceDeblendConfig]( 

71 doc="Deblend config", target=SourceDeblendTask 

72 ) 

73 

74 seed = Field[int]( 

75 doc="Random rng seed", 

76 default=42, 

77 ) 

78 

79 def setDefaults(self): 

80 super().setDefaults() 

81 

82 # defaults for measurement config 

83 self.meas.plugins.names = [ 

84 "base_SdssCentroid", 

85 "base_PsfFlux", 

86 "base_SkyCoord", 

87 ] 

88 

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 

94 

95 # goes with SdssShape above 

96 self.meas.slots.shape = None 

97 

98 # fix odd issue where it things things are near the edge 

99 self.meas.plugins['base_SdssCentroid'].binmax = 1 

100 

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() 

122 

123 # the defaults changed from from stdev to pixel_std but 

124 # we don't want that 

125 

126 self.detect.thresholdType = "stdev" 

127 # our changes from defaults 

128 self.detect.reEstimateBackground = False 

129 

130 self.detect.thresholdValue = DEFAULT_THRESH 

131 

132 # these will be ignored when finding the image standard deviation 

133 self.detect.statsMask = util.get_stats_mask() 

134 

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 

141 

142 

143class DetectAndDeblendTask(Task): 

144 ConfigClass = DetectAndDeblendConfig 

145 _DefaultName = "detect_and_deblend" 

146 

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) 

154 

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] 

162 

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) 

166 

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) 

170 

171 if result is not None: 

172 sources = result.sources 

173 self.deblend.run(detexp, sources) 

174 

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 

184 

185 source_id = source.getId() 

186 

187 with replacer.sourceInserted(source_id): 

188 self.meas.callMeasure(source, detexp) 

189 

190 else: 

191 sources = [] 

192 

193 if show: 

194 vis.show_exp(detexp, use_mpl=True, sources=sources) 

195 

196 return sources, detexp 

197 

198 

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. 

209 

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 

213 

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 

224 

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 

235 

236 config = DetectAndDeblendConfig() 

237 config.setDefaults() 

238 

239 util.override_config(config, config_override) 

240 

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) 

247 

248 

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. 

259 

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 

273 

274 Returns 

275 ------- 

276 ndarray with results or None 

277 """ 

278 

279 if len(sources) == 0: 

280 return None 

281 

282 pgauss_fitter = get_pgauss_fitter(config) 

283 

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 = [] 

298 

299 # bmasks will be different within the loop below due to the replacer 

300 bmasks = get_bmasks(sources=sources, exposure=detexp) 

301 

302 for i, source in enumerate(sources): 

303 if source.get('deblend_nChild') != 0: 

304 continue 

305 

306 bmask = bmasks[i] 

307 

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 

324 

325 if stamp_flags != 0: 

326 

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 ) 

336 

337 else: 

338 

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 ) 

351 

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 ) 

359 

360 this_pgauss_res = fit_mbobs_wavg( 

361 mbobs=mbobs, 

362 fitter=pgauss_fitter, 

363 bmask_flags=0, 

364 ) 

365 

366 this_res = _get_combined_struct(this_gauss_res, this_pgauss_res) 

367 this_res['stamp_flags'] = stamp_flags 

368 

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 ) 

377 

378 results.append(res) 

379 

380 if len(results) > 0: 

381 results = eu.numpy_util.combine_arrlist(results) 

382 else: 

383 results = None 

384 

385 return results 

386 

387 

388def get_pgauss_fitter(config): 

389 return ngmix.prepsfmom.PGaussMom(fwhm=config['pgauss']['fwhm']) 

390 

391 

392def _get_combined_struct(gauss_res, pgauss_res): 

393 

394 skip = ['pgauss_g', 'pgauss_g_cov', 'shear_bands'] 

395 keep_dt = [('stamp_flags', 'i4')] 

396 

397 for pdt in pgauss_res.dtype.descr: 

398 

399 n = pdt[0] 

400 

401 if 'psf' in n or n in skip: 

402 continue 

403 

404 keep_dt.append(pdt) 

405 

406 out = eu.numpy_util.add_fields(gauss_res, keep_dt) 

407 

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] 

411 

412 return out 

413 

414 

415def get_bmasks(sources, exposure): 

416 """ 

417 get a list of all the bmasks for the sources 

418 

419 Parameters 

420 ---------- 

421 sources: lsst.afw.table.SourceCatalog 

422 The sources 

423 exposure: lsst.afw.image.ExposureF 

424 The exposure 

425 

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 

435 

436 

437def get_bmask(source, exposure): 

438 """ 

439 get bmask based on original peak position 

440 

441 Parameters 

442 ---------- 

443 sources: lsst.afw.table.SourceRecord 

444 The sources 

445 exposure: lsst.afw.image.ExposureF 

446 The exposure 

447 

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 

456 

457 

458def extract_obs(exp, source): 

459 """ 

460 convert an image object into an ngmix.Observation, including 

461 a psf observation 

462 

463 parameters 

464 ---------- 

465 imobj: lsst.afw.image.ExposureF 

466 The exposure 

467 source: lsst.afw.table.SourceRecord 

468 The source record 

469 

470 returns 

471 -------- 

472 obs: ngmix.Observation 

473 The Observation unless all the weight are zero, in which 

474 case AllZeroWeightError is raised 

475 """ 

476 

477 im = exp.image.array 

478 

479 wt = _extract_weight(exp) 

480 if np.all(wt <= 0): 

481 raise AllZeroWeightError('all weights <= 0') 

482 

483 bmask = exp.mask.array 

484 jacob = _extract_jacobian_at_source( 

485 exp=exp, 

486 source=source, 

487 ) 

488 

489 orig_cen = source.getCentroid() 

490 

491 psf_im = extract_psf_image(exposure=exp, orig_cen=orig_cen) 

492 

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 

496 

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]) 

501 

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 

506 

507 meta = {'orig_cen': orig_cen} 

508 

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 ) 

522 

523 return obs 

524 

525 

526def _get_stamp_mbobs(mbexp, source, stamp_size, clip=False): 

527 """ 

528 Get a postage stamp MultibandExposure 

529 

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. 

542 

543 If clip is False and the bbox does not fit, a 

544 lsst.pex.exceptions.LengthError is raised 

545 

546 Only relevant if stamp_size is sent. Default False 

547 

548 Returns 

549 ------- 

550 lsst.afw.image.ExposureF 

551 """ 

552 

553 bbox = _get_bbox(mbexp, source, stamp_size, clip=clip) 

554 

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 ) 

562 

563 obslist = ngmix.ObsList(meta={'band': band}) 

564 obslist.append(obs) 

565 mbobs.append(obslist) 

566 

567 return mbobs 

568 

569 

570def _get_bbox(mbexp, source, stamp_size, clip=False): 

571 """ 

572 Get a bounding box at the location of the specified source. 

573 

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. 

586 

587 If clip is False and the bbox does not fit, a 

588 lsst.pex.exceptions.LengthError is raised 

589 

590 Only relevant if stamp_size is sent. Default False 

591 

592 Returns 

593 ------- 

594 lsst.geom.Box2I 

595 """ 

596 

597 fp = source.getFootprint() 

598 peak = fp.getPeaks()[0] 

599 

600 x_peak, y_peak = peak.getIx(), peak.getIy() 

601 

602 bbox = geom.Box2I( 

603 geom.Point2I(x_peak, y_peak), 

604 geom.Extent2I(1, 1), 

605 ) 

606 bbox.grow(stamp_size // 2) 

607 

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 ) 

619 

620 return bbox 

621 

622 

623def extract_psf_image(exposure, orig_cen): 

624 """ 

625 get the psf associated with this image. 

626 

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 

629 

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 

636 

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") 

646 

647 psfim = np.asarray(psfim, dtype='f4') 

648 

649 shape = psfim.shape 

650 assert shape[0] == shape[1], 'require square psf images' 

651 assert shape[0] % 2 != 0, 'require odd psf images' 

652 

653 return psfim 

654 

655 

656def _extract_psf_image_fix(exposure, orig_cen): 

657 """ 

658 get the psf associated with this image 

659 

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 

663 

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") 

672 

673 psfim = np.asarray(psfim, dtype='f4') 

674 

675 psfim = util.trim_odd_image(psfim) 

676 return psfim 

677 

678 

679def _extract_weight(exp): 

680 """ 

681 TODO get the estimated sky variance rather than this hack 

682 TODO should we zero out other bits? 

683 

684 extract a weight map 

685 

686 Areas with NO_DATA will get zero weight. 

687 

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 

691 

692 parameters 

693 ---------- 

694 exp: sub exposure object 

695 """ 

696 

697 # TODO implement bit checking 

698 var_image = exp.variance.array 

699 

700 weight = var_image.copy() 

701 

702 weight[:, :] = 0 

703 

704 wuse = np.where(var_image > 0) 

705 

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') 

711 

712 return weight 

713 

714 

715def _extract_jacobian_at_source(exp, source): 

716 """ 

717 extract an ngmix.Jacobian from the image object 

718 and object record 

719 

720 exp: lsst.afw.image.ExposureF 

721 An exposure object 

722 source: lsst.afw.table.SourceRecord 

723 The source record created during detection 

724 

725 returns 

726 -------- 

727 Jacobian: ngmix.Jacobian 

728 The local jacobian 

729 """ 

730 

731 orig_cen = exp.getWcs().skyToPixel(source.getCoord()) 

732 

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 ) 

742 

743 return get_jacobian(exp, orig_cen) 

744 

745 

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 ] 

767 

768 return dt 

769 

770 

771def get_output_struct(res): 

772 """ 

773 get the output struct 

774 

775 Parameters 

776 ---------- 

777 res: ndarray 

778 The result from running metadetect.fitting.fit_mbobs_wavg 

779 

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) 

787 

788 for subdt in dt: 

789 name = subdt[0] 

790 dtype = subdt[1] 

791 

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 

800 

801 return output 

802 

803 

804def get_output(wcs, source, res, bmask, stamp_size, exp_bbox): 

805 """ 

806 get the output structure, copying in results 

807 

808 The following fields are not set: 

809 psfrec_flags, psfrec_g, psfrec_T 

810 mfrac 

811 

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 

824 

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) 

831 

832 orig_cen = source.getCentroid() 

833 

834 skypos = wcs.pixelToSky(orig_cen) 

835 

836 peak = source.getFootprint().getPeaks()[0] 

837 peak_loc = peak.getI() 

838 

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 ) 

847 

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() 

855 

856 output['ra'] = skypos.getRa().asDegrees() 

857 output['dec'] = skypos.getDec().asDegrees() 

858 

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 

863 

864 return output 

865 

866 

867class MissingDataError(Exception): 

868 """ 

869 Some number was out of range 

870 """ 

871 

872 def __init__(self, value): 

873 super().__init__(value) 

874 self.value = value 

875 

876 def __str__(self): 

877 return repr(self.value) 

878 

879 

880class AllZeroWeightError(Exception): 

881 """ 

882 Some number was out of range 

883 """ 

884 

885 def __init__(self, value): 

886 super().__init__(value) 

887 self.value = value 

888 

889 def __str__(self): 

890 return repr(self.value)