Loading [MathJax]/extensions/tex2jax.js
LSST Applications g044012fb7c+bcde43e0cd,g04a91732dc+5ded33dca8,g07dc498a13+f701f15b83,g1409bbee79+f701f15b83,g1a7e361dbc+f701f15b83,g1fd858c14a+bb8b3b62bc,g2179890f76+b0e6738434,g35bb328faa+0eb18584fe,g3bd4b5ce2c+3bd6caf1d1,g4e0f332c67+976ceb6bc8,g53246c7159+0eb18584fe,g5477a8d5ce+51234355ef,g60b5630c4e+855785d0b2,g623d845a50+855785d0b2,g6f0c2978f1+63e17faf73,g78460c75b0+85633614c8,g786e29fd12+02b9b86fc9,g8852436030+a847daa3c3,g89139ef638+f701f15b83,g8b49a6ea8e+855785d0b2,g9125e01d80+0eb18584fe,g95236ca021+d4f98599f0,g974caa22f6+0eb18584fe,g97db14605f+39e2ace225,g989de1cb63+f701f15b83,g9f33ca652e+0886102515,gaaedd4e678+f701f15b83,gabe3b4be73+543c3c03c9,gb1101e3267+5be896e349,gb58c049af0+492386d360,gc99c83e5f0+a513197d39,gcf25f946ba+a847daa3c3,gd0fa69b896+dcb3cf4724,gd6cbbdb0b4+3fef02d88a,gde0f65d7ad+0b84c86a2a,ge278dab8ac+ae64226a64,geab183fbe5+855785d0b2,gfba249425e+0eb18584fe,w.2025.07
LSST Data Management Base Package
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
hips.py
Go to the documentation of this file.
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"""Tasks for making and manipulating HIPS images."""
23
24__all__ = ["HighResolutionHipsTask", "HighResolutionHipsConfig", "HighResolutionHipsConnections",
25 "HighResolutionHipsQuantumGraphBuilder",
26 "GenerateHipsTask", "GenerateHipsConfig", "GenerateColorHipsTask", "GenerateColorHipsConfig"]
27
28from collections import defaultdict
29import numpy as np
30import argparse
31import io
32import sys
33import re
34import warnings
35import math
36from datetime import datetime
37import hpgeom as hpg
38import healsparse as hsp
39from astropy.io import fits
40try:
41 from astropy.visualization.lupton_rgb import AsinhMapping
42except ImportError:
43 from ._fallback_asinhmapping import AsinhMapping
44from PIL import Image
45
46from lsst.sphgeom import RangeSet, HealpixPixelization
47from lsst.utils.timer import timeMethod
48from lsst.daf.butler import Butler
49from lsst.daf.butler.queries.expression_factory import ExpressionFactory
50import lsst.pex.config as pexConfig
51import lsst.pipe.base as pipeBase
52from lsst.pipe.base.quantum_graph_builder import QuantumGraphBuilder
53from lsst.pipe.base.quantum_graph_skeleton import QuantumGraphSkeleton
54import lsst.afw.geom as afwGeom
55import lsst.afw.math as afwMath
56import lsst.afw.image as afwImage
57import lsst.geom as geom
58from lsst.afw.geom import makeHpxWcs
59from lsst.resources import ResourcePath
60
61from .healSparseMapping import _is_power_of_two
62
63
64class HighResolutionHipsConnections(pipeBase.PipelineTaskConnections,
65 dimensions=("healpix9", "band"),
66 defaultTemplates={"coaddName": "deep"}):
67 coadd_exposure_handles = pipeBase.connectionTypes.Input(
68 doc="Coadded exposures to convert to HIPS format.",
69 name="{coaddName}Coadd_calexp",
70 storageClass="ExposureF",
71 dimensions=("tract", "patch", "skymap", "band"),
72 multiple=True,
73 deferLoad=True,
74 )
75 hips_exposures = pipeBase.connectionTypes.Output(
76 doc="HiPS-compatible HPX image.",
77 name="{coaddName}Coadd_hpx",
78 storageClass="ExposureF",
79 dimensions=("healpix11", "band"),
80 multiple=True,
81 )
82
83 def __init__(self, *, config=None):
84 super().__init__(config=config)
85
86 quantum_order = None
87 for dim in self.dimensions:
88 if "healpix" in dim:
89 if quantum_order is not None:
90 raise ValueError("Must not specify more than one quantum healpix dimension.")
91 quantum_order = int(dim.split("healpix")[1])
92 if quantum_order is None:
93 raise ValueError("Must specify a healpix dimension in quantum dimensions.")
94
95 if quantum_order > config.hips_order:
96 raise ValueError("Quantum healpix dimension order must not be greater than hips_order")
97
98 order = None
99 for dim in self.hips_exposures.dimensions:
100 if "healpix" in dim:
101 if order is not None:
102 raise ValueError("Must not specify more than one healpix dimension.")
103 order = int(dim.split("healpix")[1])
104 if order is None:
105 raise ValueError("Must specify a healpix dimension in hips_exposure dimensions.")
106
107 if order != config.hips_order:
108 raise ValueError("healpix dimension order must match config.hips_order.")
109
110
111class HighResolutionHipsConfig(pipeBase.PipelineTaskConfig,
112 pipelineConnections=HighResolutionHipsConnections):
113 """Configuration parameters for HighResolutionHipsTask.
114
115 Notes
116 -----
117 A HiPS image covers one HEALPix cell, with the HEALPix nside equal to
118 2**hips_order. Each cell is 'shift_order' orders deeper than the HEALPix
119 cell, with 2**shift_order x 2**shift_order sub-pixels on a side, which
120 defines the target resolution of the HiPS image. The IVOA recommends
121 shift_order=9, for 2**9=512 pixels on a side.
122
123 Table 5 from
124 https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
125 shows the relationship between hips_order, number of tiles (full
126 sky coverage), cell size, and sub-pixel size/image resolution (with
127 the default shift_order=9):
128 +------------+-----------------+--------------+------------------+
129 | hips_order | Number of Tiles | Cell Size | Image Resolution |
130 +============+=================+==============+==================+
131 | 0 | 12 | 58.63 deg | 6.871 arcmin |
132 | 1 | 48 | 29.32 deg | 3.435 arcmin |
133 | 2 | 192 | 14.66 deg | 1.718 arcmin |
134 | 3 | 768 | 7.329 deg | 51.53 arcsec |
135 | 4 | 3072 | 3.665 deg | 25.77 arcsec |
136 | 5 | 12288 | 1.832 deg | 12.88 arcsec |
137 | 6 | 49152 | 54.97 arcmin | 6.442 arcsec |
138 | 7 | 196608 | 27.48 arcmin | 3.221 arcsec |
139 | 8 | 786432 | 13.74 arcmin | 1.61 arcsec |
140 | 9 | 3145728 | 6.871 arcmin | 805.2mas |
141 | 10 | 12582912 | 3.435 arcmin | 402.6mas |
142 | 11 | 50331648 | 1.718 arcmin | 201.3mas |
143 | 12 | 201326592 | 51.53 arcsec | 100.6mas |
144 | 13 | 805306368 | 25.77 arcsec | 50.32mas |
145 +------------+-----------------+--------------+------------------+
146 """
147 hips_order = pexConfig.Field(
148 doc="HIPS image order.",
149 dtype=int,
150 default=11,
151 )
152 shift_order = pexConfig.Field(
153 doc="HIPS shift order (such that each tile is 2**shift_order pixels on a side)",
154 dtype=int,
155 default=9,
156 )
157 warp = pexConfig.ConfigField(
158 dtype=afwMath.Warper.ConfigClass,
159 doc="Warper configuration",
160 )
161
162 def setDefaults(self):
163 self.warp.warpingKernelName = "lanczos5"
164
165
166class HipsTaskNameDescriptor:
167 """Descriptor used create a DefaultName that matches the order of
168 the defined dimensions in the connections class.
169
170 Parameters
171 ----------
172 prefix : `str`
173 The prefix of the Default name, to which the order will be
174 appended.
175 """
176 def __init__(self, prefix):
177 # create a defaultName template
178 self._defaultName = f"{prefix}{{}}"
179 self._order = None
180
181 def __get__(self, obj, klass=None):
182 if klass is None:
183 raise RuntimeError(
184 "HipsTaskDescriptor was used in an unexpected context"
185 )
186 if self._order is None:
187 klassDimensions = klass.ConfigClass.ConnectionsClass.dimensions
188 for dim in klassDimensions:
189 if (match := re.match(r"^healpix(\d*)$", dim)) is not None:
190 self._order = int(match.group(1))
191 break
192 else:
193 raise RuntimeError(
194 "Could not find healpix dimension in connections class"
195 )
196 return self._defaultName.format(self._order)
197
198
199class HighResolutionHipsTask(pipeBase.PipelineTask):
200 """Task for making high resolution HiPS images."""
201 ConfigClass = HighResolutionHipsConfig
202 _DefaultName = HipsTaskNameDescriptor("highResolutionHips")
203
204 def __init__(self, **kwargs):
205 super().__init__(**kwargs)
206 self.warper = afwMath.Warper.fromConfig(self.config.warp)
207
208 @timeMethod
209 def runQuantum(self, butlerQC, inputRefs, outputRefs):
210 inputs = butlerQC.get(inputRefs)
211
212 healpix_dim = f"healpix{self.config.hips_order}"
213
214 pixels = [hips_exposure.dataId[healpix_dim]
215 for hips_exposure in outputRefs.hips_exposures]
216
217 outputs = self.run(pixels=pixels, coadd_exposure_handles=inputs["coadd_exposure_handles"])
218
219 hips_exposure_ref_dict = {hips_exposure_ref.dataId[healpix_dim]:
220 hips_exposure_ref for hips_exposure_ref in outputRefs.hips_exposures}
221 for pixel, hips_exposure in outputs.hips_exposures.items():
222 butlerQC.put(hips_exposure, hips_exposure_ref_dict[pixel])
223
224 def run(self, pixels, coadd_exposure_handles):
225 """Run the HighResolutionHipsTask.
226
227 Parameters
228 ----------
229 pixels : `Iterable` [ `int` ]
230 Iterable of healpix pixels (nest ordering) to warp to.
231 coadd_exposure_handles : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
232 Handles for the coadd exposures.
233
234 Returns
235 -------
236 outputs : `lsst.pipe.base.Struct`
237 ``hips_exposures`` is a dict with pixel (key) and hips_exposure (value)
238 """
239 self.log.info("Generating HPX images for %d pixels at order %d", len(pixels), self.config.hips_order)
240
241 npix = 2**self.config.shift_order
242 bbox_hpx = geom.Box2I(corner=geom.Point2I(0, 0),
243 dimensions=geom.Extent2I(npix, npix))
244
245 # For each healpix pixel we will create an empty exposure with the
246 # correct HPX WCS. We furthermore create a dict to hold each of
247 # the warps that will go into each HPX exposure.
248 exp_hpx_dict = {}
249 warp_dict = {}
250 for pixel in pixels:
251 wcs_hpx = afwGeom.makeHpxWcs(self.config.hips_order, pixel, shift_order=self.config.shift_order)
252 exp_hpx = afwImage.ExposureF(bbox_hpx, wcs_hpx)
253 exp_hpx_dict[pixel] = exp_hpx
254 warp_dict[pixel] = []
255
256 first_handle = True
257 # Loop over input coadd exposures to minimize i/o (this speeds things
258 # up by ~8x to batch together pixels that overlap a given coadd).
259 for handle in coadd_exposure_handles:
260 coadd_exp = handle.get()
261
262 # For each pixel, warp the coadd to the HPX WCS for the pixel.
263 for pixel in pixels:
264 warped = self.warper.warpExposure(exp_hpx_dict[pixel].getWcs(), coadd_exp, maxBBox=bbox_hpx)
265
266 exp = afwImage.ExposureF(exp_hpx_dict[pixel].getBBox(), exp_hpx_dict[pixel].getWcs())
267 exp.maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask("NO_DATA"), np.nan)
268
269 if first_handle:
270 # Make sure the mask planes, filter, and photocalib of the output
271 # exposure match the (first) input exposure.
272 exp_hpx_dict[pixel].mask.conformMaskPlanes(coadd_exp.mask.getMaskPlaneDict())
273 exp_hpx_dict[pixel].setFilter(coadd_exp.getFilter())
274 exp_hpx_dict[pixel].setPhotoCalib(coadd_exp.getPhotoCalib())
275
276 if warped.getBBox().getArea() == 0 or not np.any(np.isfinite(warped.image.array)):
277 # There is no overlap, skip.
278 self.log.debug(
279 "No overlap between output HPX %d and input exposure %s",
280 pixel,
281 handle.dataId
282 )
283 continue
284
285 exp.maskedImage.assign(warped.maskedImage, warped.getBBox())
286 warp_dict[pixel].append(exp.maskedImage)
287
288 first_handle = False
289
290 stats_flags = afwMath.stringToStatisticsProperty("MEAN")
291 stats_ctrl = afwMath.StatisticsControl()
292 stats_ctrl.setNanSafe(True)
293 stats_ctrl.setWeighted(True)
294 stats_ctrl.setCalcErrorFromInputVariance(True)
295
296 # Loop over pixels and combine the warps for each pixel.
297 # The combination is done with a simple mean for pixels that
298 # overlap in neighboring patches.
299 for pixel in pixels:
300 exp_hpx_dict[pixel].maskedImage.set(np.nan, afwImage.Mask.getPlaneBitMask("NO_DATA"), np.nan)
301
302 if not warp_dict[pixel]:
303 # Nothing in this pixel
304 self.log.debug("No data in HPX pixel %d", pixel)
305 # Remove the pixel from the output, no need to persist an
306 # empty exposure.
307 exp_hpx_dict.pop(pixel)
308 continue
309
310 exp_hpx_dict[pixel].maskedImage = afwMath.statisticsStack(
311 warp_dict[pixel],
312 stats_flags,
313 stats_ctrl,
314 [1.0]*len(warp_dict[pixel]),
315 clipped=0,
316 maskMap=[]
317 )
318
319 return pipeBase.Struct(hips_exposures=exp_hpx_dict)
320
321 @classmethod
322 def build_quantum_graph_cli(cls, argv):
323 """A command-line interface entry point to `build_quantum_graph`.
324 This method provides the implementation for the
325 ``build-high-resolution-hips-qg`` script.
326
327 Parameters
328 ----------
329 argv : `Sequence` [ `str` ]
330 Command-line arguments (e.g. ``sys.argv[1:]``).
331 """
332 parser = cls._make_cli_parser()
333
334 args = parser.parse_args(argv)
335
336 if args.subparser_name is None:
337 parser.print_help()
338 sys.exit(1)
339
340 pipeline = pipeBase.Pipeline.from_uri(args.pipeline)
341 pipeline_graph = pipeline.to_graph()
342
343 if len(pipeline_graph.tasks) != 1:
344 raise RuntimeError(f"Pipeline file {args.pipeline} may only contain one task.")
345
346 (task_node,) = pipeline_graph.tasks.values()
347
348 butler = Butler(args.butler_config, collections=args.input)
349
350 if args.subparser_name == "segment":
351 # Do the segmentation
352 hpix_pixelization = HealpixPixelization(level=args.hpix_build_order)
353 dataset = task_node.inputs["coadd_exposure_handles"].dataset_type_name
354 with butler.query() as q:
355 data_ids = list(q.join_dataset_search(dataset).data_ids("tract").with_dimension_records())
356 region_pixels = []
357 for data_id in data_ids:
358 region = data_id.region
359 pixel_range = hpix_pixelization.envelope(region)
360 for r in pixel_range.ranges():
361 region_pixels.extend(range(r[0], r[1]))
362 indices = np.unique(region_pixels)
363
364 print(f"Pixels to run at HEALPix order --hpix_build_order {args.hpix_build_order}:")
365 for pixel in indices:
366 print(pixel)
367
368 elif args.subparser_name == "build":
369 # Build the quantum graph.
370
371 # Figure out collection names.
372 if args.output_run is None:
373 if args.output is None:
374 raise ValueError("At least one of --output or --output-run options is required.")
375 args.output_run = "{}/{}".format(args.output, pipeBase.Instrument.makeCollectionTimestamp())
376
377 build_ranges = RangeSet(sorted(args.pixels))
378
379 # Metadata includes a subset of attributes defined in CmdLineFwk.
380 metadata = {
381 "input": args.input,
382 "butler_argument": args.butler_config,
383 "output": args.output,
384 "output_run": args.output_run,
385 "data_query": args.where,
386 "time": f"{datetime.now()}",
387 }
388
389 builder = HighResolutionHipsQuantumGraphBuilder(
390 pipeline_graph,
391 butler,
392 input_collections=args.input,
393 output_run=args.output_run,
394 constraint_order=args.hpix_build_order,
395 constraint_ranges=build_ranges,
396 where=args.where,
397 )
398 qg = builder.build(metadata, attach_datastore_records=True)
399 qg.saveUri(args.save_qgraph)
400
401 @classmethod
402 def _make_cli_parser(cls):
403 """Make the command-line parser.
404
405 Returns
406 -------
407 parser : `argparse.ArgumentParser`
408 """
409 parser = argparse.ArgumentParser(
410 description=(
411 "Build a QuantumGraph that runs HighResolutionHipsTask on existing coadd datasets."
412 ),
413 )
414 subparsers = parser.add_subparsers(help="sub-command help", dest="subparser_name")
415
416 parser_segment = subparsers.add_parser("segment",
417 help="Determine survey segments for workflow.")
418 parser_build = subparsers.add_parser("build",
419 help="Build quantum graph for HighResolutionHipsTask")
420
421 for sub in [parser_segment, parser_build]:
422 # These arguments are in common.
423 sub.add_argument(
424 "-b",
425 "--butler-config",
426 type=str,
427 help="Path to data repository or butler configuration.",
428 required=True,
429 )
430 sub.add_argument(
431 "-p",
432 "--pipeline",
433 type=str,
434 help="Pipeline file, limited to one task.",
435 required=True,
436 )
437 sub.add_argument(
438 "-i",
439 "--input",
440 type=str,
441 nargs="+",
442 help="Input collection(s) to search for coadd exposures.",
443 required=True,
444 )
445 sub.add_argument(
446 "-o",
447 "--hpix_build_order",
448 type=int,
449 default=1,
450 help="HEALPix order to segment sky for building quantum graph files.",
451 )
452 sub.add_argument(
453 "-w",
454 "--where",
455 type=str,
456 default="",
457 help="Data ID expression used when querying for input coadd datasets.",
458 )
459
460 parser_build.add_argument(
461 "--output",
462 type=str,
463 help=(
464 "Name of the output CHAINED collection. If this options is specified and "
465 "--output-run is not, then a new RUN collection will be created by appending "
466 "a timestamp to the value of this option."
467 ),
468 default=None,
469 metavar="COLL",
470 )
471 parser_build.add_argument(
472 "--output-run",
473 type=str,
474 help=(
475 "Output RUN collection to write resulting images. If not provided "
476 "then --output must be provided and a new RUN collection will be created "
477 "by appending a timestamp to the value passed with --output."
478 ),
479 default=None,
480 metavar="RUN",
481 )
482 parser_build.add_argument(
483 "-q",
484 "--save-qgraph",
485 type=str,
486 help="Output filename for QuantumGraph.",
487 required=True,
488 )
489 parser_build.add_argument(
490 "-P",
491 "--pixels",
492 type=int,
493 nargs="+",
494 help="Pixels at --hpix_build_order to generate quantum graph.",
495 required=True,
496 )
497
498 return parser
499
500
501class HighResolutionHipsQuantumGraphBuilder(QuantumGraphBuilder):
502 """A custom a `lsst.pipe.base.QuantumGraphBuilder` for running
503 `HighResolutionHipsTask` only.
504
505 This is a workaround for incomplete butler query support for HEALPix
506 dimensions.
507
508 Parameters
509 ----------
510 pipeline_graph : `lsst.pipe.base.PipelineGraph`
511 Pipeline graph with exactly one task, which must be a configuration
512 of `HighResolutionHipsTask`.
513 butler : `lsst.daf.butler.Butler`
514 Client for the butler data repository. May be read-only.
515 input_collections : `str` or `Iterable` [ `str` ], optional
516 Collection or collections to search for input datasets, in order.
517 If not provided, ``butler.collections`` will be searched.
518 output_run : `str`, optional
519 Name of the output collection. If not provided, ``butler.run`` will
520 be used.
521 constraint_order : `int`
522 HEALPix order used to constrain which quanta are generated, via
523 ``constraint_indices``. This should be a coarser grid (smaller
524 order) than the order used for the task's quantum and output data
525 IDs, and ideally something between the spatial scale of a patch or
526 the data repository's "common skypix" system (usually ``htm7``).
527 constraint_ranges : `lsst.sphgeom.RangeSet`
528 RangeSet that describes constraint pixels (HEALPix NEST, with order
529 ``constraint_order``) to constrain generated quanta.
530 where : `str`, optional
531 A boolean `str` expression of the form accepted by
532 `lsst.daf.butler.Butler` to constrain input datasets. This may
533 contain a constraint on tracts, patches, or bands, but not HEALPix
534 indices. Constraints on tracts and patches should usually be
535 unnecessary, however - existing coadds that overlap the given
536 HEALpix indices will be selected without such a constraint, and
537 providing one may reject some that should normally be included.
538 """
539
540 def __init__(
541 self,
542 pipeline_graph,
543 butler,
544 *,
545 input_collections=None,
546 output_run=None,
547 constraint_order,
548 constraint_ranges,
549 where="",
550 ):
551 super().__init__(pipeline_graph, butler, input_collections=input_collections, output_run=output_run)
552 self.constraint_order = constraint_order
553 self.constraint_ranges = constraint_ranges
554 self.where = where
555
556 def process_subgraph(self, subgraph):
557 # Docstring inherited.
558 (task_node,) = subgraph.tasks.values()
559
560 # Since we know this is the only task in the pipeline, we know there
561 # is only one overall input and one regular output.
562 (input_dataset_type_node,) = subgraph.inputs_of(task_node.label).values()
563 assert input_dataset_type_node is not None, "PipelineGraph should be resolved by base class."
564 (output_edge,) = task_node.outputs.values()
565 output_dataset_type_node = subgraph.dataset_types[output_edge.parent_dataset_type_name]
566 (hpx_output_dimension,) = (
567 self.butler.dimensions.skypix_dimensions[d]
568 for d in output_dataset_type_node.dimensions.skypix
569 )
570 constraint_hpx_pixelization = (
571 self.butler.dimensions.skypix_dimensions[f"healpix{self.constraint_order}"].pixelization
572 )
573 common_skypix_pixelization = self.butler.dimensions.commonSkyPix.pixelization
574
575 # We will need all the pixels at the quantum resolution as well.
576 # '4' appears here frequently because it's the number of pixels at
577 # level N in a single pixel at level (N-1).
578 (hpx_dimension,) = (
579 self.butler.dimensions.skypix_dimensions[d] for d in task_node.dimensions.names if d != "band"
580 )
581 hpx_pixelization = hpx_dimension.pixelization
582 if hpx_pixelization.level < self.constraint_order:
583 raise ValueError(f"Quantum order {hpx_pixelization.level} must be < {self.constraint_order}")
584 hpx_ranges = self.constraint_ranges.scaled(4**(hpx_pixelization.level - self.constraint_order))
585
586 # We can be generous in looking for pixels here, because we constrain
587 # by actual patch regions below.
588 common_skypix_ranges = RangeSet()
589 for begin, end in self.constraint_ranges:
590 for hpx_index in range(begin, end):
591 constraint_hpx_region = constraint_hpx_pixelization.pixel(hpx_index)
592 common_skypix_ranges |= common_skypix_pixelization.envelope(constraint_hpx_region)
593
594 # To keep the query from getting out of hand (and breaking) we simplify
595 # until we have fewer than 100 ranges which seems to work fine.
596 for simp in range(1, 10):
597 if len(common_skypix_ranges) < 100:
598 break
599 common_skypix_ranges.simplify(simp)
600
601 # Use that RangeSet to assemble a WHERE constraint expression. This
602 # could definitely get too big if the "constraint healpix" order is too
603 # fine. It's important to use `x IN (a..b)` rather than
604 # `(x >= a AND x <= b)` to avoid some pessimizations that are present
605 # in the butler expression handling (at least as of ~w_2025_05).
606 xf = ExpressionFactory(self.butler.dimensions)
607 common_skypix_proxy = getattr(xf, self.butler.dimensions.commonSkyPix.name)
608 where_terms = []
609 for begin, end in common_skypix_ranges:
610 if begin == end - 1:
611 where_terms.append(common_skypix_proxy == begin)
612 else:
613 where_terms.append(common_skypix_proxy.in_range(begin, end))
614 where = xf.any(*where_terms)
615 # Query for input datasets with this constraint, and ask for expanded
616 # data IDs because we want regions. Immediately group this by patch so
617 # we don't do later geometric stuff n_bands more times than we need to.
618 with self.butler.query() as query:
619 input_refs = query.datasets(
620 input_dataset_type_node.dataset_type,
621 collections=self.input_collections
622 ).where(
623 where,
624 self.where,
625 ).with_dimension_records()
626 inputs_by_patch = defaultdict(set)
627 patch_dimensions = self.butler.dimensions.conform(["patch"])
628 skeleton = QuantumGraphSkeleton([task_node.label])
629 for input_ref in input_refs:
630 dataset_key = skeleton.add_dataset_node(input_ref.datasetType.name, input_ref.dataId)
631 skeleton.set_dataset_ref(input_ref, dataset_key)
632 inputs_by_patch[input_ref.dataId.subset(patch_dimensions)].add(dataset_key)
633 if not inputs_by_patch:
634 message_body = "\n".join(input_refs.explain_no_results())
635 raise RuntimeError(f"No inputs found:\n{message_body}")
636
637 # Iterate over patches and compute the set of output healpix pixels
638 # that overlap each one. Use that to associate inputs with output
639 # pixels, but only for the output pixels we've already identified.
640 inputs_by_hpx = defaultdict(set)
641 for patch_data_id, input_keys_for_patch in inputs_by_patch.items():
642 patch_hpx_ranges = hpx_pixelization.envelope(patch_data_id.region)
643 for begin, end in patch_hpx_ranges & hpx_ranges:
644 for hpx_index in range(begin, end):
645 inputs_by_hpx[hpx_index].update(input_keys_for_patch)
646
647 # Iterate over the dict we just created and create preliminary quanta.
648 for hpx_index, input_keys_for_hpx_index in inputs_by_hpx.items():
649 # Group inputs by band.
650 input_keys_by_band = defaultdict(list)
651 for input_key in input_keys_for_hpx_index:
652 input_ref = skeleton.get_dataset_ref(input_key)
653 assert input_ref is not None, "Code above adds the same nodes to the graph with refs."
654 input_keys_by_band[input_ref.dataId["band"]].append(input_key)
655 # Iterate over bands to make quanta.
656 for band, input_keys_for_band in input_keys_by_band.items():
657 data_id = self.butler.registry.expandDataId({hpx_dimension.name: hpx_index, "band": band})
658 quantum_key = skeleton.add_quantum_node(task_node.label, data_id)
659 # Add inputs to the skelton
660 skeleton.add_input_edges(quantum_key, input_keys_for_band)
661 # Add the regular outputs.
662 hpx_pixel_ranges = RangeSet(hpx_index)
663 hpx_output_ranges = hpx_pixel_ranges.scaled(
664 4**(task_node.config.hips_order - hpx_pixelization.level)
665 )
666 for begin, end in hpx_output_ranges:
667 for hpx_output_index in range(begin, end):
668 dataset_key = skeleton.add_dataset_node(
669 output_dataset_type_node.name,
670 self.butler.registry.expandDataId(
671 {hpx_output_dimension: hpx_output_index, "band": band}
672 ),
673 )
674 skeleton.add_output_edge(quantum_key, dataset_key)
675 # Add auxiliary outputs (log, metadata).
676 for write_edge in task_node.iter_all_outputs():
677 if write_edge.connection_name == output_edge.connection_name:
678 continue
679 dataset_key = skeleton.add_dataset_node(write_edge.parent_dataset_type_name, data_id)
680 skeleton.add_output_edge(quantum_key, dataset_key)
681 return skeleton
682
683
684class HipsPropertiesSpectralTerm(pexConfig.Config):
685 lambda_min = pexConfig.Field(
686 doc="Minimum wavelength (nm)",
687 dtype=float,
688 )
689 lambda_max = pexConfig.Field(
690 doc="Maximum wavelength (nm)",
691 dtype=float,
692 )
693
694
695class HipsPropertiesConfig(pexConfig.Config):
696 """Configuration parameters for writing a HiPS properties file."""
697 creator_did_template = pexConfig.Field(
698 doc=("Unique identifier of the HiPS - Format: IVOID. "
699 "Use ``{band}`` to substitute the band name."),
700 dtype=str,
701 optional=False,
702 )
703 obs_collection = pexConfig.Field(
704 doc="Short name of original data set - Format: one word",
705 dtype=str,
706 optional=True,
707 )
708 obs_description_template = pexConfig.Field(
709 doc=("Data set description - Format: free text, longer free text "
710 "description of the dataset. Use ``{band}`` to substitute "
711 "the band name."),
712 dtype=str,
713 )
714 prov_progenitor = pexConfig.ListField(
715 doc="Provenance of the original data - Format: free text",
716 dtype=str,
717 default=[],
718 )
719 obs_title_template = pexConfig.Field(
720 doc=("Data set title format: free text, but should be short. "
721 "Use ``{band}`` to substitute the band name."),
722 dtype=str,
723 optional=False,
724 )
725 spectral_ranges = pexConfig.ConfigDictField(
726 doc=("Mapping from band to lambda_min, lamba_max (nm). May be approximate."),
727 keytype=str,
728 itemtype=HipsPropertiesSpectralTerm,
729 default={},
730 )
731 initial_ra = pexConfig.Field(
732 doc="Initial RA (deg) (default for HiPS viewer). If not set will use a point in MOC.",
733 dtype=float,
734 optional=True,
735 )
736 initial_dec = pexConfig.Field(
737 doc="Initial Declination (deg) (default for HiPS viewer). If not set will use a point in MOC.",
738 dtype=float,
739 optional=True,
740 )
741 initial_fov = pexConfig.Field(
742 doc="Initial field-of-view (deg). If not set will use ~1 healpix tile.",
743 dtype=float,
744 optional=True,
745 )
746 obs_ack = pexConfig.Field(
747 doc="Observation acknowledgements (free text).",
748 dtype=str,
749 optional=True,
750 )
751 t_min = pexConfig.Field(
752 doc="Time (MJD) of earliest observation included in HiPS",
753 dtype=float,
754 optional=True,
755 )
756 t_max = pexConfig.Field(
757 doc="Time (MJD) of latest observation included in HiPS",
758 dtype=float,
759 optional=True,
760 )
761
762 def validate(self):
763 super().validate()
764
765 if self.obs_collection is not None:
766 if re.search(r"\s", self.obs_collection):
767 raise ValueError("obs_collection cannot contain any space characters.")
768
769 def setDefaults(self):
770 # Values here taken from
771 # https://github.com/lsst-dm/dax_obscore/blob/44ac15029136e2ec15/configs/dp02.yaml#L46
772 u_term = HipsPropertiesSpectralTerm()
773 u_term.lambda_min = 330.
774 u_term.lambda_max = 400.
775 self.spectral_ranges["u"] = u_term
776 g_term = HipsPropertiesSpectralTerm()
777 g_term.lambda_min = 402.
778 g_term.lambda_max = 552.
779 self.spectral_ranges["g"] = g_term
780 r_term = HipsPropertiesSpectralTerm()
781 r_term.lambda_min = 552.
782 r_term.lambda_max = 691.
783 self.spectral_ranges["r"] = r_term
784 i_term = HipsPropertiesSpectralTerm()
785 i_term.lambda_min = 691.
786 i_term.lambda_max = 818.
787 self.spectral_ranges["i"] = i_term
788 z_term = HipsPropertiesSpectralTerm()
789 z_term.lambda_min = 818.
790 z_term.lambda_max = 922.
791 self.spectral_ranges["z"] = z_term
792 y_term = HipsPropertiesSpectralTerm()
793 y_term.lambda_min = 970.
794 y_term.lambda_max = 1060.
795 self.spectral_ranges["y"] = y_term
796
797
798class GenerateHipsConnections(pipeBase.PipelineTaskConnections,
799 dimensions=("instrument", "band"),
800 defaultTemplates={"coaddName": "deep"}):
801 hips_exposure_handles = pipeBase.connectionTypes.Input(
802 doc="HiPS-compatible HPX images.",
803 name="{coaddName}Coadd_hpx",
804 storageClass="ExposureF",
805 dimensions=("healpix11", "band"),
806 multiple=True,
807 deferLoad=True,
808 )
809
810
811class GenerateHipsConfig(pipeBase.PipelineTaskConfig,
812 pipelineConnections=GenerateHipsConnections):
813 """Configuration parameters for GenerateHipsTask."""
814 # WARNING: In general PipelineTasks are not allowed to do any outputs
815 # outside of the butler. This task has been given (temporary)
816 # Special Dispensation because of the nature of HiPS outputs until
817 # a more controlled solution can be found.
818 hips_base_uri = pexConfig.Field(
819 doc="URI to HiPS base for output.",
820 dtype=str,
821 optional=False,
822 )
823 min_order = pexConfig.Field(
824 doc="Minimum healpix order for HiPS tree.",
825 dtype=int,
826 default=3,
827 )
828 properties = pexConfig.ConfigField(
829 dtype=HipsPropertiesConfig,
830 doc="Configuration for properties file.",
831 )
832 allsky_tilesize = pexConfig.Field(
833 dtype=int,
834 doc="Allsky tile size; must be power of 2. HiPS standard recommends 64x64 tiles.",
835 default=64,
836 check=_is_power_of_two,
837 )
838 png_gray_asinh_minimum = pexConfig.Field(
839 doc="AsinhMapping intensity to be mapped to black for grayscale png scaling.",
840 dtype=float,
841 default=0.0,
842 )
843 png_gray_asinh_stretch = pexConfig.Field(
844 doc="AsinhMapping linear stretch for grayscale png scaling.",
845 dtype=float,
846 default=2.0,
847 )
848 png_gray_asinh_softening = pexConfig.Field(
849 doc="AsinhMapping softening parameter (Q) for grayscale png scaling.",
850 dtype=float,
851 default=8.0,
852 )
853
854
855class GenerateHipsTask(pipeBase.PipelineTask):
856 """Task for making a HiPS tree with FITS and grayscale PNGs."""
857 ConfigClass = GenerateHipsConfig
858 _DefaultName = "generateHips"
859 color_task = False
860
861 @timeMethod
862 def runQuantum(self, butlerQC, inputRefs, outputRefs):
863 inputs = butlerQC.get(inputRefs)
864
865 dims = inputRefs.hips_exposure_handles[0].dataId.dimensions.names
866 order = None
867 for dim in dims:
868 if "healpix" in dim:
869 order = int(dim.split("healpix")[1])
870 healpix_dim = dim
871 break
872 else:
873 raise RuntimeError("Could not determine healpix order for input exposures.")
874
875 hips_exposure_handle_dict = {
876 (hips_exposure_handle.dataId[healpix_dim],
877 hips_exposure_handle.dataId["band"]): hips_exposure_handle
878 for hips_exposure_handle in inputs["hips_exposure_handles"]
879 }
880
881 data_bands = {hips_exposure_handle.dataId["band"]
882 for hips_exposure_handle in inputs["hips_exposure_handles"]}
883 bands = self._check_data_bands(data_bands)
884
885 self.run(
886 bands=bands,
887 max_order=order,
888 hips_exposure_handle_dict=hips_exposure_handle_dict,
889 do_color=self.color_task,
890 )
891
892 def _check_data_bands(self, data_bands):
893 """Check that the data has only a single band.
894
895 Parameters
896 ----------
897 data_bands : `set` [`str`]
898 Bands from the input data.
899
900 Returns
901 -------
902 bands : `list` [`str`]
903 List of single band to process.
904
905 Raises
906 ------
907 RuntimeError if there is not exactly one band.
908 """
909 if len(data_bands) != 1:
910 raise RuntimeError("GenerateHipsTask can only use data from a single band.")
911
912 return list(data_bands)
913
914 @timeMethod
915 def run(self, bands, max_order, hips_exposure_handle_dict, do_color=False):
916 """Run the GenerateHipsTask.
917
918 Parameters
919 ----------
920 bands : `list [ `str` ]
921 List of bands to be processed (or single band).
922 max_order : `int`
923 HEALPix order of the maximum (native) HPX exposures.
924 hips_exposure_handle_dict : `dict` {`int`: `lsst.daf.butler.DeferredDatasetHandle`}
925 Dict of handles for the HiPS high-resolution exposures.
926 Key is (pixel number, ``band``).
927 do_color : `bool`, optional
928 Do color pngs instead of per-band grayscale.
929 """
930 min_order = self.config.min_order
931
932 if not do_color:
933 png_grayscale_mapping = AsinhMapping(
934 self.config.png_gray_asinh_minimum,
935 self.config.png_gray_asinh_stretch,
936 Q=self.config.png_gray_asinh_softening,
937 )
938 else:
939 png_color_mapping = AsinhMapping(
940 self.config.png_color_asinh_minimum,
941 self.config.png_color_asinh_stretch,
942 Q=self.config.png_color_asinh_softening,
943 )
944
945 bcb = self.config.blue_channel_band
946 gcb = self.config.green_channel_band
947 rcb = self.config.red_channel_band
948 colorstr = f"{bcb}{gcb}{rcb}"
949
950 # The base path is based on the hips_base_uri.
951 hips_base_path = ResourcePath(self.config.hips_base_uri, forceDirectory=True)
952
953 # We need to unique-ify the pixels because they show up for multiple bands.
954 # The output of this is a sorted array.
955 pixels = np.unique(np.array([pixel for pixel, _ in hips_exposure_handle_dict.keys()]))
956
957 # Add a "gutter" pixel at the end. Start with 0 which maps to 0 always.
958 pixels = np.append(pixels, [0])
959
960 # Convert the pixels to each order that will be generated.
961 pixels_shifted = {}
962 pixels_shifted[max_order] = pixels
963 for order in range(max_order - 1, min_order - 1, -1):
964 pixels_shifted[order] = np.right_shift(pixels_shifted[order + 1], 2)
965
966 # And set the gutter to an illegal pixel value.
967 for order in range(min_order, max_order + 1):
968 pixels_shifted[order][-1] = -1
969
970 # Read in the first pixel for determining image properties.
971 exp0 = list(hips_exposure_handle_dict.values())[0].get()
972 bbox = exp0.getBBox()
973 npix = bbox.getWidth()
974 shift_order = int(np.round(np.log2(npix)))
975
976 # Create blank exposures for each level, including the highest order.
977 # We also make sure we create blank exposures for any bands used in the color
978 # PNGs, even if they aren't available.
979 exposures = {}
980 for band in bands:
981 for order in range(min_order, max_order + 1):
982 exp = exp0.Factory(bbox=bbox)
983 exp.image.array[:, :] = np.nan
984 exposures[(band, order)] = exp
985
986 # Loop over all pixels, avoiding the gutter.
987 for pixel_counter, pixel in enumerate(pixels[:-1]):
988 self.log.debug("Working on high resolution pixel %d", pixel)
989 for band in bands:
990 # Read all the exposures here for the highest order.
991 # There will always be at least one band with a HiPS image available
992 # at the highest order. However, for color images it is possible that
993 # not all bands have coverage so we require this check.
994 if (pixel, band) in hips_exposure_handle_dict:
995 exposures[(band, max_order)] = hips_exposure_handle_dict[(pixel, band)].get()
996
997 # Go up the HiPS tree.
998 # We only write pixels and rebin to fill the parent pixel when we are
999 # done with a current pixel, which is determined if the next pixel
1000 # has a different pixel number.
1001 for order in range(max_order, min_order - 1, -1):
1002 if pixels_shifted[order][pixel_counter + 1] == pixels_shifted[order][pixel_counter]:
1003 # This order is not done, and so none of the other orders will be.
1004 break
1005
1006 # We can now write out the images for each band.
1007 # Note this will always trigger at the max order where each pixel is unique.
1008 if not do_color:
1009 for band in bands:
1010 self._write_hips_image(
1011 hips_base_path.join(f"band_{band}", forceDirectory=True),
1012 order,
1013 pixels_shifted[order][pixel_counter],
1014 exposures[(band, order)].image,
1015 png_grayscale_mapping,
1016 shift_order=shift_order,
1017 )
1018 else:
1019 # Make a color png.
1020 self._write_hips_color_png(
1021 hips_base_path.join(f"color_{colorstr}", forceDirectory=True),
1022 order,
1023 pixels_shifted[order][pixel_counter],
1024 exposures[(self.config.red_channel_band, order)].image,
1025 exposures[(self.config.green_channel_band, order)].image,
1026 exposures[(self.config.blue_channel_band, order)].image,
1027 png_color_mapping,
1028 )
1029
1030 log_level = self.log.INFO if order == (max_order - 3) else self.log.DEBUG
1031 self.log.log(
1032 log_level,
1033 "Completed HiPS generation for %s, order %d, pixel %d (%d/%d)",
1034 ",".join(bands),
1035 order,
1036 pixels_shifted[order][pixel_counter],
1037 pixel_counter,
1038 len(pixels) - 1,
1039 )
1040
1041 # When we are at the top of the tree, erase top level images and continue.
1042 if order == min_order:
1043 for band in bands:
1044 exposures[(band, order)].image.array[:, :] = np.nan
1045 continue
1046
1047 # Now average the images for each band.
1048 for band in bands:
1049 arr = exposures[(band, order)].image.array.reshape(npix//2, 2, npix//2, 2)
1050 with warnings.catch_warnings():
1051 warnings.simplefilter("ignore")
1052 binned_image_arr = np.nanmean(arr, axis=(1, 3))
1053
1054 # Fill the next level up. We figure out which of the four
1055 # sub-pixels the current pixel occupies.
1056 sub_index = (pixels_shifted[order][pixel_counter]
1057 - np.left_shift(pixels_shifted[order - 1][pixel_counter], 2))
1058
1059 # Fill exposure at the next level up.
1060 exp = exposures[(band, order - 1)]
1061
1062 # Fill the correct subregion.
1063 if sub_index == 0:
1064 exp.image.array[npix//2:, 0: npix//2] = binned_image_arr
1065 elif sub_index == 1:
1066 exp.image.array[0: npix//2, 0: npix//2] = binned_image_arr
1067 elif sub_index == 2:
1068 exp.image.array[npix//2:, npix//2:] = binned_image_arr
1069 elif sub_index == 3:
1070 exp.image.array[0: npix//2, npix//2:] = binned_image_arr
1071 else:
1072 # This should be impossible.
1073 raise ValueError("Illegal pixel sub index")
1074
1075 # Erase the previous exposure.
1076 if order < max_order:
1077 exposures[(band, order)].image.array[:, :] = np.nan
1078
1079 # Write the properties files and MOCs.
1080 if not do_color:
1081 for band in bands:
1082 band_pixels = np.array([pixel
1083 for pixel, band_ in hips_exposure_handle_dict.keys()
1084 if band_ == band])
1085 band_pixels = np.sort(band_pixels)
1086
1087 self._write_properties_and_moc(
1088 hips_base_path.join(f"band_{band}", forceDirectory=True),
1089 max_order,
1090 band_pixels,
1091 exp0,
1092 shift_order,
1093 band,
1094 False,
1095 )
1096 self._write_allsky_file(
1097 hips_base_path.join(f"band_{band}", forceDirectory=True),
1098 min_order,
1099 )
1100 else:
1101 self._write_properties_and_moc(
1102 hips_base_path.join(f"color_{colorstr}", forceDirectory=True),
1103 max_order,
1104 pixels[:-1],
1105 exp0,
1106 shift_order,
1107 colorstr,
1108 True,
1109 )
1110 self._write_allsky_file(
1111 hips_base_path.join(f"color_{colorstr}", forceDirectory=True),
1112 min_order,
1113 )
1114
1115 def _write_hips_image(self, hips_base_path, order, pixel, image, png_mapping, shift_order=9):
1116 """Write a HiPS image.
1117
1118 Parameters
1119 ----------
1120 hips_base_path : `lsst.resources.ResourcePath`
1121 Resource path to the base of the HiPS directory tree.
1122 order : `int`
1123 HEALPix order of the HiPS image to write.
1124 pixel : `int`
1125 HEALPix pixel of the HiPS image.
1126 image : `lsst.afw.image.Image`
1127 Image to write.
1128 png_mapping : `astropy.visualization.lupton_rgb.AsinhMapping`
1129 Mapping to convert image to scaled png.
1130 shift_order : `int`, optional
1131 HPX shift_order.
1132 """
1133 # WARNING: In general PipelineTasks are not allowed to do any outputs
1134 # outside of the butler. This task has been given (temporary)
1135 # Special Dispensation because of the nature of HiPS outputs until
1136 # a more controlled solution can be found.
1137
1138 dir_number = self._get_dir_number(pixel)
1139 hips_dir = hips_base_path.join(
1140 f"Norder{order}",
1141 forceDirectory=True
1142 ).join(
1143 f"Dir{dir_number}",
1144 forceDirectory=True
1145 )
1146
1147 wcs = makeHpxWcs(order, pixel, shift_order=shift_order)
1148
1149 uri = hips_dir.join(f"Npix{pixel}.fits")
1150
1151 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1152 image.writeFits(temporary_uri.ospath, metadata=wcs.getFitsMetadata())
1153
1154 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1155
1156 # And make a grayscale png as well
1157
1158 with np.errstate(invalid="ignore"):
1159 vals = 255 - png_mapping.map_intensity_to_uint8(image.array).astype(np.uint8)
1160
1161 vals[~np.isfinite(image.array) | (image.array < 0)] = 0
1162 im = Image.fromarray(vals[::-1, :], "L")
1163
1164 uri = hips_dir.join(f"Npix{pixel}.png")
1165
1166 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1167 im.save(temporary_uri.ospath)
1168
1169 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1170
1171 def _write_hips_color_png(
1172 self,
1173 hips_base_path,
1174 order,
1175 pixel,
1176 image_red,
1177 image_green,
1178 image_blue,
1179 png_mapping,
1180 ):
1181 """Write a color png HiPS image.
1182
1183 Parameters
1184 ----------
1185 hips_base_path : `lsst.resources.ResourcePath`
1186 Resource path to the base of the HiPS directory tree.
1187 order : `int`
1188 HEALPix order of the HiPS image to write.
1189 pixel : `int`
1190 HEALPix pixel of the HiPS image.
1191 image_red : `lsst.afw.image.Image`
1192 Input for red channel of output png.
1193 image_green : `lsst.afw.image.Image`
1194 Input for green channel of output png.
1195 image_blue : `lsst.afw.image.Image`
1196 Input for blue channel of output png.
1197 png_mapping : `astropy.visualization.lupton_rgb.AsinhMapping`
1198 Mapping to convert image to scaled png.
1199 """
1200 # WARNING: In general PipelineTasks are not allowed to do any outputs
1201 # outside of the butler. This task has been given (temporary)
1202 # Special Dispensation because of the nature of HiPS outputs until
1203 # a more controlled solution can be found.
1204
1205 dir_number = self._get_dir_number(pixel)
1206 hips_dir = hips_base_path.join(
1207 f"Norder{order}",
1208 forceDirectory=True
1209 ).join(
1210 f"Dir{dir_number}",
1211 forceDirectory=True
1212 )
1213
1214 # We need to convert nans to the minimum values in the mapping.
1215 arr_red = image_red.array.copy()
1216 arr_red[np.isnan(arr_red)] = png_mapping.minimum[0]
1217 arr_green = image_green.array.copy()
1218 arr_green[np.isnan(arr_green)] = png_mapping.minimum[1]
1219 arr_blue = image_blue.array.copy()
1220 arr_blue[np.isnan(arr_blue)] = png_mapping.minimum[2]
1221
1222 image_array = png_mapping.make_rgb_image(arr_red, arr_green, arr_blue)
1223
1224 im = Image.fromarray(image_array[::-1, :, :], mode="RGB")
1225
1226 uri = hips_dir.join(f"Npix{pixel}.png")
1227
1228 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1229 im.save(temporary_uri.ospath)
1230
1231 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1232
1233 def _write_properties_and_moc(
1234 self,
1235 hips_base_path,
1236 max_order,
1237 pixels,
1238 exposure,
1239 shift_order,
1240 band,
1241 multiband
1242 ):
1243 """Write HiPS properties file and MOC.
1244
1245 Parameters
1246 ----------
1247 hips_base_path : : `lsst.resources.ResourcePath`
1248 Resource path to the base of the HiPS directory tree.
1249 max_order : `int`
1250 Maximum HEALPix order.
1251 pixels : `np.ndarray` (N,)
1252 Array of pixels used.
1253 exposure : `lsst.afw.image.Exposure`
1254 Sample HPX exposure used for generating HiPS tiles.
1255 shift_order : `int`
1256 HPX shift order.
1257 band : `str`
1258 Band (or color).
1259 multiband : `bool`
1260 Is band multiband / color?
1261 """
1262 area = hpg.nside_to_pixel_area(2**max_order, degrees=True)*len(pixels)
1263
1264 initial_ra = self.config.properties.initial_ra
1265 initial_dec = self.config.properties.initial_dec
1266 initial_fov = self.config.properties.initial_fov
1267
1268 if initial_ra is None or initial_dec is None or initial_fov is None:
1269 # We want to point to an arbitrary pixel in the footprint.
1270 # Just take the median pixel value for simplicity.
1271 temp_pixels = pixels.copy()
1272 if temp_pixels.size % 2 == 0:
1273 temp_pixels = np.append(temp_pixels, [temp_pixels[0]])
1274 medpix = int(np.median(temp_pixels))
1275 _initial_ra, _initial_dec = hpg.pixel_to_angle(2**max_order, medpix)
1276 _initial_fov = hpg.nside_to_resolution(2**max_order, units='arcminutes')/60.
1277
1278 if initial_ra is None or initial_dec is None:
1279 initial_ra = _initial_ra
1280 initial_dec = _initial_dec
1281 if initial_fov is None:
1282 initial_fov = _initial_fov
1283
1284 self._write_hips_properties_file(
1285 hips_base_path,
1286 self.config.properties,
1287 band,
1288 multiband,
1289 exposure,
1290 max_order,
1291 shift_order,
1292 area,
1293 initial_ra,
1294 initial_dec,
1295 initial_fov,
1296 )
1297
1298 # Write the MOC coverage
1299 self._write_hips_moc_file(
1300 hips_base_path,
1301 max_order,
1302 pixels,
1303 )
1304
1305 def _write_hips_properties_file(
1306 self,
1307 hips_base_path,
1308 properties_config,
1309 band,
1310 multiband,
1311 exposure,
1312 max_order,
1313 shift_order,
1314 area,
1315 initial_ra,
1316 initial_dec,
1317 initial_fov
1318 ):
1319 """Write HiPS properties file.
1320
1321 Parameters
1322 ----------
1323 hips_base_path : `lsst.resources.ResourcePath`
1324 ResourcePath at top of HiPS tree. File will be written
1325 to this path as ``properties``.
1326 properties_config : `lsst.pipe.tasks.hips.HipsPropertiesConfig`
1327 Configuration for properties values.
1328 band : `str`
1329 Name of band(s) for HiPS tree.
1330 multiband : `bool`
1331 Is multiband / color?
1332 exposure : `lsst.afw.image.Exposure`
1333 Sample HPX exposure used for generating HiPS tiles.
1334 max_order : `int`
1335 Maximum HEALPix order.
1336 shift_order : `int`
1337 HPX shift order.
1338 area : `float`
1339 Coverage area in square degrees.
1340 initial_ra : `float`
1341 Initial HiPS RA position (degrees).
1342 initial_dec : `float`
1343 Initial HiPS Dec position (degrees).
1344 initial_fov : `float`
1345 Initial HiPS display size (degrees).
1346 """
1347 # WARNING: In general PipelineTasks are not allowed to do any outputs
1348 # outside of the butler. This task has been given (temporary)
1349 # Special Dispensation because of the nature of HiPS outputs until
1350 # a more controlled solution can be found.
1351 def _write_property(fh, name, value):
1352 """Write a property name/value to a file handle.
1353
1354 Parameters
1355 ----------
1356 fh : file handle (blah)
1357 Open for writing.
1358 name : `str`
1359 Name of property
1360 value : `str`
1361 Value of property
1362 """
1363 # This ensures that the name has no spaces or space-like characters,
1364 # per the HiPS standard.
1365 if re.search(r"\s", name):
1366 raise ValueError(f"``{name}`` cannot contain any space characters.")
1367 if "=" in name:
1368 raise ValueError(f"``{name}`` cannot contain an ``=``")
1369
1370 fh.write(f"{name:25}= {value}\n")
1371
1372 if exposure.image.array.dtype == np.dtype("float32"):
1373 bitpix = -32
1374 elif exposure.image.array.dtype == np.dtype("float64"):
1375 bitpix = -64
1376 elif exposure.image.array.dtype == np.dtype("int32"):
1377 bitpix = 32
1378
1379 date_iso8601 = datetime.utcnow().isoformat(timespec="seconds") + "Z"
1380 pixel_scale = hpg.nside_to_resolution(2**(max_order + shift_order), units='degrees')
1381
1382 uri = hips_base_path.join("properties")
1383 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1384 with open(temporary_uri.ospath, "w") as fh:
1385 _write_property(
1386 fh,
1387 "creator_did",
1388 properties_config.creator_did_template.format(band=band),
1389 )
1390 if properties_config.obs_collection is not None:
1391 _write_property(fh, "obs_collection", properties_config.obs_collection)
1392 _write_property(
1393 fh,
1394 "obs_title",
1395 properties_config.obs_title_template.format(band=band),
1396 )
1397 if properties_config.obs_description_template is not None:
1398 _write_property(
1399 fh,
1400 "obs_description",
1401 properties_config.obs_description_template.format(band=band),
1402 )
1403 if len(properties_config.prov_progenitor) > 0:
1404 for prov_progenitor in properties_config.prov_progenitor:
1405 _write_property(fh, "prov_progenitor", prov_progenitor)
1406 if properties_config.obs_ack is not None:
1407 _write_property(fh, "obs_ack", properties_config.obs_ack)
1408 _write_property(fh, "obs_regime", "Optical")
1409 _write_property(fh, "data_pixel_bitpix", str(bitpix))
1410 _write_property(fh, "dataproduct_type", "image")
1411 _write_property(fh, "moc_sky_fraction", str(area/41253.))
1412 _write_property(fh, "data_ucd", "phot.flux")
1413 _write_property(fh, "hips_creation_date", date_iso8601)
1414 _write_property(fh, "hips_builder", "lsst.pipe.tasks.hips.GenerateHipsTask")
1415 _write_property(fh, "hips_creator", "Vera C. Rubin Observatory")
1416 _write_property(fh, "hips_version", "1.4")
1417 _write_property(fh, "hips_release_date", date_iso8601)
1418 _write_property(fh, "hips_frame", "equatorial")
1419 _write_property(fh, "hips_order", str(max_order))
1420 _write_property(fh, "hips_tile_width", str(exposure.getBBox().getWidth()))
1421 _write_property(fh, "hips_status", "private master clonableOnce")
1422 if multiband:
1423 _write_property(fh, "hips_tile_format", "png")
1424 _write_property(fh, "dataproduct_subtype", "color")
1425 else:
1426 _write_property(fh, "hips_tile_format", "png fits")
1427 _write_property(fh, "hips_pixel_bitpix", str(bitpix))
1428 _write_property(fh, "hips_pixel_scale", str(pixel_scale))
1429 _write_property(fh, "hips_initial_ra", str(initial_ra))
1430 _write_property(fh, "hips_initial_dec", str(initial_dec))
1431 _write_property(fh, "hips_initial_fov", str(initial_fov))
1432 if multiband:
1433 if self.config.blue_channel_band in properties_config.spectral_ranges:
1434 em_min = properties_config.spectral_ranges[
1435 self.config.blue_channel_band
1436 ].lambda_min/1e9
1437 else:
1438 self.log.warning("blue band %s not in self.config.spectral_ranges.", band)
1439 em_min = 3e-7
1440 if self.config.red_channel_band in properties_config.spectral_ranges:
1441 em_max = properties_config.spectral_ranges[
1442 self.config.red_channel_band
1443 ].lambda_max/1e9
1444 else:
1445 self.log.warning("red band %s not in self.config.spectral_ranges.", band)
1446 em_max = 1e-6
1447 else:
1448 if band in properties_config.spectral_ranges:
1449 em_min = properties_config.spectral_ranges[band].lambda_min/1e9
1450 em_max = properties_config.spectral_ranges[band].lambda_max/1e9
1451 else:
1452 self.log.warning("band %s not in self.config.spectral_ranges.", band)
1453 em_min = 3e-7
1454 em_max = 1e-6
1455 _write_property(fh, "em_min", str(em_min))
1456 _write_property(fh, "em_max", str(em_max))
1457 if properties_config.t_min is not None:
1458 _write_property(fh, "t_min", properties_config.t_min)
1459 if properties_config.t_max is not None:
1460 _write_property(fh, "t_max", properties_config.t_max)
1461
1462 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1463
1464 def _write_hips_moc_file(self, hips_base_path, max_order, pixels, min_uniq_order=1):
1465 """Write HiPS MOC file.
1466
1467 Parameters
1468 ----------
1469 hips_base_path : `lsst.resources.ResourcePath`
1470 ResourcePath to top of HiPS tree. File will be written as
1471 to this path as ``Moc.fits``.
1472 max_order : `int`
1473 Maximum HEALPix order.
1474 pixels : `np.ndarray`
1475 Array of pixels covered.
1476 min_uniq_order : `int`, optional
1477 Minimum HEALPix order for looking for fully covered pixels.
1478 """
1479 # WARNING: In general PipelineTasks are not allowed to do any outputs
1480 # outside of the butler. This task has been given (temporary)
1481 # Special Dispensation because of the nature of HiPS outputs until
1482 # a more controlled solution can be found.
1483
1484 # Make the initial list of UNIQ pixels
1485 uniq = 4*(4**max_order) + pixels
1486
1487 # Make a healsparse map which provides easy degrade/comparisons.
1488 hspmap = hsp.HealSparseMap.make_empty(2**min_uniq_order, 2**max_order, dtype=np.float32)
1489 hspmap[pixels] = 1.0
1490
1491 # Loop over orders, degrade each time, and look for pixels with full coverage.
1492 for uniq_order in range(max_order - 1, min_uniq_order - 1, -1):
1493 hspmap = hspmap.degrade(2**uniq_order, reduction="sum")
1494 pix_shift = np.right_shift(pixels, 2*(max_order - uniq_order))
1495 # Check if any of the pixels at uniq_order have full coverage.
1496 covered, = np.isclose(hspmap[pix_shift], 4**(max_order - uniq_order)).nonzero()
1497 if covered.size == 0:
1498 # No pixels at uniq_order are fully covered, we're done.
1499 break
1500 # Replace the UNIQ pixels that are fully covered.
1501 uniq[covered] = 4*(4**uniq_order) + pix_shift[covered]
1502
1503 # Remove duplicate pixels.
1504 uniq = np.unique(uniq)
1505
1506 # Output to fits.
1507 tbl = np.zeros(uniq.size, dtype=[("UNIQ", "i8")])
1508 tbl["UNIQ"] = uniq
1509
1510 order = np.log2(tbl["UNIQ"]//4).astype(np.int32)//2
1511 moc_order = np.max(order)
1512
1513 hdu = fits.BinTableHDU(tbl)
1514 hdu.header["PIXTYPE"] = "HEALPIX"
1515 hdu.header["ORDERING"] = "NUNIQ"
1516 hdu.header["COORDSYS"] = "C"
1517 hdu.header["MOCORDER"] = moc_order
1518 hdu.header["MOCTOOL"] = "lsst.pipe.tasks.hips.GenerateHipsTask"
1519
1520 uri = hips_base_path.join("Moc.fits")
1521
1522 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1523 hdu.writeto(temporary_uri.ospath)
1524
1525 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1526
1527 def _write_allsky_file(self, hips_base_path, allsky_order):
1528 """Write an Allsky.png file.
1529
1530 Parameters
1531 ----------
1532 hips_base_path : `lsst.resources.ResourcePath`
1533 Resource path to the base of the HiPS directory tree.
1534 allsky_order : `int`
1535 HEALPix order of the minimum order to make allsky file.
1536 """
1537 tile_size = self.config.allsky_tilesize
1538
1539 # The Allsky file format is described in
1540 # https://www.ivoa.net/documents/HiPS/20170519/REC-HIPS-1.0-20170519.pdf
1541 # From S4.3.2:
1542 # The Allsky file is built as an array of tiles, stored side by side in
1543 # the left-to-right order. The width of this array must be the square
1544 # root of the number of the tiles of the order. For instance, the width
1545 # of this array at order 3 is 27 ( (int)sqrt(768) ). To avoid having a
1546 # too large Allsky file, the resolution of each tile may be reduced but
1547 # must stay a power of two (typically 64x64 pixels rather than 512x512).
1548
1549 n_tiles = hpg.nside_to_npixel(hpg.order_to_nside(allsky_order))
1550 n_tiles_wide = int(np.floor(np.sqrt(n_tiles)))
1551 n_tiles_high = int(np.ceil(n_tiles / n_tiles_wide))
1552
1553 allsky_image = None
1554
1555 allsky_order_uri = hips_base_path.join(f"Norder{allsky_order}", forceDirectory=True)
1556 pixel_regex = re.compile(r"Npix([0-9]+)\.png$")
1557 png_uris = list(
1558 ResourcePath.findFileResources(
1559 candidates=[allsky_order_uri],
1560 file_filter=pixel_regex,
1561 )
1562 )
1563
1564 for png_uri in png_uris:
1565 matches = re.match(pixel_regex, png_uri.basename())
1566 pix_num = int(matches.group(1))
1567 tile_image = Image.open(io.BytesIO(png_uri.read()))
1568 row = math.floor(pix_num//n_tiles_wide)
1569 column = pix_num % n_tiles_wide
1570 box = (column*tile_size, row*tile_size, (column + 1)*tile_size, (row + 1)*tile_size)
1571 tile_image_shrunk = tile_image.resize((tile_size, tile_size))
1572
1573 if allsky_image is None:
1574 allsky_image = Image.new(
1575 tile_image.mode,
1576 (n_tiles_wide*tile_size, n_tiles_high*tile_size),
1577 )
1578 allsky_image.paste(tile_image_shrunk, box)
1579
1580 uri = allsky_order_uri.join("Allsky.png")
1581
1582 with ResourcePath.temporary_uri(suffix=uri.getExtension()) as temporary_uri:
1583 allsky_image.save(temporary_uri.ospath)
1584
1585 uri.transfer_from(temporary_uri, transfer="copy", overwrite=True)
1586
1587 def _get_dir_number(self, pixel):
1588 """Compute the directory number from a pixel.
1589
1590 Parameters
1591 ----------
1592 pixel : `int`
1593 HEALPix pixel number.
1594
1595 Returns
1596 -------
1597 dir_number : `int`
1598 HiPS directory number.
1599 """
1600 return (pixel//10000)*10000
1601
1602
1603class GenerateColorHipsConnections(pipeBase.PipelineTaskConnections,
1604 dimensions=("instrument", ),
1605 defaultTemplates={"coaddName": "deep"}):
1606 hips_exposure_handles = pipeBase.connectionTypes.Input(
1607 doc="HiPS-compatible HPX images.",
1608 name="{coaddName}Coadd_hpx",
1609 storageClass="ExposureF",
1610 dimensions=("healpix11", "band"),
1611 multiple=True,
1612 deferLoad=True,
1613 )
1614
1615
1616class GenerateColorHipsConfig(GenerateHipsConfig,
1617 pipelineConnections=GenerateColorHipsConnections):
1618 """Configuration parameters for GenerateColorHipsTask."""
1619 blue_channel_band = pexConfig.Field(
1620 doc="Band to use for blue channel of color pngs.",
1621 dtype=str,
1622 default="g",
1623 )
1624 green_channel_band = pexConfig.Field(
1625 doc="Band to use for green channel of color pngs.",
1626 dtype=str,
1627 default="r",
1628 )
1629 red_channel_band = pexConfig.Field(
1630 doc="Band to use for red channel of color pngs.",
1631 dtype=str,
1632 default="i",
1633 )
1634 png_color_asinh_minimum = pexConfig.Field(
1635 doc="AsinhMapping intensity to be mapped to black for color png scaling.",
1636 dtype=float,
1637 default=0.0,
1638 )
1639 png_color_asinh_stretch = pexConfig.Field(
1640 doc="AsinhMapping linear stretch for color png scaling.",
1641 dtype=float,
1642 default=5.0,
1643 )
1644 png_color_asinh_softening = pexConfig.Field(
1645 doc="AsinhMapping softening parameter (Q) for color png scaling.",
1646 dtype=float,
1647 default=8.0,
1648 )
1649
1650
1651class GenerateColorHipsTask(GenerateHipsTask):
1652 """Task for making a HiPS tree with color pngs."""
1653 ConfigClass = GenerateColorHipsConfig
1654 _DefaultName = "generateColorHips"
1655 color_task = True
1656
1657 def _check_data_bands(self, data_bands):
1658 """Check the data for configured bands.
1659
1660 Warn if any color bands are missing data.
1661
1662 Parameters
1663 ----------
1664 data_bands : `set` [`str`]
1665 Bands from the input data.
1666
1667 Returns
1668 -------
1669 bands : `list` [`str`]
1670 List of bands in bgr color order.
1671 """
1672 if len(data_bands) == 0:
1673 raise RuntimeError("GenerateColorHipsTask must have data from at least one band.")
1674
1675 if self.config.blue_channel_band not in data_bands:
1676 self.log.warning(
1677 "Color png blue_channel_band %s not in dataset.",
1678 self.config.blue_channel_band
1679 )
1680 if self.config.green_channel_band not in data_bands:
1681 self.log.warning(
1682 "Color png green_channel_band %s not in dataset.",
1683 self.config.green_channel_band
1684 )
1685 if self.config.red_channel_band not in data_bands:
1686 self.log.warning(
1687 "Color png red_channel_band %s not in dataset.",
1688 self.config.red_channel_band
1689 )
1690
1691 bands = [
1692 self.config.blue_channel_band,
1693 self.config.green_channel_band,
1694 self.config.red_channel_band,
1695 ]
1696
1697 return bands
Pass parameters to a Statistics object.
Definition Statistics.h:83
An integer coordinate rectangle.
Definition Box.h:55
A RangeSet is a set of unsigned 64 bit integers.
Definition RangeSet.h:106
std::shared_ptr< lsst::afw::image::Image< PixelT > > statisticsStack(std::vector< std::shared_ptr< lsst::afw::image::Image< PixelT > > > &images, Property flags, StatisticsControl const &sctrl=StatisticsControl(), std::vector< lsst::afw::image::VariancePixel > const &wvector=std::vector< lsst::afw::image::VariancePixel >(0))
A function to compute some statistics of a stack of Images.
Property stringToStatisticsProperty(std::string const property)
Conversion function to switch a string to a Property (see Statistics.h)