Coverage for python / lsst / atmospec / processStar.py: 28%

1# This file is part of atmospec. 

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/>. 

21 

22__all__ = ['ProcessStarTask', 'ProcessStarTaskConfig'] 

23 

24import importlib.resources 

25import shutil 

26import numpy as np 

27import matplotlib.pyplot as plt 

28 

29import lsstDebug 

30import lsst.afw.image as afwImage 

31import lsst.geom as geom 

32from lsst.ip.isr import IsrTaskLSST 

33import lsst.pex.config as pexConfig 

34from lsst.pex.config import FieldValidationError 

35import lsst.pipe.base as pipeBase 

36import lsst.pipe.base.connectionTypes as cT 

37from lsst.pipe.base.task import TaskError 

38 

39from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask 

40from lsst.meas.algorithms import ReferenceObjectLoader 

41from lsst.meas.astrom import AstrometryTask, FitAffineWcsTask 

42 

43import lsst.afw.detection as afwDetect 

44 

45from .spectraction import SpectractorShim 

46from .utils import getLinearStagePosition, isDispersedExp, getFilterAndDisperserFromExp 

47 

48COMMISSIONING = False # allows illegal things for on the mountain usage. 

49 

50# TODO: 

51# Sort out read noise and gain 

52# remove dummy image totally 

53# talk to Jeremy about turning the image beforehand and giving new coords 

54# deal with not having ambient temp 

55# Gen3ification 

56# astropy warning for units on save 

57# but actually just remove all manual saves entirely, I think? 

58# Make SED persistable 

59# Move to QFM for star finding failover case 

60# Remove old cruft functions 

61# change spectractions run method to be ~all kwargs with *,... 

62 

63 

64class ProcessStarTaskConnections(pipeBase.PipelineTaskConnections, 

65 dimensions=("instrument", "visit", "detector")): 

66 inputExp = cT.Input( 

67 name="icExp", 

68 doc="Image-characterize output exposure", 

69 storageClass="ExposureF", 

70 dimensions=("instrument", "visit", "detector"), 

71 multiple=False, 

72 ) 

73 inputCentroid = cT.Input( 

74 name="atmospecCentroid", 

75 doc="The main star centroid in yaml format.", 

76 storageClass="StructuredDataDict", 

77 dimensions=("instrument", "visit", "detector"), 

78 multiple=False, 

79 ) 

80 spectractorSpectrum = cT.Output( 

81 name="spectractorSpectrum", 

82 doc="The Spectractor output spectrum.", 

83 storageClass="SpectractorSpectrum", 

84 dimensions=("instrument", "visit", "detector"), 

85 ) 

86 spectractorImage = cT.Output( 

87 name="spectractorImage", 

88 doc="The Spectractor output image.", 

89 storageClass="SpectractorImage", 

90 dimensions=("instrument", "visit", "detector"), 

91 ) 

92 spectrumForwardModelFitParameters = cT.Output( 

93 name="spectrumForwardModelFitParameters", 

94 doc="The full forward model fit parameters.", 

95 storageClass="SpectractorFitParameters", 

96 dimensions=("instrument", "visit", "detector"), 

97 ) 

98 spectrumLibradtranFitParameters = cT.Output( 

99 name="spectrumLibradtranFitParameters", 

100 doc="The fitted Spectractor atmospheric parameters from fitting the atmosphere with libradtran" 

101 " on the spectrum.", 

102 storageClass="SpectractorFitParameters", 

103 dimensions=("instrument", "visit", "detector"), 

104 ) 

105 spectrogramLibradtranFitParameters = cT.Output( 

106 name="spectrogramLibradtranFitParameters", 

107 doc="The fitted Spectractor atmospheric parameters from fitting the atmosphere with libradtran" 

108 " directly on the spectrogram.", 

109 storageClass="SpectractorFitParameters", 

110 dimensions=("instrument", "visit", "detector"), 

111 ) 

112 

113 def __init__(self, *, config=None): 

114 super().__init__(config=config) 

115 if not config.doFullForwardModelDeconvolution: 

116 self.outputs.remove("spectrumForwardModelFitParameters") 

117 if not config.doFitAtmosphere: 

118 self.outputs.remove("spectrumLibradtranFitParameters") 

119 if not config.doFitAtmosphereOnSpectrogram: 

120 self.outputs.remove("spectrogramLibradtranFitParameters") 

121 

122 

123class ProcessStarTaskConfig(pipeBase.PipelineTaskConfig, 

124 pipelineConnections=ProcessStarTaskConnections): 

125 """Configuration parameters for ProcessStarTask.""" 

126 # Spectractor parameters: 

127 targetCentroidMethod = pexConfig.ChoiceField( 

128 dtype=str, 

129 doc="Method to get target centroid. " 

130 "SPECTRACTOR_FIT_TARGET_CENTROID internally.", 

131 default="auto", 

132 allowed={ 

133 # note that although this config option controls 

134 # SPECTRACTOR_FIT_TARGET_CENTROID, it doesn't map there directly, 

135 # because Spectractor only has the concepts of guess, fit and wcs, 

136 # and it calls "exact" "guess" internally, so that's remapped. 

137 "auto": "If the upstream astrometric fit succeeded, and therefore" 

138 " the centroid is an exact one, use that as an ``exact`` value," 

139 " otherwise tell Spectractor to ``fit`` the centroid", 

140 "exact": "Use a given input value as source of truth.", 

141 "fit": "Fit a 2d Moffat model to the target.", 

142 "WCS": "Use the target's catalog location and the image's wcs.", 

143 } 

144 ) 

145 rotationAngleMethod = pexConfig.ChoiceField( 

146 dtype=str, 

147 doc="Method used to get the image rotation angle. " 

148 "SPECTRACTOR_COMPUTE_ROTATION_ANGLE internally.", 

149 default="disperser", 

150 allowed={ 

151 # XXX MFL: probably need to use setDefaults to set this based on 

152 # the disperser. I think Ronchi gratings want hessian and the 

153 # holograms want disperser. 

154 "False": "Do not rotate the image.", 

155 "disperser": "Use the disperser angle geometry as specified in the disperser definition.", 

156 "hessian": "Compute the angle from the image using a Hessian transform.", 

157 } 

158 ) 

159 doDeconvolveSpectrum = pexConfig.Field( 

160 dtype=bool, 

161 doc="Deconvolve the spectrogram with a simple 2D PSF analysis? " 

162 "SPECTRACTOR_DECONVOLUTION_PSF2D internally.", 

163 default=True, 

164 ) 

165 doFullForwardModelDeconvolution = pexConfig.Field( 

166 dtype=bool, 

167 doc="Deconvolve the spectrogram with full forward model? " 

168 "SPECTRACTOR_DECONVOLUTION_FFM internally.", 

169 default=True, 

170 ) 

171 deconvolutionSigmaClip = pexConfig.Field( 

172 dtype=float, 

173 doc="Sigma clipping level for the deconvolution when fitting the full forward model? " 

174 "SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP internally.", 

175 default=100, 

176 ) 

177 doSubtractBackground = pexConfig.Field( 

178 dtype=bool, 

179 doc="Subtract the background with Spectractor? " 

180 "SPECTRACTOR_BACKGROUND_SUBTRACTION internally.", 

181 default=True, 

182 ) 

183 rebin = pexConfig.Field( 

184 dtype=int, 

185 doc="Rebinning factor to use on the input image, in pixels. " 

186 "CCD_REBIN internally.", 

187 default=2, # TODO Change to 1 once speed issues are resolved 

188 ) 

189 xWindow = pexConfig.Field( 

190 dtype=int, 

191 doc="Window x size to search for the target object. Ignored if targetCentroidMethod in ('exact, wcs')" 

192 "XWINDOW internally.", 

193 default=150, 

194 ) 

195 yWindow = pexConfig.Field( 

196 dtype=int, 

197 doc="Window y size to search for the targeted object. Ignored if targetCentroidMethod in " 

198 "('exact, wcs')" 

199 "YWINDOW internally.", 

200 default=150, 

201 ) 

202 xWindowRotated = pexConfig.Field( 

203 dtype=int, 

204 doc="Window x size to search for the target object in the rotated image. " 

205 "Ignored if rotationAngleMethod=False" 

206 "XWINDOW_ROT internally.", 

207 default=50, 

208 ) 

209 yWindowRotated = pexConfig.Field( 

210 dtype=int, 

211 doc="Window y size to search for the target object in the rotated image. " 

212 "Ignored if rotationAngleMethod=False" 

213 "YWINDOW_ROT internally.", 

214 default=50, 

215 ) 

216 pixelShiftPrior = pexConfig.Field( 

217 dtype=float, 

218 doc="Prior on the reliability of the centroid estimate in pixels. " 

219 "PIXSHIFT_PRIOR internally.", 

220 default=5, 

221 check=lambda x: x > 0, 

222 ) 

223 doFilterRotatedImage = pexConfig.Field( 

224 dtype=bool, 

225 doc="Apply a filter to the rotated image? If not True, this creates residuals and correlated noise. " 

226 "ROT_PREFILTER internally.", 

227 default=True, 

228 ) 

229 imageRotationSplineOrder = pexConfig.Field( 

230 dtype=int, 

231 doc="Order of the spline used when rotating the image. " 

232 "ROT_ORDER internally.", 

233 default=5, 

234 # XXX min value of 3 for allowed range, max 5 

235 ) 

236 rotationAngleMin = pexConfig.Field( 

237 dtype=float, 

238 doc="In the Hessian analysis to compute the rotation angle, cut all angles below this, in degrees. " 

239 "ROT_ANGLE_MIN internally.", 

240 default=-10, 

241 ) 

242 rotationAngleMax = pexConfig.Field( 

243 dtype=float, 

244 doc="In the Hessian analysis to compute rotation angle, cut all angles above this, in degrees. " 

245 "ROT_ANGLE_MAX internally.", 

246 default=10, 

247 ) 

248 plotLineWidth = pexConfig.Field( 

249 dtype=float, 

250 doc="Line width parameter for plotting. " 

251 "LINEWIDTH internally.", 

252 default=2, 

253 ) 

254 verbose = pexConfig.Field( 

255 dtype=bool, 

256 doc="Set verbose mode? " 

257 "VERBOSE internally.", 

258 default=True, # sets INFO level logging in Spectractor 

259 ) 

260 spectractorDebugMode = pexConfig.Field( 

261 dtype=bool, 

262 doc="Set spectractor debug mode? " 

263 "DEBUG internally.", 

264 default=True, 

265 ) 

266 spectractorDebugLogging = pexConfig.Field( 

267 dtype=bool, 

268 doc="Set spectractor debug logging? " 

269 "DEBUG_LOGGING internally.", 

270 default=False 

271 ) 

272 doDisplay = pexConfig.Field( 

273 dtype=bool, 

274 doc="Display plots, for example when running in a notebook? " 

275 "DISPLAY internally.", 

276 default=True 

277 ) 

278 lambdaMin = pexConfig.Field( 

279 dtype=int, 

280 doc="Minimum wavelength for spectral extraction (in nm). " 

281 "LAMBDA_MIN internally.", 

282 default=300 

283 ) 

284 lambdaMax = pexConfig.Field( 

285 dtype=int, 

286 doc=" maximum wavelength for spectrum extraction (in nm). " 

287 "LAMBDA_MAX internally.", 

288 default=1100 

289 ) 

290 lambdaStep = pexConfig.Field( 

291 dtype=float, 

292 doc="Step size for the wavelength array (in nm). " 

293 "LAMBDA_STEP internally.", 

294 default=1, 

295 ) 

296 spectralOrder = pexConfig.ChoiceField( 

297 dtype=int, 

298 doc="The spectral order to extract. " 

299 "SPEC_ORDER internally.", 

300 default=1, 

301 allowed={ 

302 1: "The first order spectrum in the positive y direction", 

303 -1: "The first order spectrum in the negative y direction", 

304 2: "The second order spectrum in the positive y direction", 

305 -2: "The second order spectrum in the negative y direction", 

306 } 

307 ) 

308 signalWidth = pexConfig.Field( # TODO: change this to be set wrt the focus/seeing, i.e. FWHM from imChar 

309 dtype=int, 

310 doc="Half transverse width of the signal rectangular window in pixels. " 

311 "PIXWIDTH_SIGNAL internally.", 

312 default=40, 

313 ) 

314 backgroundDistance = pexConfig.Field( 

315 dtype=int, 

316 doc="Distance from dispersion axis to analyse the background in pixels. " 

317 "PIXDIST_BACKGROUND internally.", 

318 default=140, 

319 ) 

320 backgroundWidth = pexConfig.Field( 

321 dtype=int, 

322 doc="Transverse width of the background rectangular window in pixels. " 

323 "PIXWIDTH_BACKGROUND internally.", 

324 default=40, 

325 ) 

326 backgroundBoxSize = pexConfig.Field( 

327 dtype=int, 

328 doc="Box size for sextractor evaluation of the background. " 

329 "PIXWIDTH_BOXSIZE internally.", 

330 default=20, 

331 ) 

332 backgroundOrder = pexConfig.Field( 

333 dtype=int, 

334 doc="The order of the polynomial background to fit in the transverse direction. " 

335 "BGD_ORDER internally.", 

336 default=1, 

337 ) 

338 psfType = pexConfig.ChoiceField( 

339 dtype=str, 

340 doc="The PSF model type to use. " 

341 "PSF_TYPE internally.", 

342 default="Moffat", 

343 allowed={ 

344 "Moffat": "A Moffat function", 

345 "MoffatGauss": "A Moffat plus a Gaussian" 

346 } 

347 ) 

348 psfPolynomialOrder = pexConfig.Field( 

349 dtype=int, 

350 doc="The order of the polynomials to model wavelength dependence of the PSF shape parameters. " 

351 "PSF_POLY_ORDER internally.", 

352 default=2 

353 ) 

354 psfRegularization = pexConfig.Field( 

355 dtype=float, 

356 doc="Regularisation parameter for the chisq minimisation to extract the spectrum. " 

357 "PSF_FIT_REG_PARAM internally.", 

358 default=1, 

359 # XXX allowed range strictly positive 

360 ) 

361 psfTransverseStepSize = pexConfig.Field( 

362 dtype=int, 

363 doc="Step size in pixels for the first transverse PSF1D fit. " 

364 "PSF_PIXEL_STEP_TRANSVERSE_FIT internally.", 

365 default=10, 

366 ) 

367 psfFwhmClip = pexConfig.Field( 

368 dtype=float, 

369 doc="PSF is not evaluated outside a region larger than max(signalWidth, psfFwhmClip*fwhm) pixels. " 

370 "PSF_FWHM_CLIP internally.", 

371 default=2, 

372 ) 

373 calibBackgroundOrder = pexConfig.Field( 

374 dtype=int, 

375 doc="Order of the background polynomial to fit. " 

376 "CALIB_BGD_ORDER internally.", 

377 default=3, 

378 ) 

379 calibPeakWidth = pexConfig.Field( 

380 dtype=int, 

381 doc="Half-range to look for local extrema in pixels around tabulated line values. " 

382 "CALIB_PEAK_WIDTH internally.", 

383 default=7 

384 ) 

385 calibBackgroundWidth = pexConfig.Field( 

386 dtype=int, 

387 doc="Size of the peak sides to use to fit spectrum base line. " 

388 "CALIB_BGD_WIDTH internally.", 

389 default=15, 

390 ) 

391 calibSavgolWindow = pexConfig.Field( 

392 dtype=int, 

393 doc="Window size for the savgol filter in pixels. " 

394 "CALIB_SAVGOL_WINDOW internally.", 

395 default=5, 

396 ) 

397 calibSavgolOrder = pexConfig.Field( 

398 dtype=int, 

399 doc="Polynomial order for the savgol filter. " 

400 "CALIB_SAVGOL_ORDER internally.", 

401 default=2, 

402 ) 

403 transmissionSystematicError = pexConfig.Field( 

404 dtype=float, 

405 doc="The systematic error on the instrumental transmission. OBS_TRANSMISSION_SYSTEMATICS internally", 

406 default=0.005 

407 ) 

408 instrumentTransmissionOverride = pexConfig.Field( 

409 dtype=str, 

410 doc="File to use for the full instrumental transmission. Must be located in the" 

411 " $SPECTRACTOR_DIR/spectractor/simulation/AuxTelThroughput/ directory." 

412 " OBS_FULL_INSTRUMENT_TRANSMISSON internally.", 

413 default="multispectra_holo4_003_HD142331_AuxTel_throughput.txt" 

414 ) 

415 offsetFromMainStar = pexConfig.Field( 

416 dtype=int, 

417 doc="Number of pixels from the main star's centroid to start extraction", 

418 default=100 

419 ) 

420 spectrumLengthPixels = pexConfig.Field( 

421 dtype=int, 

422 doc="Length of the spectrum in pixels", 

423 default=5000 

424 ) 

425 # ProcessStar own parameters 

426 isr = pexConfig.ConfigurableField( 

427 target=IsrTaskLSST, 

428 doc="Task to perform instrumental signature removal", 

429 ) 

430 charImage = pexConfig.ConfigurableField( 

431 target=CharacterizeImageTask, 

432 doc="""Task to characterize a science exposure: 

433 - detect sources, usually at high S/N 

434 - estimate the background, which is subtracted from the image and returned as field "background" 

435 - estimate a PSF model, which is added to the exposure 

436 - interpolate over defects and cosmic rays, updating the image, variance and mask planes 

437 """, 

438 ) 

439 doWrite = pexConfig.Field( 

440 dtype=bool, 

441 doc="Write out the results?", 

442 default=True, 

443 ) 

444 doFlat = pexConfig.Field( 

445 dtype=bool, 

446 doc="Flatfield the image?", 

447 default=True 

448 ) 

449 doCosmics = pexConfig.Field( 

450 dtype=bool, 

451 doc="Repair cosmic rays?", 

452 default=True 

453 ) 

454 doDisplayPlots = pexConfig.Field( 

455 dtype=bool, 

456 doc="Matplotlib show() the plots, so they show up in a notebook or X window", 

457 default=False 

458 ) 

459 doSavePlots = pexConfig.Field( 

460 dtype=bool, 

461 doc="Save matplotlib plots to output rerun?", 

462 default=False 

463 ) 

464 forceObjectName = pexConfig.Field( 

465 dtype=str, 

466 doc="A supplementary name for OBJECT. Will be forced to apply to ALL visits, so this should only" 

467 " ONLY be used for immediate commissioning debug purposes. All long term fixes should be" 

468 " supplied as header fix-up yaml files.", 

469 default="" 

470 ) 

471 referenceFilterOverride = pexConfig.Field( 

472 dtype=str, 

473 doc="Which filter in the reference catalog to match to?", 

474 default="phot_g_mean" 

475 ) 

476 # This is a post-processing function in Spectractor and therefore isn't 

477 # controlled by its top-level function, and thus doesn't map to a 

478 # spectractor.parameters ALL_CAPS config option 

479 doFitAtmosphere = pexConfig.Field( 

480 dtype=bool, 

481 doc="Use uvspec to fit the atmosphere? Requires the binary to be available.", 

482 default=False 

483 ) 

484 # This is a post-processing function in Spectractor and therefore isn't 

485 # controlled by its top-level function, and thus doesn't map to a 

486 # spectractor.parameters ALL_CAPS config option 

487 doFitAtmosphereOnSpectrogram = pexConfig.Field( 

488 dtype=bool, 

489 doc="Experimental option to use uvspec to fit the atmosphere directly on the spectrogram?" 

490 " Requires the binary to be available.", 

491 default=False 

492 ) 

493 

494 def setDefaults(self): 

495 self.charImage.doWriteExposure = False 

496 

497 self.charImage.doApCorr = False 

498 self.charImage.doMeasurePsf = False 

499 self.charImage.repair.cosmicray.nCrPixelMax = 100000 

500 self.charImage.repair.doCosmicRay = False 

501 if self.charImage.doMeasurePsf: 

502 self.charImage.measurePsf.starSelector['objectSize'].signalToNoiseMin = 10.0 

503 self.charImage.measurePsf.starSelector['objectSize'].fluxMin = 5000.0 

504 self.charImage.detection.includeThresholdMultiplier = 3 

505 

506 def validate(self): 

507 super().validate() 

508 uvspecPath = shutil.which('uvspec') 

509 if uvspecPath is None and self.doFitAtmosphere is True: 

510 raise FieldValidationError(self.__class__.doFitAtmosphere, self, "uvspec is not in the path," 

511 " but doFitAtmosphere is True.") 

512 if uvspecPath is None and self.doFitAtmosphereOnSpectrogram is True: 

513 raise FieldValidationError(self.__class__.doFitAtmosphereOnSpectrogram, self, "uvspec is not in" 

514 " the path, but doFitAtmosphere is True.") 

515 

516 

517class ProcessStarTask(pipeBase.PipelineTask): 

518 """Task for the spectral extraction of single-star dispersed images. 

519 

520 For a full description of how this tasks works, see the run() method. 

521 """ 

522 

523 ConfigClass = ProcessStarTaskConfig 

524 _DefaultName = "processStar" 

525 

526 def __init__(self, **kwargs): 

527 super().__init__(**kwargs) 

528 self.makeSubtask("isr") 

529 self.makeSubtask("charImage") 

530 

531 self.debug = lsstDebug.Info(__name__) 

532 if self.debug.enabled: 

533 self.log.info("Running with debug enabled...") 

534 # If we're displaying, test it works and save displays for later. 

535 # It's worth testing here as displays are flaky and sometimes 

536 # can't be contacted, and given processing takes a while, 

537 # it's a shame to fail late due to display issues. 

538 if self.debug.display: 

539 try: 

540 import lsst.afw.display as afwDisp 

541 afwDisp.setDefaultBackend(self.debug.displayBackend) 

542 afwDisp.Display.delAllDisplays() 

543 # pick an unlikely number to be safe xxx replace this 

544 self.disp1 = afwDisp.Display(987, open=True) 

545 

546 im = afwImage.ImageF(2, 2) 

547 im.array[:] = np.ones((2, 2)) 

548 self.disp1.mtv(im) 

549 self.disp1.erase() 

550 afwDisp.setDefaultMaskTransparency(90) 

551 except NameError: 

552 self.debug.display = False 

553 self.log.warning('Failed to setup/connect to display! Debug display has been disabled') 

554 

555 if self.debug.notHeadless: 

556 pass # other backend options can go here 

557 else: # this stop windows popping up when plotting. When headless, use 'agg' backend too 

558 plt.interactive(False) 

559 

560 self.config.validate() 

561 self.config.freeze() 

562 

563 def findObjects(self, exp, nSigma=None, grow=0): 

564 """Find the objects in a postISR exposure.""" 

565 nPixMin = self.config.mainStarNpixMin 

566 if not nSigma: 

567 nSigma = self.config.mainStarNsigma 

568 if not grow: 

569 grow = self.config.mainStarGrow 

570 isotropic = self.config.mainStarGrowIsotropic 

571 

572 threshold = afwDetect.Threshold(nSigma, afwDetect.Threshold.STDEV) 

573 footPrintSet = afwDetect.FootprintSet(exp.getMaskedImage(), threshold, "DETECTED", nPixMin) 

574 if grow > 0: 

575 footPrintSet = afwDetect.FootprintSet(footPrintSet, grow, isotropic) 

576 return footPrintSet 

577 

578 def _getEllipticity(self, shape): 

579 """Calculate the ellipticity given a quadrupole shape. 

580 

581 Parameters 

582 ---------- 

583 shape : `lsst.afw.geom.ellipses.Quadrupole` 

584 The quadrupole shape 

585 

586 Returns 

587 ------- 

588 ellipticity : `float` 

589 The magnitude of the ellipticity 

590 """ 

591 ixx = shape.getIxx() 

592 iyy = shape.getIyy() 

593 ixy = shape.getIxy() 

594 ePlus = (ixx - iyy) / (ixx + iyy) 

595 eCross = 2*ixy / (ixx + iyy) 

596 return (ePlus**2 + eCross**2)**0.5 

597 

598 def getRoundestObject(self, footPrintSet, parentExp, fluxCut=1e-15): 

599 """Get the roundest object brighter than fluxCut from a footPrintSet. 

600 

601 Parameters 

602 ---------- 

603 footPrintSet : `lsst.afw.detection.FootprintSet` 

604 The set of footprints resulting from running detection on parentExp 

605 

606 parentExp : `lsst.afw.image.exposure` 

607 The parent exposure for the footprint set. 

608 

609 fluxCut : `float` 

610 The flux, below which, sources are rejected. 

611 

612 Returns 

613 ------- 

614 source : `lsst.afw.detection.Footprint` 

615 The winning footprint from the input footPrintSet 

616 """ 

617 self.log.debug("ellipticity\tflux/1e6\tcentroid") 

618 sourceDict = {} 

619 for fp in footPrintSet.getFootprints(): 

620 shape = fp.getShape() 

621 e = self._getEllipticity(shape) 

622 flux = fp.getSpans().flatten(parentExp.image.array, parentExp.image.getXY0()).sum() 

623 self.log.debug("%.4f\t%.2f\t%s"%(e, flux/1e6, str(fp.getCentroid()))) 

624 if flux > fluxCut: 

625 sourceDict[e] = fp 

626 

627 return sourceDict[sorted(sourceDict.keys())[0]] 

628 

629 def getBrightestObject(self, footPrintSet, parentExp, roundnessCut=1e9): 

630 """Get the brightest object rounder than the cut from a footPrintSet. 

631 

632 Parameters 

633 ---------- 

634 footPrintSet : `lsst.afw.detection.FootprintSet` 

635 The set of footprints resulting from running detection on parentExp 

636 

637 parentExp : `lsst.afw.image.exposure` 

638 The parent exposure for the footprint set. 

639 

640 roundnessCut : `float` 

641 The ellipticity, above which, sources are rejected. 

642 

643 Returns 

644 ------- 

645 source : `lsst.afw.detection.Footprint` 

646 The winning footprint from the input footPrintSet 

647 """ 

648 self.log.debug("ellipticity\tflux\tcentroid") 

649 sourceDict = {} 

650 for fp in footPrintSet.getFootprints(): 

651 shape = fp.getShape() 

652 e = self._getEllipticity(shape) 

653 flux = fp.getSpans().flatten(parentExp.image.array, parentExp.image.getXY0()).sum() 

654 self.log.debug("%.4f\t%.2f\t%s"%(e, flux/1e6, str(fp.getCentroid()))) 

655 if e < roundnessCut: 

656 sourceDict[flux] = fp 

657 

658 return sourceDict[sorted(sourceDict.keys())[-1]] 

659 

660 def findMainSource(self, exp): 

661 """Return the x,y of the brightest or roundest object in an exposure. 

662 

663 Given a postISR exposure, run source detection on it, and return the 

664 centroid of the main star. Depending on the task configuration, this 

665 will either be the roundest object above a certain flux cutoff, or 

666 the brightest object which is rounder than some ellipticity cutoff. 

667 

668 Parameters 

669 ---------- 

670 exp : `afw.image.Exposure` 

671 The postISR exposure in which to find the main star 

672 

673 Returns 

674 ------- 

675 x, y : `tuple` of `float` 

676 The centroid of the main star in the image 

677 

678 Notes 

679 ----- 

680 Behavior of this method is controlled by many task config params 

681 including, for the detection stage: 

682 config.mainStarNpixMin 

683 config.mainStarNsigma 

684 config.mainStarGrow 

685 config.mainStarGrowIsotropic 

686 

687 And post-detection, for selecting the main source: 

688 config.mainSourceFindingMethod 

689 config.mainStarFluxCut 

690 config.mainStarRoundnessCut 

691 """ 

692 # TODO: probably replace all this with QFM 

693 fpSet = self.findObjects(exp) 

694 if self.config.mainSourceFindingMethod == 'ROUNDEST': 

695 source = self.getRoundestObject(fpSet, exp, fluxCut=self.config.mainStarFluxCut) 

696 elif self.config.mainSourceFindingMethod == 'BRIGHTEST': 

697 source = self.getBrightestObject(fpSet, exp, 

698 roundnessCut=self.config.mainStarRoundnessCut) 

699 else: 

700 # should be impossible as this is a choice field, but still 

701 raise RuntimeError("Invalid source finding method " 

702 f"selected: {self.config.mainSourceFindingMethod}") 

703 return source.getCentroid() 

704 

705 def updateMetadata(self, exp, **kwargs): 

706 """Update an exposure's metadata with set items from the visit info. 

707 

708 Spectractor expects many items, like the hour angle and airmass, to be 

709 in the metadata, so pull them out of the visit info etc and put them 

710 into the main metadata. Also updates the metadata with any supplied 

711 kwargs. 

712 

713 Parameters 

714 ---------- 

715 exp : `lsst.afw.image.Exposure` 

716 The exposure to update. 

717 **kwargs : `dict` 

718 The items to add. 

719 """ 

720 md = exp.getMetadata() 

721 vi = exp.getInfo().getVisitInfo() 

722 

723 ha = vi.getBoresightHourAngle().asDegrees() 

724 airmass = vi.getBoresightAirmass() 

725 

726 md['HA'] = ha 

727 md.setComment('HA', 'Hour angle of observation start') 

728 

729 md['AIRMASS'] = airmass 

730 md.setComment('AIRMASS', 'Airmass at observation start') 

731 

732 if 'centroid' in kwargs: 

733 centroid = kwargs['centroid'] 

734 else: 

735 centroid = (None, None) 

736 

737 md['OBJECTX'] = centroid[0] 

738 md.setComment('OBJECTX', 'x pixel coordinate of object centroid') 

739 

740 md['OBJECTY'] = centroid[1] 

741 md.setComment('OBJECTY', 'y pixel coordinate of object centroid') 

742 

743 exp.setMetadata(md) 

744 

745 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

746 inputs = butlerQC.get(inputRefs) 

747 inputs['dataIdDict'] = inputRefs.inputExp.dataId 

748 

749 outputs = self.run(**inputs) 

750 butlerQC.put(outputs, outputRefs) 

751 

752 def getNormalizedTargetName(self, target): 

753 """Normalize the name of the target. 

754 

755 All targets which start with 'spec:' are converted to the name of the 

756 star without the leading 'spec:'. Any objects with mappings defined in 

757 data/nameMappings.txt are converted to the mapped name. 

758 

759 Parameters 

760 ---------- 

761 target : `str` 

762 The name of the target. 

763 

764 Returns 

765 ------- 

766 normalizedTarget : `str` 

767 The normalized name of the target. 

768 """ 

769 target = target.replace('spec:', '') 

770 

771 with importlib.resources.path("lsst.atmospec", "resources/data/nameMappings.txt") as nameMappingsFile: 

772 names, mappedNames = np.loadtxt(nameMappingsFile, dtype=str, unpack=True) 

773 assert len(names) == len(mappedNames) 

774 conversions = {name: mapped for name, mapped in zip(names, mappedNames)} 

775 

776 if target in conversions.keys(): 

777 converted = conversions[target] 

778 self.log.info(f"Converted target name {target} to {converted}") 

779 return converted 

780 return target 

781 

782 def _getSpectractorTargetSetting(self, inputCentroid): 

783 """Calculate the value to set SPECTRACTOR_FIT_TARGET_CENTROID to. 

784 

785 Parameters 

786 ---------- 

787 inputCentroid : `dict` 

788 The `atmospecCentroid` dict, as received in the task input data. 

789 

790 Returns 

791 ------- 

792 centroidMethod : `str` 

793 The value to set SPECTRACTOR_FIT_TARGET_CENTROID to. 

794 """ 

795 

796 # if mode is auto and the astrometry worked then it's an exact 

797 # centroid, and otherwise we fit, as per docs on this option. 

798 if self.config.targetCentroidMethod == 'auto': 

799 if inputCentroid['astrometricMatch'] is True: 

800 self.log.info("Auto centroid is using exact centroid for target from the astrometry") 

801 return 'guess' # this means exact 

802 else: 

803 self.log.info("Auto centroid is using FIT in Spectractor to get the target centroid") 

804 return 'fit' # this means exact 

805 

806 # this is just renaming the config parameter because guess sounds like 

807 # an instruction, and really we're saying to take this as given. 

808 if self.config.targetCentroidMethod == 'exact': 

809 return 'guess' 

810 

811 # all other options fall through 

812 return self.config.targetCentroidMethod 

813 

814 def run(self, *, inputExp, inputCentroid, dataIdDict): 

815 if not isDispersedExp(inputExp): 

816 raise RuntimeError(f"Exposure is not a dispersed image {dataIdDict}") 

817 starNames = self.loadStarNames() 

818 

819 overrideDict = { 

820 # normal config parameters 

821 'SPECTRACTOR_FIT_TARGET_CENTROID': self._getSpectractorTargetSetting(inputCentroid), 

822 'SPECTRACTOR_COMPUTE_ROTATION_ANGLE': self.config.rotationAngleMethod, 

823 'SPECTRACTOR_DECONVOLUTION_PSF2D': self.config.doDeconvolveSpectrum, 

824 'SPECTRACTOR_DECONVOLUTION_FFM': self.config.doFullForwardModelDeconvolution, 

825 'SPECTRACTOR_DECONVOLUTION_SIGMA_CLIP': self.config.deconvolutionSigmaClip, 

826 'SPECTRACTOR_BACKGROUND_SUBTRACTION': self.config.doSubtractBackground, 

827 'CCD_REBIN': self.config.rebin, 

828 'XWINDOW': self.config.xWindow, 

829 'YWINDOW': self.config.yWindow, 

830 'XWINDOW_ROT': self.config.xWindowRotated, 

831 'YWINDOW_ROT': self.config.yWindowRotated, 

832 'PIXSHIFT_PRIOR': self.config.pixelShiftPrior, 

833 'ROT_PREFILTER': self.config.doFilterRotatedImage, 

834 'ROT_ORDER': self.config.imageRotationSplineOrder, 

835 'ROT_ANGLE_MIN': self.config.rotationAngleMin, 

836 'ROT_ANGLE_MAX': self.config.rotationAngleMax, 

837 'LINEWIDTH': self.config.plotLineWidth, 

838 'VERBOSE': self.config.verbose, 

839 'DEBUG': self.config.spectractorDebugMode, 

840 'DEBUG_LOGGING': self.config.spectractorDebugLogging, 

841 'DISPLAY': self.config.doDisplay, 

842 'LAMBDA_MIN': self.config.lambdaMin, 

843 'LAMBDA_MAX': self.config.lambdaMax, 

844 'LAMBDA_STEP': self.config.lambdaStep, 

845 'SPEC_ORDER': self.config.spectralOrder, 

846 'PIXWIDTH_SIGNAL': self.config.signalWidth, 

847 'PIXDIST_BACKGROUND': self.config.backgroundDistance, 

848 'PIXWIDTH_BACKGROUND': self.config.backgroundWidth, 

849 'PIXWIDTH_BOXSIZE': self.config.backgroundBoxSize, 

850 'BGD_ORDER': self.config.backgroundOrder, 

851 'PSF_TYPE': self.config.psfType, 

852 'PSF_POLY_ORDER': self.config.psfPolynomialOrder, 

853 'PSF_FIT_REG_PARAM': self.config.psfRegularization, 

854 'PSF_PIXEL_STEP_TRANSVERSE_FIT': self.config.psfTransverseStepSize, 

855 'PSF_FWHM_CLIP': self.config.psfFwhmClip, 

856 'CALIB_BGD_ORDER': self.config.calibBackgroundOrder, 

857 'CALIB_PEAK_WIDTH': self.config.calibPeakWidth, 

858 'CALIB_BGD_WIDTH': self.config.calibBackgroundWidth, 

859 'CALIB_SAVGOL_WINDOW': self.config.calibSavgolWindow, 

860 'CALIB_SAVGOL_ORDER': self.config.calibSavgolOrder, 

861 'OBS_TRANSMISSION_SYSTEMATICS': self.config.transmissionSystematicError, 

862 'OBS_FULL_INSTRUMENT_TRANSMISSON': self.config.instrumentTransmissionOverride, 

863 

864 # Hard-coded parameters 

865 'OBS_NAME': 'AUXTEL', 

866 'CCD_IMSIZE': 4000, # short axis - we trim the CCD to square 

867 'CCD_MAXADU': 170000, # XXX need to set this from camera value 

868 'CCD_GAIN': 1.1, # set programatically later, this is default nominal value 

869 'OBS_NAME': 'AUXTEL', 

870 'OBS_ALTITUDE': 2.66299616375123, # XXX get this from / check with utils value 

871 'OBS_LATITUDE': -30.2446389756252, # XXX get this from / check with utils value 

872 'OBS_EPOCH': "J2000.0", 

873 'OBS_CAMERA_DEC_FLIP_SIGN': 1, 

874 'OBS_CAMERA_RA_FLIP_SIGN': 1, 

875 'OBS_SURFACE': 9636, 

876 'PAPER': False, 

877 'SAVE': False, 

878 'DISTANCE2CCD_ERR': 0.4, 

879 

880 # Parameters set programatically 

881 'LAMBDAS': np.arange(self.config.lambdaMin, 

882 self.config.lambdaMax, 

883 self.config.lambdaStep), 

884 'CALIB_BGD_NPARAMS': self.config.calibBackgroundOrder + 1, 

885 

886 # Parameters set elsewhere 

887 # OBS_CAMERA_ROTATION 

888 # DISTANCE2CCD 

889 } 

890 

891 supplementDict = {'CALLING_CODE': 'LSST_DM', 

892 'STAR_NAMES': starNames} 

893 

894 # anything that changes between dataRefs! 

895 resetParameters = {}