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

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 template = pipeBase.connectionTypes.Input( 

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

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

54 storageClass="ExposureF", 

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

56 ) 

57 difference = pipeBase.connectionTypes.Input( 

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

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

60 storageClass="ExposureF", 

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

62 ) 

63 diaSources = pipeBase.connectionTypes.Input( 

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

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

66 storageClass="ArrowAstropy", 

67 name="{fakesType}dia_source_detector", 

68 ) 

69 psfMatchingKernel = pipeBase.connectionTypes.Input( 

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

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

72 storageClass="MatchingKernel", 

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

74 ) 

75 spatiallySampledMetrics = pipeBase.connectionTypes.Output( 

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

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

78 storageClass="ArrowAstropy", 

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

80 ) 

81 

82 

83class SpatiallySampledMetricsConfig(pipeBase.PipelineTaskConfig, 

84 pipelineConnections=SpatiallySampledMetricsConnections): 

85 """Config for SpatiallySampledMetricsTask 

86 """ 

87 metricsMaskPlanes = lsst.pex.config.ListField( 

88 dtype=str, 

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

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

91 'INEXACT_PSF', 'INJECTED', 'INJECTED_TEMPLATE', 'INTRP', 'NOT_DEBLENDED', 

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

93 'UNMASKEDNAN', 

94 ), 

95 ) 

96 metricSources = pexConfig.ConfigurableField( 

97 target=SkyObjectsTask, 

98 doc="Generate QA metric sources", 

99 ) 

100 

101 def setDefaults(self): 

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

103 

104 

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

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

107 """ 

108 ConfigClass = SpatiallySampledMetricsConfig 

109 _DefaultName = "spatiallySampledMetrics" 

110 

111 def __init__(self, **kwargs): 

112 super().__init__(**kwargs) 

113 

114 self.makeSubtask("metricSources") 

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

116 self.schema.addField( 

117 "x", "F", 

118 "X location of the metric evaluation.", 

119 units="pixel") 

120 self.schema.addField( 

121 "y", "F", 

122 "Y location of the metric evaluation.", 

123 units="pixel") 

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

125 doc="Metric evaluation objects.") 

126 self.schema.addField( 

127 "source_density", "F", 

128 "Density of diaSources at location.", 

129 units="count/degree^2") 

130 self.schema.addField( 

131 "dipole_density", "F", 

132 "Density of dipoles at location.", 

133 units="count/degree^2") 

134 self.schema.addField( 

135 "dipole_direction", "F", 

136 "Mean dipole orientation.", 

137 units="radian") 

138 self.schema.addField( 

139 "dipole_separation", "F", 

140 "Mean dipole separation.", 

141 units="pixel") 

142 self.schema.addField( 

143 "template_value", "F", 

144 "Median of template at location.", 

145 units="nJy") 

146 self.schema.addField( 

147 "science_value", "F", 

148 "Median of science at location.", 

149 units="nJy") 

150 self.schema.addField( 

151 "diffim_value", "F", 

152 "Median of diffim at location.", 

153 units="nJy") 

154 self.schema.addField( 

155 "science_psfSize", "F", 

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

157 units="pixel") 

158 self.schema.addField( 

159 "template_psfSize", "F", 

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

161 units="pixel") 

162 for maskPlane in self.config.metricsMaskPlanes: 

163 self.schema.addField( 

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

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

166 ) 

167 self.schema.addField( 

168 "psfMatchingKernel_sum", "F", 

169 "PSF matching kernel sum at location.") 

170 self.schema.addField( 

171 "psfMatchingKernel_dx", "F", 

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

173 units="pixel") 

174 self.schema.addField( 

175 "psfMatchingKernel_dy", "F", 

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

177 units="pixel") 

178 self.schema.addField( 

179 "psfMatchingKernel_length", "F", 

180 "PSF matching kernel centroid offset module.", 

181 units="arcsecond") 

182 self.schema.addField( 

183 "psfMatchingKernel_position_angle", "F", 

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

185 units="radian") 

186 self.schema.addField( 

187 "psfMatchingKernel_direction", "F", 

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

189 units="radian") 

190 

191 @timeMethod 

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

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

194 

195 Parameters 

196 ---------- 

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

198 Science exposure that the template was subtracted from. 

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

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

201 the difference image. 

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

203 Result of subtracting template from the science image. 

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

205 The catalog of detected sources. 

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

207 The PSF matching kernel of the subtraction to evaluate. 

208 

209 Returns 

210 ------- 

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

212 ``spatiallySampledMetrics`` : `astropy.table.Table` 

213 Image quality metrics spatially sampled locations. 

214 """ 

215 

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

217 

218 spatiallySampledMetrics = afwTable.SourceCatalog(self.schema) 

219 spatiallySampledMetrics.getTable().setIdFactory(idFactory) 

220 

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

222 

223 metricsMaskPlanes = [] 

224 for maskPlane in self.config.metricsMaskPlanes: 

225 try: 

226 metricsMaskPlanes.append(maskPlane) 

227 except InvalidParameterError: 

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

229 

230 for src in spatiallySampledMetrics: 

231 self._evaluateLocalMetric(src, science, template, difference, diaSources, 

232 metricsMaskPlanes=metricsMaskPlanes, 

233 psfMatchingKernel=psfMatchingKernel) 

234 

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

236 

237 def _evaluateLocalMetric(self, src, science, template, difference, diaSources, 

238 metricsMaskPlanes, psfMatchingKernel): 

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

240 

241 Parameters 

242 ---------- 

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

244 The source record to be updated with metric calculations. 

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

246 The catalog of detected sources. 

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

248 The science image. 

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

250 Result of subtracting template from the science image. 

251 metricsMaskPlanes : `list` of `str` 

252 Mask planes to calculate metrics from. 

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

254 The PSF matching kernel of the subtraction to evaluate. 

255 """ 

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

257 pix = bbox.getCenter() 

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

259 try: 

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

261 except (InvalidParameterError, RangeError): 

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

263 

264 metricRegionSize = 100 

265 bbox.grow(metricRegionSize) 

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

267 nPix = bbox.getArea() 

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

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

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

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

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

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

274 selectSources = diaSources[bbox.contains(diaSources['x'], diaSources['y'])] 

275 sourceDensity = len(selectSources)/area 

276 dipoleSources = selectSources[selectSources["isDipole"]] 

277 dipoleDensity = len(dipoleSources)/area 

278 

279 if dipoleSources: 

280 meanDipoleOrientation = angleMean(dipoleSources["dipoleAngle"]) 

281 src.set('dipole_direction', meanDipoleOrientation) 

282 meanDipoleSeparation = np.mean(dipoleSources["dipoleLength"]) 

283 src.set('dipole_separation', meanDipoleSeparation) 

284 

285 templateVal = np.median(template[bbox].image.array) 

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

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

288 src.set('source_density', sourceDensity) 

289 src.set('dipole_density', dipoleDensity) 

290 src.set('template_value', templateVal) 

291 src.set('science_value', scienceVal) 

292 src.set('diffim_value', diffimVal) 

293 for maskPlane in metricsMaskPlanes: 

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

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

296 ) 

297 

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

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

300 

301 point1 = lsst.geom.SpherePoint( 

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

303 lsst.geom.radians) 

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

305 bearing = point1.bearingTo(point2) 

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

307 pa = pa_ref_angle - bearing 

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

309 pa = pa.wrapCtr() 

310 position_angle = pa.asRadians() 

311 

312 src.set('psfMatchingKernel_sum', krnlSum) 

313 src.set('psfMatchingKernel_dx', dx) 

314 src.set('psfMatchingKernel_dy', dy) 

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

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

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