Coverage for python / lsst / analysis / tools / tasks / astrometricCatalogMatch.py: 33%

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__ = ( 

23 "AstrometricCatalogMatchConfig", 

24 "AstrometricCatalogMatchTask", 

25 "AstrometricCatalogMatchVisitConfig", 

26 "AstrometricCatalogMatchVisitTask", 

27) 

28 

29import lsst.geom 

30import lsst.pex.config as pexConfig 

31import lsst.pipe.base as pipeBase 

32import numpy as np 

33from astropy.table import Table 

34from astropy.time import Time 

35from lsst.pex.config.configurableActions import ConfigurableActionStructField 

36from lsst.pipe.tasks.loadReferenceCatalog import LoadReferenceCatalogTask 

37 

38from ..actions.vector import StarSelector, VisitPlotFlagSelector 

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

40 

41 

42class AstrometricCatalogMatchConnections( 

43 CatalogMatchConnections, 

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

45 defaultTemplates={"targetCatalog": "objectTable_tract", "refCatalog": "gaia_dr3_20230707"}, 

46): 

47 matchedCatalog = pipeBase.connectionTypes.Output( 

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

49 name="{targetCatalog}_{refCatalog}_match_astrom", 

50 storageClass="ArrowAstropy", 

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

52 ) 

53 

54 

55class AstrometricCatalogMatchConfig( 

56 CatalogMatchConfig, pipelineConnections=AstrometricCatalogMatchConnections 

57): 

58 bands = pexConfig.ListField[str]( 

59 doc="The bands to persist downstream", 

60 default=["u", "g", "r", "i", "z", "y"], 

61 ) 

62 extraPerBandColumns = pexConfig.ListField[str]( 

63 doc="Other columns to load that should be loaded for each band individually.", 

64 default=["epoch"], 

65 ) 

66 extraColumns = pexConfig.ListField[str]( 

67 doc="Other catalog columns to persist to downstream tasks", 

68 default=["patch", "ebv", "refBand"], 

69 ) 

70 

71 def setDefaults(self): 

72 super().setDefaults() 

73 self.matchesRefCat = True 

74 self.referenceCatalogLoader.doApplyColorTerms = False 

75 self.referenceCatalogLoader.refObjLoader.requireProperMotion = True 

76 self.referenceCatalogLoader.refObjLoader.anyFilterMapsToThis = "phot_g_mean" 

77 self.connections.refCatalog = "gaia_dr3_20230707" 

78 

79 

80class AstrometricCatalogMatchTask(CatalogMatchTask): 

81 """Match a tract-level catalog to a reference catalog""" 

82 

83 ConfigClass = AstrometricCatalogMatchConfig 

84 _DefaultName = "analysisToolsAstrometricCatalogMatch" 

85 

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

87 # Docs inherited from base class 

88 

89 inputs = butlerQC.get(inputRefs) 

90 columns = self.prepColumns(self.config.bands) 

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

92 

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

94 

95 loaderTask = LoadReferenceCatalogTask( 

96 config=self.config.referenceCatalogLoader, 

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

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

99 refCats=inputs["refCat"], 

100 ) 

101 

102 skymap = inputs.pop("skymap") 

103 # We cannot load the reference catalog at different epochs for each 

104 # object. Instead use the median of the epochs corresponding to the 

105 # reference band for each object. 

106 refEpochs = np.array( 

107 [ 

108 obj[obj["refBand"] + "_epoch"] 

109 for obj in table 

110 if not np.ma.is_masked(obj[obj["refBand"] + "_epoch"]) 

111 ] 

112 ) 

113 epoch = Time(np.median(refEpochs), format="mjd") 

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

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

116 

117 butlerQC.put(outputs, outputRefs) 

118 

119 

120class AstrometricCatalogMatchVisitConnections( 

121 pipeBase.PipelineTaskConnections, 

122 dimensions=("visit",), 

123 defaultTemplates={"targetCatalog": "sourceTable_visit", "refCatalog": "gaia_dr3_20230707"}, 

124): 

125 catalog = pipeBase.connectionTypes.Input( 

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

127 storageClass="ArrowAstropy", 

128 name="{targetCatalog}", 

129 dimensions=("visit",), 

130 deferLoad=True, 

131 ) 

132 

133 refCat = pipeBase.connectionTypes.PrerequisiteInput( 

134 doc="The astrometry reference catalog to match to loaded input catalog sources.", 

135 name="{refCatalog}", 

136 storageClass="SimpleCatalog", 

137 dimensions=("skypix",), 

138 deferLoad=True, 

139 multiple=True, 

140 ) 

141 

142 visitSummaryTable = pipeBase.connectionTypes.Input( 

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

144 storageClass="ExposureCatalog", 

145 name="finalVisitSummary", 

146 dimensions=("visit",), 

147 ) 

148 

149 matchedCatalog = pipeBase.connectionTypes.Output( 

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

151 name="{targetCatalog}_{refCatalog}_match_astrom", 

152 storageClass="ArrowAstropy", 

153 dimensions=("visit",), 

154 ) 

155 

156 

157class AstrometricCatalogMatchVisitConfig( 

158 AstrometricCatalogMatchConfig, pipelineConnections=AstrometricCatalogMatchVisitConnections 

159): 

160 selectorActions = ConfigurableActionStructField( 

161 doc="Which selectors to use to narrow down the data for QA plotting.", 

162 default={"flagSelector": VisitPlotFlagSelector}, 

163 ) 

164 

165 extraColumns = pexConfig.ListField[str]( 

166 doc="Other catalog columns to persist to downstream tasks", 

167 default=["psfFlux", "psfFluxErr"], 

168 ) 

169 

170 bands = pexConfig.ListField[str]( 

171 doc="The bands to persist downstream", 

172 default=[], 

173 ) 

174 

175 def setDefaults(self): 

176 self.matchesRefCat = True 

177 self.idColumn = "sourceId" 

178 # sourceSelectorActions.sourceSelector is StarSelector 

179 self.sourceSelectorActions.sourceSelector = StarSelector() 

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

181 # extraColumnSelectors.selector1 is SnSelector 

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

183 # extraColumnSelectors.selector2 is GalaxySelector 

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

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

186 self.extraColumnSelectors.selector4 = VisitPlotFlagSelector 

187 self.referenceCatalogLoader.doApplyColorTerms = False 

188 self.referenceCatalogLoader.refObjLoader.requireProperMotion = False 

189 self.referenceCatalogLoader.refObjLoader.anyFilterMapsToThis = "phot_g_mean" 

190 

191 

192class AstrometricCatalogMatchVisitTask(AstrometricCatalogMatchTask): 

193 """Match a visit-level catalog to a reference catalog""" 

194 

195 ConfigClass = AstrometricCatalogMatchVisitConfig 

196 _DefaultName = "analysisToolsAstrometricCatalogMatchVisit" 

197 

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

199 # Docs inherited from base class 

200 

201 inputs = butlerQC.get(inputRefs) 

202 

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

204 columns.extend(self.config.extraColumns) 

205 for selectorAction in [ 

206 self.config.selectorActions, 

207 self.config.sourceSelectorActions, 

208 self.config.extraColumnSelectors, 

209 ]: 

210 for selector in selectorAction: 

211 selectorSchema = selector.getFormattedInputSchema() 

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

213 

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

215 

216 loaderTask = LoadReferenceCatalogTask( 

217 config=self.config.referenceCatalogLoader, 

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

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

220 refCats=inputs["refCat"], 

221 ) 

222 

223 visitSummaryTable = inputs.pop("visitSummaryTable") 

224 loadedRefCat, badDetectors = self._loadRefCat(loaderTask, visitSummaryTable) 

225 for badDetector in badDetectors: 

226 mask = table["detector"] != badDetector 

227 self.log.info( 

228 "Trimming %d sources from detector %s with incomplete or missing astrometry.", 

229 len(table) - sum(mask), 

230 badDetector, 

231 ) 

232 table = table[mask] 

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

234 

235 butlerQC.put(outputs, outputRefs) 

236 

237 def _loadRefCat(self, loaderTask, visitSummaryTable): 

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

239 

240 Parameters 

241 ---------- 

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

243 The table of visit information 

244 

245 Returns 

246 ------- 

247 refCat : `astropy.table.Table` 

248 Reference catalog. 

249 badDetectors : `list` [ `int` ] 

250 IDs of detectors with invalid visit summaries (should be skipped). 

251 """ 

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

253 corners = [] 

254 badDetectors = [] 

255 for visSum in visitSummaryTable: 

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

257 # If the coordinates are nan skip this detector. 

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

259 badDetectors.append(visSum.getId()) 

260 break 

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

262 if not corners: 

263 raise pipeBase.NoWorkFound("No visit summary corners were finite.") 

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

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

266 radius = visitBoundingCircle.getOpeningAngle() 

267 

268 # Get the observation date of the visit 

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

270 

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

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

273 # dataframe 

274 

275 filterName = self.config.referenceCatalogLoader.refObjLoader.anyFilterMapsToThis 

276 loadedRefCat = loaderTask.getSkyCircleCatalog(center, radius, filterName, epoch=epoch) 

277 

278 return Table(loadedRefCat), badDetectors