Coverage for python / lsst / pipe / tasks / objectMasks.py: 0%

1# This file is part of pipe_tasks. 

2# 

3# Developed for the LSST Data Management System. 

4# This product includes software developed by the LSST Project 

5# (https://www.lsst.org). 

6# See the COPYRIGHT file at the top-level directory of this distribution 

7# for details of code ownership. 

8# 

9# This program is free software: you can redistribute it and/or modify 

10# it under the terms of the GNU General Public License as published by 

11# the Free Software Foundation, either version 3 of the License, or 

12# (at your option) any later version. 

13# 

14# This program is distributed in the hope that it will be useful, 

15# but WITHOUT ANY WARRANTY; without even the implied warranty of 

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

17# GNU General Public License for more details. 

18# 

19# You should have received a copy of the GNU General Public License 

20# along with this program. If not, see <https://www.gnu.org/licenses/>. 

21 

22__all__ = ["ObjectMaskCatalog", "RegionFileFormatter"] 

23 

24import re 

25import logging 

26import lsst.daf.base as dafBase 

27import lsst.geom as geom 

28import lsst.afw.table as afwTable 

29from typing import Any 

30from lsst.daf.butler import FormatterV2, FormatterNotImplementedError 

31from lsst.resources import ResourcePath 

32 

33 

34class ObjectMaskCatalog: 

35 """Class to support bright object masks 

36 """ 

37 

38 def __init__(self): 

39 schema = afwTable.SimpleTable.makeMinimalSchema() 

40 schema.addField("type", str, "type of region (e.g. box, circle)", size=10) 

41 schema.addField("radius", "Angle", "radius of mask (if type == circle") 

42 schema.addField("height", "Angle", "height of mask (if type == box)") 

43 schema.addField("width", "Angle", "width of mask (if type == box)") 

44 schema.addField("angle", "Angle", "rotation of mask (if type == box)") 

45 schema.addField("mag", float, "object's magnitude") 

46 

47 self._catalog = afwTable.SimpleCatalog(schema) 

48 self._catalog.table.setMetadata(dafBase.PropertyList()) 

49 

50 self.table = self._catalog.table 

51 self.addNew = self._catalog.addNew 

52 

53 def __len__(self): 

54 return len(self._catalog) 

55 

56 def __iter__(self): 

57 return iter(self._catalog) 

58 

59 def __getitem__(self, i): 

60 return self._catalog.__getitem__(i) 

61 

62 def __setitem__(self, i, v): 

63 return self._catalog.__setitem__(i, v) 

64 

65 @classmethod 

66 def read(cls, fileName): 

67 """Read a ds9 region file, returning a ObjectMaskCatalog object 

68 

69 The files should be structured as follows: 

70 

71 # Description of catalogue as a comment 

72 # CATALOG: catalog-id-string 

73 # TRACT: 0 

74 # PATCH: 5,4 

75 # FILTER: HSC-I 

76 

77 wcs; fk5 

78 

79 circle(RA, DEC, RADIUS) # ID: 1, mag: 12.34 

80 box(RA, DEC, XSIZE, YSIZE, THETA) # ID: 2, mag: 23.45 

81 ... 

82 

83 The ", mag: XX.YY" is optional 

84 

85 The commented lines must be present, with the relevant fields such as 

86 tract patch and filter filled in. The coordinate system must be listed 

87 as above. Each patch is specified as a box or circle, with RA, DEC, 

88 and dimensions specified in decimal degrees (with or without an 

89 explicit "d"). 

90 

91 Only (axis-aligned) boxes and circles are currently supported as 

92 region definitions. 

93 """ 

94 

95 log = logging.getLogger("lsst.ObjectMaskCatalog") 

96 

97 brightObjects = cls() 

98 checkedWcsIsFk5 = False 

99 NaN = float("NaN")*geom.degrees 

100 

101 nFormatError = 0 # number of format errors seen 

102 with open(fileName) as fd: 

103 for lineNo, line in enumerate(fd.readlines(), 1): 

104 line = line.rstrip() 

105 

106 if re.search(r"^\s*#", line): 

107 # 

108 # Parse any line of the form "# key : value" and put them into the metadata. 

109 # 

110 # The medatdata values must be defined as outlined in the above docstring 

111 # 

112 # The value of these three keys will be checked, 

113 # so get them right! 

114 # 

115 mat = re.search(r"^\s*#\s*([a-zA-Z][a-zA-Z0-9_]+)\s*:\s*(.*)", line) 

116 if mat: 

117 key, value = mat.group(1).lower(), mat.group(2) 

118 if key == "tract": 

119 value = int(value) 

120 

121 brightObjects.table.getMetadata().set(key, value) 

122 

123 line = re.sub(r"^\s*#.*", "", line) 

124 if not line: 

125 continue 

126 

127 if re.search(r"^\s*wcs\s*;\s*fk5\s*$", line, re.IGNORECASE): 

128 checkedWcsIsFk5 = True 

129 continue 

130 

131 # This regular expression parses the regions file for each region to be masked, 

132 # with the format as specified in the above docstring. 

133 mat = re.search(r"^\s*(box|circle)" 

134 r"(?:\s+|\s*\(\s*)" # open paren or space 

135 r"([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)([d]*)" # ra + units 

136 r"(?:\s+|\s*,\s*)" # sep 

137 r"([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)([d]*)" # dec + units 

138 r"(?:\s+|\s*,\s*)" # sep 

139 r"([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)([d]*)" # param1 + units 

140 r"(?:" # start optional 1 

141 r"(?:\s+|\s*,\s*)" # sep 

142 r"([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)([d]*)" # param2 + units 

143 r"(?:" # start optional 2 

144 r"(?:\s+|\s*,\s*)" # sep 

145 r"([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?)([d]*)" # param3 + units 

146 ")?" # end optional 2 

147 ")?" # end optional 1 

148 r"(?:\s*|\s*\)\s*)" # close paren or space 

149 r"#\s*ID:[\w\s]*(\d+)" # start comment, ID 

150 r"(?:\s*,?\s*mag:\s*([+\-]?(?:0|[1-9]\d*)(?:\.\d*)?(?:[eE][+\-]?\d+)?))?" 

151 r"\s*$", line) 

152 if mat: 

153 _type, ra, raUnit, dec, decUnit, \ 

154 param1, param1Unit, param2, param2Unit, param3, param3Unit, \ 

155 _id, mag = mat.groups() 

156 

157 _id = int(_id) 

158 if mag is None: 

159 mag = NaN 

160 else: 

161 mag = float(mag) 

162 

163 ra = convertToAngle(ra, raUnit, "ra", fileName, lineNo) 

164 dec = convertToAngle(dec, decUnit, "dec", fileName, lineNo) 

165 

166 radius = NaN 

167 width = NaN 

168 height = NaN 

169 angle = 0.0*geom.degrees 

170 

171 if _type == "box": 

172 width = convertToAngle(param1, param1Unit, "width", fileName, lineNo) 

173 height = convertToAngle(param2, param2Unit, "height", fileName, lineNo) 

174 if param3 is not None: 

175 angle = convertToAngle(param3, param3Unit, "angle", fileName, lineNo) 

176 

177 if angle != 0.0: 

178 log.warning("Rotated boxes are not supported: \"%s\" at %s:%d", 

179 line, fileName, lineNo) 

180 nFormatError += 1 

181 elif _type == "circle": 

182 radius = convertToAngle(param1, param1Unit, "radius", fileName, lineNo) 

183 

184 if not (param2 is None and param3 is None): 

185 log.warning("Extra parameters for circle: \"%s\" at %s:%d", 

186 line, fileName, lineNo) 

187 nFormatError += 1 

188 

189 rec = brightObjects.addNew() 

190 # N.b. rec["coord"] = Coord is not supported, so we have to use the setter 

191 rec["type"] = _type 

192 rec["id"] = _id 

193 rec["mag"] = mag 

194 rec.setCoord(geom.SpherePoint(ra, dec)) 

195 

196 rec["angle"] = angle 

197 rec["height"] = height 

198 rec["width"] = width 

199 rec["radius"] = radius 

200 else: 

201 log.warning("Unexpected line \"%s\" at %s:%d", line, fileName, lineNo) 

202 nFormatError += 1 

203 

204 if nFormatError > 0: 

205 raise RuntimeError("Saw %d formatting errors in %s" % (nFormatError, fileName)) 

206 

207 if not checkedWcsIsFk5: 

208 raise RuntimeError("Expected to see a line specifying an fk5 wcs in %s" % fileName) 

209 

210 # This makes the deep copy contiguous in memory so that a ColumnView can be exposed to Numpy 

211 brightObjects._catalog = brightObjects._catalog.copy(True) 

212 

213 return brightObjects 

214 

215 

216def convertToAngle(var, varUnit, what, fileName, lineNo): 

217 """Given a variable and its units, return an geom.Angle 

218 

219 what, fileName, and lineNo are used to generate helpful error messages 

220 """ 

221 var = float(var) 

222 

223 if varUnit in ("d", "", None): 

224 pass 

225 elif varUnit == "'": 

226 var /= 60.0 

227 elif varUnit == '"': 

228 var /= 3600.0 

229 else: 

230 raise RuntimeError("unsupported unit \"%s\" for %s at %s:%d" % 

231 (varUnit, what, fileName, lineNo)) 

232 

233 return var*geom.degrees 

234 

235 

236class RegionFileFormatter(FormatterV2): 

237 """Plugin for reading DS9 region file catalogs with Gen3 Butler. 

238 """ 

239 default_extension = ".reg" 

240 can_read_from_local_file = True 

241 

242 def read_from_local_file( 

243 self, path: str, component: str | None = None, expected_size: int = -1 

244 ) -> Any: 

245 # Docstring inherited. 

246 pytype = self.file_descriptor.storageClass.pytype 

247 return pytype.read(path) 

248 

249 def write_local_file(self, in_memory_dataset: Any, uri: ResourcePath) -> None: 

250 # Docstring inherited. 

251 raise FormatterNotImplementedError("Write not implemented.")