Coverage for python/lsst/pipe/tasks/healSparseMapping.py: 41%

1# This file is part of pipe_tasks. 

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__ = ["HealSparseInputMapTask", "HealSparseInputMapConfig", 

23 "HealSparseMapFormatter", "HealSparsePropertyMapConnections", 

24 "HealSparsePropertyMapConfig", "HealSparsePropertyMapTask", 

25 "ConsolidateHealSparsePropertyMapConnections", 

26 "ConsolidateHealSparsePropertyMapConfig", 

27 "ConsolidateHealSparsePropertyMapTask"] 

28 

29from collections import defaultdict 

30import astropy.units 

31import esutil 

32import warnings 

33import numbers 

34import numpy as np 

35import hpgeom as hpg 

36import healsparse as hsp 

37 

38import lsst.pex.config as pexConfig 

39import lsst.pipe.base as pipeBase 

40import lsst.geom 

41import lsst.afw.geom as afwGeom 

42from lsst.daf.butler import Formatter 

43from lsst.skymap import BaseSkyMap 

44from lsst.utils.timer import timeMethod 

45from .healSparseMappingProperties import (BasePropertyMap, BasePropertyMapConfig, 

46 PropertyMapMap, compute_approx_psf_size_and_shape) 

47 

48 

49class HealSparseMapFormatter(Formatter): 

50 """Interface for reading and writing healsparse.HealSparseMap files.""" 

51 unsupportedParameters = frozenset() 

52 supportedExtensions = frozenset({".hsp", ".fit", ".fits"}) 

53 extension = '.hsp' 

54 

55 def read(self, component=None): 

56 # Docstring inherited from Formatter.read. 

57 path = self.fileDescriptor.location.path 

58 

59 if component == 'coverage': 

60 try: 

61 data = hsp.HealSparseCoverage.read(path) 

62 except (OSError, RuntimeError): 

63 raise ValueError(f"Unable to read healsparse map with URI {self.fileDescriptor.location.uri}") 

64 

65 return data 

66 

67 if self.fileDescriptor.parameters is None: 

68 pixels = None 

69 degrade_nside = None 

70 else: 

71 pixels = self.fileDescriptor.parameters.get('pixels', None) 

72 degrade_nside = self.fileDescriptor.parameters.get('degrade_nside', None) 

73 try: 

74 data = hsp.HealSparseMap.read(path, pixels=pixels, degrade_nside=degrade_nside) 

75 except (OSError, RuntimeError): 

76 raise ValueError(f"Unable to read healsparse map with URI {self.fileDescriptor.location.uri}") 

77 

78 return data 

79 

80 def write(self, inMemoryDataset): 

81 # Docstring inherited from Formatter.write. 

82 # Update the location with the formatter-preferred file extension 

83 self.fileDescriptor.location.updateExtension(self.extension) 

84 inMemoryDataset.write(self.fileDescriptor.location.path, clobber=True) 

85 

86 

87def _is_power_of_two(value): 

88 """Check that value is a power of two. 

89 

90 Parameters 

91 ---------- 

92 value : `int` 

93 Value to check. 

94 

95 Returns 

96 ------- 

97 is_power_of_two : `bool` 

98 True if value is a power of two; False otherwise, or 

99 if value is not an integer. 

100 """ 

101 if not isinstance(value, numbers.Integral): 101 ↛ 102line 101 didn't jump to line 102 because the condition on line 101 was never true

102 return False 

103 

104 # See https://stackoverflow.com/questions/57025836 

105 # Every power of 2 has exactly 1 bit set to 1; subtracting 

106 # 1 therefore flips every preceding bit. If you and that 

107 # together with the original value it must be 0. 

108 return (value & (value - 1) == 0) and value != 0 

109 

110 

111class HealSparseInputMapConfig(pexConfig.Config): 

112 """Configuration parameters for HealSparseInputMapTask""" 

113 nside = pexConfig.Field( 

114 doc="Mapping healpix nside. Must be power of 2.", 

115 dtype=int, 

116 default=32768, 

117 check=_is_power_of_two, 

118 ) 

119 nside_coverage = pexConfig.Field( 

120 doc="HealSparse coverage map nside. Must be power of 2.", 

121 dtype=int, 

122 default=256, 

123 check=_is_power_of_two, 

124 ) 

125 bad_mask_min_coverage = pexConfig.Field( 

126 doc=("Minimum area fraction of a map healpixel pixel that must be " 

127 "covered by bad pixels to be removed from the input map. " 

128 "This is approximate."), 

129 dtype=float, 

130 default=0.5, 

131 ) 

132 

133 

134class HealSparseInputMapTask(pipeBase.Task): 

135 """Task for making a HealSparse input map.""" 

136 

137 ConfigClass = HealSparseInputMapConfig 

138 _DefaultName = "healSparseInputMap" 

139 

140 def __init__(self, **kwargs): 

141 pipeBase.Task.__init__(self, **kwargs) 

142 

143 self.ccd_input_map = None 

144 self.cell_input_map = None 

145 

146 def build_ccd_input_map(self, bbox, wcs, ccds): 

147 """Build a map from ccd valid polygons or bounding boxes. 

148 

149 Parameters 

150 ---------- 

151 bbox : `lsst.geom.Box2I` 

152 Bounding box for region to build input map. 

153 wcs : `lsst.afw.geom.SkyWcs` 

154 WCS object for region to build input map. 

155 ccds : `lsst.afw.table.ExposureCatalog` 

156 Exposure catalog with ccd data from coadd inputs. 

157 """ 

158 with warnings.catch_warnings(): 

159 # Healsparse will emit a warning if nside coverage is greater than 

160 # 128. In the case of generating patch input maps, and not global 

161 # maps, high nside coverage works fine, so we can suppress this 

162 # warning. 

163 warnings.simplefilter("ignore") 

164 self.ccd_input_map = hsp.HealSparseMap.make_empty(nside_coverage=self.config.nside_coverage, 

165 nside_sparse=self.config.nside, 

166 dtype=hsp.WIDE_MASK, 

167 wide_mask_maxbits=len(ccds)) 

168 self._wcs = wcs 

169 self._bbox = bbox 

170 self._ccds = ccds 

171 

172 pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds() 

173 hpix_area_arcsec2 = hpg.nside_to_pixel_area(self.config.nside, degrees=True)*(3600.**2.) 

174 self._min_bad = self.config.bad_mask_min_coverage*hpix_area_arcsec2/(pixel_scale**2.) 

175 

176 metadata = {} 

177 self._bits_per_visit_ccd = {} 

178 self._bits_per_visit = defaultdict(list) 

179 for bit, ccd_row in enumerate(ccds): 

180 metadata[f"B{bit:04d}CCD"] = ccd_row["ccd"] 

181 metadata[f"B{bit:04d}VIS"] = ccd_row["visit"] 

182 metadata[f"B{bit:04d}WT"] = ccd_row["weight"] 

183 

184 self._bits_per_visit_ccd[(ccd_row["visit"], ccd_row["ccd"])] = bit 

185 self._bits_per_visit[ccd_row["visit"]].append(bit) 

186 

187 ccd_poly = ccd_row.getValidPolygon() 

188 if ccd_poly is None: 

189 ccd_poly = afwGeom.Polygon(lsst.geom.Box2D(ccd_row.getBBox())) 

190 # Detectors need to be rendered with their own wcs. 

191 ccd_poly_radec = self._pixels_to_radec(ccd_row.getWcs(), ccd_poly.convexHull().getVertices()) 

192 

193 # Create a ccd healsparse polygon 

194 poly = hsp.Polygon(ra=ccd_poly_radec[: -1, 0], 

195 dec=ccd_poly_radec[: -1, 1], 

196 value=[bit]) 

197 self.ccd_input_map.set_bits_pix(poly.get_pixels(nside=self.ccd_input_map.nside_sparse), 

198 [bit]) 

199 

200 # Cut down to the overall bounding box with associated wcs. 

201 bbox_afw_poly = afwGeom.Polygon(lsst.geom.Box2D(bbox)) 

202 bbox_poly_radec = self._pixels_to_radec(self._wcs, 

203 bbox_afw_poly.convexHull().getVertices()) 

204 bbox_poly = hsp.Polygon(ra=bbox_poly_radec[: -1, 0], dec=bbox_poly_radec[: -1, 1], 

205 value=np.arange(self.ccd_input_map.wide_mask_maxbits)) 

206 with warnings.catch_warnings(): 

207 warnings.simplefilter("ignore") 

208 bbox_poly_map = bbox_poly.get_map_like(self.ccd_input_map) 

209 self.ccd_input_map = hsp.and_intersection([self.ccd_input_map, bbox_poly_map]) 

210 self.ccd_input_map.metadata = metadata 

211 

212 # Create a temporary map to hold the count of bad pixels in each healpix pixel 

213 dtype = [(f"v{visit}", np.int64) for visit in self._bits_per_visit.keys()] 

214 

215 with warnings.catch_warnings(): 

216 # Healsparse will emit a warning if nside coverage is greater than 

217 # 128. In the case of generating patch input maps, and not global 

218 # maps, high nside coverage works fine, so we can suppress this 

219 # warning. 

220 warnings.simplefilter("ignore") 

221 self._ccd_input_bad_count_map = hsp.HealSparseMap.make_empty( 

222 nside_coverage=self.config.nside_coverage, 

223 nside_sparse=self.config.nside, 

224 dtype=dtype, 

225 primary=dtype[0][0]) 

226 

227 self._ccd_input_pixels = self.ccd_input_map.valid_pixels 

228 

229 # Don't set input bad map if there are no ccds which overlap the bbox. 

230 if len(self._ccd_input_pixels) > 0: 

231 # Ensure these are sorted. 

232 self._ccd_input_pixels = np.sort(self._ccd_input_pixels) 

233 

234 self._ccd_input_bad_count_map[self._ccd_input_pixels] = np.zeros(1, dtype=dtype) 

235 

236 def mask_warp_bbox(self, bbox, visit, mask, bit_mask_value): 

237 """Mask a subregion from a visit. 

238 This must be run after build_ccd_input_map initializes 

239 the overall map. 

240 

241 Parameters 

242 ---------- 

243 bbox : `lsst.geom.Box2I` 

244 Bounding box from region to mask. 

245 visit : `int` 

246 Visit number corresponding to warp with mask. 

247 mask : `lsst.afw.image.MaskX` 

248 Mask plane from warp exposure. 

249 bit_mask_value : `int` 

250 Bit mask to check for bad pixels. 

251 

252 Raises 

253 ------ 

254 RuntimeError : Raised if build_ccd_input_map was not run first. 

255 """ 

256 if self.ccd_input_map is None: 

257 raise RuntimeError("Must run build_ccd_input_map before mask_warp_bbox") 

258 

259 if len(self._ccd_input_pixels) == 0: 

260 # This tract has no coverage, so there is nothing to do. 

261 return 

262 

263 # Find the bad pixels and convert to healpix 

264 bad_pixels = np.where(mask.array & bit_mask_value) 

265 if len(bad_pixels[0]) == 0: 

266 # No bad pixels 

267 return 

268 

269 # Bad pixels come from warps which use the overall wcs. 

270 bad_ra, bad_dec = self._wcs.pixelToSkyArray( 

271 bad_pixels[1].astype(np.float64) + bbox.getMinX(), 

272 bad_pixels[0].astype(np.float64) + bbox.getMinY(), 

273 degrees=True, 

274 ) 

275 bad_hpix = hpg.angle_to_pixel(self.config.nside, bad_ra, bad_dec) 

276 

277 # Check if any of these "bad" pixels are in the valid footprint. 

278 match_input, match_bad = esutil.numpy_util.match(self._ccd_input_pixels, bad_hpix, presorted=True) 

279 if len(match_bad) == 0: 

280 return 

281 

282 bad_hpix = bad_hpix[match_bad] 

283 

284 # Create a view of the column we need to add to. 

285 count_map_visit = self._ccd_input_bad_count_map[f"v{visit}"] 

286 # Add the bad pixels to the accumulator. Note that the view 

287 # cannot append pixels, but the match above ensures we are 

288 # only adding to pixels that are already in the coverage 

289 # map and initialized. 

290 count_map_visit.update_values_pix(bad_hpix, 1, operation="add") 

291 

292 def finalize_ccd_input_map_mask(self): 

293 """Use accumulated mask information to finalize the masking of 

294 ccd_input_map. 

295 

296 Raises 

297 ------ 

298 RuntimeError : Raised if build_ccd_input_map was not run first. 

299 """ 

300 if self.ccd_input_map is None: 

301 raise RuntimeError("Must run build_ccd_input_map before finalize_ccd_input_map_mask.") 

302 

303 count_map_arr = self._ccd_input_bad_count_map[self._ccd_input_pixels] 

304 for visit in self._bits_per_visit: 

305 to_mask, = np.where(count_map_arr[f"v{visit}"] > self._min_bad) 

306 if to_mask.size == 0: 

307 continue 

308 self.ccd_input_map.clear_bits_pix(self._ccd_input_pixels[to_mask], 

309 self._bits_per_visit[visit]) 

310 

311 # Clear memory 

312 self._ccd_input_bad_count_map = None 

313 

314 def initialize_cell_input_map(self, bbox, wcs, visit_detectors): 

315 """Initialize the cell input map. 

316 

317 Parameters 

318 ---------- 

319 bbox : `lsst.geom.Box2I` 

320 Bounding box for region to build input map. 

321 wcs : `lsst.afw.geom.SkyWcs` 

322 WCS object for region to build input map. 

323 visit_detectors : `list` [`tuple`] 

324 List of visit/detector tuples. 

325 """ 

326 with warnings.catch_warnings(): 

327 # Healsparse will emit a warning if nside coverage is greater than 

328 # 128. In the case of generating patch input maps, and not global 

329 # maps, high nside coverage works fine, so we can suppress this 

330 # warning. 

331 warnings.simplefilter("ignore") 

332 self.cell_input_map = hsp.HealSparseMap.make_empty( 

333 nside_coverage=self.config.nside_coverage, 

334 nside_sparse=self.config.nside, 

335 dtype=hsp.WIDE_MASK, 

336 wide_mask_maxbits=len(visit_detectors), 

337 ) 

338 

339 self._wcs = wcs 

340 self._bbox = bbox 

341 self._visit_detectors = visit_detectors 

342 

343 metadata = {} 

344 self._bits_per_visit_detector = {} 

345 self._bits_per_visit = defaultdict(list) 

346 for bit, (visit, detector) in enumerate(visit_detectors): 

347 metadata[f"B{bit:04d}CCD"] = detector 

348 metadata[f"B{bit:04d}VIS"] = visit 

349 # Weight will be filled later. 

350 metadata[f"B{bit:04d}WT"] = 0.0 

351 

352 self._bits_per_visit_detector[(visit, detector)] = bit 

353 self._bits_per_visit[visit].append(bit) 

354 

355 self.cell_input_map.metadata = metadata 

356 

357 self._cell_pixels = {} 

358 self._visit_detector_cache = None 

359 self._detector_map_cache = None 

360 

361 def build_cell_input_map(self, cell): 

362 """Add a cell to the input map. 

363 

364 Parameters 

365 ---------- 

366 cell : `lsst.skymap.cellInfo.CellInfo` 

367 Cell to initialize. 

368 """ 

369 if self.cell_input_map is None: 

370 raise RuntimeError("Must run initialize_cell_input_map() before build_cell_input_map()") 

371 

372 # The input map only needs the *inner* sky polygon. 

373 cell_poly = cell.getInnerSkyPolygon() 

374 vertices = np.asarray([[v.x(), v.y(), v.z()] for v in cell_poly.getVertices()]) 

375 pixels = hpg.query_polygon_vec(self.config.nside, vertices) 

376 

377 self._cell_pixels[cell.sequential_index] = pixels 

378 

379 def add_warp_to_cell_input_map(self, ccd_row, weight, cell): 

380 """Add a warp to the input map for a given cell. 

381 

382 Parameters 

383 ---------- 

384 ccd_row : `lsst.afw.table.ExposureRecord` 

385 Row from the ccd table. 

386 weight : `float` 

387 Weight to use for this detector. 

388 cell : `lsst.skymap.cellInfo.CellInfo` 

389 Cell that overlaps the ccd_table_row. 

390 """ 

391 visit = int(ccd_row["visit"]) 

392 detector = int(ccd_row["ccd"]) 

393 

394 if (bit := self._bits_per_visit_detector.get((visit, detector), None)) is None: 

395 raise RuntimeError(f"Visit {visit} / detector {detector} not expected in map.") 

396 

397 if (cell_pixels := self._cell_pixels.get(cell.sequential_index, None)) is None: 

398 raise RuntimeError(f"Cell {cell.sequential_index} not expected in map.") 

399 

400 if (visit, detector) != self._visit_detector_cache: 

401 self._visit_detector_cache = (visit, detector) 

402 

403 ccd_poly = ccd_row.validPolygon 

404 if ccd_poly is None: 

405 ccd_poly = afwGeom.Polygon(lsst.geom.Box2D(ccd_row.getBBox())) 

406 # Detectors need to be rendered with their own wcs. 

407 ccd_poly_radec = self._pixels_to_radec(ccd_row.getWcs(), ccd_poly.convexHull().getVertices()) 

408 

409 poly = hsp.Polygon( 

410 ra=ccd_poly_radec[: -1, 0], 

411 dec=ccd_poly_radec[: -1, 1], 

412 value=True, 

413 ) 

414 with warnings.catch_warnings(): 

415 # Healsparse will emit a warning if nside coverage is greater 

416 # than 128. In the case of generating patch input maps, and 

417 # not global maps, high nside coverage works fine, so we can 

418 # suppress this warning. 

419 warnings.simplefilter("ignore") 

420 

421 self._detector_map_cache = poly.get_map( 

422 nside_coverage=self.config.nside_coverage, 

423 nside_sparse=self.config.nside, 

424 dtype=np.bool_, 

425 ) 

426 

427 self.cell_input_map.metadata[f"B{bit:04d}WT"] = weight 

428 

429 overlap = self._detector_map_cache[cell_pixels] 

430 self.cell_input_map.set_bits_pix(cell_pixels[overlap], bit) 

431 

432 def _pixels_to_radec(self, wcs, pixels): 

433 """Convert pixels to ra/dec positions using a wcs. 

434 

435 Parameters 

436 ---------- 

437 wcs : `lsst.afw.geom.SkyWcs` 

438 WCS object. 

439 pixels : `list` [`lsst.geom.Point2D`] 

440 List of pixels to convert. 

441 

442 Returns 

443 ------- 

444 radec : `numpy.ndarray` 

445 Nx2 array of ra/dec positions associated with pixels. 

446 """ 

447 sph_pts = wcs.pixelToSky(pixels) 

448 return np.array([(sph.getRa().asDegrees(), sph.getDec().asDegrees()) 

449 for sph in sph_pts]) 

450 

451 

452class HealSparsePropertyMapConnections(pipeBase.PipelineTaskConnections, 

453 dimensions=("tract", "band", "skymap",), 

454 defaultTemplates={"coaddName": "deep", 

455 "calexpType": ""}): 

456 input_maps = pipeBase.connectionTypes.Input( 

457 doc="Healsparse bit-wise coadd input maps", 

458 name="{coaddName}Coadd_inputMap", 

459 storageClass="HealSparseMap", 

460 dimensions=("tract", "patch", "skymap", "band"), 

461 multiple=True, 

462 deferLoad=True, 

463 ) 

464 coadd_exposures = pipeBase.connectionTypes.Input( 

465 doc="Coadded exposures associated with input_maps", 

466 name="{coaddName}Coadd_calexp", 

467 storageClass="ExposureF", 

468 dimensions=("tract", "patch", "skymap", "band"), 

469 multiple=True, 

470 deferLoad=True, 

471 ) 

472 visit_summaries = pipeBase.connectionTypes.Input( 

473 doc="Visit summary tables with aggregated statistics", 

474 name="finalVisitSummary", 

475 storageClass="ExposureCatalog", 

476 dimensions=("instrument", "visit"), 

477 multiple=True, 

478 deferLoad=True, 

479 ) 

480 sky_map = pipeBase.connectionTypes.Input( 

481 doc="Input definition of geometry/bbox and projection/wcs for coadded exposures", 

482 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

483 storageClass="SkyMap", 

484 dimensions=("skymap",), 

485 ) 

486 

487 # Create output connections for all possible maps defined in the 

488 # registry. The vars() trick used here allows us to set class attributes 

489 # programatically. Taken from 

490 # https://stackoverflow.com/questions/2519807/ 

491 # setting-a-class-attribute-with-a-given-name-in-python-while-defining-the-class 

492 for name in BasePropertyMap.registry: 

493 vars()[f"{name}_map_min"] = pipeBase.connectionTypes.Output( 

494 doc=f"Minimum-value map of {name}", 

495 name=f"{{coaddName}}Coadd_{name}_map_min", 

496 storageClass="HealSparseMap", 

497 dimensions=("tract", "skymap", "band"), 

498 ) 

499 vars()[f"{name}_map_max"] = pipeBase.connectionTypes.Output( 

500 doc=f"Maximum-value map of {name}", 

501 name=f"{{coaddName}}Coadd_{name}_map_max", 

502 storageClass="HealSparseMap", 

503 dimensions=("tract", "skymap", "band"), 

504 ) 

505 vars()[f"{name}_map_mean"] = pipeBase.connectionTypes.Output( 

506 doc=f"Mean-value map of {name}", 

507 name=f"{{coaddName}}Coadd_{name}_map_mean", 

508 storageClass="HealSparseMap", 

509 dimensions=("tract", "skymap", "band"), 

510 ) 

511 vars()[f"{name}_map_weighted_mean"] = pipeBase.connectionTypes.Output( 

512 doc=f"Weighted mean-value map of {name}", 

513 name=f"{{coaddName}}Coadd_{name}_map_weighted_mean", 

514 storageClass="HealSparseMap", 

515 dimensions=("tract", "skymap", "band"), 

516 ) 

517 vars()[f"{name}_map_sum"] = pipeBase.connectionTypes.Output( 

518 doc=f"Sum-value map of {name}", 

519 name=f"{{coaddName}}Coadd_{name}_map_sum", 

520 storageClass="HealSparseMap", 

521 dimensions=("tract", "skymap", "band"), 

522 ) 

523 

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

525 super().__init__(config=config) 

526 

527 # Not all possible maps in the registry will be configured to run. 

528 # Here we remove the unused connections. 

529 for name in BasePropertyMap.registry: 

530 if name not in config.property_maps: 

531 prop_config = BasePropertyMapConfig() 

532 prop_config.do_min = False 

533 prop_config.do_max = False 

534 prop_config.do_mean = False 

535 prop_config.do_weighted_mean = False 

536 prop_config.do_sum = False 

537 else: 

538 prop_config = config.property_maps[name] 

539 

540 if not prop_config.do_min: 

541 self.outputs.remove(f"{name}_map_min") 

542 if not prop_config.do_max: 

543 self.outputs.remove(f"{name}_map_max") 

544 if not prop_config.do_mean: 

545 self.outputs.remove(f"{name}_map_mean") 

546 if not prop_config.do_weighted_mean: 

547 self.outputs.remove(f"{name}_map_weighted_mean") 

548 if not prop_config.do_sum: 

549 self.outputs.remove(f"{name}_map_sum") 

550 

551 

552class HealSparsePropertyMapConfig(pipeBase.PipelineTaskConfig, 

553 pipelineConnections=HealSparsePropertyMapConnections): 

554 """Configuration parameters for HealSparsePropertyMapTask""" 

555 property_maps = BasePropertyMap.registry.makeField( 

556 multi=True, 

557 default=["exposure_time", 

558 "psf_size", 

559 "psf_e1", 

560 "psf_e2", 

561 "psf_maglim", 

562 "sky_noise", 

563 "sky_background", 

564 "dcr_dra", 

565 "dcr_ddec", 

566 "dcr_e1", 

567 "dcr_e2", 

568 "epoch"], 

569 doc="Property map computation objects", 

570 ) 

571 

572 def setDefaults(self): 

573 self.property_maps["exposure_time"].do_sum = True 

574 self.property_maps["psf_size"].do_weighted_mean = True 

575 self.property_maps["psf_e1"].do_weighted_mean = True 

576 self.property_maps["psf_e2"].do_weighted_mean = True 

577 self.property_maps["psf_maglim"].do_weighted_mean = True 

578 self.property_maps["sky_noise"].do_weighted_mean = True 

579 self.property_maps["sky_background"].do_weighted_mean = True 

580 self.property_maps["dcr_dra"].do_weighted_mean = True 

581 self.property_maps["dcr_ddec"].do_weighted_mean = True 

582 self.property_maps["dcr_e1"].do_weighted_mean = True 

583 self.property_maps["dcr_e2"].do_weighted_mean = True 

584 self.property_maps["epoch"].do_mean = True 

585 self.property_maps["epoch"].do_min = True 

586 self.property_maps["epoch"].do_max = True 

587 

588 

589class HealSparsePropertyMapTask(pipeBase.PipelineTask): 

590 """Task to compute Healsparse property maps. 

591 

592 This task will compute individual property maps (per tract, per 

593 map type, per band). These maps cover the full coadd tract, and 

594 are not truncated to the inner tract region. 

595 """ 

596 ConfigClass = HealSparsePropertyMapConfig 

597 _DefaultName = "healSparsePropertyMapTask" 

598 

599 def __init__(self, **kwargs): 

600 super().__init__(**kwargs) 

601 self.property_maps = PropertyMapMap() 

602 for name, config, PropertyMapClass in self.config.property_maps.apply(): 

603 self.property_maps[name] = PropertyMapClass(config, name) 

604 

605 @timeMethod 

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

607 inputs = butlerQC.get(inputRefs) 

608 

609 sky_map = inputs.pop("sky_map") 

610 

611 tract = butlerQC.quantum.dataId["tract"] 

612 band = butlerQC.quantum.dataId["band"] 

613 

614 input_map_dict = {ref.dataId["patch"]: ref for ref in inputs["input_maps"]} 

615 coadd_dict = {ref.dataId["patch"]: ref for ref in inputs["coadd_exposures"]} 

616 

617 visit_summary_dict = {ref.dataId["visit"]: ref.get() 

618 for ref in inputs["visit_summaries"]} 

619 

620 self.run(sky_map, tract, band, coadd_dict, input_map_dict, visit_summary_dict) 

621 

622 # Write the outputs 

623 for name, property_map in self.property_maps.items(): 

624 if property_map.config.do_min: 

625 butlerQC.put(property_map.min_map, 

626 getattr(outputRefs, f"{name}_map_min")) 

627 if property_map.config.do_max: 

628 butlerQC.put(property_map.max_map, 

629 getattr(outputRefs, f"{name}_map_max")) 

630 if property_map.config.do_mean: 

631 butlerQC.put(property_map.mean_map, 

632 getattr(outputRefs, f"{name}_map_mean")) 

633 if property_map.config.do_weighted_mean: 

634 butlerQC.put(property_map.weighted_mean_map, 

635 getattr(outputRefs, f"{name}_map_weighted_mean")) 

636 if property_map.config.do_sum: 

637 butlerQC.put(property_map.sum_map, 

638 getattr(outputRefs, f"{name}_map_sum")) 

639 

640 def run(self, sky_map, tract, band, coadd_dict, input_map_dict, visit_summary_dict): 

641 """Run the healsparse property task. 

642 

643 Parameters 

644 ---------- 

645 sky_map : Sky map object 

646 tract : `int` 

647 Tract number. 

648 band : `str` 

649 Band name for logging. 

650 coadd_dict : `dict` [`int`: `lsst.daf.butler.DeferredDatasetHandle`] 

651 Dictionary of coadd exposure references. Keys are patch numbers. 

652 input_map_dict : `dict` [`int`: `lsst.daf.butler.DeferredDatasetHandle`] 

653 Dictionary of input map references. Keys are patch numbers. 

654 visit_summary_dict : `dict` [`int`: `lsst.afw.table.ExposureCatalog`] 

655 Dictionary of visit summary tables. Keys are visit numbers. 

656 

657 Raises 

658 ------ 

659 RepeatableQuantumError 

660 If visit_summary_dict is missing any visits or detectors found in an 

661 input map. This leads to an inconsistency between what is in the coadd 

662 (via the input map) and the visit summary tables which contain data 

663 to compute the maps. 

664 """ 

665 tract_info = sky_map[tract] 

666 

667 tract_maps_initialized = False 

668 

669 for patch in input_map_dict.keys(): 

670 self.log.info("Making maps for band %s, tract %d, patch %d.", 

671 band, tract, patch) 

672 

673 patch_info = tract_info[patch] 

674 

675 input_map = input_map_dict[patch].get() 

676 

677 # Extract input map metadata. 

678 input_bit_weight_dict = {} 

679 for bit in range(input_map.wide_mask_maxbits): 

680 # Not all bits may be listed because maxbits must be a multiple 

681 # of 8. 

682 if f"B{bit:04d}CCD" in input_map.metadata: 

683 visit = input_map.metadata[f"B{bit:04d}VIS"] 

684 detector = input_map.metadata[f"B{bit:04d}CCD"] 

685 weight = input_map.metadata[f"B{bit:04d}WT"] 

686 input_bit_weight_dict[(visit, detector)] = (bit, weight) 

687 

688 # Initialize the tract maps as soon as we have the first input 

689 # map for getting nside information. 

690 if not tract_maps_initialized: 690 ↛ 705line 690 didn't jump to line 705 because the condition on line 690 was always true

691 # We use the first input map nside information to initialize 

692 # the tract maps 

693 nside_coverage = self._compute_nside_coverage_tract(tract_info) 

694 nside = input_map.nside_sparse 

695 

696 do_compute_approx_psf = False 

697 # Initialize the tract maps 

698 for property_map in self.property_maps: 

699 property_map.initialize_tract_maps(nside_coverage, nside) 

700 if property_map.requires_psf: 

701 do_compute_approx_psf = True 

702 

703 tract_maps_initialized = True 

704 

705 if input_map.valid_pixels.size == 0: 

706 self.log.warning("No valid pixels for band %s, tract %d, patch %d; skipping.", 

707 band, tract, patch) 

708 continue 

709 

710 # LSST coadds are now always in nJy, and the lsst.images formats 

711 # don't even have a PhotoCalib anymore. 

712 coadd_zeropoint = float((1.0 * astropy.units.nJy).to_value(astropy.units.ABmag)) 

713 

714 # Crop input_map to the inner polygon of the patch 

715 poly_vertices = patch_info.getInnerSkyPolygon(tract_info.getWcs()).getVertices() 

716 patch_radec = self._vertices_to_radec(poly_vertices) 

717 patch_poly = hsp.Polygon(ra=patch_radec[:, 0], dec=patch_radec[:, 1], 

718 value=np.arange(input_map.wide_mask_maxbits)) 

719 with warnings.catch_warnings(): 

720 # Healsparse will emit a warning if nside coverage is greater than 

721 # 128. In the case of generating patch input maps, and not global 

722 # maps, high nside coverage works fine, so we can suppress this 

723 # warning. 

724 warnings.simplefilter("ignore") 

725 patch_poly_map = patch_poly.get_map_like(input_map) 

726 input_map = hsp.and_intersection([input_map, patch_poly_map]) 

727 

728 valid_pixels, vpix_ra, vpix_dec = input_map.valid_pixels_pos(return_pixels=True) 

729 

730 # Check if there are no valid pixels for the inner (unique) patch region 

731 if valid_pixels.size == 0: 731 ↛ 732line 731 didn't jump to line 732 because the condition on line 731 was never true

732 continue 

733 

734 # Initialize the value accumulators 

735 for property_map in self.property_maps: 

736 property_map.initialize_values(valid_pixels.size) 

737 property_map.zeropoint = coadd_zeropoint 

738 

739 # Initialize the weight and counter accumulators 

740 total_weights = np.zeros(valid_pixels.size) 

741 total_inputs = np.zeros(valid_pixels.size, dtype=np.int32) 

742 

743 for (visit, detector), (bit, weight) in input_bit_weight_dict.items(): 

744 # Which pixels in the map are used by this visit/detector 

745 inmap, = np.where(input_map.check_bits_pix(valid_pixels, [bit])) 

746 

747 # Check if there are any valid pixels in the map from this deteector. 

748 if inmap.size == 0: 748 ↛ 749line 748 didn't jump to line 749 because the condition on line 748 was never true

749 continue 

750 

751 # Retrieve the correct visitSummary row 

752 if visit not in visit_summary_dict: 752 ↛ 753line 752 didn't jump to line 753 because the condition on line 752 was never true

753 msg = f"Visit {visit} not found in visit_summaries." 

754 raise pipeBase.RepeatableQuantumError(msg) 

755 row = visit_summary_dict[visit].find(detector) 

756 if row is None: 756 ↛ 757line 756 didn't jump to line 757 because the condition on line 756 was never true

757 msg = f"Visit {visit} / detector {detector} not found in visit_summaries." 

758 raise pipeBase.RepeatableQuantumError(msg) 

759 

760 x, y = row.wcs.skyToPixelArray(vpix_ra[inmap], vpix_dec[inmap], degrees=True) 

761 scalings = self._compute_calib_scale(row, x, y) 

762 

763 if do_compute_approx_psf: 763 ↛ 766line 763 didn't jump to line 766 because the condition on line 763 was always true

764 psf_array = compute_approx_psf_size_and_shape(row, vpix_ra[inmap], vpix_dec[inmap]) 

765 else: 

766 psf_array = None 

767 

768 total_weights[inmap] += weight 

769 total_inputs[inmap] += 1 

770 

771 # Accumulate the values 

772 for property_map in self.property_maps: 

773 property_map.accumulate_values(inmap, 

774 vpix_ra[inmap], 

775 vpix_dec[inmap], 

776 weight, 

777 scalings, 

778 row, 

779 psf_array=psf_array) 

780 

781 # Finalize the mean values and set the tract maps 

782 for property_map in self.property_maps: 

783 property_map.finalize_mean_values(total_weights, total_inputs) 

784 property_map.set_map_values(valid_pixels) 

785 

786 def _compute_calib_scale(self, ccd_row, x, y): 

787 """Compute calibration scaling values. 

788 

789 Parameters 

790 ---------- 

791 ccd_row : `lsst.afw.table.ExposureRecord` 

792 Exposure metadata for a given detector exposure.