Coverage for python / lsst / dax / apdb / pixelization.py: 21%

1# This file is part of dax_apdb. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (http://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 <http://www.gnu.org/licenses/>. 

21 

22from __future__ import annotations 

23 

24__all__ = ["Pixelization"] 

25 

26import logging 

27 

28from lsst import sphgeom 

29 

30_LOG = logging.getLogger(__name__) 

31 

32 

33class Pixelization: 

34 """Wrapper for pixelization classes from `sphgeom` with configurable 

35 pixelization type and parameters. 

36 

37 Parameters 

38 ---------- 

39 pixelization : `str` 

40 Name of a pixelization type, one of ""htm", "q3c", "mq3c", or 

41 "healpix". 

42 pix_level : `int` 

43 Pixelization level. 

44 pix_max_ranges : `int` 

45 Maximum number of ranges returned from `envelope()` method. 

46 """ 

47 

48 def __init__(self, pixelization: str, pix_level: int, pix_max_ranges: int): 

49 self._pix_max_ranges = pix_max_ranges 

50 self._is_healpix = False 

51 

52 if pixelization == "htm": 

53 self.pixelator = sphgeom.HtmPixelization(pix_level) 

54 elif pixelization == "q3c": 

55 self.pixelator = sphgeom.Q3cPixelization(pix_level) 

56 elif pixelization == "mq3c": 

57 self.pixelator = sphgeom.Mq3cPixelization(pix_level) 

58 elif pixelization == "healpix": 

59 # Healpix does not support maxRanges. 

60 self._pix_max_ranges = 0 

61 self._is_healpix = True 

62 self.pixelator = sphgeom.HealpixPixelization(pix_level) 

63 else: 

64 raise ValueError(f"unknown pixelization: {pixelization}") 

65 

66 def pixels(self, region: sphgeom.Region) -> list[int]: 

67 """Compute set of the pixel indices for given region. 

68 

69 Parameters 

70 ---------- 

71 region : `lsst.sphgeom.Region` 

72 """ 

73 # We want finest set of pixels, so ask as many pixel as reasonable, but 

74 # healpix does not support non-zero maxRanges. 

75 ranges = self.pixelator.envelope(region, 0 if self._is_healpix else 1_000_000) 

76 indices = [] 

77 for lower, upper in ranges: 

78 indices += list(range(lower, upper)) 

79 return indices 

80 

81 def circle_pixels(self, ra: float, dec: float, pad_arcsec: float) -> list[int]: 

82 """Make a list of spatial partitions that a small circle touches. 

83 

84 Parameters 

85 ---------- 

86 ra, dec : `float` 

87 Center of a circle, degrees. 

88 pad_arcsec : `float` 

89 Radius of a circle in arcseconds. 

90 

91 Returns 

92 ------- 

93 pixels : `list` [`int`] 

94 All pixels that envelop the circle. 

95 """ 

96 lon_lat = sphgeom.LonLat.fromDegrees(ra, dec) 

97 center = sphgeom.UnitVector3d(lon_lat) 

98 region = sphgeom.Circle(center, sphgeom.Angle.fromDegrees(pad_arcsec / 3600.0)) 

99 return self.pixels(region) 

100 

101 def pixel(self, direction: sphgeom.UnitVector3d) -> int: 

102 """Compute the index of the pixel for given direction. 

103 

104 Parameters 

105 ---------- 

106 direction : `lsst.sphgeom.UnitVector3d` 

107 """ 

108 index = self.pixelator.index(direction) 

109 return index 

110 

111 def region(self, pixel: int) -> sphgeom.Region: 

112 """Return region corresponding to a pixel index. 

113 

114 Parameters 

115 ---------- 

116 pixel : `int` 

117 Pixel index. 

118 

119 Returns 

120 ------- 

121 region : `lsst.sphgeom.Region` 

122 Region for a given pixel index. 

123 """ 

124 region = self.pixelator.pixel(pixel) 

125 return region 

126 

127 def envelope(self, region: sphgeom.Region) -> list[tuple[int, int]]: 

128 """Generate a set of HTM indices covering specified region. 

129 

130 Parameters 

131 ---------- 

132 region: `sphgeom.Region` 

133 Region that needs to be indexed. 

134 

135 Returns 

136 ------- 

137 ranges : `list` of `tuple` 

138 Sequence of ranges, range is a tuple (minHtmID, maxHtmID). 

139 """ 

140 _LOG.debug("region: %s", region) 

141 indices = self.pixelator.envelope(region, self._pix_max_ranges) 

142 

143 if _LOG.isEnabledFor(logging.DEBUG): 

144 for irange in indices.ranges(): 

145 _LOG.debug( 

146 "range: %s %s", 

147 self.pixelator.toString(irange[0]), 

148 self.pixelator.toString(irange[1]), 

149 ) 

150 

151 return indices.ranges()