Coverage for python / lsst / ip / diffim / computeSpatiallySampledMetrics.py: 24%

1# This file is part of ip_diffim. 

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 

22import numpy as np 

23 

24import lsst.geom 

25 

26import lsst.afw.table as afwTable 

27import lsst.pipe.base as pipeBase 

28import lsst.pex.config as pexConfig 

29 

30from lsst.ip.diffim.utils import getPsfFwhm, angleMean, evaluateMaskFraction, getKernelCenterDisplacement 

31from lsst.meas.algorithms import SkyObjectsTask 

32from lsst.pex.exceptions import InvalidParameterError, RangeError 

33from lsst.utils.timer import timeMethod 

34 

35import lsst.utils 

36 

37__all__ = ["SpatiallySampledMetricsConfig", "SpatiallySampledMetricsTask"] 

38 

39 

40class SpatiallySampledMetricsConnections(pipeBase.PipelineTaskConnections, 

41 dimensions=("instrument", "visit", "detector"), 

42 defaultTemplates={"coaddName": "deep", 

43 "warpTypeSuffix": "", 

44 "fakesType": ""}): 

45 science = pipeBase.connectionTypes.Input( 

46 doc="Input science exposure.", 

47 dimensions=("instrument", "visit", "detector"), 

48 storageClass="ExposureF", 

49 name="{fakesType}calexp" 

50 ) 

51 matchedTemplate = pipeBase.connectionTypes.Input( 

52 doc="Warped and PSF-matched template used to create the difference image.", 

53 dimensions=("instrument", "visit", "detector"), 

54 storageClass="ExposureF", 

55 name="{fakesType}{coaddName}Diff_matchedExp", 

56 ) 

57 template = pipeBase.connectionTypes.Input( 

58 doc="Warped and not PSF-matched template used to create the difference image.", 

59 dimensions=("instrument", "visit", "detector"), 

60 storageClass="ExposureF", 

61 name="{fakesType}{coaddName}Diff_templateExp", 

62 ) 

63 difference = pipeBase.connectionTypes.Input( 

64 doc="Difference image with detection mask plane filled in.", 

65 dimensions=("instrument", "visit", "detector"), 

66 storageClass="ExposureF", 

67 name="{fakesType}{coaddName}Diff_differenceExp", 

68 ) 

69 diaSources = pipeBase.connectionTypes.Input( 

70 doc="Filtered diaSources on the difference image.", 

71 dimensions=("instrument", "visit", "detector"), 

72 storageClass="SourceCatalog", 

73 name="{fakesType}{coaddName}Diff_candidateDiaSrc", 

74 ) 

75 psfMatchingKernel = pipeBase.connectionTypes.Input( 

76 doc="Kernel used to PSF match the science and template images.", 

77 dimensions=("instrument", "visit", "detector"), 

78 storageClass="MatchingKernel", 

79 name="{fakesType}{coaddName}Diff_psfMatchKernel", 

80 ) 

81 spatiallySampledMetrics = pipeBase.connectionTypes.Output( 

82 doc="Summary metrics computed at randomized locations.", 

83 dimensions=("instrument", "visit", "detector"), 

84 storageClass="ArrowAstropy", 

85 name="{fakesType}{coaddName}Diff_spatiallySampledMetrics", 

86 ) 

87 

88 

89class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig, 

90 pipelineConnections=SpatiallySampledMetricsConnections): 

91 """Config for SpatiallySampledMetricsTask 

92 """ 

93 metricsMaskPlanes = lsst.pex.config.ListField( 

94 dtype=str, 

95 doc="List of mask planes to include in metrics", 

96 default=('BAD', 'CLIPPED', 'CR', 'DETECTED', 'DETECTED_NEGATIVE', 'EDGE', 

97 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED', 

98 'NO_DATA', 'REJECTED', 'SAT', 'SAT_TEMPLATE', 'SENSOR_EDGE', 'STREAK', 'SUSPECT', 

99 'UNMASKEDNAN', 

100 ), 

101 ) 

102 metricSources = pexConfig.ConfigurableField( 

103 target=SkyObjectsTask, 

104 doc="Generate QA metric sources", 

105 ) 

106 

107 def setDefaults(self): 

108 self.metricSources.avoidMask = ["NO_DATA", "EDGE"] 

109 

110 

111class SpatiallySampledMetricsTask(lsst.pipe.base.PipelineTask): 

112 """Detect and measure sources on a difference image. 

113 """ 

114 ConfigClass = SpatiallySampledMetricsConfig 

115 _DefaultName = "spatiallySampledMetrics" 

116 

117 def __init__(self, **kwargs): 

118 super().__init__(**kwargs) 

119 

120 self.makeSubtask("metricSources") 

121 self.schema = afwTable.SourceTable.makeMinimalSchema() 

122 self.schema.addField( 

123 "x", "F", 

124 "X location of the metric evaluation.", 

125 units="pixel") 

126 self.schema.addField( 

127 "y", "F", 

128 "Y location of the metric evaluation.", 

129 units="pixel") 

130 self.metricSources.skySourceKey = self.schema.addField("sky_source", type="Flag", 

131 doc="Metric evaluation objects.") 

132 self.schema.addField( 

133 "source_density", "F", 

134 "Density of diaSources at location.", 

135 units="count/degree^2") 

136 self.schema.addField( 

137 "dipole_density", "F", 

138 "Density of dipoles at location.", 

139 units="count/degree^2") 

140 self.schema.addField( 

141 "dipole_direction", "F", 

142 "Mean dipole orientation.", 

143 units="radian") 

144 self.schema.addField( 

145 "dipole_separation", "F", 

146 "Mean dipole separation.", 

147 units="pixel") 

148 self.schema.addField( 

149 "template_value", "F", 

150 "Median of template at location.", 

151 units="nJy") 

152 self.schema.addField( 

153 "science_value", "F", 

154 "Median of science at location.", 

155 units="nJy") 

156 self.schema.addField( 

157 "diffim_value", "F", 

158 "Median of diffim at location.", 

159 units="nJy") 

160 self.schema.addField( 

161 "science_psfSize", "F", 

162 "Width of the science image PSF at location.", 

163 units="pixel") 

164 self.schema.addField( 

165 "template_psfSize", "F", 

166 "Width of the template image PSF at location.", 

167 units="pixel") 

168 for maskPlane in self.config.metricsMaskPlanes: 

169 self.schema.addField( 

170 "%s_mask_fraction"%maskPlane.lower(), "F", 

171 "Fraction of pixels with %s mask"%maskPlane 

172 ) 

173 self.schema.addField( 

174 "psfMatchingKernel_sum", "F", 

175 "PSF matching kernel sum at location.") 

176 self.schema.addField( 

177 "psfMatchingKernel_dx", "F", 

178 "PSF matching kernel centroid offset in x at location.", 

179 units="pixel") 

180 self.schema.addField( 

181 "psfMatchingKernel_dy", "F", 

182 "PSF matching kernel centroid offset in y at location.", 

183 units="pixel") 

184 self.schema.addField( 

185 "psfMatchingKernel_length", "F", 

186 "PSF matching kernel centroid offset module.", 

187 units="arcsecond") 

188 self.schema.addField( 

189 "psfMatchingKernel_position_angle", "F", 

190 "PSF matching kernel centroid offset position angle.", 

191 units="radian") 

192 self.schema.addField( 

193 "psfMatchingKernel_direction", "F", 

194 "PSF matching kernel centroid offset direction in detector plane.", 

195 units="radian") 

196 

197 @timeMethod 

198 def run(self, science, matchedTemplate, template, difference, diaSources, psfMatchingKernel): 

199 """Calculate difference image metrics on specific locations across the images 

200 

201 Parameters 

202 ---------- 

203 science : `lsst.afw.image.ExposureF` 

204 Science exposure that the template was subtracted from. 

205 matchedTemplate : `lsst.afw.image.ExposureF` 

206 Warped and PSF-matched template that was used produce the 

207 difference image. 

208 template : `lsst.afw.image.ExposureF` 

209 Warped and non PSF-matched template that was used produce 

210 the difference image. 

211 difference : `lsst.afw.image.ExposureF` 

212 Result of subtracting template from the science image. 

213 diaSources : `lsst.afw.table.SourceCatalog` 

214 The catalog of detected sources. 

215 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel` 

216 The PSF matching kernel of the subtraction to evaluate. 

217 

218 Returns 

219 ------- 

220 results : `lsst.pipe.base.Struct` 

221 ``spatiallySampledMetrics`` : `astropy.table.Table` 

222 Image quality metrics spatially sampled locations. 

223 """ 

224 

225 idFactory = lsst.meas.base.IdGenerator().make_table_id_factory() 

226 

227 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema) 

228 spatiallySampledMetrics.getTable().setIdFactory(idFactory) 

229 

230 self.metricSources.run(mask=science.mask, seed=difference.info.id, catalog=spatiallySampledMetrics) 

231 

232 metricsMaskPlanes = [] 

233 for maskPlane in self.config.metricsMaskPlanes: 

234 try: 

235 metricsMaskPlanes.append(maskPlane) 

236 except InvalidParameterError: 

237 self.log.info("Unable to calculate metrics for mask plane %s: not in image"%maskPlane) 

238 

239 for src in spatiallySampledMetrics: 

240 self._evaluateLocalMetric(src, science, matchedTemplate, template, difference, diaSources, 

241 metricsMaskPlanes=metricsMaskPlanes, 

242 psfMatchingKernel=psfMatchingKernel) 

243 

244 return pipeBase.Struct(spatiallySampledMetrics=spatiallySampledMetrics.asAstropy()) 

245 

246 def _evaluateLocalMetric(self, src, science, matchedTemplate, template, difference, diaSources, 

247 metricsMaskPlanes, psfMatchingKernel): 

248 """Calculate image quality metrics at spatially sampled locations. 

249 

250 Parameters 

251 ---------- 

252 src : `lsst.afw.table.SourceRecord` 

253 The source record to be updated with metric calculations. 

254 diaSources : `lsst.afw.table.SourceCatalog` 

255 The catalog of detected sources. 

256 science : `lsst.afw.image.Exposure` 

257 The science image. 

258 matchedTemplate : `lsst.afw.image.Exposure` 

259 The reference image, warped and psf-matched to the science image. 

260 difference : `lsst.afw.image.Exposure` 

261 Result of subtracting template from the science image. 

262 metricsMaskPlanes : `list` of `str` 

263 Mask planes to calculate metrics from. 

264 psfMatchingKernel : `~lsst.afw.math.LinearCombinationKernel` 

265 The PSF matching kernel of the subtraction to evaluate. 

266 """ 

267 bbox = src.getFootprint().getBBox() 

268 pix = bbox.getCenter() 

269 src.set('science_psfSize', getPsfFwhm(science.psf, position=pix)) 

270 try: 

271 src.set('template_psfSize', getPsfFwhm(template.psf, position=pix)) 

272 except (InvalidParameterError, RangeError): 

273 src.set('template_psfSize', np.nan) 

274 

275 metricRegionSize = 100 

276 bbox.grow(metricRegionSize) 

277 bbox = bbox.clippedTo(science.getBBox()) 

278 nPix = bbox.getArea() 

279 pixScale = science.wcs.getPixelScale(bbox.getCenter()) 

280 area = nPix*pixScale.asDegrees()**2 

281 peak = src.getFootprint().getPeaks()[0] 

282 src.set('x', peak['i_x']) 

283 src.set('y', peak['i_y']) 

284 src.setCoord(science.wcs.pixelToSky(peak['i_x'], peak['i_y'])) 

285 selectSources = diaSources[bbox.contains(diaSources.getX(), diaSources.getY())] 

286 sourceDensity = len(selectSources)/area 

287 dipoleSources = selectSources[selectSources["ip_diffim_DipoleFit_classification"]] 

288 dipoleDensity = len(dipoleSources)/area 

289 

290 if dipoleSources: 

291 meanDipoleOrientation = angleMean(dipoleSources["ip_diffim_DipoleFit_orientation"]) 

292 src.set('dipole_direction', meanDipoleOrientation) 

293 meanDipoleSeparation = np.mean(dipoleSources["ip_diffim_DipoleFit_separation"]) 

294 src.set('dipole_separation', meanDipoleSeparation) 

295 

296 templateVal = np.median(matchedTemplate[bbox].image.array) 

297 scienceVal = np.median(science[bbox].image.array) 

298 diffimVal = np.median(difference[bbox].image.array) 

299 src.set('source_density', sourceDensity) 

300 src.set('dipole_density', dipoleDensity) 

301 src.set('template_value', templateVal) 

302 src.set('science_value', scienceVal) 

303 src.set('diffim_value', diffimVal) 

304 for maskPlane in metricsMaskPlanes: 

305 src.set("%s_mask_fraction"%maskPlane.lower(), 

306 evaluateMaskFraction(difference.mask[bbox], maskPlane) 

307 ) 

308 

309 krnlSum, dx, dy, direction, length = getKernelCenterDisplacement( 

310 psfMatchingKernel, src.get('x'), src.get('y')) 

311 

312 point1 = lsst.geom.SpherePoint( 

313 src.get('coord_ra'), src.get('coord_dec'), 

314 lsst.geom.radians) 

315 point2 = science.wcs.pixelToSky(src.get('x') + dx, src.get('y') + dy) 

316 bearing = point1.bearingTo(point2) 

317 pa_ref_angle = lsst.geom.Angle(np.pi/2, lsst.geom.radians) 

318 pa = pa_ref_angle - bearing 

319 # Wrap around to get Delta_RA from -pi to +pi 

320 pa = pa.wrapCtr() 

321 position_angle = pa.asRadians() 

322 

323 src.set('psfMatchingKernel_sum', krnlSum) 

324 src.set('psfMatchingKernel_dx', dx) 

325 src.set('psfMatchingKernel_dy', dy) 

326 src.set('psfMatchingKernel_length', length*pixScale.asArcseconds()) 

327 src.set('psfMatchingKernel_position_angle', position_angle) # in E of N position angle 

328 src.set('psfMatchingKernel_direction', direction) # direction offset in detector