Coverage for python / lsst / analysis / tools / tasks / photometricCatalogMatch.py: 23%

99 statements  

« prev     ^ index     » next       coverage.py v7.13.5, created at 2026-04-15 00:23 +0000

1# This file is part of analysis_tools. 

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__ = ("PhotometricCatalogMatchConfig", "PhotometricCatalogMatchTask") 

23 

24 

25import lsst.geom 

26import lsst.pipe.base as pipeBase 

27import numpy as np 

28from astropy.table import Table 

29from lsst.pipe.tasks.loadReferenceCatalog import LoadReferenceCatalogTask 

30 

31from ..actions.vector import StarSelector, VisitPlotFlagSelector 

32from ..tasks.catalogMatch import CatalogMatchConfig, CatalogMatchConnections, CatalogMatchTask 

33 

34 

35class PhotometricCatalogMatchConnections( 

36 CatalogMatchConnections, 

37 dimensions=("tract", "skymap"), 

38 defaultTemplates={"targetCatalog": "objectTable_tract", "refCatalog": "ps1_pv3_3pi_20170110"}, 

39): 

40 matchedCatalog = pipeBase.connectionTypes.Output( 

41 doc="Catalog with matched target and reference objects with separations", 

42 name="{targetCatalog}_{refCatalog}_match_photom", 

43 storageClass="ArrowAstropy", 

44 dimensions=("tract", "skymap"), 

45 ) 

46 

47 

48class PhotometricCatalogMatchConfig( 

49 CatalogMatchConfig, 

50 pipelineConnections=PhotometricCatalogMatchConnections, 

51): 

52 def setDefaults(self): 

53 super().setDefaults() 

54 self.matchesRefCat = True 

55 self.referenceCatalogLoader.doReferenceSelection = False 

56 self.referenceCatalogLoader.doApplyColorTerms = True 

57 

58 

59class PhotometricCatalogMatchTask(CatalogMatchTask): 

60 """A wrapper task to provide the information that 

61 is specific to the photometric reference catalog. 

62 """ 

63 

64 ConfigClass = PhotometricCatalogMatchConfig 

65 _DefaultName = "analysisToolsPhotometricCatalogMatch" 

66 

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

68 """Run the matching to the photometric reference catalog. 

69 

70 Parameters 

71 ---------- 

72 butlerQC : `lsst.pipe.base.QuantumContext` 

73 inputRefs : `lsst.pipe.base.InputQuantizedConnection` 

74 outputRefs : `lsst.pipe.base.OutputQuantizedConnection` 

75 """ 

76 

77 inputs = butlerQC.get(inputRefs) 

78 bands = [] 

79 for filterName in self.config.filterNames: 

80 bands.append(self.config.referenceCatalogLoader.refObjLoader.filterMap[filterName]) 

81 

82 # For some reason the imsim filterMaps don't work the same way 

83 # as the HSC ones do, this is a bit hacky but fixes this. 

84 # This code makes the assumption that filterMap is a dict 

85 # mapping observed filter names: band, but filterMap is 

86 # currently defined as a mapping of: 

87 # observed filter names -> reference catalog filter names. 

88 # Therefore, when the reference catalog filter name is 

89 # not the band name we need this work around. 

90 # TODO: workaround for DM-46728 

91 if bands[0].startswith("lsst") or "sim" in bands[0] or "smeared" in bands[0]: 

92 bands = self.config.filterNames 

93 elif bands[0].startswith("monster"): 

94 # for the_monster_* the reference catalog filter name is 

95 # "monster_{system}_{band}" the last character is the band 

96 bands = [band[-1] for band in bands] 

97 

98 columns = self.prepColumns(bands) 

99 table = inputs["catalog"].get(parameters={"columns": columns}) 

100 

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

102 

103 loaderTask = LoadReferenceCatalogTask( 

104 config=self.config.referenceCatalogLoader, 

105 dataIds=[ref.dataId for ref in inputRefs.refCat], 

106 name=inputs["refCat"][0].ref.datasetType.name, 

107 refCats=inputs["refCat"], 

108 ) 

109 

110 skymap = inputs.pop("skymap") 

111 loadedRefCat = self._loadRefCat(loaderTask, skymap[tract]) 

112 outputs = self.run(targetCatalog=table, refCatalog=loadedRefCat, bands=bands) 

113 

114 butlerQC.put(outputs, outputRefs) 

115 

116 

117class PhotometricCatalogMatchVisitConnections( 

118 pipeBase.PipelineTaskConnections, 

119 dimensions=("visit",), 

120 defaultTemplates={"targetCatalog": "sourceTable_visit", "refCatalog": "ps1_pv3_3pi_20170110"}, 

121): 

122 catalog = pipeBase.connectionTypes.Input( 

123 doc="The visit-wide catalog to make plots from.", 

124 storageClass="ArrowAstropy", 

125 name="{targetCatalog}", 

126 dimensions=("visit",), 

127 deferLoad=True, 

128 ) 

129 

130 refCat = pipeBase.connectionTypes.PrerequisiteInput( 

131 doc="The photometric reference catalog to match to.", 

132 name="{refCatalog}", 

133 storageClass="SimpleCatalog", 

134 dimensions=("skypix",), 

135 deferLoad=True, 

136 multiple=True, 

137 ) 

138 

139 visitSummaryTable = pipeBase.connectionTypes.Input( 

140 doc="A summary table of the ccds in the visit", 

141 storageClass="ExposureCatalog", 

142 name="finalVisitSummary", 

143 dimensions=("visit",), 

144 ) 

145 

146 matchedCatalog = pipeBase.connectionTypes.Output( 

147 doc="Catalog with matched target and reference objects with separations", 

148 name="{targetCatalog}_{refCatalog}_match_photom", 

149 storageClass="ArrowAstropy", 

150 dimensions=("visit",), 

151 ) 

152 

153 

154class PhotometricCatalogMatchVisitConfig( 

155 PhotometricCatalogMatchConfig, pipelineConnections=PhotometricCatalogMatchVisitConnections 

156): 

157 def setDefaults(self): 

158 super().setDefaults() 

159 self.idColumn = "sourceId" 

160 self.matchesRefCat = True 

161 self.filterNames = [] 

162 self.extraPerBandColumns = [] 

163 self.patchColumn = "" 

164 self.selectorBands = [] 

165 self.extraColumns = ["x", "y", "ap09Flux", "ap09FluxErr"] 

166 self.selectorActions.flagSelector = VisitPlotFlagSelector 

167 self.sourceSelectorActions.sourceSelector = StarSelector() 

168 self.sourceSelectorActions.sourceSelector.vectorKey = "extendedness" 

169 self.extraColumnSelectors.selector1.fluxType = "psfFlux" 

170 self.extraColumnSelectors.selector2.vectorKey = "extendedness" 

171 self.extraColumnSelectors.selector3.vectorKey = "extendedness" 

172 self.extraColumnSelectors.selector4 = VisitPlotFlagSelector 

173 

174 

175class PhotometricCatalogMatchVisitTask(PhotometricCatalogMatchTask): 

176 """A wrapper task to provide the information that 

177 is specific to the photometric reference catalog. 

178 """ 

179 

180 ConfigClass = PhotometricCatalogMatchVisitConfig 

181 _DefaultName = "analysisToolsPhotometricCatalogMatchVisit" 

182 

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

184 """Run the matching to the photometric reference catalog. 

185 

186 Parameters 

187 ---------- 

188 butlerQC : `lsst.pipe.base.QuantumContext` 

189 inputRefs : `lsst.pipe.base.InputQuantizedConnection` 

190 outputRefs : `lsst.pipe.base.OutputQuantizedConnection` 

191 """ 

192 

193 inputs = butlerQC.get(inputRefs) 

194 physicalFilter = inputs["catalog"].dataId["physical_filter"] 

195 

196 # For some reason the imsim filterMaps don't work the same way as 

197 # the HSC ones do, this is a bit hacky but fixes this 

198 # This code makes the assumption that filterMap is a dict 

199 # mapping observed filter names: band, but filterMap is 

200 # currently defined as a mapping of: 

201 # observed filter names -> reference catalog filter names. 

202 # Therefore, when the reference catalog filter name is 

203 # not the band name we need this work around. 

204 # TODO: workaround for DM-46728 

205 if "sim" in physicalFilter: 

206 physicalFilter = physicalFilter[0] 

207 bands = [physicalFilter] 

208 elif self.config.referenceCatalogLoader.refObjLoader.filterMap[physicalFilter].startswith("monster"): 

209 # filtermap gives the column name for the refcat 

210 # which is monster_{system}_{band} 

211 bands = [self.config.referenceCatalogLoader.refObjLoader.filterMap[physicalFilter][-1]] 

212 else: 

213 bands = [self.config.referenceCatalogLoader.refObjLoader.filterMap[physicalFilter]] 

214 

215 # No bands needed for visit tables 

216 # but we do need them later for the matching 

217 columns = [self.config.idColumn, "coord_ra", "coord_dec", "detector"] 

218 columns.extend(self.config.extraColumns) 

219 for selectorAction in [ 

220 self.config.selectorActions, 

221 self.config.sourceSelectorActions, 

222 self.config.extraColumnSelectors, 

223 ]: 

224 for selector in selectorAction: 

225 selectorSchema = selector.getFormattedInputSchema() 

226 columns += [s[0] for s in selectorSchema] 

227 

228 table = inputs["catalog"].get(parameters={"columns": columns}) 

229 

230 loaderTask = LoadReferenceCatalogTask( 

231 config=self.config.referenceCatalogLoader, 

232 dataIds=[ref.dataId for ref in inputRefs.refCat], 

233 name=inputs["refCat"][0].ref.datasetType.name, 

234 refCats=inputs["refCat"], 

235 ) 

236 

237 visitSummaryTable = inputs.pop("visitSummaryTable") 

238 loadedRefCat = self._loadRefCat(loaderTask, visitSummaryTable, physicalFilter) 

239 outputs = self.run(targetCatalog=table, refCatalog=loadedRefCat, bands=bands) 

240 

241 # The matcher adds the band to the front of the columns 

242 # but the visit plots aren't expecting it 

243 cols = list(outputs.matchedCatalog.columns) 

244 for col in cols: 

245 if col[:2] == bands[0] + "_": 

246 outputs.matchedCatalog.rename_column(col, col[2:]) 

247 

248 butlerQC.put(outputs, outputRefs) 

249 

250 def _loadRefCat(self, loaderTask, visitSummaryTable, physicalFilter): 

251 """Make a reference catalog with coordinates in degrees 

252 

253 Parameters 

254 ---------- 

255 visitSummaryTable : `lsst.afw.table.ExposureCatalog` 

256 The table of visit information 

257 """ 

258 # Get convex hull around the detectors, then get its center and radius 

259 corners = [] 

260 for visSum in visitSummaryTable: 

261 for ra, dec in zip(visSum["raCorners"], visSum["decCorners"]): 

262 # If the coordinates are nan then don't keep going 

263 # because it crashes later 

264 if not np.isfinite(ra) or not np.isfinite(dec): 

265 raise pipeBase.NoWorkFound("Visit summary corners not finite") 

266 corners.append(lsst.geom.SpherePoint(ra, dec, units=lsst.geom.degrees).getVector()) 

267 visitBoundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(corners).getBoundingCircle() 

268 center = lsst.geom.SpherePoint(visitBoundingCircle.getCenter()) 

269 radius = visitBoundingCircle.getOpeningAngle() 

270 

271 # Get the observation date of the visit 

272 epoch = visSum.getVisitInfo().getDate().toAstropy() 

273 

274 # Load the reference catalog in the skyCircle of the detectors, then 

275 # convert the coordinates to degrees and convert the catalog to a 

276 # dataframe 

277 

278 loadedRefCat = loaderTask.getSkyCircleCatalog(center, radius, [physicalFilter], epoch=epoch) 

279 

280 return Table(loadedRefCat)