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

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 numpy as np

31import os

32import requests

33import lsst.geom as geom

34import lsst.afw.detection as afwDetection

35import lsst.afw.image as afwImage

36import lsst.afw.math as afwMath

37from lsst.pex.exceptions import InvalidParameterError, RangeError

38import lsst.pipe.base

39from lsst.utils.logging import getLogger

42_LOG = getLogger(__name__)

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

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

47 position.

49 Parameters

50 ----------

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

52 The PSF matching kernel to evaluate.

53 x : `float`

54 The x position on the detector to evaluate the kernel

55 y : `float`

56 The y position on the detector to evaluate the kernel

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

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

60 Returns

61 -------

62 kernel_sum : `float`

63 The sum of the kernel on the desired location

64 dx : `float`

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

66 center of the extraction of the kernel

67 dy : `float`

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

69 center of the extraction of the kernel

70 pos_angle: `float`

71 The position angle in detector coordinates of the displacement

72 length : `float`

73 The displacement module of the kernel centroid in pixel units

74 """

76 if image is None:

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

79 # obtain the kernel image

80 hsize = kernel.getWidth()//2

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

83 data = image.array

84 h, w = data.shape

85 xx = np.arange(w)

86 yy = np.arange(h)

88 # create sum vectors and estimate weighted average

89 vx = data.sum(axis=0)

90 vx /= vx.sum()

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

93 vy = data.sum(axis=1)

94 vy /= vy.sum()

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

97 # obtain position angle and norm of displacement

98 pos_angle = np.arctan2(dy, dx)

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

100

101 return kernel_sum, dx, dy, pos_angle, length

102

103

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

105 """Directly calculate the horizontal and vertical widths

106 of a PSF at half its maximum value.

107

108 Parameters

109 ----------

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

111 Point spread function (PSF) to evaluate.

112 average : `bool`, optional

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

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

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

116 average position is used.

117

118 Returns

119 -------

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

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

122 Returns the widths along the Y and X axes,

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

124

125 See Also

126 --------

127 evaluateMeanPsfFwhm

128 """

129 if position is None:

130 position = psf.getAveragePosition()

131 shape = psf.computeShape(position)

132 sigmaToFwhm = 2*np.log(2*np.sqrt(2))

133

134 if average:

135 return sigmaToFwhm*shape.getTraceRadius()

136 else:

137 return [sigmaToFwhm*np.sqrt(shape.getIxx()), sigmaToFwhm*np.sqrt(shape.getIyy())]

138

139

140def evaluateMeanPsfFwhm(exposure: afwImage.Exposure,

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

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

143

144 Parameters

145 ----------

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

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

148 The exposure must contain a `psf` attribute.

149 fwhmExposureBuffer : `float`

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

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

152 exposure.

153 fwhmExposureGrid : `int`

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

155

156 Returns

157 -------

158 meanFwhm : `float`

159 The mean PSF FWHM on the exposure.

160

161 Raises

162 ------

163 ValueError

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

165

166 See Also

167 --------

168 `getPsfFwhm`

169 `computeAveragePsf`

170 """

171

172 psf = exposure.psf

173

174 bbox = exposure.getBBox()

175 xmax, ymax = bbox.getMax()

176 xmin, ymin = bbox.getMin()

177

178 xbuffer = fwhmExposureBuffer*(xmax-xmin)

179 ybuffer = fwhmExposureBuffer*(ymax-ymin)

180

181 width = []

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

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

184 ):

185 pos = geom.Point2D(x, y)

186 try:

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

188 except (InvalidParameterError, RangeError):

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

190 continue

191

192 width.append(fwhm)

193

194 if not width:

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

196

197 return np.nanmean(width)

198

199

200def computeAveragePsf(exposure: afwImage.Exposure,

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

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

203

204 Parameters

205 ----------

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

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

208 The exposure must contain a `psf` attribute.

209 psfExposureBuffer : `float`

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

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

212 exposure.

213 psfExposureGrid : `int`

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

215

216 Returns

217 -------

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

219 The average PSF across the exposure.

220

221 Raises

222 ------

223 ValueError

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

225

226 See Also

227 --------

228 `evaluateMeanPsfFwhm`

229 """

230

231 psf = exposure.psf

232

233 bbox = exposure.getBBox()

234 xmax, ymax = bbox.getMax()

235 xmin, ymin = bbox.getMin()

236

237 xbuffer = psfExposureBuffer*(xmax-xmin)

238 ybuffer = psfExposureBuffer*(ymax-ymin)

239

240 nImg = 0

241 psfArray = None

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

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

244 ):

245 pos = geom.Point2D(x, y)

246 try:

247 singleImage = psf.computeKernelImage(pos)

248 except InvalidParameterError:

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

250 continue

251

252 if psfArray is None:

253 psfArray = singleImage.array

254 else:

255 psfArray += singleImage.array

256 nImg += 1

257

258 if psfArray is None:

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

260

261 psfImage = afwImage.ImageD(psfArray/nImg)

262 return psfImage

263

264

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

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

267

268 Parameters

269 ----------

270 image : `lsst.afw.image.Image`

271 Image or variance plane of an exposure to evaluate.

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

273 Mask plane to use for excluding pixels.

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

275 Statistics control object for configuring the calculation.

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

277 The type of statistic to compute. Typical values are

278 ``afwMath.MEANCLIP`` or ``afwMath.STDEVCLIP``.

279

280 Returns

281 -------

282 value : `float`

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

284 """

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

286 return statObj.getValue(statistic)

287

288

289def computePSFNoiseEquivalentArea(psf):

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

291

292 Parameters

293 ----------

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

295

296 Returns

297 -------

298 nea : `float`

299 """

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

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

302 return nea

303

304

305def angleMean(angles):

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

307

308 Parameters

309 ----------

310 angles : `ndarray`

311 An array of angles, in radians

312

313 Returns

314 -------

315 `lsst.geom.Angle`

316 The mean angle

317 """

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

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

320

321

322def evaluateMaskFraction(mask, maskPlane):

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

324

325 Parameters

326 ----------

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

328 The mask to evaluate the fraction on

329 maskPlane : `str`

330 The particular mask plane to evaluate the fraction

331

332 Returns

333 -------

334 value : `float`

335 The calculated fraction of masked pixels

336 """

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

338 return nMaskSet/mask.array.size

339

340

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

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

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

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

345 well-subtracted.

346

347 Parameters

348 ----------

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

350 Science exposure that was subtracted.

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

352 Result of subtracting template and science.

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

354 Good calibration sources detected on science image; these

355 footprints are what the metrics are computed on.

356

357 Returns

358 -------

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

360

361 ``differenceFootprintRatioMean`` : `float`

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

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

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

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

366 pixels in each star or sky footprint):

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

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

369 ``differenceFootprintRatioStdev`` : `float`

370 Standard Deviation across footprints of the above ratio.

371 ``differenceFootprintSkyRatioMean`` : `float`

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

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

374 pixels in each sky source footprint):

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

376 ``differenceFootprintSkyRatioStdev`` : `float`

377 Standard Deivation across footprints of the above sky ratio.

378 """

379 def footprint_mean(sources, sky=0):

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

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

382 diffim values if provided.

383 """

384 n = len(sources)

385 science_footprints = np.zeros(n)

386 difference_footprints = np.zeros(n)

387 ratio = np.zeros(n)

388 for i, record in enumerate(sources):

389 footprint = record.getFootprint()

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

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

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

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

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

395 return science_footprints, difference_footprints, ratio

396

397 if "sky_source" in stars.schema:

398 sky = stars["sky_source"]

399 selectStars = stars[~sky]

400 if sky_sources is None:

401 sky_sources = stars[sky]

402 else:

403 selectStars = stars

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

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

406 sky_science, sky_difference, sky_ratio = footprint_mean(sky_sources)

407 sky_mean = sky_ratio.mean()

408 sky_std = sky_ratio.std()

409 sky_difference = sky_difference.mean()

410 else:

411 sky_mean = np.nan

412 sky_std = np.nan

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

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

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

416 differenceFootprintRatioStdev=ratio.std(),

417 differenceFootprintSkyRatioMean=sky_mean,

418 differenceFootprintSkyRatioStdev=sky_std,

419 )

420

421

422def populate_sattle_visit_cache(visit_info, historical=False):

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

424

425 Parameters

426 ----------

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

428 Visit info for the science exposure being processed.

429 historical: `bool`

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

431

432 Raises

433 ------

434 requests.HTTPError

435 Raised if sattle call does not return success.

436 """

437

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

439

440 exposure_time_days = visit_info.getExposureTime() / 86400.0

441 exposure_end_mjd = visit_mjd + exposure_time_days / 2.0

442 exposure_start_mjd = visit_mjd - exposure_time_days / 2.0

443

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

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

446

447 r = requests.put(

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

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

450 "exposure_start_mjd": exposure_start_mjd,

451 "exposure_end_mjd": exposure_end_mjd,

452 "boresight_ra": boresight_ra,

453 "boresight_dec": boresight_dec,

454 "historical": historical})

455

456 r.raise_for_status()

457

458

459def checkMask(mask, sources, excludeMaskPlanes):

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

461

462 Parameters

463 ----------

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

465 The image mask plane to use to reject sources

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

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

468 The source catalog to evaluate.

469 excludeMaskPlanes : `list` of `str`

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

471

472 Returns

473 -------

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

475 Array indicating whether each source in the catalog should be

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

477 mask plane at its location.

478 """

479 setExcludeMaskPlanes = [

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

481 ]

482

483 excludePixelMask = mask.getPlaneBitMask(setExcludeMaskPlanes)

484

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

486 for i, source in enumerate(sources):

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

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

489 good[i] = False

490 continue

491

492 # Reject footprints with any bad mask bits set.

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

494 good[i] = False

495 continue

496 return good

497

498

499def setSourceFootprints(sources, kernelSize):

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

501

502 Parameters

503 ----------

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

505 The source catalog to add footprints to.

506 kernelSize : `int`

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

508

509 Returns

510 -------

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

512 The modified source catalog

513 """

514 size = 2*kernelSize + 1

515 for source in sources:

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

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

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

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

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

521 source.setFootprint(boxFootprint)

522 return sources

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

156 statements