22__all__ = [
"HealSparseInputMapTask",
"HealSparseInputMapConfig",
23 "HealSparseMapFormatter",
"HealSparsePropertyMapConnections",
24 "HealSparsePropertyMapConfig",
"HealSparsePropertyMapTask",
25 "ConsolidateHealSparsePropertyMapConnections",
26 "ConsolidateHealSparsePropertyMapConfig",
27 "ConsolidateHealSparsePropertyMapTask"]
29from collections
import defaultdict
35import healsparse
as hsp
41from lsst.daf.butler
import Formatter
43from lsst.utils.timer
import timeMethod
44from .healSparseMappingProperties
import (BasePropertyMap, BasePropertyMapConfig,
45 PropertyMapMap, compute_approx_psf_size_and_shape)
49 """Interface for reading and writing healsparse.HealSparseMap files."""
50 unsupportedParameters = frozenset()
51 supportedExtensions = frozenset({
".hsp",
".fit",
".fits"})
54 def read(self, component=None):
56 path = self.fileDescriptor.location.path
58 if component ==
'coverage':
60 data = hsp.HealSparseCoverage.read(path)
61 except (OSError, RuntimeError):
62 raise ValueError(f
"Unable to read healsparse map with URI {self.fileDescriptor.location.uri}")
66 if self.fileDescriptor.parameters
is None:
70 pixels = self.fileDescriptor.parameters.get(
'pixels',
None)
71 degrade_nside = self.fileDescriptor.parameters.get(
'degrade_nside',
None)
73 data = hsp.HealSparseMap.read(path, pixels=pixels, degrade_nside=degrade_nside)
74 except (OSError, RuntimeError):
75 raise ValueError(f
"Unable to read healsparse map with URI {self.fileDescriptor.location.uri}")
79 def write(self, inMemoryDataset):
82 self.fileDescriptor.location.updateExtension(self.
extension)
83 inMemoryDataset.write(self.fileDescriptor.location.path, clobber=
True)
87 """Check that value is a power of two.
96 is_power_of_two : `bool`
97 True if value is a power of two; False otherwise, or
98 if value is not an integer.
100 if not isinstance(value, numbers.Integral):
107 return (value & (value - 1) == 0)
and value != 0
111 """Configuration parameters for HealSparseInputMapTask"""
112 nside = pexConfig.Field(
113 doc=
"Mapping healpix nside. Must be power of 2.",
116 check=_is_power_of_two,
118 nside_coverage = pexConfig.Field(
119 doc=
"HealSparse coverage map nside. Must be power of 2.",
122 check=_is_power_of_two,
124 bad_mask_min_coverage = pexConfig.Field(
125 doc=(
"Minimum area fraction of a map healpixel pixel that must be "
126 "covered by bad pixels to be removed from the input map. "
127 "This is approximate."),
134 """Task for making a HealSparse input map."""
136 ConfigClass = HealSparseInputMapConfig
137 _DefaultName =
"healSparseInputMap"
140 pipeBase.Task.__init__(self, **kwargs)
145 """Build a map from ccd valid polygons or bounding boxes.
149 bbox : `lsst.geom.Box2I`
150 Bounding box for region to build input map.
151 wcs : `lsst.afw.geom.SkyWcs`
152 WCS object for region to build input map.
153 ccds : `lsst.afw.table.ExposureCatalog`
154 Exposure catalog with ccd data from coadd inputs.
156 with warnings.catch_warnings():
161 warnings.simplefilter(
"ignore")
162 self.
ccd_input_map = hsp.HealSparseMap.make_empty(nside_coverage=self.config.nside_coverage,
163 nside_sparse=self.config.nside,
165 wide_mask_maxbits=len(ccds))
170 pixel_scale = wcs.getPixelScale(bbox.getCenter()).asArcseconds()
171 hpix_area_arcsec2 = hpg.nside_to_pixel_area(self.config.nside, degrees=
True)*(3600.**2.)
172 self.
_min_bad = self.config.bad_mask_min_coverage*hpix_area_arcsec2/(pixel_scale**2.)
177 for bit, ccd_row
in enumerate(ccds):
178 metadata[f
"B{bit:04d}CCD"] = ccd_row[
"ccd"]
179 metadata[f
"B{bit:04d}VIS"] = ccd_row[
"visit"]
180 metadata[f
"B{bit:04d}WT"] = ccd_row[
"weight"]
185 ccd_poly = ccd_row.getValidPolygon()
189 ccd_poly_radec = self.
_pixels_to_radec(ccd_row.getWcs(), ccd_poly.convexHull().getVertices())
192 poly = hsp.Polygon(ra=ccd_poly_radec[: -1, 0],
193 dec=ccd_poly_radec[: -1, 1],
201 bbox_afw_poly.convexHull().getVertices())
202 bbox_poly = hsp.Polygon(ra=bbox_poly_radec[: -1, 0], dec=bbox_poly_radec[: -1, 1],
204 with warnings.catch_warnings():
205 warnings.simplefilter(
"ignore")
211 dtype = [(f
"v{visit}", np.int64)
for visit
in self.
_bits_per_visit.keys()]
213 with warnings.catch_warnings():
218 warnings.simplefilter(
"ignore")
220 nside_coverage=self.config.nside_coverage,
221 nside_sparse=self.config.nside,
235 """Mask a subregion from a visit.
236 This must be run after build_ccd_input_map initializes
241 bbox : `lsst.geom.Box2I`
242 Bounding box from region to mask.
244 Visit number corresponding to warp with mask.
245 mask : `lsst.afw.image.MaskX`
246 Mask plane from warp exposure.
247 bit_mask_value : `int`
248 Bit mask to check for bad pixels.
252 RuntimeError : Raised if build_ccd_input_map was not run first.
255 raise RuntimeError(
"Must run build_ccd_input_map before mask_warp_bbox")
262 bad_pixels = np.where(mask.array & bit_mask_value)
263 if len(bad_pixels[0]) == 0:
268 bad_ra, bad_dec = self.
_wcs.pixelToSkyArray(
269 bad_pixels[1].astype(np.float64) + bbox.getMinX(),
270 bad_pixels[0].astype(np.float64) + bbox.getMinY(),
273 bad_hpix = hpg.angle_to_pixel(self.config.nside, bad_ra, bad_dec)
276 match_input, match_bad = esutil.numpy_util.match(self.
_ccd_input_pixels, bad_hpix, presorted=
True)
277 if len(match_bad) == 0:
280 bad_hpix = bad_hpix[match_bad]
288 count_map_visit.update_values_pix(bad_hpix, 1, operation=
"add")
291 """Use accumulated mask information to finalize the masking of
296 RuntimeError : Raised if build_ccd_input_map was not run first.
299 raise RuntimeError(
"Must run build_ccd_input_map before finalize_ccd_input_map_mask.")
303 to_mask, = np.where(count_map_arr[f
"v{visit}"] > self.
_min_bad)
304 if to_mask.size == 0:
313 """Convert pixels to ra/dec positions using a wcs.
317 wcs : `lsst.afw.geom.SkyWcs`
319 pixels : `list` [`lsst.geom.Point2D`]
320 List of pixels to convert.
324 radec : `numpy.ndarray`
325 Nx2 array of ra/dec positions associated with pixels.
327 sph_pts = wcs.pixelToSky(pixels)
328 return np.array([(sph.getRa().asDegrees(), sph.getDec().asDegrees())
333 dimensions=(
"tract",
"band",
"skymap",),
334 defaultTemplates={
"coaddName":
"deep",
336 input_maps = pipeBase.connectionTypes.Input(
337 doc=
"Healsparse bit-wise coadd input maps",
338 name=
"{coaddName}Coadd_inputMap",
339 storageClass=
"HealSparseMap",
340 dimensions=(
"tract",
"patch",
"skymap",
"band"),
344 coadd_exposures = pipeBase.connectionTypes.Input(
345 doc=
"Coadded exposures associated with input_maps",
346 name=
"{coaddName}Coadd",
347 storageClass=
"ExposureF",
348 dimensions=(
"tract",
"patch",
"skymap",
"band"),
352 visit_summaries = pipeBase.connectionTypes.Input(
353 doc=
"Visit summary tables with aggregated statistics",
354 name=
"finalVisitSummary",
355 storageClass=
"ExposureCatalog",
356 dimensions=(
"instrument",
"visit"),
360 sky_map = pipeBase.connectionTypes.Input(
361 doc=
"Input definition of geometry/bbox and projection/wcs for coadded exposures",
362 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
363 storageClass=
"SkyMap",
364 dimensions=(
"skymap",),
372 for name
in BasePropertyMap.registry:
373 vars()[f
"{name}_map_min"] = pipeBase.connectionTypes.Output(
374 doc=f
"Minimum-value map of {name}",
375 name=f
"{{coaddName}}Coadd_{name}_map_min",
376 storageClass=
"HealSparseMap",
377 dimensions=(
"tract",
"skymap",
"band"),
379 vars()[f
"{name}_map_max"] = pipeBase.connectionTypes.Output(
380 doc=f
"Maximum-value map of {name}",
381 name=f
"{{coaddName}}Coadd_{name}_map_max",
382 storageClass=
"HealSparseMap",
383 dimensions=(
"tract",
"skymap",
"band"),
385 vars()[f
"{name}_map_mean"] = pipeBase.connectionTypes.Output(
386 doc=f
"Mean-value map of {name}",
387 name=f
"{{coaddName}}Coadd_{name}_map_mean",
388 storageClass=
"HealSparseMap",
389 dimensions=(
"tract",
"skymap",
"band"),
391 vars()[f
"{name}_map_weighted_mean"] = pipeBase.connectionTypes.Output(
392 doc=f
"Weighted mean-value map of {name}",
393 name=f
"{{coaddName}}Coadd_{name}_map_weighted_mean",
394 storageClass=
"HealSparseMap",
395 dimensions=(
"tract",
"skymap",
"band"),
397 vars()[f
"{name}_map_sum"] = pipeBase.connectionTypes.Output(
398 doc=f
"Sum-value map of {name}",
399 name=f
"{{coaddName}}Coadd_{name}_map_sum",
400 storageClass=
"HealSparseMap",
401 dimensions=(
"tract",
"skymap",
"band"),
404 def __init__(self, *, config=None):
405 super().__init__(config=config)
409 for name
in BasePropertyMap.registry:
410 if name
not in config.property_maps:
412 prop_config.do_min =
False
413 prop_config.do_max =
False
414 prop_config.do_mean =
False
415 prop_config.do_weighted_mean =
False
416 prop_config.do_sum =
False
418 prop_config = config.property_maps[name]
420 if not prop_config.do_min:
421 self.outputs.remove(f
"{name}_map_min")
422 if not prop_config.do_max:
423 self.outputs.remove(f
"{name}_map_max")
424 if not prop_config.do_mean:
425 self.outputs.remove(f
"{name}_map_mean")
426 if not prop_config.do_weighted_mean:
427 self.outputs.remove(f
"{name}_map_weighted_mean")
428 if not prop_config.do_sum:
429 self.outputs.remove(f
"{name}_map_sum")
432class HealSparsePropertyMapConfig(pipeBase.PipelineTaskConfig,
433 pipelineConnections=HealSparsePropertyMapConnections):
434 """Configuration parameters for HealSparsePropertyMapTask"""
435 property_maps = BasePropertyMap.registry.makeField(
437 default=[
"exposure_time",
449 doc=
"Property map computation objects",
452 def setDefaults(self):
453 self.property_maps[
"exposure_time"].do_sum =
True
454 self.property_maps[
"psf_size"].do_weighted_mean =
True
455 self.property_maps[
"psf_e1"].do_weighted_mean =
True
456 self.property_maps[
"psf_e2"].do_weighted_mean =
True
457 self.property_maps[
"psf_maglim"].do_weighted_mean =
True
458 self.property_maps[
"sky_noise"].do_weighted_mean =
True
459 self.property_maps[
"sky_background"].do_weighted_mean =
True
460 self.property_maps[
"dcr_dra"].do_weighted_mean =
True
461 self.property_maps[
"dcr_ddec"].do_weighted_mean =
True
462 self.property_maps[
"dcr_e1"].do_weighted_mean =
True
463 self.property_maps[
"dcr_e2"].do_weighted_mean =
True
464 self.property_maps[
"epoch"].do_mean =
True
465 self.property_maps[
"epoch"].do_min =
True
466 self.property_maps[
"epoch"].do_max =
True
469class HealSparsePropertyMapTask(pipeBase.PipelineTask):
470 """Task to compute Healsparse property maps.
472 This task will compute individual property maps (per tract, per
473 map type, per band). These maps cover the full coadd tract, and
474 are not truncated to the inner tract region.
476 ConfigClass = HealSparsePropertyMapConfig
477 _DefaultName =
"healSparsePropertyMapTask"
479 def __init__(self, **kwargs):
480 super().__init__(**kwargs)
482 for name, config, PropertyMapClass
in self.config.property_maps.apply():
483 self.property_maps[name] = PropertyMapClass(config, name)
486 def runQuantum(self, butlerQC, inputRefs, outputRefs):
487 inputs = butlerQC.get(inputRefs)
489 sky_map = inputs.pop(
"sky_map")
491 tract = butlerQC.quantum.dataId[
"tract"]
492 band = butlerQC.quantum.dataId[
"band"]
494 input_map_dict = {ref.dataId[
"patch"]: ref
for ref
in inputs[
"input_maps"]}
495 coadd_dict = {ref.dataId[
"patch"]: ref
for ref
in inputs[
"coadd_exposures"]}
497 visit_summary_dict = {ref.dataId[
"visit"]: ref.get()
498 for ref
in inputs[
"visit_summaries"]}
500 self.run(sky_map, tract, band, coadd_dict, input_map_dict, visit_summary_dict)
503 for name, property_map
in self.property_maps.items():
504 if property_map.config.do_min:
505 butlerQC.put(property_map.min_map,
506 getattr(outputRefs, f
"{name}_map_min"))
507 if property_map.config.do_max:
508 butlerQC.put(property_map.max_map,
509 getattr(outputRefs, f
"{name}_map_max"))
510 if property_map.config.do_mean:
511 butlerQC.put(property_map.mean_map,
512 getattr(outputRefs, f
"{name}_map_mean"))
513 if property_map.config.do_weighted_mean:
514 butlerQC.put(property_map.weighted_mean_map,
515 getattr(outputRefs, f
"{name}_map_weighted_mean"))
516 if property_map.config.do_sum:
517 butlerQC.put(property_map.sum_map,
518 getattr(outputRefs, f
"{name}_map_sum"))
520 def run(self, sky_map, tract, band, coadd_dict, input_map_dict, visit_summary_dict):
521 """Run the healsparse property task.
525 sky_map : Sky map object
529 Band name for logging.
530 coadd_dict : `dict` [`int`: `lsst.daf.butler.DeferredDatasetHandle`]
531 Dictionary of coadd exposure references. Keys are patch numbers.
532 input_map_dict : `dict` [`int`: `lsst.daf.butler.DeferredDatasetHandle`]
533 Dictionary of input map references. Keys are patch numbers.
534 visit_summary_dict : `dict` [`int`: `lsst.afw.table.ExposureCatalog`]
535 Dictionary of visit summary tables. Keys are visit numbers.
539 RepeatableQuantumError
540 If visit_summary_dict is missing any visits or detectors found in an
541 input map. This leads to an inconsistency between what is in the coadd
542 (via the input map) and the visit summary tables which contain data
545 tract_info = sky_map[tract]
547 tract_maps_initialized =
False
549 for patch
in input_map_dict.keys():
550 self.log.info(
"Making maps for band %s, tract %d, patch %d.",
553 patch_info = tract_info[patch]
555 input_map = input_map_dict[patch].get()
559 if not tract_maps_initialized:
562 nside_coverage = self._compute_nside_coverage_tract(tract_info)
563 nside = input_map.nside_sparse
565 do_compute_approx_psf =
False
567 for property_map
in self.property_maps:
568 property_map.initialize_tract_maps(nside_coverage, nside)
569 if property_map.requires_psf:
570 do_compute_approx_psf =
True
572 tract_maps_initialized =
True
574 if input_map.valid_pixels.size == 0:
575 self.log.warning(
"No valid pixels for band %s, tract %d, patch %d; skipping.",
579 coadd_photo_calib = coadd_dict[patch].get(component=
"photoCalib")
580 coadd_inputs = coadd_dict[patch].get(component=
"coaddInputs")
582 coadd_zeropoint = 2.5*np.log10(coadd_photo_calib.getInstFluxAtZeroMagnitude())
585 poly_vertices = patch_info.getInnerSkyPolygon(tract_info.getWcs()).getVertices()
586 patch_radec = self._vertices_to_radec(poly_vertices)
587 patch_poly = hsp.Polygon(ra=patch_radec[:, 0], dec=patch_radec[:, 1],
588 value=np.arange(input_map.wide_mask_maxbits))
589 with warnings.catch_warnings():
594 warnings.simplefilter(
"ignore")
595 patch_poly_map = patch_poly.get_map_like(input_map)
596 input_map = hsp.and_intersection([input_map, patch_poly_map])
598 valid_pixels, vpix_ra, vpix_dec = input_map.valid_pixels_pos(return_pixels=
True)
601 if valid_pixels.size == 0:
605 for property_map
in self.property_maps:
606 property_map.initialize_values(valid_pixels.size)
607 property_map.zeropoint = coadd_zeropoint
610 total_weights = np.zeros(valid_pixels.size)
611 total_inputs = np.zeros(valid_pixels.size, dtype=np.int32)
613 for bit, ccd_row
in enumerate(coadd_inputs.ccds):
615 inmap, = np.where(input_map.check_bits_pix(valid_pixels, [bit]))
622 visit = ccd_row[
"visit"]
623 detector_id = ccd_row[
"ccd"]
624 weight = ccd_row[
"weight"]
626 x, y = ccd_row.getWcs().skyToPixelArray(vpix_ra[inmap], vpix_dec[inmap], degrees=
True)
627 scalings = self._compute_calib_scale(ccd_row, x, y)
629 if do_compute_approx_psf:
630 psf_array = compute_approx_psf_size_and_shape(ccd_row, vpix_ra[inmap], vpix_dec[inmap])
634 total_weights[inmap] += weight
635 total_inputs[inmap] += 1
638 if visit
not in visit_summary_dict:
639 msg = f
"Visit {visit} not found in visit_summaries."
640 raise pipeBase.RepeatableQuantumError(msg)
641 row = visit_summary_dict[visit].find(detector_id)
643 msg = f
"Visit {visit} / detector_id {detector_id} not found in visit_summaries."
644 raise pipeBase.RepeatableQuantumError(msg)
647 for property_map
in self.property_maps:
648 property_map.accumulate_values(inmap,
657 for property_map
in self.property_maps:
658 property_map.finalize_mean_values(total_weights, total_inputs)
659 property_map.set_map_values(valid_pixels)
661 def _compute_calib_scale(self, ccd_row, x, y):
662 """Compute calibration scaling values.
666 ccd_row : `lsst.afw.table.ExposureRecord`
667 Exposure metadata for a given detector exposure.
669 Array of x positions.
671 Array of y positions.
675 calib_scale : `np.ndarray`
676 Array of calibration scale values.
678 photo_calib = ccd_row.getPhotoCalib()
679 bf = photo_calib.computeScaledCalibration()
680 if bf.getBBox() == ccd_row.getBBox():
682 calib_scale = photo_calib.getCalibrationMean()*bf.evaluate(x, y)
685 calib_scale = photo_calib.getCalibrationMean()
689 def _vertices_to_radec(self, vertices):
690 """Convert polygon vertices to ra/dec.
694 vertices : `list` [ `lsst.sphgeom.UnitVector3d` ]
695 Vertices for bounding polygon.
699 radec : `numpy.ndarray`
700 Nx2 array of ra/dec positions (in degrees) associated with vertices.
703 radec = np.array([(x.getLon().asDegrees(), x.getLat().asDegrees())
for
707 def _compute_nside_coverage_tract(self, tract_info):
708 """Compute the optimal coverage nside for a tract.
712 tract_info : `lsst.skymap.tractInfo.ExplicitTractInfo`
713 Tract information object.
717 nside_coverage : `int`
718 Optimal coverage nside for a tract map.
720 num_patches = tract_info.getNumPatches()
723 patch_info = tract_info.getPatchInfo(0)
724 vertices = patch_info.getInnerSkyPolygon(tract_info.getWcs()).getVertices()
725 radec = self._vertices_to_radec(vertices)
726 delta_ra = np.max(radec[:, 0]) - np.min(radec[:, 0])
727 delta_dec = np.max(radec[:, 1]) - np.min(radec[:, 1])
728 patch_area = delta_ra*delta_dec*np.cos(np.deg2rad(np.mean(radec[:, 1])))
730 tract_area = num_patches[0]*num_patches[1]*patch_area
732 nside_coverage_tract = 32
733 while hpg.nside_to_pixel_area(nside_coverage_tract, degrees=
True) > tract_area:
734 nside_coverage_tract = 2*nside_coverage_tract
737 nside_coverage_tract = int(np.clip(nside_coverage_tract/2, 32, 128))
739 return nside_coverage_tract
742class ConsolidateHealSparsePropertyMapConnections(pipeBase.PipelineTaskConnections,
743 dimensions=(
"band",
"skymap",),
744 defaultTemplates={
"coaddName":
"deep"}):
745 sky_map = pipeBase.connectionTypes.Input(
746 doc=
"Input definition of geometry/bbox and projection/wcs for coadded exposures",
747 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
748 storageClass=
"SkyMap",
749 dimensions=(
"skymap",),
757 for name
in BasePropertyMap.registry:
758 vars()[f
"{name}_map_min"] = pipeBase.connectionTypes.Input(
759 doc=f
"Minimum-value map of {name}",
760 name=f
"{{coaddName}}Coadd_{name}_map_min",
761 storageClass=
"HealSparseMap",
762 dimensions=(
"tract",
"skymap",
"band"),
766 vars()[f
"{name}_consolidated_map_min"] = pipeBase.connectionTypes.Output(
767 doc=f
"Minumum-value map of {name}",
768 name=f
"{{coaddName}}Coadd_{name}_consolidated_map_min",
769 storageClass=
"HealSparseMap",
770 dimensions=(
"skymap",
"band"),
772 vars()[f
"{name}_map_max"] = pipeBase.connectionTypes.Input(
773 doc=f
"Maximum-value map of {name}",
774 name=f
"{{coaddName}}Coadd_{name}_map_max",
775 storageClass=
"HealSparseMap",
776 dimensions=(
"tract",
"skymap",
"band"),
780 vars()[f
"{name}_consolidated_map_max"] = pipeBase.connectionTypes.Output(
781 doc=f
"Minumum-value map of {name}",
782 name=f
"{{coaddName}}Coadd_{name}_consolidated_map_max",
783 storageClass=
"HealSparseMap",
784 dimensions=(
"skymap",
"band"),
786 vars()[f
"{name}_map_mean"] = pipeBase.connectionTypes.Input(
787 doc=f
"Mean-value map of {name}",
788 name=f
"{{coaddName}}Coadd_{name}_map_mean",
789 storageClass=
"HealSparseMap",
790 dimensions=(
"tract",
"skymap",
"band"),
794 vars()[f
"{name}_consolidated_map_mean"] = pipeBase.connectionTypes.Output(
795 doc=f
"Minumum-value map of {name}",
796 name=f
"{{coaddName}}Coadd_{name}_consolidated_map_mean",
797 storageClass=
"HealSparseMap",
798 dimensions=(
"skymap",
"band"),
800 vars()[f
"{name}_map_weighted_mean"] = pipeBase.connectionTypes.Input(
801 doc=f
"Weighted mean-value map of {name}",
802 name=f
"{{coaddName}}Coadd_{name}_map_weighted_mean",
803 storageClass=
"HealSparseMap",
804 dimensions=(
"tract",
"skymap",
"band"),
808 vars()[f
"{name}_consolidated_map_weighted_mean"] = pipeBase.connectionTypes.Output(
809 doc=f
"Minumum-value map of {name}",
810 name=f
"{{coaddName}}Coadd_{name}_consolidated_map_weighted_mean",
811 storageClass=
"HealSparseMap",
812 dimensions=(
"skymap",
"band"),
814 vars()[f
"{name}_map_sum"] = pipeBase.connectionTypes.Input(
815 doc=f
"Sum-value map of {name}",
816 name=f
"{{coaddName}}Coadd_{name}_map_sum",
817 storageClass=
"HealSparseMap",
818 dimensions=(
"tract",
"skymap",
"band"),
822 vars()[f
"{name}_consolidated_map_sum"] = pipeBase.connectionTypes.Output(
823 doc=f
"Minumum-value map of {name}",
824 name=f
"{{coaddName}}Coadd_{name}_consolidated_map_sum",
825 storageClass=
"HealSparseMap",
826 dimensions=(
"skymap",
"band"),
829 def __init__(self, *, config=None):
830 super().__init__(config=config)
834 for name
in BasePropertyMap.registry:
835 if name
not in config.property_maps:
837 prop_config.do_min =
False
838 prop_config.do_max =
False
839 prop_config.do_mean =
False
840 prop_config.do_weighted_mean =
False
841 prop_config.do_sum =
False
843 prop_config = config.property_maps[name]
845 if not prop_config.do_min:
846 self.inputs.remove(f
"{name}_map_min")
847 self.outputs.remove(f
"{name}_consolidated_map_min")
848 if not prop_config.do_max:
849 self.inputs.remove(f
"{name}_map_max")
850 self.outputs.remove(f
"{name}_consolidated_map_max")
851 if not prop_config.do_mean:
852 self.inputs.remove(f
"{name}_map_mean")
853 self.outputs.remove(f
"{name}_consolidated_map_mean")
854 if not prop_config.do_weighted_mean:
855 self.inputs.remove(f
"{name}_map_weighted_mean")
856 self.outputs.remove(f
"{name}_consolidated_map_weighted_mean")
857 if not prop_config.do_sum:
858 self.inputs.remove(f
"{name}_map_sum")
859 self.outputs.remove(f
"{name}_consolidated_map_sum")
862class ConsolidateHealSparsePropertyMapConfig(pipeBase.PipelineTaskConfig,
863 pipelineConnections=ConsolidateHealSparsePropertyMapConnections):
864 """Configuration parameters for ConsolidateHealSparsePropertyMapTask"""
865 property_maps = BasePropertyMap.registry.makeField(
867 default=[
"exposure_time",
879 doc=
"Property map computation objects",
881 nside_coverage = pexConfig.Field(
882 doc=
"Consolidated HealSparse coverage map nside. Must be power of 2.",
885 check=_is_power_of_two,
888 def setDefaults(self):
889 self.property_maps[
"exposure_time"].do_sum =
True
890 self.property_maps[
"psf_size"].do_weighted_mean =
True
891 self.property_maps[
"psf_e1"].do_weighted_mean =
True
892 self.property_maps[
"psf_e2"].do_weighted_mean =
True
893 self.property_maps[
"psf_maglim"].do_weighted_mean =
True
894 self.property_maps[
"sky_noise"].do_weighted_mean =
True
895 self.property_maps[
"sky_background"].do_weighted_mean =
True
896 self.property_maps[
"dcr_dra"].do_weighted_mean =
True
897 self.property_maps[
"dcr_ddec"].do_weighted_mean =
True
898 self.property_maps[
"dcr_e1"].do_weighted_mean =
True
899 self.property_maps[
"dcr_e2"].do_weighted_mean =
True
900 self.property_maps[
"epoch"].do_mean =
True
901 self.property_maps[
"epoch"].do_min =
True
902 self.property_maps[
"epoch"].do_max =
True
905class ConsolidateHealSparsePropertyMapTask(pipeBase.PipelineTask):
906 """Task to consolidate HealSparse property maps.
908 This task will take all the individual tract-based maps (per map type,
909 per band) and consolidate them into one survey-wide map (per map type,
910 per band). Each tract map is truncated to its inner region before
913 ConfigClass = ConsolidateHealSparsePropertyMapConfig
914 _DefaultName =
"consolidateHealSparsePropertyMapTask"
916 def __init__(self, **kwargs):
917 super().__init__(**kwargs)
919 for name, config, PropertyMapClass
in self.config.property_maps.apply():
920 self.property_maps[name] = PropertyMapClass(config, name)
923 def runQuantum(self, butlerQC, inputRefs, outputRefs):
924 inputs = butlerQC.get(inputRefs)
926 sky_map = inputs.pop(
"sky_map")
929 for name
in self.config.property_maps.names:
930 for type_
in [
'min',
'max',
'mean',
'weighted_mean',
'sum']:
931 map_type = f
"{name}_map_{type_}"
932 if map_type
in inputs:
933 input_refs = {ref.dataId[
'tract']: ref
934 for ref
in inputs[map_type]}
935 consolidated_map = self.consolidate_map(sky_map, input_refs)
936 butlerQC.put(consolidated_map,
937 getattr(outputRefs, f
"{name}_consolidated_map_{type_}"))
939 def consolidate_map(self, sky_map, input_refs):
940 """Consolidate the healsparse property maps.
944 sky_map : Sky map object
945 input_refs : `dict` [`int`: `lsst.daf.butler.DeferredDatasetHandle`]
946 Dictionary of tract_id mapping to dataref.
950 consolidated_map : `healsparse.HealSparseMap`
951 Consolidated HealSparse map.
956 nside_coverage_inputs =
None
957 for tract_id
in input_refs:
958 cov = input_refs[tract_id].get(component=
'coverage')
960 cov_mask = cov.coverage_mask
961 nside_coverage_inputs = cov.nside_coverage
963 cov_mask |= cov.coverage_mask
965 cov_pix_inputs, = np.where(cov_mask)
968 if nside_coverage_inputs == self.config.nside_coverage:
969 cov_pix = cov_pix_inputs
970 elif nside_coverage_inputs > self.config.nside_coverage:
973 bit_shift = hsp.utils._compute_bitshift(self.config.nside_coverage,
974 nside_coverage_inputs)
975 cov_pix = np.right_shift(cov_pix_inputs, bit_shift)
979 bit_shift = hsp.utils._compute_bitshift(nside_coverage_inputs,
980 self.config.nside_coverage)
981 cov_pix = np.left_shift(cov_pix_inputs, bit_shift)
984 consolidated_map =
None
985 for tract_id
in input_refs:
986 input_map = input_refs[tract_id].get()
987 if consolidated_map
is None:
988 consolidated_map = hsp.HealSparseMap.make_empty(
989 self.config.nside_coverage,
990 input_map.nside_sparse,
992 sentinel=input_map._sentinel,
994 metadata=input_map.metadata,
998 vpix, ra, dec = input_map.valid_pixels_pos(return_pixels=
True)
999 vpix_tract_ids = sky_map.findTractIdArray(ra, dec, degrees=
True)
1001 in_tract = (vpix_tract_ids == tract_id)
1003 consolidated_map[vpix[in_tract]] = input_map[vpix[in_tract]]
1005 return consolidated_map
A floating-point coordinate rectangle geometry.
LonLat represents a spherical coordinate (longitude/latitude angle) pair.