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 numpy as np 

30from astropy.table import Table 

31from astropy.time import Time 

32 

33import lsst.geom 

34import lsst.pex.config as pexConfig 

35import lsst.pipe.base as pipeBase 

36from lsst.pex.config.configurableActions import ConfigurableActionStructField 

37from lsst.pipe.tasks.loadReferenceCatalog import LoadReferenceCatalogTask 

38 

39from ..actions.vector import StarSelector, VisitPlotFlagSelector 

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

41 

42 

43class AstrometricCatalogMatchConnections( 

44 CatalogMatchConnections, 

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

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

47): 

48 matchedCatalog = pipeBase.connectionTypes.Output( 

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

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

51 storageClass="ArrowAstropy", 

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

53 ) 

54 

55 

56class AstrometricCatalogMatchConfig( 

57 CatalogMatchConfig, pipelineConnections=AstrometricCatalogMatchConnections 

58): 

59 bands = pexConfig.ListField[str]( 

60 doc="The bands to persist downstream", 

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

62 ) 

63 extraPerBandColumns = pexConfig.ListField[str]( 

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

65 default=["epoch"], 

66 ) 

67 extraColumns = pexConfig.ListField[str]( 

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

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

70 ) 

71 

72 def setDefaults(self): 

73 super().setDefaults() 

74 self.matchesRefCat = True 

75 self.referenceCatalogLoader.doApplyColorTerms = False 

76 self.referenceCatalogLoader.refObjLoader.requireProperMotion = True 

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

78 self.connections.refCatalog = "gaia_dr3_20230707" 

79 

80 

81class AstrometricCatalogMatchTask(CatalogMatchTask): 

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

83 

84 ConfigClass = AstrometricCatalogMatchConfig 

85 _DefaultName = "analysisToolsAstrometricCatalogMatch" 

86 

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

88 # Docs inherited from base class 

89 

90 inputs = butlerQC.get(inputRefs) 

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

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

93 

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

95 

96 loaderTask = LoadReferenceCatalogTask( 

97 config=self.config.referenceCatalogLoader, 

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

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

100 refCats=inputs["refCat"], 

101 ) 

102 

103 skymap = inputs.pop("skymap") 

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

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

106 # reference band for each object. 

107 refEpochs = np.array( 

108 [ 

109 obj[obj["refBand"] + "_epoch"] 

110 for obj in table 

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

112 ] 

113 ) 

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

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

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

117 

118 butlerQC.put(outputs, outputRefs) 

119 

120 

121class AstrometricCatalogMatchVisitConnections( 

122 pipeBase.PipelineTaskConnections, 

123 dimensions=("visit",), 

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

125): 

126 catalog = pipeBase.connectionTypes.Input( 

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

128 storageClass="ArrowAstropy", 

129 name="{targetCatalog}", 

130 dimensions=("visit",), 

131 deferLoad=True, 

132 ) 

133 

134 refCat = pipeBase.connectionTypes.PrerequisiteInput( 

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

136 name="{refCatalog}", 

137 storageClass="SimpleCatalog", 

138 dimensions=("skypix",), 

139 deferLoad=True, 

140 multiple=True, 

141 ) 

142 

143 visitSummaryTable = pipeBase.connectionTypes.Input( 

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

145 storageClass="ExposureCatalog", 

146 name="finalVisitSummary", 

147 dimensions=("visit",), 

148 ) 

149 

150 matchedCatalog = pipeBase.connectionTypes.Output( 

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

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

153 storageClass="ArrowAstropy", 

154 dimensions=("visit",), 

155 ) 

156 

157 

158class AstrometricCatalogMatchVisitConfig( 

159 AstrometricCatalogMatchConfig, pipelineConnections=AstrometricCatalogMatchVisitConnections 

160): 

161 selectorActions = ConfigurableActionStructField( 

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

163 default={"flagSelector": VisitPlotFlagSelector}, 

164 ) 

165 

166 extraColumns = pexConfig.ListField[str]( 

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

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

169 ) 

170 

171 bands = pexConfig.ListField[str]( 

172 doc="The bands to persist downstream", 

173 default=[], 

174 ) 

175 

176 def setDefaults(self): 

177 self.matchesRefCat = True 

178 self.idColumn = "sourceId" 

179 # sourceSelectorActions.sourceSelector is StarSelector 

180 self.sourceSelectorActions.sourceSelector = StarSelector() 

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

182 # extraColumnSelectors.selector1 is SnSelector 

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

184 # extraColumnSelectors.selector2 is GalaxySelector 

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

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

187 self.extraColumnSelectors.selector4 = VisitPlotFlagSelector 

188 self.referenceCatalogLoader.doApplyColorTerms = False 

189 self.referenceCatalogLoader.refObjLoader.requireProperMotion = False 

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

191 

192 

193class AstrometricCatalogMatchVisitTask(AstrometricCatalogMatchTask): 

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

195 

196 ConfigClass = AstrometricCatalogMatchVisitConfig 

197 _DefaultName = "analysisToolsAstrometricCatalogMatchVisit" 

198 

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

200 # Docs inherited from base class 

201 

202 inputs = butlerQC.get(inputRefs) 

203 

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

205 columns.extend(self.config.extraColumns) 

206 for selectorAction in [ 

207 self.config.selectorActions, 

208 self.config.sourceSelectorActions, 

209 self.config.extraColumnSelectors, 

210 ]: 

211 for selector in selectorAction: 

212 selectorSchema = selector.getFormattedInputSchema() 

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

214 

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

216 

217 loaderTask = LoadReferenceCatalogTask( 

218 config=self.config.referenceCatalogLoader, 

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

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

221 refCats=inputs["refCat"], 

222 ) 

223 

224 visitSummaryTable = inputs.pop("visitSummaryTable") 

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

226 for badDetector in badDetectors: 

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

228 self.log.info( 

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

230 len(table) - sum(mask), 

231 badDetector, 

232 ) 

233 table = table[mask] 

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

235 

236 butlerQC.put(outputs, outputRefs) 

237 

238 def _loadRefCat(self, loaderTask, visitSummaryTable): 

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

240 

241 Parameters 

242 ---------- 

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

244 The table of visit information 

245 

246 Returns 

247 ------- 

248 refCat : `astropy.table.Table` 

249 Reference catalog. 

250 badDetectors : `list` [ `int` ] 

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

252 """ 

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

254 corners = [] 

255 badDetectors = [] 

256 for visSum in visitSummaryTable: 

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

258 # If the coordinates are nan skip this detector. 

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

260 badDetectors.append(visSum.getId()) 

261 break 

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

263 if not corners: 

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

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

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

267 radius = visitBoundingCircle.getOpeningAngle() 

268 

269 # Get the observation date of the visit 

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

271 

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

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

274 # dataframe 

275 

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

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

278 

279 return Table(loadedRefCat), badDetectors