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

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 ]

29import itertools

30import os

31import requests

33from astropy.stats import gaussian_sigma_to_fwhm

34import numpy as np

36import lsst.afw.detection as afwDetection

37import lsst.afw.image as afwImage

38import lsst.afw.math as afwMath

39import lsst.geom as geom

40from lsst.pex.exceptions import InvalidParameterError, RangeError

41import lsst.pipe.base

42from lsst.utils.logging import getLogger

45_LOG = getLogger(__name__)

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

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

50 position.

52 Parameters

53 ----------

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

55 The PSF matching kernel to evaluate.

56 x : `float`

57 The x position on the detector to evaluate the kernel

58 y : `float`

59 The y position on the detector to evaluate the kernel

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

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

63 Returns

64 -------

65 kernel_sum : `float`

66 The sum of the kernel on the desired location

67 dx : `float`

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

69 center of the extraction of the kernel

70 dy : `float`

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

72 center of the extraction of the kernel

73 pos_angle: `float`

74 The position angle in detector coordinates of the displacement

75 length : `float`

76 The displacement module of the kernel centroid in pixel units

77 """

79 if image is None:

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

82 # obtain the kernel image

83 hsize = kernel.getWidth()//2

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

86 data = image.array

87 h, w = data.shape

88 xx = np.arange(w)

89 yy = np.arange(h)

91 # create sum vectors and estimate weighted average

92 vx = data.sum(axis=0)

93 vx /= vx.sum()

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

96 vy = data.sum(axis=1)

97 vy /= vy.sum()

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

100 # obtain position angle and norm of displacement

101 pos_angle = np.arctan2(dy, dx)

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

103

104 return kernel_sum, dx, dy, pos_angle, length

105

106

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

108 """Directly calculate the horizontal and vertical widths

109 of a PSF at half its maximum value.

110

111 Parameters

112 ----------

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

114 Point spread function (PSF) to evaluate.

115 average : `bool`, optional

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

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

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

119 average position is used.

120

121 Returns

122 -------

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

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

125 Returns the widths along the Y and X axes,

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

127

128 See Also

129 --------

130 evaluateMeanPsfFwhm

131 """

132 if position is None:

133 position = psf.getAveragePosition()

134 shape = psf.computeShape(position)

135

136 if average:

137 return gaussian_sigma_to_fwhm*shape.getTraceRadius()

138 else:

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

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

141

142

143def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,

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

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

146

147 Parameters

148 ----------

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

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

151 The exposure must contain a `psf` attribute.

152 fwhmExposureBuffer : `float`

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

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

155 exposure.

156 fwhmExposureGrid : `int`

157 Grid size to compute the mean FWHM in an exposure.

158

159 Returns

160 -------

161 meanFwhm : `float`

162 The mean PSF FWHM on the exposure.

163

164 Raises

165 ------

166 ValueError

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

168

169 See Also

170 --------

171 `getPsfFwhm`

172 `computeAveragePsf`

173 """

174

175 psf = exposure.psf

176

177 bbox = exposure.getBBox()

178 xmax, ymax = bbox.getMax()

179 xmin, ymin = bbox.getMin()

180

181 xbuffer = fwhmExposureBuffer*(xmax-xmin)

182 ybuffer = fwhmExposureBuffer*(ymax-ymin)

183

184 width = []

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

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

187 ):

188 pos = geom.Point2D(x, y)

189 try:

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

191 except (InvalidParameterError, RangeError):

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

193 continue

194

195 width.append(fwhm)

196

197 if not width:

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

199

200 return np.nanmean(width)

201

202

203def computeAveragePsf(exposure: afwImage.Exposure,

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

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

206

207 Parameters

208 ----------

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

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

211 The exposure must contain a `psf` attribute.

212 psfExposureBuffer : `float`

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

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

215 exposure.

216 psfExposureGrid : `int`

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

218

219 Returns

220 -------

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

222 The average PSF across the exposure.

223

224 Raises

225 ------

226 ValueError

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

228

229 See Also

230 --------

231 `evaluateMeanPsfFwhm`

232 """

233

234 psf = exposure.psf

235

236 bbox = exposure.getBBox()

237 xmax, ymax = bbox.getMax()

238 xmin, ymin = bbox.getMin()

239

240 xbuffer = psfExposureBuffer*(xmax-xmin)

241 ybuffer = psfExposureBuffer*(ymax-ymin)

242

243 nImg = 0

244 psfArray = None

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

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

247 ):

248 pos = geom.Point2D(x, y)

249 try:

250 singleImage = psf.computeKernelImage(pos)

251 except InvalidParameterError:

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

253 continue

254

255 if psfArray is None:

256 psfArray = singleImage.array

257 else:

258 psfArray += singleImage.array

259 nImg += 1

260

261 if psfArray is None:

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

263

264 psfImage = afwImage.ImageD(psfArray/nImg)

265 return psfImage

266

267

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

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

270

271 Parameters

272 ----------

273 image : `lsst.afw.image.Image`

274 Image or variance plane of an exposure to evaluate.

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

276 Mask plane to use for excluding pixels.

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

278 Statistics control object for configuring the calculation.

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

280 The type of statistic to compute. Typical values are

281 ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.

282

283 Returns

284 -------

285 value : `float`

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

287 """

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

289 return statObj.getValue(statistic)

290

291

292def computePSFNoiseEquivalentArea(psf):

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

294

295 Parameters

296 ----------

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

298

299 Returns

300 -------

301 nea : `float`

302 """

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

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

305 return nea

306

307

308def angleMean(angles):

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

310

311 Parameters

312 ----------

313 angles : `ndarray`

314 An array of angles, in radians

315

316 Returns

317 -------

318 `lsst.geom.Angle`

319 The mean angle

320 """

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

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

323

324

325def evaluateMaskFraction(mask, maskPlane):

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

327

328 Parameters

329 ----------

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

331 The mask to evaluate the fraction on

332 maskPlane : `str`

333 The particular mask plane to evaluate the fraction

334

335 Returns

336 -------

337 value : `float`

338 The calculated fraction of masked pixels

339 """

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

341 return nMaskSet/mask.array.size

342

343

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

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

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

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

348 well-subtracted.

349

350 Parameters

351 ----------

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

353 Science exposure that was subtracted.

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

355 Result of subtracting template and science.

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

357 Good calibration sources detected on science image; these

358 footprints are what the metrics are computed on.

359

360 Returns

361 -------

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

363

364 ``differenceFootprintRatioMean`` : `float`

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

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

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

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

369 pixels in each star or sky footprint):

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

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

372 ``differenceFootprintRatioStdev`` : `float`

373 Standard Deviation across footprints of the above ratio.

374 ``differenceFootprintSkyRatioMean`` : `float`

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

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

377 pixels in each sky source footprint):

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

379 ``differenceFootprintSkyRatioStdev`` : `float`

380 Standard Deivation across footprints of the above sky ratio.

381 """

382 def footprint_mean(sources, sky=0):

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

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

385 diffim values if provided.

386 """

387 n = len(sources)

388 science_footprints = np.zeros(n)

389 difference_footprints = np.zeros(n)

390 ratio = np.zeros(n)

391 for i, record in enumerate(sources):

392 footprint = record.getFootprint()

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

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

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

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

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

398 return science_footprints, difference_footprints, ratio

399

400 if "sky_source" in stars.schema:

401 sky = stars["sky_source"]

402 selectStars = stars[~sky]

403 if sky_sources is None:

404 sky_sources = stars[sky]

405 else:

406 selectStars = stars

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

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

409 sky_science, sky_difference, sky_ratio = footprint_mean(sky_sources)

410 sky_mean = sky_ratio.mean()

411 sky_std = sky_ratio.std()

412 sky_difference = sky_difference.mean()

413 else:

414 sky_mean = np.nan

415 sky_std = np.nan

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

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

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

419 differenceFootprintRatioStdev=ratio.std(),

420 differenceFootprintSkyRatioMean=sky_mean,

421 differenceFootprintSkyRatioStdev=sky_std,

422 )

423

424

425def populate_sattle_visit_cache(visit_info, historical=False):

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

427

428 Parameters

429 ----------

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

431 Visit info for the science exposure being processed.

432 historical: `bool`

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

434

435 Raises

436 ------

437 requests.HTTPError

438 Raised if sattle call does not return success.

439 """

440

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

442

443 exposure_time_days = visit_info.getExposureTime() / 86400.0

444 exposure_end_mjd = visit_mjd + exposure_time_days / 2.0

445 exposure_start_mjd = visit_mjd - exposure_time_days / 2.0

446

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

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

449

450 r = requests.put(

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

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

453 "exposure_start_mjd": exposure_start_mjd,

454 "exposure_end_mjd": exposure_end_mjd,

455 "boresight_ra": boresight_ra,

456 "boresight_dec": boresight_dec,

457 "historical": historical})

458

459 r.raise_for_status()

460

461

462def checkMask(mask, sources, excludeMaskPlanes):

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

464

465 Parameters

466 ----------

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

468 The image mask plane to use to reject sources

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

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

471 The source catalog to evaluate.

472 excludeMaskPlanes : `list` of `str`

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

474

475 Returns

476 -------

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

478 Array indicating whether each source in the catalog should be

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

480 mask plane at its location.

481 """

482 setExcludeMaskPlanes = [

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

484 ]

485

486 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)

487

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

489 for i, source in enumerate(sources):

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

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

492 good[i] = False

493 continue

494

495 # Reject footprints with any bad mask bits set.

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

497 good[i] = False

498 continue

499 return good

500

501

502def setSourceFootprints(sources, kernelSize):

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

504

505 Parameters

506 ----------

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

508 The source catalog to add footprints to.

509 kernelSize : `int`

510 The "radius" of the footprint, i.e half the size of the bounding box.

511

512 Returns

513 -------

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

515 The modified source catalog

516 """

517 size = 2*kernelSize + 1

518 for source in sources:

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

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

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

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

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

524 source.setFootprint(boxFootprint)

525 return sources

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

156 statements