Coverage for python/lsst/ip/diffim/utils.py: 82%

1# This file is part of ip_diffim.

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.

9# This program is free software: you can redistribute it and/or modify

10# it under the terms of the GNU General Public License as published by

11# the Free Software Foundation, either version 3 of the License, or

12# (at your option) any later version.

13#

14# This program is distributed in the hope that it will be useful,

15# but WITHOUT ANY WARRANTY; without even the implied warranty of

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

17# GNU General Public License for more details.

18#

19# You should have received a copy of the GNU General Public License

20# along with this program. If not, see <https://www.gnu.org/licenses/>.

22"""Support utilities for Measuring sources"""

25__all__ = ["evaluateMeanPsfFwhm", "getPsfFwhm", "getKernelCenterDisplacement",

26 "computeDifferenceImageMetrics", "checkMask", "setSourceFootprints",

27 "record_from_obs_info", "get_difference_image_provenance",

28 ]

30import itertools

31import os

32import requests

34from astropy.stats import gaussian_sigma_to_fwhm

35import astropy.units as u

36import numpy as np

38import lsst.afw.detection as afwDetection

39import lsst.afw.image as afwImage

40import lsst.afw.math as afwMath

41from lsst.daf.butler import Timespan

42import lsst.geom as geom

43from lsst.pex.exceptions import InvalidParameterError, RangeError

44import lsst.pipe.base

45from lsst.utils.logging import getLogger

48_LOG = getLogger(__name__)

51def record_from_obs_info(obs_info, instrument, visit, universe):

52 """Build an ``exposure`` dimension record from an observation info.

54 Parameters

55 ----------

56 obs_info : `astro_metadata_translator.ObservationInfo`

57 Observation info of an image that the image to be converted is

58 derived from.

59 instrument : `str`

60 Name of the instrument, usually from the quantum data ID.

61 visit : `int`

62 Id of the visit, usually from the quantum data ID. It is also used

63 as the id of the exposure.

64 universe : `lsst.daf.butler.DimensionUniverse`

65 Dimension universe that defines the record's schema, usually

66 ``butlerQC.quantum.dataId.universe``.

68 Returns

69 -------

70 record : `lsst.daf.butler.DimensionRecord`

71 An ``exposure`` record holding the fields that `lsst.images` reads

72 when it converts a legacy image. Fields it does not read, such as

73 ``dark_time`` and ``target_name``, are `None`.

75 Raises

76 ------

77 ValueError

78 Raised if ``obs_info`` is `None`, which means the image it came from

79 was not in an `lsst.images` format.

81 Notes

82 -----

83 The instrument and the visit are passed separately because an

84 observation info does not carry them in the form the conversion needs:

85 its ``instrument`` is read from the legacy FITS header rather than from

86 the record, and its ``exposure_id`` is the id that the visit info held.

87 """

88 if obs_info is None:

89 raise ValueError("Cannot rebuild an exposure record without an observation info.")

90 azimuth = zenithAngle = None

91 if obs_info.altaz_begin is not None: 91 ↛ 94line 91 didn't jump to line 94 because the condition on line 91 was always true

92 azimuth = obs_info.altaz_begin.az.deg

93 zenithAngle = 90.0 - obs_info.altaz_begin.alt.deg

94 exposureTime = obs_info.exposure_time_requested

95 return universe["exposure"].RecordClass(

96 instrument=instrument,

97 id=visit,

98 obs_id=obs_info.observation_id,

99 group=obs_info.exposure_group,

100 day_obs=obs_info.observing_day,

101 physical_filter=obs_info.physical_filter,

102 exposure_time=None if exposureTime is None else exposureTime.to_value(u.s),

103 seq_num=obs_info.observation_counter,

104 seq_start=obs_info.group_counter_start,

105 seq_end=obs_info.group_counter_end,

106 can_see_sky=obs_info.can_see_sky,

107 azimuth=azimuth,

108 zenith_angle=zenithAngle,

109 timespan=Timespan(begin=obs_info.datetime_begin, end=obs_info.datetime_end),

110 )

111

112

113def get_difference_image_provenance(image, attribute, connection, dataId):

114 """Get the PSF matching kernel or the template records of an input

115 difference image.

116

117 Parameters

118 ----------

119 image : `lsst.images.DifferenceImage`

120 Input image read by a task in future mode.

121 attribute : {"kernel", "templates"}

122 Name of the attribute to get.

123 connection : `str`

124 Name of the input connection the image was read from.

125 dataId : `lsst.daf.butler.DataCoordinate`

126 Data ID of the quantum, for the error message.

127

128 Returns

129 -------

130 value : `lsst.images.convolution_kernels.ConvolutionKernel` or \

131 `list` [`lsst.images.DifferenceImageTemplateInfo`]

132 Value of the attribute.

133

134 Raises

135 ------

136 lsst.pipe.base.InvalidQuantumError

137 Raised if the image does not have the attribute, which happens when

138 the input was written as a legacy Exposure and converted on read.

139 """

140 try:

141 return getattr(image, attribute)

142 except AttributeError as e:

143 what = "PSF matching kernel" if attribute == "kernel" else "template coadd records"

144 raise lsst.pipe.base.InvalidQuantumError(

145 f"Input {connection!r} for {dataId} has no {what}, which is required for image_type='future'."

146 " This is likely due to reading a legacy Exposure, which lacks this component. To fix,"

147 " reprocess this image starting from `getTemplate`."

148 ) from e

149

150

151def getKernelCenterDisplacement(kernel, x, y, image=None):

152 """Calculate the PSF matching kernel peak offset from the nominal

153 position.

154

155 Parameters

156 ----------

157 kernel : `~lsst.afw.math.LinearCombinationKernel`

158 The PSF matching kernel to evaluate.

159 x : `float`

160 The x position on the detector to evaluate the kernel

161 y : `float`

162 The y position on the detector to evaluate the kernel

163 image : `~lsst.afw.image.ImageD`

164 The image to use as base for computing kernel pixel values

165

166 Returns

167 -------

168 kernel_sum : `float`

169 The sum of the kernel on the desired location

170 dx : `float`

171 The displacement of the kernel averaged peak, with respect to the

172 center of the extraction of the kernel

173 dy : `float`

174 The displacement of the kernel averaged peak, with respect to the

175 center of the extraction of the kernel

176 pos_angle: `float`

177 The position angle in detector coordinates of the displacement

178 length : `float`

179 The displacement module of the kernel centroid in pixel units

180 """

181

182 if image is None: 182 ↛ 186line 182 didn't jump to line 186 because the condition on line 182 was always true

183 image = afwImage.ImageD(kernel.getDimensions())

184

185 # obtain the kernel image

186 hsize = kernel.getWidth()//2

187 kernel_sum = kernel.computeImage(image, doNormalize=False, x=x, y=y)

188

189 data = image.array

190 h, w = data.shape

191 xx = np.arange(w)

192 yy = np.arange(h)

193

194 # create sum vectors and estimate weighted average

195 vx = data.sum(axis=0)

196 vx /= vx.sum()

197 dx = np.dot(vx, xx) - hsize

198

199 vy = data.sum(axis=1)

200 vy /= vy.sum()

201 dy = np.dot(vy, yy) - hsize

202

203 # obtain position angle and norm of displacement

204 pos_angle = np.arctan2(dy, dx)

205 length = np.sqrt(dx**2 + dy**2)

206

207 return kernel_sum, dx, dy, pos_angle, length

208

209

210def getPsfFwhm(psf, average=True, position=None):

211 """Directly calculate the horizontal and vertical widths

212 of a PSF at half its maximum value.

213

214 Parameters

215 ----------

216 psf : `~lsst.afw.detection.Psf`

217 Point spread function (PSF) to evaluate.

218 average : `bool`, optional

219 Set to return the average width over Y and X axes.

220 position : `~lsst.geom.Point2D`, optional

221 The position at which to evaluate the PSF. If `None`, then the

222 average position is used.

223

224 Returns

225 -------

226 psfSize : `float` | `tuple` [`float`]

227 The FWHM of the PSF computed at its average position.

228 Returns the widths along the Y and X axes,

229 or the average of the two if `average` is set.

230

231 See Also

232 --------

233 evaluateMeanPsfFwhm

234 """

235 if position is None:

236 position = psf.getAveragePosition()

237 shape = psf.computeShape(position)

238

239 if average:

240 return gaussian_sigma_to_fwhm*shape.getTraceRadius()

241 else:

242 return [gaussian_sigma_to_fwhm*np.sqrt(shape.getIxx()),

243 gaussian_sigma_to_fwhm*np.sqrt(shape.getIyy())]

244

245

246def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,

247 fwhmExposureBuffer: float, fwhmExposureGrid: int) -> float:

248 """Get the mean PSF FWHM by evaluating it on a grid within an exposure.

249

250 Parameters

251 ----------

252 exposure : `~lsst.afw.image.Exposure`

253 The exposure for which the mean FWHM of the PSF is to be computed.

254 The exposure must contain a `psf` attribute.

255 fwhmExposureBuffer : `float`

256 Fractional buffer margin to be left out of all sides of the image

257 during the construction of the grid to compute mean PSF FWHM in an

258 exposure.

259 fwhmExposureGrid : `int`

260 Grid size to compute the average PSF FWHM in an exposure (pixels).

261

262 Returns

263 -------

264 meanFwhm : `float`

265 The mean PSF FWHM on the exposure.

266

267 Raises

268 ------

269 ValueError

270 Raised if the PSF cannot be computed at any of the grid points.

271

272 See Also

273 --------

274 `getPsfFwhm`

275 `computeAveragePsf`

276 """

277

278 psf = exposure.psf

279

280 bbox = exposure.getBBox()

281 xmax, ymax = bbox.getMax()

282 xmin, ymin = bbox.getMin()

283

284 xbuffer = fwhmExposureBuffer*(xmax-xmin)

285 ybuffer = fwhmExposureBuffer*(ymax-ymin)

286

287 width = []

288 for (x, y) in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, fwhmExposureGrid),

289 np.linspace(ymin+ybuffer, ymax-ybuffer, fwhmExposureGrid)

290 ):

291 pos = geom.Point2D(x, y)

292 try:

293 fwhm = getPsfFwhm(psf, average=True, position=pos)

294 except (InvalidParameterError, RangeError):

295 _LOG.debug("Unable to compute PSF FWHM at position (%f, %f).", x, y)

296 continue

297

298 width.append(fwhm)

299

300 if not width: 300 ↛ 301line 300 didn't jump to line 301 because the condition on line 300 was never true

301 raise ValueError("Unable to compute PSF FWHM at any position on the exposure.")

302

303 return np.nanmean(width)

304

305

306def computeAveragePsf(exposure: afwImage.Exposure,

307 psfExposureBuffer: float, psfExposureGrid: int) -> afwImage.ImageD:

308 """Get the average PSF by evaluating it on a grid within an exposure.

309

310 Parameters

311 ----------

312 exposure : `~lsst.afw.image.Exposure`

313 The exposure for which the average PSF is to be computed.

314 The exposure must contain a `psf` attribute.

315 psfExposureBuffer : `float`

316 Fractional buffer margin to be left out of all sides of the image

317 during the construction of the grid to compute average PSF in an

318 exposure.

319 psfExposureGrid : `int`

320 Grid size to compute the average PSF in an exposure.

321

322 Returns

323 -------

324 psfImage : `~lsst.afw.image.Image`

325 The average PSF across the exposure.

326

327 Raises

328 ------

329 ValueError

330 Raised if the PSF cannot be computed at any of the grid points.

331

332 See Also

333 --------

334 `evaluateMeanPsfFwhm`

335 """

336

337 psf = exposure.psf

338

339 bbox = exposure.getBBox()

340 xmax, ymax = bbox.getMax()

341 xmin, ymin = bbox.getMin()

342

343 xbuffer = psfExposureBuffer*(xmax-xmin)

344 ybuffer = psfExposureBuffer*(ymax-ymin)

345

346 nImg = 0

347 psfArray = None

348 for (x, y) in itertools.product(np.linspace(xmin+xbuffer, xmax-xbuffer, psfExposureGrid),

349 np.linspace(ymin+ybuffer, ymax-ybuffer, psfExposureGrid)

350 ):

351 pos = geom.Point2D(x, y)

352 try:

353 singleImage = psf.computeKernelImage(pos)

354 except InvalidParameterError:

355 _LOG.debug("Unable to compute PSF image at position (%f, %f).", x, y)

356 continue

357

358 if psfArray is None:

359 psfArray = singleImage.array

360 else:

361 psfArray += singleImage.array

362 nImg += 1

363

364 if psfArray is None:

365 raise ValueError("Unable to compute PSF image at any position on the exposure.")

366

367 psfImage = afwImage.ImageD(psfArray/nImg)

368 return psfImage

369

370

371def computeRobustStatistics(image, mask, statsCtrl, statistic=afwMath.MEANCLIP):

372 """Calculate a robust mean of the variance plane of an exposure.

373

374 Parameters

375 ----------

376 image : `lsst.afw.image.Image`

377 Image or variance plane of an exposure to evaluate.

378 mask : `lsst.afw.image.Mask`

379 Mask plane to use for excluding pixels.

380 statsCtrl : `lsst.afw.math.StatisticsControl`

381 Statistics control object for configuring the calculation.

382 statistic : `lsst.afw.math.Property`, optional

383 The type of statistic to compute. Typical values are

384 ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.

385

386 Returns

387 -------

388 value : `float`

389 The result of the statistic calculated from the unflagged pixels.

390 """

391 statObj = afwMath.makeStatistics(image, mask, statistic, statsCtrl)

392 return statObj.getValue(statistic)

393

394

395def computePSFNoiseEquivalentArea(psf):

396 """Compute the noise equivalent area for an image psf

397

398 Parameters

399 ----------

400 psf : `lsst.afw.detection.Psf`

401

402 Returns

403 -------

404 nea : `float`

405 """

406 psfImg = psf.computeImage(psf.getAveragePosition())

407 nea = 1./np.sum(psfImg.array**2)

408 return nea

409

410

411def angleMean(angles):

412 """Calculate the mean of an array of angles.

413

414 Parameters

415 ----------

416 angles : `ndarray`

417 An array of angles, in radians

418

419 Returns

420 -------

421 `lsst.geom.Angle`

422 The mean angle

423 """

424 complexArray = [complex(np.cos(angle), np.sin(angle)) for angle in angles]

425 return (geom.Angle(np.angle(np.mean(complexArray))))

426

427

428def evaluateMaskFraction(mask, maskPlane):

429 """Evaluate the fraction of masked pixels in a mask plane.

430

431 Parameters

432 ----------

433 mask : `lsst.afw.image.Mask`

434 The mask to evaluate the fraction on

435 maskPlane : `str`

436 The particular mask plane to evaluate the fraction

437

438 Returns

439 -------

440 value : `float`

441 The calculated fraction of masked pixels

442 """

443 nMaskSet = np.count_nonzero((mask.array & mask.getPlaneBitMask(maskPlane)))

444 return nMaskSet/mask.array.size

445

446

447def computeDifferenceImageMetrics(science, difference, stars, sky_sources=None):

448 r"""Compute quality metrics (saved to the task metadata) on the

449 difference image, at the locations of detected stars on the science

450 image. This restricts the metric to locations that should be

451 well-subtracted.

452

453 Parameters

454 ----------

455 science : `lsst.afw.image.ExposureF`

456 Science exposure that was subtracted.

457 difference : `lsst.afw.image.ExposureF`

458 Result of subtracting template and science.

459 stars : `lsst.afw.table.SourceCatalog`

460 Good calibration sources detected on science image; these

461 footprints are what the metrics are computed on.

462

463 Returns

464 -------

465 metrics : `lsst.pipe.base.Struct`

466

467 ``differenceFootprintRatioMean`` : `float`

468 Mean of the ratio of the absolute value of the difference image

469 (with the mean absolute value of the sky regions on the difference

470 image removed) to the science image, computed in the footprints

471 of stars detected on the science image (the sums below are of the

472 pixels in each star or sky footprint):

473 :math:`\mathrm{mean}_{footprints}((\sum |difference| -

474 \mathrm{mean}(\sum |difference_{sky}|)) / \sum science)`

475 ``differenceFootprintRatioStdev`` : `float`

476 Standard Deviation across footprints of the above ratio.

477 ``differenceFootprintSkyRatioMean`` : `float`

478 Mean of the ratio of the absolute value of sky source regions on

479 the difference image to the science image (the sum below is of the

480 pixels in each sky source footprint):

481 :math:`\mathrm{mean}_{footprints}(\sum |difference_{sky}| / \sum science_{sky})`

482 ``differenceFootprintSkyRatioStdev`` : `float`

483 Standard Deivation across footprints of the above sky ratio.

484 """

485 def footprint_mean(sources, sky=0):

486 """Compute ratio of the absolute value of the diffim to the science

487 image, within each source footprint, subtracting the sky from the

488 diffim values if provided.

489 """

490 n = len(sources)

491 science_footprints = np.zeros(n)

492 difference_footprints = np.zeros(n)

493 ratio = np.zeros(n)

494 for i, record in enumerate(sources):

495 footprint = record.getFootprint()

496 heavy = afwDetection.makeHeavyFootprint(footprint, science.maskedImage)

497 heavy_diff = afwDetection.makeHeavyFootprint(footprint, difference.maskedImage)

498 science_footprints[i] = abs(heavy.getImageArray()).mean()

499 difference_footprints[i] = abs(heavy_diff.getImageArray()).mean()

500 ratio[i] = abs((difference_footprints[i] - sky)/science_footprints[i])

501 return science_footprints, difference_footprints, ratio

502

503 if "sky_source" in stars.schema: 503 ↛ 509line 503 didn't jump to line 509 because the condition on line 503 was always true

504 sky = stars["sky_source"]

505 selectStars = stars[~sky]

506 if sky_sources is None:

507 sky_sources = stars[sky]

508 else:

509 selectStars = stars

510 # Note that the len() below is only evaluated if sky_sources is not None

511 if sky_sources is not None and len(sky_sources) > 0:

512 sky_science, sky_difference, sky_ratio = footprint_mean(sky_sources)

513 sky_mean = sky_ratio.mean()

514 sky_std = sky_ratio.std()

515 sky_difference = sky_difference.mean()

516 else:

517 sky_mean = np.nan

518 sky_std = np.nan

519 sky_difference = np.nanmedian(np.abs(difference.image.array))

520 science_footprints, difference_footprints, ratio = footprint_mean(selectStars, sky_difference)

521 return lsst.pipe.base.Struct(differenceFootprintRatioMean=ratio.mean(),

522 differenceFootprintRatioStdev=ratio.std(),

523 differenceFootprintSkyRatioMean=sky_mean,

524 differenceFootprintSkyRatioStdev=sky_std,

525 )

526

527

528def populate_sattle_visit_cache(visit_info, historical=False):

529 """Populate a cache of predicted satellite positions in the sattle service.

530

531 Parameters

532 ----------

533 visit_info: `lsst.afw.table.ExposureRecord.visitInfo`

534 Visit info for the science exposure being processed.

535 historical: `bool`

536 Set to True if observations are older than the current day.

537

538 Raises

539 ------

540 requests.HTTPError

541 Raised if sattle call does not return success.

542 """

543

544 visit_mjd = visit_info.getDate().toAstropy().mjd

545

546 exposure_time_days = visit_info.getExposureTime() / 86400.0

547 exposure_end_mjd = visit_mjd + exposure_time_days / 2.0

548 exposure_start_mjd = visit_mjd - exposure_time_days / 2.0

549

550 boresight_ra = visit_info.boresightRaDec.getRa().asDegrees()

551 boresight_dec = visit_info.boresightRaDec.getDec().asDegrees()

552

553 r = requests.put(

554 f'{os.getenv("SATTLE_URI_BASE")}/visit_cache',

555 json={"visit_id": visit_info.getId(),

556 "exposure_start_mjd": exposure_start_mjd,

557 "exposure_end_mjd": exposure_end_mjd,

558 "boresight_ra": boresight_ra,

559 "boresight_dec": boresight_dec,

560 "historical": historical})

561

562 r.raise_for_status()

563

564

565def checkMask(mask, sources, excludeMaskPlanes):

566 """Exclude sources that have masked pixels in their footprints.

567

568 Parameters

569 ----------

570 mask : `lsst.afw.image.Mask`

571 The image mask plane to use to reject sources

572 based on the location of their centroid on the ccd.

573 sources : `lsst.afw.table.SourceCatalog`

574 The source catalog to evaluate.

575 excludeMaskPlanes : `list` of `str`

576 List of the names of the mask planes to exclude.

577

578 Returns

579 -------

580 good : `numpy.ndarray` of `bool`

581 Array indicating whether each source in the catalog should be

582 kept (True) or rejected (False) based on the value of the

583 mask plane at its location.

584 """

585 setExcludeMaskPlanes = [

586 maskPlane for maskPlane in excludeMaskPlanes if maskPlane in mask.getMaskPlaneDict()

587 ]

588

589 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)

590

591 good = np.ones(len(sources), dtype=bool)

592 for i, source in enumerate(sources):

593 bbox = source.getFootprint().getBBox()

594 if not mask.getBBox().contains(bbox):

595 good[i] = False

596 continue

597

598 # Reject footprints with any bad mask bits set.

599 if (mask.subset(bbox).array & excludePixelMask).any():

600 good[i] = False

601 continue

602 return good

603

604

605def setSourceFootprints(sources, kernelSize):

606 """Add footprints of fixed size to a source catalog

607

608 Each source's footprint is replaced by a centered box of side

609 ``2*kernelSize + 1`` pixels (always odd), preserving the brightest peak

610 of the original footprint.

611

612 Parameters

613 ----------

614 sources : `lsst.afw.table.SourceCatalog`

615 The source catalog to add footprints to. Modified in place.

616 kernelSize : `int`

617 Width (in pixels) of the PSF-matching kernel. The footprint size is

618 ``2*kernelSize + 1`` so that the kernel can be convolved across the

619 source without falling off the stamp.

620

621 Returns

622 -------

623 sources : `lsst.afw.table.SourceCatalog`

624 The input `sources` catalog that was modified in place,

625 returned for convenience.

626 """

627 size = 2*kernelSize + 1

628 for source in sources:

629 bbox = lsst.geom.Box2I.makeCenteredBox(source.getCentroid(),

630 lsst.geom.Extent2I(size, size))

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

632 boxFootprint = lsst.afw.detection.Footprint(lsst.afw.geom.SpanSet(bbox))

633 boxFootprint.addPeak(peak.getFx(), peak.getFy(), peak.getPeakValue())

634 source.setFootprint(boxFootprint)

635 return sources