Coverage for python/lsst/analysis/ap/imageUtils.py: 98%

1# This file is part of analysis_ap. 

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"""Read AP images as `lsst.images` types, converting legacy afw exposures, 

23and cut them out on the sky. 

24""" 

25 

26__all__ = ["as_masked_image", "as_visit_image", "get_visit_image", "cutout", "cutout_box", "pixel_scale"] 

27 

28import warnings 

29 

30import astropy.units as u 

31import numpy as np 

32from astropy.coordinates import SkyCoord 

33 

34import lsst.afw.image as afwImage 

35import lsst.geom 

36from lsst.daf.butler import DatasetNotFoundError, DatasetRef 

37from lsst.images import ( 

38 Box, 

39 DifferenceImage, 

40 GeneralFrame, 

41 Image, 

42 MaskedImage, 

43 SkyProjection, 

44 VisitImage, 

45 get_legacy_difference_image_mask_planes, 

46 get_legacy_template_mask_planes, 

47 get_legacy_visit_image_mask_planes, 

48) 

49from lsst.images.fits import parse_legacy_bunit 

50 

51_IMAGE_TYPES = {"science": VisitImage, "template": DifferenceImage, "difference": DifferenceImage} 

52"""`lsst.images` type of each AP image role.""" 

53 

54_PLANE_MAPS = { 

55 "science": get_legacy_visit_image_mask_planes, 

56 "template": get_legacy_template_mask_planes, 

57 "difference": get_legacy_difference_image_mask_planes, 

58} 

59"""Legacy mask-plane map factory of each AP image role.""" 

60 

61 

62def _legacy_unit(exposure): 

63 """Return the pixel unit of a legacy exposure. 

64 

65 ``BUNIT`` is used if set. Otherwise an exposure whose photometric 

66 calibration is the identity is in nJy, and any other is instrumental, 

67 which ``from_legacy`` also calls electrons. 

68 """ 

69 if bunit := exposure.metadata.get("BUNIT"): 

70 return parse_legacy_bunit(bunit, instrumental_unit=u.electron) 

71 photo_calib = exposure.photoCalib 

72 if photo_calib is not None and photo_calib.getCalibrationMean() == 1.0: 

73 return u.nJy 

74 return u.electron 

75 

76 

77def as_visit_image(image, *, image_type="science", exposure_record=None, detector=None): 

78 """Return an AP image as an `lsst.images` type. 

79 

80 Parameters 

81 ---------- 

82 image : `lsst.images.VisitImage` or `lsst.afw.image.Exposure` 

83 Image to convert. An `lsst.images.VisitImage` (which includes 

84 `~lsst.images.DifferenceImage`) is returned unchanged. 

85 image_type : {"science", "template", "difference"}, optional 

86 Role of a legacy image, which sets the type it becomes (a 

87 `~lsst.images.VisitImage` for a science image, a 

88 `~lsst.images.DifferenceImage` otherwise) and the mask planes it 

89 keeps. 

90 exposure_record : `lsst.daf.butler.DimensionRecord`, optional 

91 ``exposure`` dimension record of a legacy image. Without it the 

92 converted observation information lacks the fields the record 

93 supplies. 

94 detector : `lsst.afw.cameraGeom.Detector` or \ 

95 `lsst.images.cameras.Detector`, optional 

96 Detector to attach to a legacy image that has none, such as a 

97 legacy template; the science image's detector is the right one. 

98 

99 Returns 

100 ------- 

101 image : `lsst.images.VisitImage` 

102 The image in `lsst.images` form. A converted image shares its image 

103 and variance pixels with ``image``. 

104 

105 Raises 

106 ------ 

107 TypeError 

108 Raised if ``image`` is neither an `lsst.images.VisitImage` nor an 

109 `lsst.afw.image.Exposure`. 

110 ValueError 

111 Raised if a legacy image has no detector and none is given. 

112 """ 

113 if isinstance(image, VisitImage): 

114 return image 

115 if not isinstance(image, afwImage.Exposure): 

116 raise TypeError(f"Cannot convert a {type(image).__name__} to an lsst.images VisitImage.") 

117 # Conversion edits the exposure's metadata, so convert a shallow copy 

118 # with its own metadata. 

119 exposure = type(image)(image, deep=False) 

120 exposure.setMetadata(image.getMetadata().deepCopy()) 

121 if exposure.getDetector() is None: 

122 if detector is None: 

123 raise ValueError(f"The legacy {image_type} image has no detector; pass the science image's " 

124 "detector to convert it.") 

125 exposure.setDetector(detector if hasattr(detector, "getBBox") else detector.to_legacy()) 

126 with warnings.catch_warnings(): 

127 # The analysis tools do not use the fields a missing record leaves 

128 # empty. 

129 warnings.filterwarnings("ignore", message="Converting without an exposure dimension record") 

130 return _IMAGE_TYPES[image_type].from_legacy(exposure, exposure_record=exposure_record, 

131 unit=_legacy_unit(image), 

132 plane_map=_PLANE_MAPS[image_type]()) 

133 

134 

135def as_masked_image(image, *, image_type="science"): 

136 """Return the pixels of an AP image as an `lsst.images.MaskedImage`. 

137 

138 Use this for images without the detector, PSF and other components a 

139 `~lsst.images.VisitImage` needs, such as a post-ISR image. 

140 

141 Parameters 

142 ---------- 

143 image : `lsst.images.MaskedImage` or `lsst.afw.image.Exposure` 

144 Image to convert. An `lsst.images.MaskedImage` (which includes 

145 `~lsst.images.VisitImage`) is returned unchanged. 

146 image_type : {"science", "template", "difference"}, optional 

147 Role of a legacy image, which sets the mask planes it keeps. 

148 

149 Returns 

150 ------- 

151 image : `lsst.images.MaskedImage` 

152 The image, mask and variance, with a sky projection if ``image`` has 

153 a WCS. A converted image shares its image and variance pixels with 

154 ``image``. 

155 

156 Raises 

157 ------ 

158 TypeError 

159 Raised if ``image`` is neither an `lsst.images.MaskedImage` nor an 

160 `lsst.afw.image.Exposure`. 

161 """ 

162 if isinstance(image, MaskedImage): 

163 return image 

164 if not isinstance(image, afwImage.Exposure): 

165 raise TypeError(f"Cannot convert a {type(image).__name__} to an lsst.images MaskedImage.") 

166 sky_projection = None 

167 if image.wcs is not None: 

168 sky_projection = SkyProjection.from_legacy(image.wcs, GeneralFrame(unit=u.pix)) 

169 return MaskedImage.from_legacy(image.maskedImage, unit=_legacy_unit(image), 

170 plane_map=_PLANE_MAPS[image_type](), sky_projection=sky_projection) 

171 

172 

173def _exposure_record(butler, data_id): 

174 """Return the ``exposure`` dimension record of a visit, or `None`. 

175 

176 Parameters 

177 ---------- 

178 butler : `lsst.daf.butler.Butler` 

179 Butler to query. 

180 data_id : `lsst.daf.butler.DataCoordinate` 

181 Data ID that includes the visit. 

182 """ 

183 records = list(butler.query_dimension_records("exposure", data_id=data_id, explain=False)) 

184 return records[0] if records else None 

185 

186 

187def get_visit_image(butler, dataset, data_id=None, *, image_type="science", detector=None, **kwargs): 

188 """Read an AP image as an `lsst.images` type, whatever its storage 

189 class. 

190 

191 Parameters 

192 ---------- 

193 butler : `lsst.daf.butler.Butler` 

194 Butler to read from. 

195 dataset : `str` or `lsst.daf.butler.DatasetRef` 

196 Dataset type name, or a resolved reference. 

197 data_id : `dict`, optional 

198 Data ID of the image; ignored when ``dataset`` is a reference. 

199 image_type : {"science", "template", "difference"}, optional 

200 Role of the image; see `as_visit_image`. 

201 detector : `lsst.afw.cameraGeom.Detector` or \ 

202 `lsst.images.cameras.Detector`, optional 

203 Detector for a legacy image that has none; see `as_visit_image`. 

204 **kwargs 

205 Additional data ID keys, passed to 

206 `lsst.daf.butler.Butler.find_dataset`. 

207 

208 Returns 

209 ------- 

210 image : `lsst.images.VisitImage` 

211 The image, converted with the ``exposure`` record of its visit if 

212 it is stored as a legacy exposure. 

213 

214 Raises 

215 ------ 

216 lsst.daf.butler.DatasetNotFoundError 

217 Raised if the butler has no such dataset. 

218 """ 

219 if isinstance(dataset, DatasetRef): 219 ↛ 220line 219 didn't jump to line 220 because the condition on line 219 was never true

220 ref = dataset 

221 else: 

222 ref = butler.find_dataset(dataset, data_id, **kwargs) 

223 if ref is None: 

224 raise DatasetNotFoundError(f"No {dataset!r} dataset with data ID {data_id}.") 

225 image = butler.get(ref) 

226 if isinstance(image, VisitImage): 

227 return image 

228 return as_visit_image(image, image_type=image_type, exposure_record=_exposure_record(butler, ref.dataId), 

229 detector=detector) 

230 

231 

232def _to_sky_coord(center): 

233 """Return a sky position as an `astropy.coordinates.SkyCoord`.""" 

234 if isinstance(center, lsst.geom.SpherePoint): 

235 return SkyCoord(center.getRa().asDegrees(), center.getDec().asDegrees(), unit="deg") 

236 return center 

237 

238 

239def cutout_box(image, bbox): 

240 """Copy the pixels of an image within a box. 

241 

242 Parameters 

243 ---------- 

244 image : `lsst.images.MaskedImage` 

245 Image to cut out of. 

246 bbox : `lsst.images.Box` 

247 Box to copy, which may extend past the edges of ``image``. 

248 

249 Returns 

250 ------- 

251 cutout : `lsst.images.Image` 

252 Copy of the image plane within ``bbox``; pixels outside ``image`` 

253 are NaN. 

254 

255 Raises 

256 ------ 

257 lsst.images.NoOverlapError 

258 Raised if ``bbox`` does not overlap ``image``. 

259 """ 

260 result = Image(np.nan, bbox=bbox, dtype=image.image.array.dtype, unit=image.unit, 

261 sky_projection=image.sky_projection) 

262 overlap = bbox.intersection(image.bbox) 

263 result[overlap].array[...] = image.image[overlap].array 

264 return result 

265 

266 

267def cutout(image, center, size): 

268 """Copy a square of pixels centered on a sky position. 

269 

270 Parameters 

271 ---------- 

272 image : `lsst.images.MaskedImage` 

273 Image to cut out of. 

274 center : `astropy.coordinates.SkyCoord` or `lsst.geom.SpherePoint` 

275 Sky position to center the cutout on. 

276 size : `int` 

277 Width and height of the cutout, in pixels. 

278 

279 Returns 

280 ------- 

281 cutout : `lsst.images.Image` 

282 Copy of the image plane in the square; pixels outside ``image`` are 

283 NaN. The square is placed as by `lsst.afw.image.Exposure.getCutout`. 

284 

285 Raises 

286 ------ 

287 lsst.images.NoOverlapError 

288 Raised if the square does not overlap ``image``. 

289 """ 

290 pixel = image.sky_projection.sky_to_pixel(_to_sky_coord(center)) 

291 legacy_box = lsst.geom.Box2I.makeCenteredBox(lsst.geom.Point2D(float(pixel.x), float(pixel.y)), 

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

293 return cutout_box(image, Box.from_legacy(legacy_box)) 

294 

295 

296def pixel_scale(image): 

297 """Return the pixel scale at the center of an image. 

298 

299 Parameters 

300 ---------- 

301 image : `lsst.images.MaskedImage` 

302 Image with a sky projection. 

303 

304 Returns 

305 ------- 

306 scale : `float` 

307 Square root of the pixel area on the sky, in arcseconds. 

308 """ 

309 x = 0.5*(image.bbox.x.start + image.bbox.x.stop - 1) 

310 y = 0.5*(image.bbox.y.start + image.bbox.y.stop - 1) 

311 sky = image.sky_projection.pixel_to_sky(x=np.array([x, x + 1.0, x]), y=np.array([y, y, y + 1.0])) 

312 dx = sky[0].spherical_offsets_to(sky[1]) 

313 dy = sky[0].spherical_offsets_to(sky[2]) 

314 area = abs(dx[0]*dy[1] - dx[1]*dy[0]) 

315 return float(np.sqrt(area.to_value(u.arcsec**2))) 

316 

317 

318def _sky_to_pixel(wcs, ra, dec): 

319 """Transform sky positions in degrees to pixels. 

320 

321 Parameters 

322 ---------- 

323 wcs : `lsst.images.SkyProjection` or `lsst.afw.geom.SkyWcs` 

324 Transform to apply. 

325 ra, dec : array-like 

326 Sky positions, in degrees. 

327 

328 Returns 

329 ------- 

330 x, y : `numpy.ndarray` 

331 Pixel positions. 

332 """ 

333 ra = np.asarray(ra, dtype=float) 

334 dec = np.asarray(dec, dtype=float) 

335 if hasattr(wcs, "sky_to_pixel"): 

336 pixels = wcs.sky_to_pixel(SkyCoord(ra, dec, unit="deg")) 

337 return np.asarray(pixels.x, dtype=float), np.asarray(pixels.y, dtype=float) 

338 return wcs.skyToPixelArray(ra, dec, degrees=True) 

339 

340 

341def _pixel_to_sky(wcs, x, y): 

342 """Transform pixel positions to sky positions in degrees. 

343 

344 Parameters 

345 ---------- 

346 wcs : `lsst.images.SkyProjection` or `lsst.afw.geom.SkyWcs` 

347 Transform to apply. 

348 x, y : array-like 

349 Pixel positions. 

350 

351 Returns 

352 ------- 

353 ra, dec : `numpy.ndarray` 

354 Sky positions, in degrees. 

355 """ 

356 x = np.asarray(x, dtype=float) 

357 y = np.asarray(y, dtype=float) 

358 if hasattr(wcs, "pixel_to_sky"): 

359 sky = wcs.pixel_to_sky(x=x, y=y) 

360 return np.asarray(sky.ra.deg), np.asarray(sky.dec.deg) 

361 return wcs.pixelToSkyArray(x, y, degrees=True) 

362 

363 

364def _box_edges(bbox): 

365 """Return the outer pixel edges of a box, ``(x_min, x_max, y_min, 

366 y_max)``. 

367 

368 Parameters 

369 ---------- 

370 bbox : `lsst.images.Box` or `lsst.geom.Box2I` 

371 Integer pixel box. 

372 """ 

373 if isinstance(bbox, lsst.geom.Box2I): 373 ↛ 375line 373 didn't jump to line 375 because the condition on line 373 was always true

374 bbox = Box.from_legacy(bbox) 

375 return (bbox.x.start - 0.5, bbox.x.stop - 0.5, bbox.y.start - 0.5, bbox.y.stop - 0.5)